diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 9dd1ba6fc6e68a292bc5c6e307aee1c9be1bd2d6..2c8a45840c9d4b56f9e9f903b8d37b79525e0fe1 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -17,9 +17,12 @@ extern bool freq_map;
 //
 //****************************************************************************************
 int main(int argc, char *argv[]) {
-  const long seed = 1121;
-  iniranf(seed);
-  setrancut(2.0);
+  
+  /*set the random value for random generator*/
+  const long seed = 1121; //the default random seed number
+  iniranf(seed); //initialize the seed
+  setrancut(2.0); //default value of the normal cut for the normal distribution
+
 
   // turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
   // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
@@ -41,16 +44,24 @@ int main(int argc, char *argv[]) {
   long lastpos = -1L;
   char str1[S_SIZE];
   
+  // paramters to read the command line in user input script
+  long CommandNo = -1L;  //the number of commands, since this value is also
+                         // the index of array UserCommandFlag[], so this
+			 // value is always less than 1 in the really case. 
+  UserCommand UserCommandFlag[500];
  
 
   // for time handling
   uint32_t start, stop;
 
+ 
   /************************************************************************
    read in files and flags
    *************************************************************************/
+
+   
   if (argc > 1) {
-    read_script(argv[1], true);
+    read_script(argv[1], true,CommandNo, UserCommandFlag);
   } else {
     fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
     exit_(1);
@@ -76,9 +87,9 @@ int main(int argc, char *argv[]) {
   
   
    // Flag factory
-  for(i=0;i<CommandNo;i++){//read user defined command by sequence
+  for(i=0; i<=CommandNo; i++){//read user defined command by sequence
    //assign user defined command
-    strcpy(CommandStr,UserCommand[i]); 
+    strcpy(CommandStr,UserCommandFlag[i].CommandStr); 
     
     //print the twiss paramters in a file defined by the name
     if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
@@ -95,18 +106,16 @@ int main(int argc, char *argv[]) {
    }
     
   //set RF voltage  
-  //if (RFvoltageFlag == true) {
   else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
     printf("\nSetting RF voltage:\n");
     printf("    Old RF voltage is: %lf [MV]\n", get_RFVoltage(ElemIndex("cav"))
         / 1e6);
-    set_RFVoltage(ElemIndex("cav"), RFvolt);
+    set_RFVoltage(ElemIndex("cav"), UserCommandFlag[i].RFvolt);
     printf("    New RF voltage is: %lf [MV]\n", get_RFVoltage(ElemIndex("cav"))
         / 1e6);
   }
 
   // Chamber factory
- // if (ReadChamberFlag == true)
   else if(strcmp(CommandStr,"ReadChamberFlag") == 0) {
     ReadCh(chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
   PrintCh(); // print chamber into chamber.out
@@ -135,23 +144,24 @@ int main(int argc, char *argv[]) {
     double qf_bn_fit = 0.0, qd_bn_fit = 0.0;
 
     /* quadrupole field strength before fitting*/
-    qf_bn = Cell[Elem_GetPos(ElemIndex(qf), 1)].Elem.M->PB[HOMmax + 2];
-    qd_bn = Cell[Elem_GetPos(ElemIndex(qd), 1)].Elem.M->PB[HOMmax + 2];
+    qf_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
+    qd_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
 
     /* fitting tunes*/
     fprintf(stdout,
-        "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n", qf, qd,
-        targetnux, targetnuz);
-    FitTune(ElemIndex(qf), ElemIndex(qd), targetnux, targetnuz);
+        "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n", UserCommandFlag[i].qf, 
+	UserCommandFlag[i].qd,UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
+    FitTune(ElemIndex(UserCommandFlag[i].qf), ElemIndex(UserCommandFlag[i].qd), 
+            UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
 
     /* integrated field strength after fitting*/
-    qf_bn_fit = Cell[Elem_GetPos(ElemIndex(qf), 1)].Elem.M->PB[HOMmax + 2];
-    qd_bn_fit = Cell[Elem_GetPos(ElemIndex(qd), 1)].Elem.M->PB[HOMmax + 2];
+    qf_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
+    qd_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
     /* print out the quadrupole strengths before and after the fitting*/
     printf("Before the fitting, the quadrupole field strengths are: \n");
-    printf(" %s = %f,    %s = %f\n", qf, qf_bn, qd, qd_bn);
+    printf(" %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn);
     printf("After the fitting, the quadrupole field strengths are: \n");
-    printf(" %s = %f,    %s = %f\n", qf, qf_bn_fit, qd, qd_bn_fit);
+    printf(" %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit);
 
     /* Compute and get Twiss parameters */
     Ring_GetTwiss(chroma = true, 0.0);
@@ -166,31 +176,34 @@ int main(int argc, char *argv[]) {
         0.0;
 
     /* quadrupole field strength before fitting*/
-    qf1_bn = Cell[Elem_GetPos(ElemIndex(qf1), 1)].Elem.M->PB[HOMmax + 2];
-    qf2_bn = Cell[Elem_GetPos(ElemIndex(qf2), 1)].Elem.M->PB[HOMmax + 2];
-    qd1_bn = Cell[Elem_GetPos(ElemIndex(qd1), 1)].Elem.M->PB[HOMmax + 2];
-    qd2_bn = Cell[Elem_GetPos(ElemIndex(qd2), 1)].Elem.M->PB[HOMmax + 2];
+    qf1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
+    qf2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
+    qd1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
+    qd2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
 
     /* fitting tunes*/
     fprintf(
         stdout,
         "\n Fitting tunes for Soleil ring: %s %s %s %s, targetnux = %f, targetnuz = %f \n",
-        qf1, qf2, qd1, qd2, targetnux, targetnuz);
-    FitTune4(ElemIndex(qf1), ElemIndex(qf2), ElemIndex(qd1), ElemIndex(qd2),
-        targetnux, targetnuz);
+        UserCommandFlag[i].qf1, UserCommandFlag[i].qf2, UserCommandFlag[i].qd1, UserCommandFlag[i].qd2, 
+	UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
+    FitTune4(ElemIndex(UserCommandFlag[i].qf1), ElemIndex(UserCommandFlag[i].qf2),
+             ElemIndex(UserCommandFlag[i].qd1), ElemIndex(UserCommandFlag[i].qd2),
+             UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
 
     /* integrated field strength after fitting*/
-    qf1_bn_fit = Cell[Elem_GetPos(ElemIndex(qf1), 1)].Elem.M->PB[HOMmax + 2];
-    qf2_bn_fit = Cell[Elem_GetPos(ElemIndex(qf2), 1)].Elem.M->PB[HOMmax + 2];
-    qd1_bn_fit = Cell[Elem_GetPos(ElemIndex(qd1), 1)].Elem.M->PB[HOMmax + 2];
-    qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(qd2), 1)].Elem.M->PB[HOMmax + 2];
+    qf1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
+    qf2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
+    qd1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
+    qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
     /* print out the quadrupole strengths before and after the fitting*/
     printf("Before the fitting, the quadrupole field strengths are: \n");
-    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", qf1, qf1_bn,
-        qf2, qf2_bn, qd1, qd1_bn, qd2, qd2_bn);
+    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", UserCommandFlag[i].qf1, qf1_bn,
+        UserCommandFlag[i].qf2, qf2_bn, UserCommandFlag[i].qd1, qd1_bn, UserCommandFlag[i].qd2, qd2_bn);
     printf("After the fitting, the quadrupole field strengths are: \n");
-    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", qf1,
-        qf1_bn_fit, qf2, qf2_bn_fit, qd1, qd1_bn_fit, qd2, qd2_bn_fit);
+    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", UserCommandFlag[i].qf1,
+        qf1_bn_fit, UserCommandFlag[i].qf2, qf2_bn_fit, UserCommandFlag[i].qd1, qd1_bn_fit, 
+	UserCommandFlag[i].qd2, qd2_bn_fit);
 
     /* Compute and get Twiss parameters */
     Ring_GetTwiss(chroma = true, 0.0);
@@ -198,35 +211,33 @@ int main(int argc, char *argv[]) {
   }
 
   /* fit the chromaticities*/
- // if (FitChromFlag == true) {
   else if(strcmp(CommandStr,"FitChromFlag") == 0) {
   
     double sxm1_bn = 0.0, sxm2_bn = 0.0;
     double sxm1_bn_fit = 0.0, sxm2_bn_fit = 0.0;
-    sxm1_bn = Cell[Elem_GetPos(ElemIndex(sxm1), 1)].Elem.M->PB[HOMmax + 3];
-    sxm2_bn = Cell[Elem_GetPos(ElemIndex(sxm2), 1)].Elem.M->PB[HOMmax + 3];
-
-    /* fit the chromaticities*/
-    fprintf(
-        stdout,
-        "\n Fitting chromaticities: %s %s, targetksix = %f,  targetksiz = %f\n",
-        sxm1, sxm2, targetksix, targetksiz);
-    FitChrom(ElemIndex(sxm1), ElemIndex(sxm2), targetksix, targetksiz);
-
-    sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(sxm1), 1)].Elem.M->PB[HOMmax + 3];
-    sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(sxm2), 1)].Elem.M->PB[HOMmax + 3];
+    sxm1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
+    sxm2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
+
+    fprintf(stdout,"\n Fitting chromaticities: %s %s, targetksix = %f,  targetksiz = %f\n",
+        UserCommandFlag[i].sxm1, UserCommandFlag[i].sxm2, UserCommandFlag[i].targetksix,
+	 UserCommandFlag[i].targetksiz);
+	 
+    FitChrom(ElemIndex(UserCommandFlag[i].sxm1), ElemIndex(UserCommandFlag[i].sxm2), 
+             UserCommandFlag[i].targetksix, UserCommandFlag[i].targetksiz);
+
+    sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
+    sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
     /* print out the sextupole strengths before and after the fitting*/
     printf("Before the fitting, the sextupole field strengths are \n");
-    printf("    %s = %f,    %s = %f\n", sxm1, sxm1_bn, sxm2, sxm2_bn);
+    printf("    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn);
     printf("After the fitting, the sextupole field strengths are: \n");
-    printf("    %s = %f,    %s = %f\n", sxm1, sxm1_bn_fit, sxm2, sxm2_bn_fit);
+    printf("    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit);
 
     Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
     printglob(); /* print parameter list */
   }
 
   // coupling calculation
- // if (CouplingFlag == true) {
   else if(strcmp(CommandStr,"CouplingFlag") == 0) {
     Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
     printlatt("linlat_coupling.out"); /* dump linear lattice functions into "linlat.dat" */
@@ -240,7 +251,7 @@ int main(int argc, char *argv[]) {
   //if (ErrorCouplingFlag == true) {
   else if(strcmp(CommandStr,"ErrorCouplingFlag") == 0) {
     prtmfile("flat_file.dat"); //print the elements without rotation errors
-    SetErr(err_seed, err_rms);
+    SetErr(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
     prtmfile("flat_file_errcoupling_full.dat"); //print the elements with rotation errors
     Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
     printlatt("linlat_errcoupling.out"); /* dump linear lattice functions into "linlat.dat" */
@@ -251,11 +262,10 @@ int main(int argc, char *argv[]) {
     printglob(); /* print parameter list */
   }
   
-  // add coupling by random rotating of the half quadrupole magnets, delicated for soleil
-  //if (ErrorCoupling2Flag == true) {
+//  add coupling by random rotating of the half quadrupole magnets, delicated for soleil
   else if(strcmp(CommandStr,"ErrorCoupling2Flag") == 0) {
     prtmfile("flat_file.dat"); //print the elements without rotation errors
-    SetErr2(err_seed, err_rms);
+    SetErr2(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
     prtmfile("flat_file_errcoupling_half.dat"); //print the elements with rotation errors
     Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
     printlatt("linlat_errcoupling2.out"); /* dump linear lattice functions into "linlat.dat" */
@@ -265,13 +275,10 @@ int main(int argc, char *argv[]) {
     printglob(); /* print parameter list */
   }
 
-  // WARNING Fit tunes and chromaticities before applying errors !!!!
   //set multipoles in all magnets
   // read multipole error from a file
- // if (ReadMultipoleFlag == true) {
   else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
-    fprintf(stdout,
-        "\n Read Multipoles file for lattice with thick sextupoles \n");
+    fprintf(stdout,"\n Read Multipoles file for lattice with thick sextupoles \n");
     ReadFieldErr(multipole_file);
     //first print the full lattice with error as a flat file
     prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors  /* very important file for debug*/
@@ -279,67 +286,57 @@ int main(int argc, char *argv[]) {
     printglob();
   }
 
-  //Old version to multipole errors in soleil lattice, is obsoleted
- // if (MultipoleFlag == true) {
-  //  if (ThinsextFlag == true) {
-   //   fprintf(stdout,
-   //       "\n Setting Multipoles for lattice with thin sextupoles \n");
-   //   Multipole_thinsext(fic_hcorr, fic_vcorr, fic_skew); /* for thin sextupoles */
-   //   prtmfile("flat_file_errmultipole.dat"); // writes flat file   /* very important file for debug*/
-   //   Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
-   //   printglob();
-   // } else {
-   //   fprintf(stdout,
-   //       "\n Setting Multipoles for lattice with thick sextupoles \n");
-   //   Multipole_thicksext(fic_hcorr, fic_vcorr, fic_skew); /* for thick sextupoles */
-   //   prtmfile("flat_file_errmultipole.dat"); // writes flat file   /* very important file for debug*/
-   //   Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
-   //   printglob();
-   // }
- // }
+ 
 
   /******************************************************************************************/
   // COMPUTATION PART after setting the model
   /******************************************************************************************/
-
   // computes TuneShift with amplitudes
- // if (AmplitudeTuneShiftFlag == true) {
   else if(strcmp(CommandStr,"AmplitudeTuneShiftFlag") == 0) {
-      TunesShiftWithAmplitude(_AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint,
-          _AmplitudeTuneShift_nturn, _AmplitudeTuneShift_xmax,
-          _AmplitudeTuneShift_ymax, _AmplitudeTuneShift_delta);
+      TunesShiftWithAmplitude(UserCommandFlag[i]._AmplitudeTuneShift_nxpoint,
+                              UserCommandFlag[i]._AmplitudeTuneShift_nypoint,
+                              UserCommandFlag[i]._AmplitudeTuneShift_nturn, 
+			      UserCommandFlag[i]._AmplitudeTuneShift_xmax,
+                              UserCommandFlag[i]._AmplitudeTuneShift_ymax, 
+			      UserCommandFlag[i]._AmplitudeTuneShift_delta);
     } 
-   //compute tuneshift with energy
-  //if (EnergyTuneShiftFlag == true) {
-  else if(strcmp(CommandStr,"EnergyTuneShiftFlag") == 0) {
-      TunesShiftWithEnergy(_EnergyTuneShift_npoint, _EnergyTuneShift_nturn,
-          _EnergyTuneShift_deltamax);
-    }
+  // compute tuneshift with energy
+ else if(strcmp(CommandStr,"EnergyTuneShiftFlag") == 0) {
+     TunesShiftWithEnergy(UserCommandFlag[i]._EnergyTuneShift_npoint, 
+                          UserCommandFlag[i]._EnergyTuneShift_nturn,
+                          UserCommandFlag[i]._EnergyTuneShift_deltamax);
+   }
 
   // Computes FMA
-  //if (FmapFlag == true) {
   else if(strcmp(CommandStr,"FmapFlag") == 0) {
     printf("\n begin Fmap calculation for on momentum particles: \n");
-    fmap(_FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn, _FmapFlag_xmax,
-        _FmapFlag_ymax, _FmapFlag_delta, _FmapFlag_diffusion);
+    fmap(UserCommandFlag[i]._FmapFlag_nxpoint, 
+         UserCommandFlag[i]._FmapFlag_nypoint, 
+	 UserCommandFlag[i]._FmapFlag_nturn, 
+	 UserCommandFlag[i]._FmapFlag_xmax,
+         UserCommandFlag[i]._FmapFlag_ymax, 
+	 UserCommandFlag[i]._FmapFlag_delta, 
+	 UserCommandFlag[i]._FmapFlag_diffusion);
   }
 
   // Compute FMA dp
- // if (FmapdpFlag == true) {
   else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
     printf("\n begin Fmap calculation for off momentum particles: \n");
-    fmapdp(_FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn,
-        _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z,
-        _FmapdpFlag_diffusion);
+    fmapdp(UserCommandFlag[i]._FmapdpFlag_nxpoint, 
+           UserCommandFlag[i]._FmapdpFlag_nepoint, 
+	   UserCommandFlag[i]._FmapdpFlag_nturn,
+           UserCommandFlag[i]._FmapdpFlag_xmax, 
+	   UserCommandFlag[i]._FmapdpFlag_emax, 
+	   UserCommandFlag[i]._FmapdpFlag_z,
+           UserCommandFlag[i]._FmapdpFlag_diffusion);
   }
 
- // if (CodeComparaisonFlag) {
-  else if(strcmp(CommandStr,"CodeComparaisonFlag") == 0) {
-    fmap(200, 100, 1026, -32e-3, 7e-3, 0.0, true);
-  }
+//  // if (CodeComparaisonFlag) {
+//   else if(strcmp(CommandStr,"CodeComparaisonFlag") == 0) {
+//     fmap(200, 100, 1026, -32e-3, 7e-3, 0.0, true);
+//   }
 
   // MOMENTUM ACCEPTANCE
- // if (MomentumAccFlag == true) {
   else if(strcmp(CommandStr,"MomentumAccFlag") == 0) {
     bool cavityflag, radiationflag;
     /* record the initial values*/
@@ -347,10 +344,10 @@ int main(int argc, char *argv[]) {
     radiationflag = globval.radiation;
 
     /* set the dimension for the momentum tracking*/
-    if (strncmp("6D", TrackDim, 2) == 0) {
+    if (strncmp("6D", UserCommandFlag[i].TrackDim, 2) == 0) {
       globval.Cavity_on = true;
       globval.radiation = true;
-    } else if (strncmp("4D", TrackDim, 2) == 0) {
+    } else if (strncmp("4D", UserCommandFlag[i].TrackDim, 2) == 0) {
       globval.Cavity_on = false;
       globval.radiation = false;
     } else {
@@ -360,10 +357,14 @@ int main(int argc, char *argv[]) {
 
     /* calculate momentum acceptance*/
     printf("\n Calculate momentum acceptance: \n");
-    MomentumAcceptance(_MomentumAccFlag_istart, _MomentumAccFlag_istop,
-        _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp,
-        _MomentumAccFlag_nstepp, _MomentumAccFlag_deltaminn,
-        _MomentumAccFlag_deltamaxn, _MomentumAccFlag_nstepn);
+    MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_istart, 
+                       UserCommandFlag[i]._MomentumAccFlag_istop,
+                       UserCommandFlag[i]._MomentumAccFlag_deltaminp, 
+		       UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
+                       UserCommandFlag[i]._MomentumAccFlag_nstepp, 
+		       UserCommandFlag[i]._MomentumAccFlag_deltaminn,
+                       UserCommandFlag[i]._MomentumAccFlag_deltamaxn, 
+		       UserCommandFlag[i]._MomentumAccFlag_nstepn);
 
     /* restore the initial values*/
     globval.Cavity_on = cavityflag;
@@ -371,13 +372,12 @@ int main(int argc, char *argv[]) {
   }
 
   // induced amplitude 
- // if (InducedAmplitudeFlag == true) {
   else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) {
    printf("\n Calculate induced amplitude: \n");
     InducedAmplitude(193L);
   }
 
- // if (EtaFlag == true) {
+
   else if(strcmp(CommandStr,"EtaFlag") == 0) {
     // compute cod and Twiss parameters for different energy offsets
     for (int ii = 0; ii <= 40; ii++) {
@@ -395,7 +395,7 @@ int main(int argc, char *argv[]) {
     }
   }
 
-  //if (PhaseSpaceFlag == true) {
+// track to get phase space
   else if(strcmp(CommandStr,"PhaseSpaceFlag") == 0) {
     bool cavityflag, radiationflag;
     /* record the initial values*/
@@ -403,23 +403,29 @@ int main(int argc, char *argv[]) {
     radiationflag = globval.radiation;
 
     /* set the dimension for the momentum tracking*/
-    if (strncmp("6D", _Phase_Dim, 2) == 0) {
+    if (strncmp("6D", UserCommandFlag[i]._Phase_Dim, 2) == 0) {
       globval.Cavity_on = true;
-    } else if (strncmp("4D", _Phase_Dim, 2) == 0) {
+    } else if (strncmp("4D", UserCommandFlag[i]._Phase_Dim, 2) == 0) {
       globval.Cavity_on = false;
     } else {
       printf("MomentumAccFlag: Error!!! _Phase_Dim must be '4D' or '6D'\n");
       exit_(1);
     };
     /* setting damping */
-    if (_Phase_Damping == true) {
+    if ( UserCommandFlag[i]._Phase_Damping == true) {
       globval.radiation = true;
     } else  {
       globval.radiation = false;
     }
 
     start = stampstart();
-    Phase(_Phase_X, _Phase_Px, _Phase_Y, _Phase_Py, _Phase_delta, _Phase_ctau, _Phase_nturn);
+    Phase(UserCommandFlag[i]._Phase_X, 
+          UserCommandFlag[i]._Phase_Px, 
+	  UserCommandFlag[i]._Phase_Y, 
+	  UserCommandFlag[i]._Phase_Py, 
+	  UserCommandFlag[i]._Phase_delta, 
+	  UserCommandFlag[i]._Phase_ctau, 
+	  UserCommandFlag[i]._Phase_nturn);
     printf("the simulation time for phase space in tracy 3 is \n");
     stop = stampstop(start);
 
@@ -428,11 +434,11 @@ int main(int argc, char *argv[]) {
     globval.radiation = radiationflag;
   }
   
- else
-    printf("Wrong!!!!!");
- }//end of looking for user defined flag
+ 
   
 
+ else if (strcmp(CommandStr,"TouschekFlag") == 0 ||strcmp(CommandStr,"IBSFlag") == 0 ||
+          strcmp(CommandStr,"TousTrackFlag") == 0 ){ 
   //  ????????????? NSRL version, Check with Soleil version "MomentumAcceptance"
   // IBS & TOUSCHEK
   int k, n_turns;
@@ -440,7 +446,7 @@ int main(int argc, char *argv[]) {
   FILE *outf;
   const double Qb = 5e-9; // e charge in one bunch
 
-  if (TouschekFlag == true) {
+  
 //  else if(strcmp(CommandStr,"TouschekFlag") == 0) {
     double sum_delta[globval.Cell_nLoc + 1][2];
     double sum2_delta[globval.Cell_nLoc + 1][2];
@@ -485,7 +491,7 @@ int main(int argc, char *argv[]) {
         sigma_delta, sigma_s);
 
  // Intra Beam Scattering(IBS)
-    if (IBSFlag == true) {
+    if (strcmp(CommandStr,"IBSFlag") == 0) {
       // initialize eps_IBS with eps_SR
       for (k = 0; k < 3; k++)
         eps[k] = globval.eps[k];
@@ -500,7 +506,7 @@ int main(int argc, char *argv[]) {
     // at each element location which is tracked over n turns, 
     //the vacuum chamber is read from the file "chamber_file"  
     // finally, write the momentum acceptance to file "mom_aper.out".
-    if (TousTrackFlag == true) {
+    if (strcmp(CommandStr,"TousTrackFlag") == 0) {
       globval.Aperture_on = true;
       ReadCh(chamber_file);
       //      LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
@@ -518,6 +524,9 @@ int main(int argc, char *argv[]) {
     }
 
   }
+  else
+    printf("Wrong!!!!!");
+ }//end of looking for user defined flag
 
   
   return 0;
diff --git a/tracy/tracy/inc/read_script.h b/tracy/tracy/inc/read_script.h
index 432d321be4e0bf60f264e1908d376198578266e4..58533977508bf67d789fb4c916477e2208d5a151 100644
--- a/tracy/tracy/inc/read_script.h
+++ b/tracy/tracy/inc/read_script.h
@@ -1,80 +1,118 @@
-//extern long max_line;
-extern char UserCommand[500][max_str];  //string array with bool flag
-extern long CommandNo;
-
-/* global flag, used in script_read() and main() for read the input from script*/
-
-//extern bool QuadFringeOnFlag;
-
-//extern bool RFvoltageFlag; 
-extern double RFvolt;
-
-//extern bool TuneTracFlag;
-//extern bool ChromTracFlag ;
-
-//extern bool FmapFlag;
-extern long _FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn;
-extern double _FmapFlag_xmax, _FmapFlag_ymax, _FmapFlag_delta;
-extern bool _FmapFlag_diffusion;
-
-//extern bool FmapdpFlag;
-extern long _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn;
-extern double _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z;
-extern bool _FmapdpFlag_diffusion;	
-
-//extern bool AmplitudeTuneShiftFlag;
-extern long _AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint;
-extern long _AmplitudeTuneShift_nturn;
-extern double _AmplitudeTuneShift_xmax, _AmplitudeTuneShift_ymax, _AmplitudeTuneShift_delta;
-
-//extern bool EnergyTuneShiftFlag;
-extern long _EnergyTuneShift_npoint, _EnergyTuneShift_nturn;
-extern double _EnergyTuneShift_deltamax;
-
-
-//extern bool MomentumAccFlag;
-extern char TrackDim[3];
-extern long  _MomentumAccFlag_istart, _MomentumAccFlag_istop,
-            _MomentumAccFlag_nstepn, _MomentumAccFlag_nstepp;
-extern double _MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn;
-extern double _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp;
+/* class to read the parameters of the bool flags in the user input script */
+ class UserCommand
+ {
+   public:
+   char  CommandStr[max_str];   // bool flag name
+   
+   /* parameters */
+   double RFvolt;   // RF voltage
+    // FitTuneFlag ; 
+   char qf[max_str],qd[max_str]; 
+   double targetnux, targetnuz;
+   // FitTune4Flag  ; 
+    char qf1[max_str],qf2[max_str],qd1[max_str],qd2[max_str]; 
+ // FitChromFlag ; 
+   char sxm1[max_str],sxm2[max_str]; 
+   double targetksix, targetksiz ;
+  //add coupling error to full/half quadrupoles
+   long err_seed; 
+   double err_rms;
+  //AmplitudeTuneShiftFlag;
+   long _AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint;
+   long _AmplitudeTuneShift_nturn;
+   double _AmplitudeTuneShift_xmax, _AmplitudeTuneShift_ymax, _AmplitudeTuneShift_delta;
+ 
+ //EnergyTuneShiftFlag;
+   long _EnergyTuneShift_npoint, _EnergyTuneShift_nturn;
+   double _EnergyTuneShift_deltamax;
+ 
+ //extern bool FmapFlag;
+   long _FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn;
+   double _FmapFlag_xmax, _FmapFlag_ymax, _FmapFlag_delta;
+   bool _FmapFlag_diffusion; 
+ 
+ //extern bool FmapdpFlag;
+   long _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn;
+   double _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z;
+   bool _FmapdpFlag_diffusion;	
+
+   
+  //MomentumAccFlag;
+   char TrackDim[3];
+   long  _MomentumAccFlag_istart, _MomentumAccFlag_istop,
+         _MomentumAccFlag_nstepn, _MomentumAccFlag_nstepp;
+   double _MomentumAccFlag_deltaminn, _MomentumAccFlag_deltamaxn;
+   double _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp;
+ 
+  // /* Phase space */
+   double _Phase_X, _Phase_Px, _Phase_Y, _Phase_Py,_Phase_delta, _Phase_ctau;
+   long _Phase_nturn;
+   char _Phase_Dim[3];
+   bool _Phase_Damping;
+ 
+   //Touschek lifetime
+   bool TouschekFlag, IBSFlag, TousTrackFlag;
+   
+    
+    //set default values
+  UserCommand(void)   //constructor
+ {  
+  
+  // /* fmap for on momentum particle*/
+   _FmapFlag_nxpoint=31L, _FmapFlag_nypoint=21L, _FmapFlag_nturn=516L;
+   _FmapFlag_xmax=0.025, _FmapFlag_ymax=0.005, _FmapFlag_delta=0.0;
+   _FmapFlag_diffusion = true;
+
+/*fmap for off momentum particle*/
+ _FmapdpFlag_nxpoint=31L, _FmapdpFlag_nepoint=21L, _FmapdpFlag_nturn=516L;
+ _FmapdpFlag_xmax=0.025, _FmapdpFlag_emax=0.005, _FmapdpFlag_z=0.0;
+ _FmapdpFlag_diffusion = true;			
+
+/* tune shift with amplitude*/
+ _AmplitudeTuneShift_nxpoint=31L;  _AmplitudeTuneShift_nypoint=21L;
+ _AmplitudeTuneShift_nturn=516L;   _AmplitudeTuneShift_xmax=0.025;
+ _AmplitudeTuneShift_ymax=0.005,   _AmplitudeTuneShift_delta=0.0;
+
+/* tune shift with energy*/
+ _EnergyTuneShift_npoint=31L;  _EnergyTuneShift_nturn=516L;
+ _EnergyTuneShift_deltamax=0.06;
+
+/* random rotation coupling error*/
+ err_seed=0L;  err_rms=0.0;  
+
+
+
+/* momentum acceptance */
+ TrackDim[3] = '6D';
+ _MomentumAccFlag_istart=1L, _MomentumAccFlag_istop=108L,
+ _MomentumAccFlag_nstepn=100L, _MomentumAccFlag_nstepp=100L;
+ _MomentumAccFlag_deltaminn=-0.01, _MomentumAccFlag_deltamaxn=-0.05;
+ _MomentumAccFlag_deltaminp=0.01, _MomentumAccFlag_deltamaxp=0.05;
 
 /* Phase space */
-extern double _Phase_X, _Phase_Px, _Phase_Y, _Phase_Py,
-    _Phase_delta, _Phase_ctau;
-extern long _Phase_nturn;
-extern char _Phase_Dim[3];
-extern bool _Phase_Damping;
-
-//extern bool ErrorCouplingFlag ; 
-extern long err_seed; extern double err_rms;
-//extern bool ErrorCoupling2Flag ;
-//extern bool CouplingFlag ;
-
-extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
-extern char multipole_file[max_str];
-//extern bool ReadMultipoleFlag;
-//extern bool MultipoleFlag , ThinsextFlag ;
-
-//extern bool FitTuneFlag ; 
-extern char qf[max_str],qd[max_str]; extern double targetnux, targetnuz;
-//extern bool FitTune4Flag  ; 
-extern char qf1[max_str],qf2[max_str],qd1[max_str],qd2[max_str]; 
-//extern bool FitChromFlag ; 
-extern char sxm1[max_str],sxm2[max_str]; extern double targetksix, targetksiz ;
-
-//extern bool ChamberFlag , 
-//extern bool ReadChamberFlag;
-//extern bool GirderErrorFlag ;
-//extern bool SigmaFlag ;
-//extern bool PX2Flag;
-//extern bool InducedAmplitudeFlag;
-//extern bool CodeComparaisonFlag;
-//extern bool EtaFlag;
-//extern bool PhaseSpaceFlag;  
-extern bool TouschekFlag, IBSFlag, TousTrackFlag;
-extern char chamber_file[max_str];
-extern char twiss_file[max_str];
-
-
-void read_script(const char *param_file_name, bool rd_lat);
+ _Phase_X=0.0, _Phase_Px=0.0, _Phase_Y=0.0, _Phase_Py=0.0,
+ _Phase_delta=0.0, _Phase_ctau=0.0;
+ _Phase_nturn=512L;
+ _Phase_Dim[3]='4D';
+ _Phase_Damping = false;
+
+  
+/* fit tunes for full quadrupole*/
+ targetnux = 0.0, targetnuz = 0.0;
+
+/* fit charomaticities*/
+ targetksix = 0.0, targetksiz = 0.0;
+};  
+ 
+ };
+
+// file names
+ extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
+ extern char multipole_file[max_str];
+ 
+
+ extern char chamber_file[max_str];
+ extern char twiss_file[max_str];
+
+// function
+ void read_script(const char *param_file_name, bool rd_lat,long& CommNo, UserCommand UserCommandFlag[]);
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index af85b6375f4acd5e828c81042e28d3e07a48d234..1ee1369cb134c9b3416c22726425eabcefc82ca3 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -23,99 +23,26 @@
        Read_Lattice()
 
    Comments:
-       Jianfeng Zhang 07/2010 soleil
+       Jianfeng Zhang 03/2011 soleil
 
 ****************************************************************************/
-/* The string array stored the command which are read from user input script*/
-//long max_line = 500L;
-//#define max_line 500;
-char UserCommand[500][max_str];  //string array which store user command flag
-long CommandNo; //number of bool flags in the user script
-
-
-/* global flag, used in script_read() and main() for read the input from script*/
-//bool QuadFringeOnFlag = false;
-
-//bool RFvoltageFlag =false; 
-double RFvolt =0.0;
-
-//char TuneTracFlag;
-//char ChromTracFlag;
-
-//bool FmapFlag = false;
-long _FmapFlag_nxpoint=31L, _FmapFlag_nypoint=21L, _FmapFlag_nturn=516L;
-double _FmapFlag_xmax=0.025, _FmapFlag_ymax=0.005, _FmapFlag_delta=0.0;
-bool _FmapFlag_diffusion = true;
-
-//bool FmapdpFlag = false;
-long _FmapdpFlag_nxpoint=31L, _FmapdpFlag_nepoint=21L, _FmapdpFlag_nturn=516L;
-double _FmapdpFlag_xmax=0.025, _FmapdpFlag_emax=0.005, _FmapdpFlag_z=0.0;
-bool _FmapdpFlag_diffusion = true;			
-
-//bool AmplitudeTuneShiftFlag = false;
-long _AmplitudeTuneShift_nxpoint=31L; long _AmplitudeTuneShift_nypoint=21L;
-long _AmplitudeTuneShift_nturn=516L; double _AmplitudeTuneShift_xmax=0.025;
-double _AmplitudeTuneShift_ymax=0.005, _AmplitudeTuneShift_delta=0.0;
-
-//bool EnergyTuneShiftFlag = false;
-long _EnergyTuneShift_npoint=31L; long _EnergyTuneShift_nturn=516L;
-double _EnergyTuneShift_deltamax=0.06;
-
-//bool ErrorCouplingFlag = false; 
-long err_seed=0L; double err_rms=0.0;  
-//bool ErrorCoupling2Flag = false;
-//bool CouplingFlag = false;
-
-
-//bool MomentumAccFlag = false;
-char TrackDim[3] = "6D";
-long  _MomentumAccFlag_istart=1L, _MomentumAccFlag_istop=108L,
-     _MomentumAccFlag_nstepn=100L, _MomentumAccFlag_nstepp=100L;
-double _MomentumAccFlag_deltaminn=-0.01, _MomentumAccFlag_deltamaxn=-0.05;
-double _MomentumAccFlag_deltaminp=0.01, _MomentumAccFlag_deltamaxp=0.05;
-
-/* Phase space */
-double _Phase_X=0.0, _Phase_Px=0.0, _Phase_Y=0.0, _Phase_Py=0.0,
-    _Phase_delta=0.0, _Phase_ctau=0.0;
-long _Phase_nturn=512L;
-char _Phase_Dim[3]="4D";
-bool _Phase_Damping = false;
-
-//bool ReadMultipoleFlag = false;
-//bool MultipoleFlag = false, ThinsextFlag = false;
 
+/* files */ 
+ //twiss file
+char twiss_file[max_str];
+ //chamber file
+char chamber_file[max_str];
+ // multipole files
+char multipole_file[max_str];
+  // multipole file for soleil lattice 
 char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
 
-/*generic function to fit tunes using 1 family of quadrupoles*/
-//bool FitTuneFlag  = false; 
-char qf[max_str],qd[max_str]; double targetnux = 0.0, targetnuz = 0.0;
-/* specific for soleil lattice in which the quadrupole is cut into 2 parts*/
-//bool FitTune4Flag  = false; 
-char qf1[max_str],qf2[max_str],qd1[max_str],qd2[max_str];
-
-//bool FitChromFlag = false; 
-char sxm1[max_str],sxm2[max_str]; double targetksix = 0.0, targetksiz = 0.0;
-//bool ChamberFlag = false, 
-//bool ReadChamberFlag = false;
-//bool GirderErrorFlag = false;
-//bool SigmaFlag = false;
-//bool PX2Flag = false;
-//bool InducedAmplitudeFlag = false;
-//bool CodeComparaisonFlag = false;
-//bool EtaFlag = false;           
-//bool PhaseSpaceFlag = false;
-bool TouschekFlag = false, IBSFlag = false, TousTrackFlag = false;
 
-char chamber_file[max_str];
-char multipole_file[max_str];
-char twiss_file[max_str];
 #define OLD_LATTICE
 
 /*  Read script   */
-void read_script(const char *param_file_name, bool rd_lat)
+void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCommand UserCommandFlag[])
  {
-  //initialize bool flag number in the user script
-  CommandNo = 0L; 
  
   
   char    str[max_str]="voidstring", dummy[max_str]="voidstring";
@@ -134,9 +61,7 @@ void read_script(const char *param_file_name, bool rd_lat)
  //bool TuneTracFlag;
   char *pch;
   
-  
-  
-
+    
   // Manipulation of the parameter file
   strcpy(dummy, param_file_name);
   pch = strstr (dummy,".prm"); // search for extension .prm
@@ -174,6 +99,7 @@ void read_script(const char *param_file_name, bool rd_lat)
       sscanf(line, "%s", name);
       
       
+      
       //find the sequence of the bool flag in user input script 
       NameLen = strlen(name);
       EndName[0] = name[NameLen-4];
@@ -183,8 +109,7 @@ void read_script(const char *param_file_name, bool rd_lat)
       
       //find the bool flag whose last 4 character are 'Flag' 
       if(strcmp(EndName,"Flag")==0)
-        CommandNo++;
-      
+           CommNo++;
       
       
       
@@ -225,149 +150,130 @@ void read_script(const char *param_file_name, bool rd_lat)
       } 
       
       /* read in bool flags */
-      //print twiss parameters flag flag
+       //print twiss parameters flag flag
       else if (strcmp("PrintTwissFlag", name) == 0){
           sscanf(line, "%*s %s", twiss_file);
-	  strcpy(UserCommand[CommandNo-1],name);
+	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
       } 
       //read chamber file flat
-     else if (strcmp("ReadChamberFlag", name) == 0){
-	  strcpy(UserCommand[CommandNo-1],name);
-      } 
-    // else if (strcmp("ReadChamberFlag", name) == 0){
-     //    ReadChamberFlag = true;
-	
-    // }
-     // else if (strcmp("QuadFringeOnFlag", name) == 0){
-      //  	globval.quad_fringe = true;
-     // }
+      else if (strcmp("ReadChamberFlag", name) == 0){
+ 	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
+       } 
+   
       else if (strcmp("QuadFringeOnFlag", name) == 0){
-	  strcpy(UserCommand[CommandNo-1],name);
+	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
       
-    //  else if (strcmp("RFvoltageFlag", name) == 0){
-     //     RFvoltageFlag = true;
-    //	  sscanf(line, "%*s %lf", &RFvolt);
-     // }
+    
       else if (strcmp("RFvoltageFlag", name) == 0){
-	  sscanf(line, "%*s %lf", &RFvolt);
-	  strcpy(UserCommand[CommandNo-1],name);
+	  sscanf(line, "%*s %lf", &(UserCommandFlag[CommNo].RFvolt));
+	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
       
       else if (strcmp("TuneTracFlag", name) == 0){
-        //  TuneTracFlag = true;
-	  strcpy(UserCommand[CommandNo-1],name);
+	 strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
       else if (strcmp("ChromTracFlag", name) == 0){
-        //  ChromTracFlag = true;
-	   strcpy(UserCommand[CommandNo-1],name);
-	   cout << UserCommand[CommandNo-1] << "\n";
+	   strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
       // FMA
-      else if (strcmp("FmapFlag", name) == 0){
-      //else if (strcmp("FmapFlag", name) == 0){
-          
+      else if (strcmp("FmapFlag", name) == 0){        
     	  strcpy(dummy, "");
           sscanf(line, "%*s %ld %ld %ld %lf %lf %lf %s", 
-          		&_FmapFlag_nxpoint, &_FmapFlag_nypoint,
-          		&_FmapFlag_nturn, &_FmapFlag_xmax, &_FmapFlag_ymax,
-          		&_FmapFlag_delta, dummy);
+          		&(UserCommandFlag[CommNo]._FmapFlag_nxpoint), 
+			&(UserCommandFlag[CommNo]._FmapFlag_nypoint),
+          		&(UserCommandFlag[CommNo]._FmapFlag_nturn), 
+			&(UserCommandFlag[CommNo]._FmapFlag_xmax), 
+			&(UserCommandFlag[CommNo]._FmapFlag_ymax),
+          		&(UserCommandFlag[CommNo]._FmapFlag_delta), 
+			dummy);
        
           // FmapFlag = true;
-         strcpy(UserCommand[CommandNo-1],name);
+         strcpy(UserCommandFlag[CommNo].CommandStr,name);
 	  
         if(strcmp(dummy, "true") == 0)
-          _FmapFlag_diffusion = true;
+          UserCommandFlag[CommNo]._FmapFlag_diffusion = true;
         else if(strcmp(dummy, "false") == 0)
-        	_FmapFlag_diffusion = false;
+          UserCommandFlag[CommNo]._FmapFlag_diffusion = false;
         else {
           printf("set boolean flag true or false for FMA Diffusion (value is %s) \n", dummy);
           exit_(1);
         }
       }
       // FMA dp
-      else if (strcmp("FmapdpFlag", name) == 0){
-     // else if (strcmp("FmapdpFlag", name) == 0){
-          
+      else if (strcmp("FmapdpFlag", name) == 0){         
     	  strcpy(dummy, "");
           sscanf(line, "%*s %ld %ld %ld %lf %lf %lf %s", 
-          		&_FmapdpFlag_nxpoint, &_FmapdpFlag_nepoint,
-          		&_FmapdpFlag_nturn, &_FmapdpFlag_xmax, &_FmapdpFlag_emax,
-          		&_FmapdpFlag_z, dummy);
+          		&(UserCommandFlag[CommNo]._FmapdpFlag_nxpoint),
+			&(UserCommandFlag[CommNo]._FmapdpFlag_nepoint),
+          		&(UserCommandFlag[CommNo]._FmapdpFlag_nturn), 
+			&(UserCommandFlag[CommNo]._FmapdpFlag_xmax), 
+			&(UserCommandFlag[CommNo]._FmapdpFlag_emax),
+          		&(UserCommandFlag[CommNo]._FmapdpFlag_z), 
+			dummy);
       
          //  FmapdpFlag = true;
-         strcpy(UserCommand[CommandNo-1],name);
+         strcpy(UserCommandFlag[CommNo].CommandStr,name);
 	 
         if(strcmp(dummy, "true") == 0)
-          _FmapdpFlag_diffusion = true;
+          UserCommandFlag[CommNo]._FmapdpFlag_diffusion = true;
         else if(strcmp(dummy, "false") == 0)
-        	_FmapdpFlag_diffusion = false;
+          UserCommandFlag[CommNo]._FmapdpFlag_diffusion = false;
         else {
           printf("set boolean flag true or false for FMAdp Diffusion (value is %s) \n", dummy);
           exit_(1);
         }
       }
-      else if (strcmp("AmplitudeTuneShiftFlag", name) == 0){
-      strcpy(UserCommand[CommandNo-1],name);
-        sscanf(line, "%*s %ld %ld %ld %lf %lf %lf", &_AmplitudeTuneShift_nxpoint,
-        		&_AmplitudeTuneShift_nypoint, &_AmplitudeTuneShift_nturn,
-        		&_AmplitudeTuneShift_xmax, &_AmplitudeTuneShift_ymax,
-        		&_AmplitudeTuneShift_delta);
-   
-        //	AmplitudeTuneShiftFlag = true;
-       
+      else if (strcmp("AmplitudeTuneShiftFlag", name) == 0){   
+        sscanf(line, "%*s %ld %ld %ld %lf %lf %lf", 
+	       &(UserCommandFlag[CommNo]._AmplitudeTuneShift_nxpoint),
+	       &(UserCommandFlag[CommNo]._AmplitudeTuneShift_nypoint), 
+	       &(UserCommandFlag[CommNo]._AmplitudeTuneShift_nturn),
+               &(UserCommandFlag[CommNo]._AmplitudeTuneShift_xmax),
+	       &(UserCommandFlag[CommNo]._AmplitudeTuneShift_ymax),
+               &(UserCommandFlag[CommNo]._AmplitudeTuneShift_delta));
+       strcpy(UserCommandFlag[CommNo].CommandStr,name);    
       }
       else if (strcmp("EnergyTuneShiftFlag", name) == 0){
-      strcpy(UserCommand[CommandNo-1],name);
-        sscanf(line, "%*s  %ld %ld %lf",  &_EnergyTuneShift_npoint,
-        		&_EnergyTuneShift_nturn, &_EnergyTuneShift_deltamax);
-      
-        	//EnergyTuneShiftFlag = true;
-        
+        sscanf(line, "%*s  %ld %ld %lf",  
+	       &(UserCommandFlag[CommNo]._EnergyTuneShift_npoint),
+               &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), 
+	       &(UserCommandFlag[CommNo]._EnergyTuneShift_deltamax));
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);    
       }
       else if (strcmp("ErrorCouplingFlag", name) == 0){
-      strcpy(UserCommand[CommandNo-1],name);
-        sscanf(line, "%*s  %ld %lf",&err_seed,&err_rms);
+        sscanf(line, "%*s  %ld %lf",&(UserCommandFlag[CommNo].err_seed),&(UserCommandFlag[CommNo].err_rms));
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);
       
-        //  ErrorCouplingFlag = true;
-        
       }
         else if (strcmp("ErrorCoupling2Flag", name) == 0){
-        sscanf(line, "%*s  %ld %lf",&err_seed,&err_rms);
-        strcpy(UserCommand[CommandNo-1],name);
-        //  ErrorCoupling2Flag = true;
+        sscanf(line, "%*s  %ld %lf",&(UserCommandFlag[CommNo].err_seed),&(UserCommandFlag[CommNo].err_rms));
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);
         
       }
       else if (strcmp("CouplingFlag", name) == 0){
-        strcpy(UserCommand[CommandNo-1],name);
-         // CouplingFlag = true;
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);
         
       }
        else if (strcmp("MomentumAccFlag", name) == 0){
-        sscanf(line, "%*s  %s %ld %ld %lf %lf %ld %lf %lf %ld",TrackDim,
-        		&_MomentumAccFlag_istart, &_MomentumAccFlag_istop,
-        		&_MomentumAccFlag_deltaminp, &_MomentumAccFlag_deltamaxp, &_MomentumAccFlag_nstepp,
-        		&_MomentumAccFlag_deltaminn, &_MomentumAccFlag_deltamaxn, &_MomentumAccFlag_nstepn
+        sscanf(line, "%*s  %s %ld %ld %lf %lf %ld %lf %lf %ld",
+	              UserCommandFlag[CommNo].TrackDim,
+		      &(UserCommandFlag[CommNo]._MomentumAccFlag_istart),
+		      &(UserCommandFlag[CommNo]._MomentumAccFlag_istop),
+        	      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminp),
+		      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxp),
+		      &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepp),
+        	      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminn),
+		      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxn),
+		      &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepn)
         		);
-        strcpy(UserCommand[CommandNo-1],name);
-        //  MomentumAccFlag = true;
-         
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);   
       }
         else if (strcmp("ReadMultipoleFlag", name) == 0){
-            strcpy(UserCommand[CommandNo-1],name);
-        //  ReadMultipoleFlag = true;
+            strcpy(UserCommandFlag[CommNo].CommandStr,name);
         
       }
-       else if (strcmp("MultipoleFlag", name) == 0){
-         strcpy(UserCommand[CommandNo-1],name);
-          //MultipoleFlag = true;
-      
-      }
-      else if (strcmp("ThinsextFlag", name) == 0){
-           strcpy(UserCommand[CommandNo-1],name);
-        //  ThinsextFlag = true;
-       
-      }
+    
       else if (strcmp("fic_hcorr", name) == 0){// read of h-correctors for multipoles
         sscanf(line, "%*s %s", str);
         sprintf(fic_hcorr,"%s%s", in_dir, str);
@@ -379,85 +285,82 @@ void read_script(const char *param_file_name, bool rd_lat)
       else if (strcmp("fic_skew", name) == 0){// read of skew quads for multipoles
         sscanf(line, "%*s %s", str);
         sprintf(fic_skew,"%s%s", in_dir, str);
-      } // generic one, fit for 1 family of quadrupoles
+      } 
+//       generic one, fit for 1 family of quadrupoles
       else if (strcmp("FitTuneFlag", name) == 0){
-        sscanf(line, "%*s %s %s %lf %lf",qf,qd,&targetnux,&targetnuz);
-        strcpy(UserCommand[CommandNo-1],name);
-        //  FitTuneFlag = true;
-        
-      }// fit for 2 families,specific for the soleil lattice in which the quadrupole is cut into 2 parts
+        sscanf(line, "%*s %s %s %lf %lf",UserCommandFlag[CommNo].qf,UserCommandFlag[CommNo].qd,
+	&(UserCommandFlag[CommNo].targetnux),&(UserCommandFlag[CommNo].targetnuz));
+       strcpy(UserCommandFlag[CommNo].CommandStr,name);
+      }
+      // fit for 2 families,specific for the soleil lattice in which the quadrupole is cut into 2 parts
        else if (strcmp("FitTune4Flag", name) == 0){
-        sscanf(line, "%*s %s %s %s %s %lf %lf",qf1,qf2,qd1,qd2,&targetnux,&targetnuz);
-        strcpy(UserCommand[CommandNo-1],name);
-        //  FitTune4Flag = true;
-        
+        sscanf(line, "%*s %s %s %s %s %lf %lf",UserCommandFlag[CommNo].qf1,UserCommandFlag[CommNo].qf2,
+	UserCommandFlag[CommNo].qd1,UserCommandFlag[CommNo].qd2,
+	&(UserCommandFlag[CommNo].targetnux),&(UserCommandFlag[CommNo].targetnuz));
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
        else if (strcmp("FitChromFlag", name) == 0){
-        sscanf(line, "%*s  %s %s %lf %lf",sxm1,sxm2,&targetksix,&targetksiz);
-        strcpy(UserCommand[CommandNo-1],name);
-        //  FitChromFlag = true;
-        
-      }
-      else if (strcmp("GirderErrorFlag", name) == 0){
-       strcpy(UserCommand[CommandNo-1],name);
-         // GirderErrorFlag = true;
-       
-      }
-       else if (strcmp("SigmaFlag", name) == 0){
-        strcpy(UserCommand[CommandNo-1],name);
-         // SigmaFlag = true;
-        
-      }
-       else if (strcmp("PX2Flag", name) == 0){
-       strcpy(UserCommand[CommandNo-1],name);
-         // PX2Flag = true;
-        
+        sscanf(line, "%*s  %s %s %lf %lf",UserCommandFlag[CommNo].sxm1,UserCommandFlag[CommNo].sxm2,
+	&(UserCommandFlag[CommNo].targetksix),&(UserCommandFlag[CommNo].targetksiz));
+        strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
+//       else if (strcmp("GirderErrorFlag", name) == 0){
+//        strcpy(UserCommand[CommandNo-1],name);
+//          // GirderErrorFlag = true;
+//        
+//       }
+//        else if (strcmp("SigmaFlag", name) == 0){
+//         strcpy(UserCommand[CommandNo-1],name);
+//          // SigmaFlag = true;
+//         
+//       }
+//        else if (strcmp("PX2Flag", name) == 0){
+//        strcpy(UserCommand[CommandNo-1],name);
+//          // PX2Flag = true;
+//         
+//       }
        else if (strcmp("InducedAmplitudeFlag", name) == 0){
-       strcpy(UserCommand[CommandNo-1],name);
-        //  InducedAmplitudeFlag = true;
-        
-      }
-       else if (strcmp("CodeComparaisonFlag", name) == 0){
-       strcpy(UserCommand[CommandNo-1],name);
-        //  CodeComparaisonFlag = true;
-       
+        strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       }
+//        else if (strcmp("CodeComparaisonFlag", name) == 0){
+//        strcpy(UserCommand[CommandNo-1],name);
+//         //  CodeComparaisonFlag = true;
+//        
+//       }
        else if (strcmp("EtaFlag", name) == 0){
-       strcpy(UserCommand[CommandNo-1],name);
-         // EtaFlag = true;
-        
+         strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       }
        else if (strcmp("PhaseSpaceFlag", name) == 0) {
-        sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld %s", _Phase_Dim,
-            &_Phase_X, &_Phase_Px, &_Phase_Y, &_Phase_Py, &_Phase_delta,
-            &_Phase_ctau, &_Phase_nturn, str);
+        sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld %s", 
+	             UserCommandFlag[CommNo]._Phase_Dim,
+                     &(UserCommandFlag[CommNo]._Phase_X), 
+		     &(UserCommandFlag[CommNo]._Phase_Px), 
+		     &(UserCommandFlag[CommNo]._Phase_Y), 
+		     &(UserCommandFlag[CommNo]._Phase_Py), 
+		     &(UserCommandFlag[CommNo]._Phase_delta),
+                     &(UserCommandFlag[CommNo]._Phase_ctau), 
+		     &(UserCommandFlag[CommNo]._Phase_nturn), 
+		     str);
         if (strcmp(str, "true") == 0)
-          _Phase_Damping = true;
+          UserCommandFlag[CommNo]._Phase_Damping = true;
         else if (strcmp(str, "false") == 0)
-          _Phase_Damping = false;
+          UserCommandFlag[CommNo]._Phase_Damping = false;
         else {
           printf("set boolean flag true or false for PhaseSpaceFlag \n");
           exit_(1);
         }
-	strcpy(UserCommand[CommandNo-1],name);
-        //PhaseSpaceFlag = true;
+	 strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       }
       //calculate Touschek lifetime
        else if (strcmp("TouschekFlag", name) == 0){
-        //strcpy(UserCommand[CommandNo-1],name);
-          TouschekFlag = true;
+	  strcpy(UserCommandFlag[CommNo].CommandStr,name); 
        
-      }
+      }//intra beam scattering
        else if (strcmp("IBSFlag", name) == 0){
-       //strcpy(UserCommand[CommandNo-1],name);
-          IBSFlag = true;
-        
-      }
+        strcpy(UserCommandFlag[CommNo].CommandStr,name); 
+      }//momentum acceptance is got by tracking, then cal touschek lifetime
        else if (strcmp("TousTrackFlag", name) == 0){
-      // strcpy(UserCommand[CommandNo-1],name);
-          TousTrackFlag = true;
-        
+        strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       }
       
       else {