diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 3e97134287d5582761345c8d728d5770d69bfcd8..956e15303e0f295e95c9552b673d8db60d07eadc 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -466,7 +466,9 @@ int main(int argc, char *argv[]) {
   /******************************************************************************************/
   // computes TuneShift with amplitudes
   else if(strcmp(CommandStr,"AmplitudeTuneShiftFlag") == 0) {
-      TunesShiftWithAmplitude(UserCommandFlag[i]._AmplitudeTuneShift_nxpoint,
+      TunesShiftWithAmplitude(UserCommandFlag[i]._AmplitudeTuneShift_nudx_file,
+		              UserCommandFlag[i]._AmplitudeTuneShift_nudz_file,
+                              UserCommandFlag[i]._AmplitudeTuneShift_nxpoint,
                               UserCommandFlag[i]._AmplitudeTuneShift_nypoint,
                               UserCommandFlag[i]._AmplitudeTuneShift_nturn, 
 			      UserCommandFlag[i]._AmplitudeTuneShift_xmax,
@@ -475,7 +477,8 @@ int main(int argc, char *argv[]) {
     } 
   // compute tuneshift with energy
  else if(strcmp(CommandStr,"EnergyTuneShiftFlag") == 0) {
-     TunesShiftWithEnergy(UserCommandFlag[i]._EnergyTuneShift_npoint, 
+     TunesShiftWithEnergy(UserCommandFlag[i]._EnergyTuneShift_nudp_file,
+                          UserCommandFlag[i]._EnergyTuneShift_npoint, 
                           UserCommandFlag[i]._EnergyTuneShift_nturn,
                           UserCommandFlag[i]._EnergyTuneShift_deltamax);
    }
@@ -483,7 +486,8 @@ int main(int argc, char *argv[]) {
   // Computes FMA
   else if(strcmp(CommandStr,"FmapFlag") == 0) {
     printf("\n begin Fmap calculation for on momentum particles: \n");
-    fmap(UserCommandFlag[i]._FmapFlag_nxpoint, 
+    fmap(UserCommandFlag[i]._FmapFlag_fmap_file,
+         UserCommandFlag[i]._FmapFlag_nxpoint, 
          UserCommandFlag[i]._FmapFlag_nypoint, 
 	 UserCommandFlag[i]._FmapFlag_nturn, 
 	 UserCommandFlag[i]._FmapFlag_xmax,
@@ -495,7 +499,8 @@ int main(int argc, char *argv[]) {
   // Compute FMA dp
   else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
     printf("\n begin Fmap calculation for off momentum particles: \n");
-    fmapdp(UserCommandFlag[i]._FmapdpFlag_nxpoint, 
+    fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
+           UserCommandFlag[i]._FmapdpFlag_nxpoint, 
            UserCommandFlag[i]._FmapdpFlag_nepoint, 
 	   UserCommandFlag[i]._FmapdpFlag_nturn,
            UserCommandFlag[i]._FmapdpFlag_xmax, 
@@ -512,15 +517,15 @@ int main(int argc, char *argv[]) {
   // MOMENTUM ACCEPTANCE
   else if(strcmp(CommandStr,"MomentumAccFlag") == 0) {
     bool cavityflag, radiationflag;
+     
     /* record the initial values*/
     cavityflag = globval.Cavity_on;
-    radiationflag = globval.radiation;
-
-    /* set the dimension for the momentum tracking*/
-    if (strncmp("6D", UserCommandFlag[i].TrackDim, 2) == 0) {
+    radiationflag = globval.radiation; 
+    
+    if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"6D") == 0) {
       globval.Cavity_on = true;
       globval.radiation = true;
-    } else if (strncmp("4D", UserCommandFlag[i].TrackDim, 2) == 0) {
+    }else if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"4D") == 0) {
       globval.Cavity_on = false;
       globval.radiation = false;
     } else {
@@ -530,15 +535,16 @@ int main(int argc, char *argv[]) {
 
     /* calculate momentum acceptance*/
     printf("\n Calculate momentum acceptance: \n");
-    MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_istart, 
+    MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
+                       UserCommandFlag[i]._MomentumAccFlag_istart, 
                        UserCommandFlag[i]._MomentumAccFlag_istop,
                        UserCommandFlag[i]._MomentumAccFlag_deltaminp, 
-		               UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
+		       UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
                        UserCommandFlag[i]._MomentumAccFlag_nstepp, 
-		               UserCommandFlag[i]._MomentumAccFlag_deltaminn,
+		       UserCommandFlag[i]._MomentumAccFlag_deltaminn,
                        UserCommandFlag[i]._MomentumAccFlag_deltamaxn, 
-		               UserCommandFlag[i]._MomentumAccFlag_nstepn,
-		               UserCommandFlag[i]._MomentumAccFlag_nturn);
+		       UserCommandFlag[i]._MomentumAccFlag_nstepn,
+		       UserCommandFlag[i]._MomentumAccFlag_nturn);
 
     /* restore the initial values*/
     globval.Cavity_on = cavityflag;
@@ -577,9 +583,9 @@ int main(int argc, char *argv[]) {
     radiationflag = globval.radiation;
 
     /* set the dimension for the momentum tracking*/
-    if (strncmp("6D", UserCommandFlag[i]._Phase_Dim, 2) == 0) {
+    if (strcmp("6D", UserCommandFlag[i]._Phase_Dim) == 0) {
       globval.Cavity_on = true;
-    } else if (strncmp("4D", UserCommandFlag[i]._Phase_Dim, 2) == 0) {
+    } else if (strcmp("4D", UserCommandFlag[i]._Phase_Dim) == 0) {
       globval.Cavity_on = false;
     } else {
       printf("MomentumAccFlag: Error!!! _Phase_Dim must be '4D' or '6D'\n");
@@ -593,7 +599,8 @@ int main(int argc, char *argv[]) {
     }
 
     start = stampstart();
-    Phase(UserCommandFlag[i]._Phase_X, 
+    Phase(UserCommandFlag[i]._Phase_phase_file,
+          UserCommandFlag[i]._Phase_X, 
           UserCommandFlag[i]._Phase_Px, 
 	  UserCommandFlag[i]._Phase_Y, 
 	  UserCommandFlag[i]._Phase_Py, 
diff --git a/tracy/tracy/inc/read_script.h b/tracy/tracy/inc/read_script.h
index 8a2050f17ce0a594cb33cfe98a7efc5f4e883d21..202a66befae38b8cbe68a2a10bef9daecd188e37 100644
--- a/tracy/tracy/inc/read_script.h
+++ b/tracy/tracy/inc/read_script.h
@@ -32,27 +32,33 @@
    long err_seed; 
    double err_rms;
   //AmplitudeTuneShiftFlag;
+   char _AmplitudeTuneShift_nudx_file[max_str];
+   char _AmplitudeTuneShift_nudz_file[max_str];
    long _AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint;
    long _AmplitudeTuneShift_nturn;
    double _AmplitudeTuneShift_xmax, _AmplitudeTuneShift_ymax, _AmplitudeTuneShift_delta;
  
  //EnergyTuneShiftFlag;
+   char _EnergyTuneShift_nudp_file[max_str];
    long _EnergyTuneShift_npoint, _EnergyTuneShift_nturn;
    double _EnergyTuneShift_deltamax;
  
  //extern bool FmapFlag;
+   char _FmapFlag_fmap_file[max_str];
    long _FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn;
    double _FmapFlag_xmax, _FmapFlag_ymax, _FmapFlag_delta;
    bool _FmapFlag_diffusion; 
  
  //extern bool FmapdpFlag;
+   char _FmapdpFlag_fmapdp_file[max_str];
    long _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn;
    double _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z;
    bool _FmapdpFlag_diffusion;	
 
    
   //MomentumAccFlag;
-   char TrackDim[3];
+   char _MomentumAccFlag_momacc_file[max_str];  //  file to save tracked momentum accpetance
+   char _MomentumAccFlag_TrackDim[3];
    long  _MomentumAccFlag_istart, _MomentumAccFlag_istop,
          _MomentumAccFlag_nstepn, _MomentumAccFlag_nstepp,
          _MomentumAccFlag_nturn;
@@ -60,35 +66,42 @@
    double _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp;
  
   // /* Phase space */
+   //char *_Phase_phase_file;
+   char _Phase_phase_file[max_str];
+   char _Phase_Dim[3];
    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
+  //set default values
   UserCommand(void)   //constructor
  {  
   
   // /* fmap for on momentum particle*/
+   strcpy(_FmapFlag_fmap_file,"fmap.out");
    _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*/
+strcpy(_FmapdpFlag_fmapdp_file,"fmapdp.out");
  _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*/
+ strcpy(_AmplitudeTuneShift_nudx_file,"nudx.out");
+ strcpy(_AmplitudeTuneShift_nudz_file,"nudz.out");
  _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*/
+strcpy(_EnergyTuneShift_nudp_file,"nudp.out");
  _EnergyTuneShift_npoint=31L;  _EnergyTuneShift_nturn=516L;
  _EnergyTuneShift_deltamax=0.06;
 
@@ -96,7 +109,8 @@
  err_seed=0L;  err_rms=0.0;  
 
 /* momentum acceptance */
- TrackDim[3] = '6D';
+ strcpy(_MomentumAccFlag_momacc_file,"momentumacceptance.out");//default file name
+ strcpy(_MomentumAccFlag_TrackDim,"6D"); //default track dimension
  _MomentumAccFlag_istart=1L, _MomentumAccFlag_istop=108L,
  _MomentumAccFlag_nstepn=100L, _MomentumAccFlag_nstepp=100L;
  _MomentumAccFlag_deltaminn=-0.01, _MomentumAccFlag_deltamaxn=-0.05;
@@ -104,10 +118,11 @@
  _MomentumAccFlag_nturn = 1026;
 
 /* Phase space */
+ strcpy(_Phase_phase_file,"phase.out");  //default phase file
+ strcpy( _Phase_Dim,"4D"); //default phase dimension
  _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*/
diff --git a/tracy/tracy/inc/soleillib.h b/tracy/tracy/inc/soleillib.h
index d721ef54a0cbc2cd5f9bc28fae3468886ab6a0fc..79629a6c7b0ec67b874745ba777ac0a47737a937 100644
--- a/tracy/tracy/inc/soleillib.h
+++ b/tracy/tracy/inc/soleillib.h
@@ -21,7 +21,7 @@ void set_vectorcod(Vector codvector[], double dP);
 void SetDecapole(void);
 
 /* Tracking */
-void Phase(double x,double xp,double y, double yp,double energy, double ctau, long Nbtour);
+void Phase(const char *phasefile,double x,double xp,double y, double yp,double energy, double ctau, long Nbtour);
 void Phase2(long pos, double x,double xp,double y, double yp,double energy, double ctau,
             long Nbtour);
 void PhasePoly(long pos, double x0,double px0, double z0, double pz0, double delta0,
@@ -31,10 +31,11 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0, do
                           double end, long Nb, long Nbtour, int num);
 void PhasePortrait2(long pos,double x0,double px0,double z0, double pz0, double delta0, double ctau,
                           double end, long Nb, long Nbtour, int num);
-void Multipole_thicksext(char const *fic_hcorr, char const *fic_vcorr, char const *fic_skew);
-void Multipole_thinsext(char const *fic_hcorr, char const *fic_vcorr, char const *fic_skew);
-void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max, long nstepp,
-                        double em_min, double em_max, long nstepm, long nturn);
+void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const char *fic_skew);
+void Multipole_thinsext(const char  *fic_hcorr, const char *fic_vcorr, const char *fic_skew);
+void MomentumAcceptance(char *MomAccFile,long deb, long fin, double ep_min,  
+                        double ep_max, long nstepp, double em_min, double em_max, 
+			long nstepm, long nturn);
 
 void Trac_Tab(double x, double px, double y, double py, double dp,
             long nmax, long pos, long *lastn, long *lastpos, FILE *outf1, double Tx[][NTURN]);
@@ -45,18 +46,19 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
                double energy, bool diffusion);
                             
 /* Frequency map analysis */
-void TunesShiftWithEnergy(long Nb, long Nbtour, double emax);
+void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nbx, 
+                             long Nbz, long Nbtour, double xmax, double zmax, double energy);
+void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax);
 //void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
 //                 double energy, bool diffusion, bool matlab);
 //void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
 //              double z, bool diffusion, bool matlab);
-void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-                 double energy, bool diffusion);
-void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
+void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+          double energy, bool diffusion);	  
+void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
               double z, bool diffusion);
 void Nu_Naff(void);
-void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, double ymax,
-                 double energy);
+
 
 /* Vacuum chamber */
 void ReadCh(const char *AperFile);
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index f54966dab3e18a82720134e40e8c0f354b9e6bb4..d30d5ccd341ccb73b91935bf755aca37caf8a19c 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -36,7 +36,8 @@ char twiss_file[max_str];
 char cod_file[max_str];
 // girder file
 char girder_file[max_str];
- //chamber file
+
+
 
  // multipole files
 char multipole_file[max_str];
@@ -64,7 +65,7 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
  {
  
   
-  char    str[max_str]="voidstring", dummy[max_str]="voidstring";
+  char    str[max_str]="voidstring", dummy[max_str]="voidstring", nextpara[max_str]="voidpara";
   char    in[max_str];  //temporary line with preceding white space
   char    *line; //line to store the command without preceding white space
   char    name[max_str]="voidname";
@@ -115,8 +116,9 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
     if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0 &&
         strcmp(line,"\r") != 0 && strcmp(line,"\r\n") != 0) {
       // get initial command token
-      sscanf(line, "%s", name);
-      
+      sscanf(line, "%s", name); 
+      //initialize the nextpara for each new command line
+      strcpy(nextpara,"voidpara");
       
       
       //find the sequence of the bool flag in user input script 
@@ -130,7 +132,7 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
       if(strcmp(EndName,"Flag")==0)
            CommNo++;
         
-      
+         
       /*read file names................. */
       /* set file directory*/
       if (strcmp("in_dir", name) == 0)
@@ -244,67 +246,80 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
       // FMA
       else if (strcmp("FmapFlag", name) == 0){        
     	  strcpy(dummy, "");
-          sscanf(line, "%*s %ld %ld %ld %lf %lf %lf %s", 
-          		&(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(UserCommandFlag[CommNo].CommandStr,name);
 	  
+	  sscanf(line, "%*s %s",nextpara);
+	 
+	  if(strcmp(nextpara,"voidpara")!=0)
+            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s", 
+              		  UserCommandFlag[CommNo]._FmapFlag_fmap_file,
+			  &(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);
+         
         if(strcmp(dummy, "true") == 0)
           UserCommandFlag[CommNo]._FmapFlag_diffusion = true;
         else if(strcmp(dummy, "false") == 0)
           UserCommandFlag[CommNo]._FmapFlag_diffusion = false;
-        else {
-          printf("set boolean flag true or false for FMA Diffusion (value is %s) \n", dummy);
-          exit_(1);
-        }
+       
+       // FmapFlag = true;
+         strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
       // FMA dp
       else if (strcmp("FmapdpFlag", name) == 0){         
     	  strcpy(dummy, "");
-          sscanf(line, "%*s %ld %ld %ld %lf %lf %lf %s", 
-          		&(UserCommandFlag[CommNo]._FmapdpFlag_nxpoint),
+	  sscanf(line, "%*s %s",nextpara);
+	 
+	  if(strcmp(nextpara,"voidpara")!=0)
+            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s", 
+          		UserCommandFlag[CommNo]._FmapdpFlag_fmapdp_file,
+			&(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(UserCommandFlag[CommNo].CommandStr,name);
 	 
         if(strcmp(dummy, "true") == 0)
           UserCommandFlag[CommNo]._FmapdpFlag_diffusion = true;
         else if(strcmp(dummy, "false") == 0)
           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){   
-        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));
+
+	  //  FmapdpFlag = true;
+         strcpy(UserCommandFlag[CommNo].CommandStr,name);
+      }
+      else if (strcmp("AmplitudeTuneShiftFlag", name) == 0){
+        sscanf(line, "%*s %s",nextpara);
+	 
+	if(strcmp(nextpara,"voidpara")!=0)   
+          sscanf(line, "%*s %s %s %ld %ld %ld %lf %lf %lf", 
+	          UserCommandFlag[CommNo]._AmplitudeTuneShift_nudx_file,
+		  UserCommandFlag[CommNo]._AmplitudeTuneShift_nudz_file,
+	         &(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){
-        sscanf(line, "%*s  %ld %ld %lf",  
-	       &(UserCommandFlag[CommNo]._EnergyTuneShift_npoint),
-               &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), 
-	       &(UserCommandFlag[CommNo]._EnergyTuneShift_deltamax));
-        strcpy(UserCommandFlag[CommNo].CommandStr,name);    
+        sscanf(line, "%*s %s",nextpara);
+	 
+	 if(strcmp(nextpara,"voidpara")!=0)
+	   sscanf(line, "%*s %s %ld %ld %lf",  
+	          UserCommandFlag[CommNo]._EnergyTuneShift_nudp_file,
+	         &(UserCommandFlag[CommNo]._EnergyTuneShift_npoint),
+                 &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), 
+	         &(UserCommandFlag[CommNo]._EnergyTuneShift_deltamax));
+        
+	strcpy(UserCommandFlag[CommNo].CommandStr,name);
+	    
       }
       else if (strcmp("ErrorCouplingFlag", name) == 0){
         sscanf(line, "%*s  %ld %lf",&(UserCommandFlag[CommNo].err_seed),&(UserCommandFlag[CommNo].err_rms));
@@ -319,10 +334,14 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
       else if (strcmp("CouplingFlag", name) == 0){
         strcpy(UserCommandFlag[CommNo].CommandStr,name);
         
-      }
+      }//momentum compact factor
        else if (strcmp("MomentumAccFlag", name) == 0){
-        sscanf(line, "%*s  %s %ld %ld %lf %lf %ld %lf %lf %ld %ld",
-	              UserCommandFlag[CommNo].TrackDim,
+         sscanf(line, "%*s %s",nextpara);
+	 
+	 if(strcmp(nextpara,"voidpara")!=0){    
+	   sscanf(line, "%*s  %s %s %ld %ld %lf %lf %ld %lf %lf %ld %ld",
+		      UserCommandFlag[CommNo]._MomentumAccFlag_momacc_file,
+		      UserCommandFlag[CommNo]._MomentumAccFlag_TrackDim,
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_istart),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_istop),
         	      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminp),
@@ -331,9 +350,9 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
         	      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminn),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxn),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepn),
-              &(UserCommandFlag[CommNo]._MomentumAccFlag_nturn)
-        		);
-        strcpy(UserCommandFlag[CommNo].CommandStr,name);   
+                      &(UserCommandFlag[CommNo]._MomentumAccFlag_nturn));
+	 } 
+	 strcpy(UserCommandFlag[CommNo].CommandStr,name);   
       }
    
 //       generic one, fit for 1 family of quadrupoles
@@ -379,26 +398,29 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
 //       }
        else if (strcmp("EtaFlag", name) == 0){
          strcpy(UserCommandFlag[CommNo].CommandStr,name); 
-      }
+      }//phase space
        else if (strcmp("PhaseSpaceFlag", name) == 0) {
-        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)
-          UserCommandFlag[CommNo]._Phase_Damping = true;
-        else if (strcmp(str, "false") == 0)
-          UserCommandFlag[CommNo]._Phase_Damping = false;
-        else {
-          printf("set boolean flag true or false for PhaseSpaceFlag \n");
-          exit_(1);
-        }
+         sscanf(line, "%*s %s",nextpara);
+	 //if the next para is not empty, then use the user set value, otherwise use default value
+	 if(strcmp(nextpara,"voidpara") != 0){
+	   sscanf(line, "%*s %s %s %lf %lf %lf %lf %lf %lf %ld %s", 
+	               UserCommandFlag[CommNo]._Phase_phase_file,
+	               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);
+	 }
+	 //radiation damping
+          if (strcmp(str, "true") == 0)
+            UserCommandFlag[CommNo]._Phase_Damping = true;
+          else if (strcmp(str, "false") == 0)
+            UserCommandFlag[CommNo]._Phase_Damping = false;
+	
 	 strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       }
       //calculate Touschek lifetime
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index 5579dc553d735926e962995cbef5210fd2c9cb3c..218e7148e177c47c52e33db5de6b4b30b30d3e80 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -685,15 +685,13 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
    Comments:
        16/01/03 add test for non zero frequency
                 add variation around the closed orbit
-
+       19/07/11 add feature to save tune shift with amplitude in the user defined file
 ****************************************************************************/
 #define nterm  4
-void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-               double energy)
+void TunesShiftWithAmplitude(const char *NudxFile, const char *NudzFile, long Nbx, 
+                             long Nbz, long Nbtour, double xmax, double zmax, double energy)
 {
   FILE * outf;
-  const char ficx[] = "nudx.out";
-  const char ficz[] = "nudz.out";
   int i = 0;
   double Tab[6][NTURN], fx[nterm], fz[nterm];
   double x = 0.0 , xp = 0.0 , z = 0.0 , zp = 0.0;
@@ -716,17 +714,17 @@ void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, doubl
   if (fabs(xmax) > 0.0){
     
     /* Opening file */
-    if ((outf = fopen(ficx, "w")) == NULL) {
-      fprintf(stdout, "NuDx: error while opening file %s\n", ficx);
+    if ((outf = fopen(NudxFile, "w")) == NULL) {
+      fprintf(stdout, "TunesShiftWithAmplitude: error while opening file %s\n", NudxFile);
       exit_(1);
     }
 
-    fprintf(outf,"# Tracy III -- %s -- %s \n", ficx, asctime2(newtime));
+    fprintf(outf,"# Tracy III -- %s -- %s \n", NudxFile, asctime2(newtime));
     fprintf(outf,"# nu = f(x) \n");
     fprintf(outf,"#    x[m]          z[m]           fx            fz \n");
 
     if ((Nbx <= 1) || (Nbz <= 1))
-      fprintf(stdout,"NuDx: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
+      fprintf(stdout,"TunesShiftWithAmplitude: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
 
     xstep = xmax/Nbx*2.0;
     x0 = 1e-6 - xmax;
@@ -761,12 +759,12 @@ void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, doubl
   if (fabs(zmax) > 0.0)
   {
     /* Opening file */
-    if ((outf = fopen(ficz, "w")) == NULL) {
-      fprintf(stdout, "NuDx: error while opening file %s\n", ficz);
+    if ((outf = fopen(NudzFile, "w")) == NULL) {
+      fprintf(stdout, "TunesShiftWithAmplitude: error while opening file %s\n", NudzFile);
       exit_(1);
     }
     
-    fprintf(outf,"# tracy III -- %s -- %s \n", ficz, asctime2(newtime));
+    fprintf(outf,"# tracy III -- %s -- %s \n", NudzFile, asctime2(newtime));
     fprintf(outf,"# nu = f(z) \n");
     fprintf(outf,"#    x[mm]         z[mm]          fx            fz \n");
 
@@ -813,8 +811,8 @@ double get_D(const double df_x, const double df_y)
 
 
 /****************************************************************************/
-/* void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-   double energy, bool diffusion, bool matlab)
+/* void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+          double energy, bool diffusion)
 
    Purpose:
        Compute a frequency map of Nbx x Nbz points
@@ -830,6 +828,7 @@ double get_D(const double df_x, const double df_y)
        Results in fmap.out
 
    Input:
+       FmapFile file to save calculated frequency map analysis
        Nbx    horizontal step number
        Nby    vertical step number
        xmax   horizontal maximum amplitude
@@ -854,14 +853,14 @@ double get_D(const double df_x, const double df_y)
    Comments:
        15/10/03 run for the diffusion: nasty patch for retrieving the closed orbit
        16/02/03 patch removed
-       
+       19/07/11 add interface of file defined by user which is used to save calculated
+                frequency map analysis
 ****************************************************************************/
 #define NTERM2  10
-void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
           double energy, bool diffusion)	  
 {
  FILE * outf;
- const char fic[] = "fmap.out";
  long i = 0L, j = 0L;
  double Tab[DIM][NTURN], Tab0[DIM][NTURN];
  double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
@@ -880,15 +879,15 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
  time(&aclock);                 /* Get time in seconds */
  newtime = localtime(&aclock);  /* Convert time to struct */
 
- if (trace) printf("Entering fmap ... results in %s\n\n",fic);
+ if (trace) printf("Entering fmap ... results in %s\n\n",FmapFile);
 
  /* Opening file */
- if ((outf = fopen(fic, "w")) == NULL) {
-   fprintf(stdout, "fmap: error while opening file %s\n", fic);
+ if ((outf = fopen(FmapFile, "w")) == NULL) {
+   fprintf(stdout, "fmap: error while opening file %s\n", FmapFile);
    exit_(1);
  }
 
- fprintf(outf,"# TRACY III SYNCHROTRON SOLEIL-- %s -- %s \n", fic, asctime2(newtime));
+ fprintf(outf,"# TRACY III -- %s -- %s \n", FmapFile, asctime2(newtime));
  fprintf(outf,"# nu = f(x) \n");
 // fprintf(outf,"#    x[mm]          z[mm]           fx             fz"
 //	 "            dfx            dfz      D=log_10(sqrt(df_x^2+df_y^2))\n");
@@ -971,8 +970,8 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
 #undef NTERM2
 
 /****************************************************************************/
-/* void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
-   double z, bool *status)
+/* void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
+              double z, bool diffusion)
 
    Purpose:
        Compute a frequency map of Nbx x Nbz points
@@ -987,6 +986,7 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
        Results in fmapdp.out
 
    Input:
+       FmapdpFile       file to save calculated frequency map analysis
        Nbx              horizontal step number
        Nbe              energy step number
        Nbtour           number of turns for tracking
@@ -1012,14 +1012,13 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
        15/10/03 run for the diffusion: nasty patch for retrieving the closed orbit
        23/10/04 for 6D turn off diffusion automatically and horizontal amplitude
        is negative for negative enrgy offset since this is true for the cod
-
+       19/07/11  add features to save calculated fmapdp in the user defined file
 ****************************************************************************/
 #define NTERM2  10
-void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
+void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
               double z, bool diffusion)
 {
  FILE * outf;
- const char fic[] = "fmapdp.out";
  long i = 0L, j = 0L;
  double Tab[DIM][NTURN], Tab0[DIM][NTURN];
  double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
@@ -1040,21 +1039,21 @@ void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
 
  if (diffusion && globval.Cavity_on == false) nturn = 2*Nbtour;
 
- if (trace) printf("Entering fmap ... results in %s\n\n",fic);
+ if (trace) printf("Entering fmap ... results in %s\n\n",FmapdpFile);
 
  /* Opening file */
- if ((outf = fopen(fic, "w")) == NULL) {
-   fprintf(stdout, "fmap: error while opening file %s\n", fic);
+ if ((outf = fopen(FmapdpFile, "w")) == NULL) {
+   fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpFile);
    exit_(1);
  }
 
- fprintf(outf,"# TRACY III SYNCHROTRON SOLEIL-- %s -- %s \n", fic, asctime2(newtime));
+ fprintf(outf,"# TRACY III -- %s -- %s \n", FmapdpFile, asctime2(newtime));
  fprintf(outf,"# nu = f(x) \n");
 // fprintf(outf,"#    dp[%%]         x[mm]          fx            fz           dfx           dfz\n");
 fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx           dfz\n");
  
  if ((Nbx <= 1) || (Nbe <= 1))
-   fprintf(stdout,"fmap: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
+   fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
 
  xp = xp0;
  zp = zp0;
@@ -1149,10 +1148,10 @@ fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx
 
 ****************************************************************************/
 #define NTERM  4
-void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
+void TunesShiftWithEnergy(const char *NudpFile,long Nb, long Nbtour, double emax)
 {
   FILE * outf;
-  const char fic[] = "nudp.out";
+ 
   long i = 0L;
 //  long lastpos = 0L;
   double Tab[DIM][NTURN];
@@ -1170,12 +1169,12 @@ void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
   if (!trace) printf("\n Entering TunesShiftWithEnergy ...\n\n");
 
   /* Opening file */
-  if ((outf = fopen(fic, "w")) == NULL) {
-    fprintf(stdout, "NuDp: error while opening file %s\n", fic);
+  if ((outf = fopen(NudpFile, "w")) == NULL) {
+    fprintf(stdout, "NuDp: error while opening file %s\n", NudpFile);
     exit_(1);
   }
 
-  fprintf(outf,"# TRACY III -- %s -- %s \n", fic, asctime2(newtime));
+  fprintf(outf,"# TRACY III -- %s -- %s \n", NudpFile, asctime2(newtime));
   fprintf(outf,"#    dP/P           fx            fz          xcod         pxcod          zcod         pzcod\n");
   if (trace) fprintf(stdout,"#    dP/P           fx            fz          xcod         pxcod          zcod         pzcod\n");
 
@@ -1250,18 +1249,20 @@ void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
        1 December 2010, Call to a Tracking round around the 6D and not 4D closed orbit
 
 ****************************************************************************/
-void Phase(double x,double xp,double y, double yp,double energy, double ctau, long Nbtour)
+void Phase(const char *phasefile, double x,double xp,double y, double yp,double energy, double ctau, long Nbtour)
 {
   double Tab[6][NTURN];
   FILE *outf;
-  const char fic[] = "phase.out";
+  const char *fic = phasefile; 
   int i;
   bool status;
   struct tm *newtime;
 
   /* Get time and date */
   newtime = GetTime();
-
+   
+ 
+  
   if (Nbtour > NTURN) {
     fprintf(stdout, "Phase: error Nbtour=%ld > NTURN=%d\n",Nbtour,NTURN);
     exit_(1);
@@ -1620,7 +1621,8 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn
 
 
 /****************************************************************************/
-/* void Multipole_thicksext(void)
+/* void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, 
+                            const char *fic_skew)
 
    Purpose:
        Set multipole in dipoles, quadrupoles, thick sextupoles, skew quadrupole,
@@ -1650,7 +1652,7 @@ void Enveloppe(double x, double px, double y, double py, double dp, double nturn
        Copy from Tracy II.
 ****************************************************************************/
 
-void Multipole_thicksext(char const *fic_hcorr, char const *fic_vcorr, char const *fic_skew)
+void Multipole_thicksext(const char *fic_hcorr, const char *fic_vcorr, const char *fic_skew)
 {
   int i = 0;
   int ndip  = 0,  /* Number of dipoles */
@@ -2297,7 +2299,8 @@ d file */
 }
 
 /****************************************************************************/
-/* void Multipole_thinsext(void)
+/* void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, 
+                           const char *fic_skew)
 
    Purpose:
        Set multipole in dipoles, quadrupoles, thin sextupoles, skew quadrupole,
@@ -2326,7 +2329,7 @@ d file */
 
 ****************************************************************************/
 
-void Multipole_thinsext(char const *fic_hcorr, char const *fic_vcorr, char const *fic_skew)
+void Multipole_thinsext(const char *fic_hcorr, const char *fic_vcorr, const char *fic_skew)
 {
   int i = 0;
   int ndip  = 0,  /* Number of dipoles */
@@ -2947,9 +2950,9 @@ void SetSkewQuad(void)
 // }
 
 /****************************************************************************/
-/* void MomentumAcceptance(long deb, long fin,
-                           double ep_min, double ep_max, long nstepp,
-                           double em_min, double em_max, long nstepm)
+/* void MomentumAcceptance(char *MomAccFile, long deb, long fin, 
+                        double ep_min, double ep_max, long nstepp, double em_min, 
+			double em_max, long nstepm, long nturn)	
    Purpose:
         Compute momemtum acceptance along the ring, track the particle with 
 	different energy, momentum acceptance is the energy when the particle
@@ -2958,16 +2961,17 @@ void SetSkewQuad(void)
          Based on the version in tracy 2.
 	
    Input:
-       deb   first element for momentum acceptance,"debut" is beginning in French
-       fin   last element for momentum acceptance,"fin"   is end in French
+       MomAccFile   file to save calculated momentum compact factor
+       deb          first element for momentum acceptance,"debut" is beginning in French
+       fin          last element for momentum acceptance,"fin"   is end in French
 
-       ep_min   minimum energy deviation for positive momentum acceptance
-       ep_max   maximum energy deviation for positive momentum acceptance
-       nstepp   number of energy steps for positive momentum acceptance
-
-       em_min minimum energy deviation for negative momentum acceptance
-       em_max maximum energy deviation for negative momentum acceptance
-       nstepm number of energy steps for negative momentum acceptance
+       ep_min       minimum energy deviation for positive momentum acceptance
+       ep_max       maximum energy deviation for positive momentum acceptance
+       nstepp       number of energy steps for positive momentum acceptance
+ 
+       em_min       minimum energy deviation for negative momentum acceptance
+       em_max       maximum energy deviation for negative momentum acceptance
+       nstepm       number of energy steps for negative momentum acceptance
 
 
        * 1 grande section droite
@@ -3002,10 +3006,14 @@ void SetSkewQuad(void)
 		       the Cell_Pass( ) is tracking from element i0+1L to i1(tracy 2).
 	   17/04/11 add number of turn
         27/06/11  fix the bug of the index in the tabz and tabpz when calling Trac( ); 
-                 fix the bug in the vertical closed orbit when calling Trac( ).                                                     		
+                 fix the bug in the vertical closed orbit when calling Trac( ). 
+        19/07/11  add the interface to save calculated momentum compact factor in the 
+	          user defined file.
+		                                                      		
 ****************************************************************************/
-void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max,
-			long nstepp, double em_min, double em_max, long nstepm, long nturn)
+void MomentumAcceptance(char *MomAccFile, long deb, long fin, 
+                        double ep_min, double ep_max, long nstepp, double em_min, 
+			double em_max, long nstepm, long nturn)		
 {
   double        dP = 0.0, dp1 = 0.0, dp2 = 0.0;
   long          lastpos = 0L,lastn = 0L;
@@ -3020,7 +3028,7 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max,
   struct tm     *newtime;  // for time
   Vector        codvector[Cell_nLocMax];
   bool          cavityflag, radiationflag;
-  bool          trace=true;
+  bool          trace=true; 
   
   x0.zero();
 
@@ -3033,12 +3041,12 @@ void MomentumAcceptance(long deb, long fin, double ep_min, double ep_max,
   /* File opening for writing */
 
   outf1 = fopen("phase.out", "w");
-  outf2 = fopen("soleil.out", "w");
+  outf2 = fopen(MomAccFile, "w");
 
-  fprintf(outf2,"# TRACY III v. SYNCHROTRON SOLEIL -- %s \n", asctime2(newtime));
+  fprintf(outf2,"# TRACY III -- %s \n", asctime2(newtime));
   fprintf(outf2,"#  i        s         dp      s_lost  name_lost \n#\n");
 
-  fprintf(outf1,"# TRACY III v. SYNCHROTRON SOLEIL -- %s \n", asctime2(newtime));
+  fprintf(outf1,"# TRACY III  -- %s \n", asctime2(newtime));
   fprintf(outf1,"#  i        x           xp            z           zp           dp          ctau\n#\n");