From 350e16da19a9312842c06f528390db770cae510b Mon Sep 17 00:00:00 2001
From: nadolski <nadolski@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Sun, 30 Jun 2013 05:12:48 +0000
Subject: [PATCH] Adding loss feature v2

---
 tracy/tools/soltracy.cc        | 553 ++++++++++++++++++---------------
 tracy/tracy/src/read_script.cc | 178 ++++++-----
 tracy/tracy/src/soleillib.cc   |  25 +-
 3 files changed, 396 insertions(+), 360 deletions(-)

diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 6bfc9e1..33f8777 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -34,8 +34,8 @@ int main(int argc, char *argv[]) {
   
   /*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
+  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 synchro radiation damping
   // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
   globval.H_exact = false;
@@ -67,8 +67,9 @@ int main(int argc, char *argv[]) {
 
    
   if (argc > 1) {
-    read_script(argv[1], true,CommandNo, UserCommandFlag);
-  } else {
+    read_script(argv[1], true, CommandNo, UserCommandFlag);
+  } 
+	else {
     fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
     exit_(1);
   }
@@ -148,11 +149,11 @@ int main(int argc, char *argv[]) {
    
   //set RF voltage  
   else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
-    printf("\nSetting RF voltage:\n");
-    printf("    Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
+    fprintf(stdout, "\nSetting RF voltage:\n");
+    fprintf(stdout, "    Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
         / 1e6);
     set_RFVoltage(globval.cav, UserCommandFlag[i].RFvolt);
-    printf("    New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
+    fprintf(stdout, "    New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
         / 1e6);
   }
 
@@ -293,10 +294,10 @@ int main(int argc, char *argv[]) {
     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", 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", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit);
+    fprintf(stdout, "Before the fitting, the quadrupole field strengths are: \n");
+    fprintf(stdout, " %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn);
+    fprintf(stdout, "After the fitting, the quadrupole field strengths are: \n");
+    fprintf(stdout, " %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);
@@ -333,10 +334,10 @@ int main(int argc, char *argv[]) {
     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", UserCommandFlag[i].qf1, qf1_bn,
+    fprintf(stdout, "    %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", UserCommandFlag[i].qf1,
+    fprintf(stdout, "    %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);
 
@@ -363,10 +364,10 @@ int main(int argc, char *argv[]) {
     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", 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", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit);
+    fprintf(stdout, "Before the fitting, the sextupole field strengths are \n");
+    fprintf(stdout, "    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn);
+    fprintf(stdout, "After the fitting, the sextupole field strengths are: \n");
+    fprintf(stdout, "    %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 */
@@ -438,90 +439,119 @@ int main(int argc, char *argv[]) {
 
   // Computes FMA
   else if(strcmp(CommandStr,"FmapFlag") == 0) {
-    printf("\n begin Fmap calculation for on momentum particles: \n");
+    fprintf(stdout, "\n begin Fmap calculation for on momentum particles: \n");
     
    #if MPI_EXEC
 
     //initialization for parallel computing
-    int myid=0, numprocs=0;
-    int namelen=0;
+    int myid = 0, numprocs = 0;
+    int namelen = 0;
     char processor_name[MPI_MAX_PROCESSOR_NAME]="";
-    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
+    
+		MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
     MPI_Comm_rank(MPI_COMM_WORLD, &myid);
     MPI_Get_processor_name(processor_name, &namelen);
-    printf("\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  
+    fprintf(stdout, "\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  
+    
+    fprintf(stdout, "\n Begin Fmap calculation for on momentum particles: \n");
     
-    printf("\n begin Fmap calculation for on momentum particles: \n");
-    //Each core or process, characterized by myid, will track different region of fmap 
+		//Each core or process, characterized by myid, will track different region of fmap 
     fmap_p(UserCommandFlag[i]._FmapFlag_fmap_file,
-	   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,
-       UserCommandFlag[i]._FmapFlag_printloss,
-       numprocs,myid);
+	    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,
+		  UserCommandFlag[i]._FmapFlag_printloss,
+		  numprocs, myid);
 
     //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
     MPI_Barrier(MPI_COMM_WORLD);
    
     //Collecting data from all files generated by cores, then write to fmap_p.out
-    if(myid==0)
-      {
-	//  system("cp 0fmap.out fmap_p.out");
+    if(myid == 0){
 	
-	FILE *outf;
-	if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL)
-	  {
-	    fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
-	    exit_(1);
-	  }
-
-	FILE *fp;
-	char buffer[81];
-	char FmapFile[50];
-	for(int j=0; j<numprocs; j++)
-	  {
-	    FmapFile[0]='\0';
-	    sprintf(FmapFile,"%d",j);
-	    strcat(FmapFile,UserCommandFlag[i]._FmapFlag_fmap_file);
-      
-	    if((fp = fopen(FmapFile, "r")) == NULL)  
-	      {
-		fprintf(stdout, "%s: error while opening file.\n", FmapFile);
-		exit_(1);
-	      }
-
-	    while(fgets(buffer, 80, fp) != NULL)
-	      {
-		fputs(buffer, outf);
-		buffer[0]='\0';
-	      }
-	    fclose(fp);
-	  }
-	fclose(outf);
-      }
+			FILE *outf;
+			 
+			// open frequency map file for final results
+			if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL){
+					fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
+					exit_(1);
+				}
+
+			// open frequency map loss file for final results
+			if (UserCommandFlag[i]._FmapFlag_printloss){
+				FILE *outf_loss;
+				char FmapLossFile[max_str+5]=" ";
+				strcpy(FmapLossFile, UserCommandFlag[i]._FmapFlag_fmap_file);
+				strcat(FmapLossFile, ".loss");
+				if ((outf_loss = fopen(FmapLossFile, "w")) == NULL){
+					fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
+					exit_(1);
+				}
+			}
+
+			FILE *fp, *fp_loss;
+			char buffer[81];
+			char FmapFile[max_str];
+		   
+			for(int j = 0; j < numprocs; j++){
+				FmapFile[0]='\0';
+				sprintf(FmapFile,"%d",j);
+				strcat(FmapFile, UserCommandFlag[i]._FmapFlag_fmap_file);
+				
+				if((fp = fopen(FmapFile, "r")) == NULL)  {
+					fprintf(stdout, "%s: error while opening file.\n", FmapFile);
+					exit_(1);
+					while(fgets(buffer, 80, fp) != NULL){
+						fputs(buffer, outf);
+						buffer[0]='\0';
+					}
+					fclose(fp);
+				}
+
+				if (UserCommandFlag[i]._FmapFlag_printloss){
+					strcpy(FmapLossFile,FmapFile);
+					strcat(FmapLossFile,".loss");
+				
+					if((fp_loss = fopen(FmapFile, "r")) == NULL)  {
+						fprintf(stdout, "%s: error while opening file.\n", FmapFile);
+						exit_(1);
+						while(fgets(buffer, 80, fp_loss) != NULL){
+							fputs(buffer, outf_loss);
+							buffer[0]='\0';
+						}
+					}
+					fclose(fp_loss);
+				}
+			}
+			fclose(outf);
+				
+			if (UserCommandFlag[i]._FmapFlag_printloss){
+				fclose(outf_loss);
+			}
+		}
     MPI_Barrier(MPI_COMM_WORLD); 
     printf("Process %d finished for on momentum fmap.\n", myid);
 
- #else
-      fmap(UserCommandFlag[i]._FmapFlag_fmap_file,
-	   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,
-	   UserCommandFlag[i]._FmapFlag_printloss);
-#endif
+ #else // Single processor
+		fmap(UserCommandFlag[i]._FmapFlag_fmap_file,
+	    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,
+	    UserCommandFlag[i]._FmapFlag_printloss);
+#endif // MPI
   }
 
   // Compute FMA dp
   else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
-    printf("\n begin Fmap calculation for off momentum particles: \n");
+    fprintf(stdout, "\n begin Fmap calculation for off momentum particles: \n");
     
    #if MPI_EXEC
 
@@ -534,70 +564,94 @@ int main(int argc, char *argv[]) {
     MPI_Get_processor_name(processor_name, &namelen);
     printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  
 
-    printf("\n begin Fmap calculation for off momentum particles: \n");
+    fprintf(stdout, "\n begin Fmap calculation for off momentum particles: \n");
     fmapdp_p(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
-	     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,
-         UserCommandFlag[i]._FmapFlag_printloss,
-         numprocs,myid);
+			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,
+			UserCommandFlag[i]._FmapspFlag_printloss, numprocs, myid);
 
     //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
     MPI_Barrier(MPI_COMM_WORLD);
    
     //Collecting data from all files generated by cores, then write to fmap_p.out
-    if(myid==0)
-     {
-       //  system("cp 0fmap.out fmap_p.out");
-
-       FILE *outf;
-       if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL)
-	 {
-	   fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
-	   exit_(1);
-	 }
-
-       FILE *fp;
-       char buffer[81];
-       char FmapFile[50];
-       for(int j=0; j<numprocs; j++)
-	 {
-	   FmapFile[0]='\0';
-	   sprintf(FmapFile,"%d",j);
-	   strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
-      
-	   if((fp = fopen(FmapFile, "r")) == NULL)  
-	     {
-	       fprintf(stdout, "%s: error while opening file.\n", FmapFile);
-	       exit_(1);
-	     }
-
-	   while(fgets(buffer, 80, fp) != NULL)
-	     {
-	       fputs(buffer, outf);
-	       buffer[0]='\0';
-	     }
-	   fclose(fp);
-	 }
-       fclose(outf);
-     }
-    MPI_Barrier(MPI_COMM_WORLD); 
-    printf("Process %d finished off momentum fmap.\n", myid);
-
- #else
-       fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
-	      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,
-	       UserCommandFlag[i]._FmapdpFlag_printloss);
+    if(myid == 0){
+			FILE *outf;
+			if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL){
+				fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
+				exit_(1);
+			}
+
+		  // open frequency map loss file for final results
+			if (UserCommandFlag[i]._FmapdpFlag_printloss){
+				FILE *outf_loss;
+				char FmapdpLossFile[max_str+5]=" ";
+				strcpy(FmapdpLossFile, UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
+				strcat(FmapdpLossFile, ".loss");
+				if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL){
+					fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmap_file);
+					exit_(1);
+				}
+			}
+
+			FILE *fp, *fp_loss;
+			char buffer[81];
+			char FmapFile[50];
+			for(int j = 0; j < numprocs; j++){
+			  FmapFile[0]='\0';
+			  sprintf(FmapFile,"%d",j);
+			  strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
+				
+			  if((fp = fopen(FmapFile, "r")) == NULL){
+					fprintf(stdout, "%s: error while opening file.\n", FmapFile);
+					exit_(1);
+				} 
+
+				while(fgets(buffer, 80, fp) != NULL){
+					fputs(buffer, outf);
+					buffer[0]='\0';
+				}
+				fclose(fp);
+				
+				if (UserCommandFlag[i]._FmapdpFlag_printloss){
+					strcpy(FmapdpLossFile, FmapdpFile);
+					strcat(FmapdpLossFile,".loss");
+				
+					if((fp_loss = fopen(FmapdpFile, "r")) == NULL)  {
+						fprintf(stdout, "%s: error while opening file.\n", FmapdpFile);
+						exit_(1);
+						while(fgets(buffer, 80, fp_loss) != NULL){
+							fputs(buffer, outf_loss);
+							buffer[0]='\0';
+						}
+					}
+					fclose(fp_loss);
+				} // if printloss
+			} // for loop
+			fclose(outf);
+				
+			if (UserCommandFlag[i]._FmapFlag_printloss){
+				fclose(outf_loss);
+			}
+		} // if myid
+
+		MPI_Barrier(MPI_COMM_WORLD); 
+    fprintf(stdout, "Process %d finished off momentum fmap.\n", myid);
+
+ #else //single processor
+	 fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
+		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,
+		UserCommandFlag[i]._FmapdpFlag_printloss);
  #endif
   }
 
@@ -626,7 +680,7 @@ int main(int argc, char *argv[]) {
     };
 
  /* calculate momentum acceptance*/
-    printf("\n Calculate momentum acceptance: \n");
+    fprintf(stdout, "\n Calculate momentum acceptance: \n");
 
  #if MPI_EXEC
     
@@ -640,126 +694,112 @@ int main(int argc, char *argv[]) {
     MPI_Get_processor_name(processor_name, &namelen);
     printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name);  
 
-    printf("\n Calculate momentum acceptance: \n");
+    fprintf(stdout, "\n Calculate momentum acceptance: \n");
     MomentumAcceptance_p(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
-                       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,
-		       UserCommandFlag[i]._MomentumAccFlag_nturn,
-		       UserCommandFlag[i]._MomentumAccFlag_zmax,
-                       numprocs,myid);
+												 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,
+												 UserCommandFlag[i]._MomentumAccFlag_nturn,
+												 UserCommandFlag[i]._MomentumAccFlag_zmax,
+												 numprocs,myid);
   //merge the files 
     //Synchronize all cores, till finish cal of different region, 
     //then files from the cores will be processed.
     MPI_Barrier(MPI_COMM_WORLD);
 
     //Collect data from all files generated by different cores, then write to user defined file
-    if(myid==0)
-      {
-	//merge the phase.out from different cpus
-	FILE *outf1;  //phase.out, for debugging
-	FILE *fp1;
-	char buffer1[81];
-	char PhaseFile[50];
-	if ((outf1 = fopen("phase_p.out", "w")) == NULL)
-	  {
-	    fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out");
-	    exit_(1);
-	  }
-
-	for(int j=0; j<numprocs; j++)
-	  {
-	    PhaseFile[0]='\0';
-	    sprintf(PhaseFile,"%d",j);
-	    strcat(PhaseFile,"phase.out");
-	    
-	    if((fp1 = fopen(PhaseFile, "r")) == NULL){
-		fprintf(stdout, "%s: error while opening file.\n", PhaseFile);
-		exit_(1);
-	      }
-
-	    while(fgets(buffer1, 80, fp1) != NULL)
-	      {
-		fputs(buffer1, outf1);
-		buffer1[0]='\0';
-	      }
-	    fclose(fp1);
-	  }
-
-	//collect data for the momentum acceptance
-	FILE *outf2; //momentum acceptance
-	FILE *fp2;
-	char buffer2[81];
-	char MAFile[50];
-	if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL)
-	  {
-	    fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file);
-	    exit_(1);
-	  }
-
-	//Collect data in Positive region
-	for(int j=0; j<numprocs; j++)
-	  {
-	    MAFile[0]='\0';
-	    sprintf(MAFile,"%d",j);
-	    strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
-
-	    if((fp2 = fopen(MAFile, "r")) == NULL)  
-	      {
-		fprintf(stdout, "%s: error while opening file.\n", MAFile);
-		exit_(1);
-	      }
+    if(myid==0){
+			//merge the phase.out from different cpus
+			FILE *outf1;  //phase.out, for debugging
+			FILE *fp1;
+			char buffer1[81];
+			char PhaseFile[50];
+			if ((outf1 = fopen("phase_p.out", "w")) == NULL){
+					fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out");
+					exit_(1);
+				}
+
+			for(int j=0; j<numprocs; j++){
+				PhaseFile[0]='\0';
+				sprintf(PhaseFile,"%d",j);
+				strcat(PhaseFile,"phase.out");
+				if((fp1 = fopen(PhaseFile, "r")) == NULL){
+					fprintf(stdout, "%s: error while opening file.\n", PhaseFile);
+					exit_(1);
+					}
+
+				while(fgets(buffer1, 80, fp1) != NULL){
+					fputs(buffer1, outf1);
+					buffer1[0]='\0';
+				}
+				fclose(fp1);
+			} // for loop
+
+		//collect data for the momentum acceptance
+		FILE *outf2; //momentum acceptance
+		FILE *fp2;
+		char buffer2[81];
+		char MAFile[50];
+		if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL){
+			fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file);
+			exit_(1);
+			}
+
+		//Collect data in Positive region
+		for(int j=0; j<numprocs; j++){
+		 MAFile[0]='\0';
+		 sprintf(MAFile,"%d",j);
+			strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
+
+			if((fp2 = fopen(MAFile, "r")) == NULL){
+				fprintf(stdout, "%s: error while opening file.\n", MAFile);
+				exit_(1);
+			}
 	 
-	    while(fgets(buffer2, 80, fp2) != NULL)
-	      {
-		if(strstr(buffer2,"Negative")!=0) break;
-
-		fputs(buffer2, outf2);
-		buffer2[0]='\0'; //reset buffer to NULL
-	      }
-	    buffer2[0]='\0';
-	    fclose(fp2);
-	  }
-
-	fprintf(outf2,"\n"); // A void line
-
-	//Collect data in Nagative region
-	for(int j=0; j<numprocs; j++)
-	  {
-	    MAFile[0]='\0';
-	    sprintf(MAFile,"%d",j);
-	    strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
-	    
-	    if((fp2 = fopen(MAFile, "r")) == NULL)  
-	      {
-		fprintf(stdout, "%s: error while opening file.\n", MAFile);
-		exit_(1);
-	      }
-
-	    bool found=false;
-	    while(fgets(buffer2, 80, fp2) != NULL)
-	      {
-		if( (strstr(buffer2,"Negative")==0) && (!found)) {  buffer2[0]='\0'; continue;}
-		if( (strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;}
-		
-		fputs(buffer2, outf2);
-		buffer2[0]='\0';
-	      }
-	    buffer2[0]='\0';
-	    fclose(fp2);
-	  }
-	fclose(outf2);
-	
-      }
+			while(fgets(buffer2, 80, fp2) != NULL){
+				if(strstr(buffer2,"Negative")!=0) break;
+				fputs(buffer2, outf2);
+				buffer2[0]='\0'; //reset buffer to NULL
+			}
+			buffer2[0]='\0';
+			fclose(fp2);
+		}
+
+		fprintf(outf2,"\n"); // A void line
+
+		//Collect data in Nagative region
+		for(int j=0; j<numprocs; j++){
+			MAFile[0]='\0';
+			sprintf(MAFile,"%d",j);
+			strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
+			
+			if((fp2 = fopen(MAFile, "r")) == NULL){
+				fprintf(stdout, "%s: error while opening file.\n", MAFile);
+				exit_(1);
+			}
+
+			bool found=false;
+			while(fgets(buffer2, 80, fp2) != NULL){
+				if((strstr(buffer2,"Negative")==0) && (!found)) {  buffer2[0]='\0'; continue;}
+				if((strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;}
+				
+				fputs(buffer2, outf2);
+				buffer2[0]='\0';
+			} // while
+				buffer2[0]='\0';
+				fclose(fp2);
+		} // for
+		fclose(outf2);
+	}
     MPI_Barrier(MPI_COMM_WORLD); 
     printf("process %d finished momentum acceptance.\n", myid);
 
-#else
+#else // single processor coputation
      MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
 			UserCommandFlag[i]._MomentumAccFlag_istart, 
 			UserCommandFlag[i]._MomentumAccFlag_istop,
@@ -780,11 +820,10 @@ int main(int argc, char *argv[]) {
 
   // induced amplitude 
   else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) {
-   printf("\n Calculate induced amplitude: \n");
-    InducedAmplitude(193L);
+		fprintf(stdout, "\n Calculate induced amplitude: \n");
+		InducedAmplitude(193L);
   }
 
-
   else if(strcmp(CommandStr,"EtaFlag") == 0) {
     // compute cod and Twiss parameters for different energy offsets
     for (int ii = 0; ii <= 40; ii++) {
@@ -829,14 +868,14 @@ int main(int argc, char *argv[]) {
     Phase(UserCommandFlag[i]._Phase_phase_file,
           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);
+					UserCommandFlag[i]._Phase_Y, 
+					UserCommandFlag[i]._Phase_Py, 
+					UserCommandFlag[i]._Phase_delta, 
+					UserCommandFlag[i]._Phase_ctau, 
+					UserCommandFlag[i]._Phase_nturn);
     printf("6D phase space tracking: \n the simulation time for phase space in tracy 3 is \n");
     stop = stampstop(start);
-
+		
     /* restore the initial values*/
     globval.Cavity_on = cavityflag;
     globval.radiation = radiationflag;
@@ -966,7 +1005,7 @@ Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc
       cout << "Number of quadrupole families used for ID correction:   " << N_Fam << endl; 
       cout << "Quadrupoles used for ID correction: " << endl;
       for (int k = 0; k < N_Fam; k++)
-        printf("%d\n",Q_Fam[k]);
+        fprintf(stdout, "%d\n",Q_Fam[k]);
     } 
     
     //initialization
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index ef7d461..dc50e7b 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -35,30 +35,30 @@ char girder_file[max_str];
 
 
 
- // multipole files; for soleil lattice
+// multipole files; for soleil lattice
 char multipole_file[max_str];
-  // multipole file for soleil lattice 
+// multipole file for soleil lattice 
 char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
- // files to set sources of coupling; for SOLEIL lattice
+// files to set sources of coupling; for SOLEIL lattice
 char virtualskewquad_file[max_str];
 
 //errors
 //char fe_file[max_str]; //the same as multipole_file[max_str]????
 
 /* COD correction */
-int nwh = 56, nwv = 56; //number of singular values for SVD correction
+int nwh = 57, nwv = 57; //number of singular values for SVD correction
 
 /* ID compensation */  
-char   IDCq_name[max_str][11];
+char IDCq_name[max_str][11];
 
 
 //multipole field error correction (SOLEIL based)
 char hcorr_file[max_str], vcorr_file[max_str]; //files with the status of hcorr/vcorr status, 
                                                   //to choose which correctors are used for orbit correction
 //default name of special elements in lattice
-char hcorr_name[max_str] = "", vcorr_name[max_str] = "";
-char skew_quad_name[max_str] = "", bpm_name[max_str] = "";
-char gs_name[max_str] = "", ge_name[max_str] = "";
+char hcorr_name[max_str]     = "", vcorr_name[max_str] = "";
+char skew_quad_name[max_str] = "", bpm_name[max_str]   = "";
+char gs_name[max_str]        = "", ge_name[max_str]    = "";
 						  
 #define OLD_LATTICE
 
@@ -226,78 +226,76 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
        //print close orbit(COD) flag
       else if (strcmp("PrintCODFlag", name) == 0){
         sscanf(line, "%*s %s",nextpara);
-	 
-	if(strcmp(nextpara,"voidpara")!=0)
-           sscanf(line, "%*s %s", UserCommandFlag[CommNo]._PrintCOD_cod_file);
-          
-	strcpy(UserCommandFlag[CommNo].CommandStr,name);  
+	      if(strcmp(nextpara,"voidpara")!=0)
+					sscanf(line, "%*s %s", UserCommandFlag[CommNo]._PrintCOD_cod_file);
+          strcpy(UserCommandFlag[CommNo].CommandStr,name);  
       }
       //print the cooridinates using tracking at each element
       else if (strcmp("PrintTrackFlag", name) == 0){
-	  sscanf(line, "%*s %s",nextpara);
+	      sscanf(line, "%*s %s",nextpara);
 	 
-	  if(strcmp(nextpara,"voidpara")!=0)
-	    sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld", 
-	                UserCommandFlag[CommNo]._PrintTrack_track_file,
-          		&(UserCommandFlag[CommNo]._PrintTrack_x), 
-			&(UserCommandFlag[CommNo]._PrintTrack_px),
-          		&(UserCommandFlag[CommNo]._PrintTrack_y), 
-			&(UserCommandFlag[CommNo]._PrintTrack_py), 
-			&(UserCommandFlag[CommNo]._PrintTrack_delta), 
-			&(UserCommandFlag[CommNo]._PrintTrack_ctau),
-			&(UserCommandFlag[CommNo]._PrintTrack_nmax));
-	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
-      }  
-      //print close orbit(COD) flag
+			if(strcmp(nextpara,"voidpara")!=0)
+				sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld", 
+					UserCommandFlag[CommNo]._PrintTrack_track_file,
+					&(UserCommandFlag[CommNo]._PrintTrack_x), 
+				  &(UserCommandFlag[CommNo]._PrintTrack_px),
+					&(UserCommandFlag[CommNo]._PrintTrack_y), 
+				  &(UserCommandFlag[CommNo]._PrintTrack_py), 
+					&(UserCommandFlag[CommNo]._PrintTrack_delta), 
+					&(UserCommandFlag[CommNo]._PrintTrack_ctau),
+					&(UserCommandFlag[CommNo]._PrintTrack_nmax));
+				strcpy(UserCommandFlag[CommNo].CommandStr,name);
+			}
+      
+			//print close orbit(COD) flag
       else if (strcmp("PrintGirderFlag", name) == 0){
-          sscanf(line, "%*s %s", girder_file);
-	  strcpy(UserCommandFlag[CommNo].CommandStr,name);  
+				sscanf(line, "%*s %s", girder_file);
+				strcpy(UserCommandFlag[CommNo].CommandStr,name);  
       }
       
       
       else if (strcmp("TuneTracFlag", name) == 0){
-	 strcpy(UserCommandFlag[CommNo].CommandStr,name);
+				strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
       else if (strcmp("ChromTracFlag", name) == 0){
-	   strcpy(UserCommandFlag[CommNo].CommandStr,name);
+				strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }
-      // FMA
-      else if (strcmp("FmapFlag", name) == 0){        
-    	  strcpy(dummy, "");
-	  strcpy(dummy2,"");
-	  sscanf(line, "%*s %s",nextpara);
+      //FMA
+      else if (strcmp("FmapFlag", name) == 0){
+    	  strcpy(dummy , "");
+				strcpy(dummy2, "");
+				sscanf(line, "%*s %s",nextpara);
 	 
-	  //if no definition of loss flag in the lattice file, then "dummy2" is empty string.
-	  if(strcmp(nextpara,"voidpara")!=0)
-            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %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,dummy2);
-         
-        if(strcmp(dummy, "true") == 0)
-          UserCommandFlag[CommNo]._FmapFlag_diffusion = true;
-        else if(strcmp(dummy, "false") == 0)
-          UserCommandFlag[CommNo]._FmapFlag_diffusion = false;
-       
-	if(strcmp(dummy2, "true") == 0)
-          UserCommandFlag[CommNo]._FmapFlag_printloss = true;
-        else if(strcmp(dummy2, "false") == 0)
-          UserCommandFlag[CommNo]._FmapFlag_printloss = false;
-	else
-	  UserCommandFlag[CommNo]._FmapFlag_printloss = false;
-	
-	
-	//cout << "debug:      " << line << endl;
-	//cout << "dummy2 =  " << dummy2 << "     lossflag =  " << UserCommandFlag[CommNo]._FmapFlag_printloss << endl;
-	
-       // FmapFlag = true;
-         strcpy(UserCommandFlag[CommNo].CommandStr,name);
-      }
+				//if no definition of loss flag in the lattice file, then "dummy2" is empty string.
+				if(strcmp(nextpara,"voidpara")!= 0)
+					sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %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,dummy2);
+						 
+				if(strcmp(dummy, "true") == 0)
+					UserCommandFlag[CommNo]._FmapFlag_diffusion = true;
+				else if(strcmp(dummy, "false") == 0)
+					UserCommandFlag[CommNo]._FmapFlag_diffusion = false;
+					 
+				if(strcmp(dummy2, "true") == 0)
+						UserCommandFlag[CommNo]._FmapFlag_printloss = true;
+				else if(strcmp(dummy2, "false") == 0)
+						UserCommandFlag[CommNo]._FmapFlag_printloss = false;
+				else
+					UserCommandFlag[CommNo]._FmapFlag_printloss = false;
+			
+				//cout << "debug:      " << line << endl;
+				//cout << "dummy2 =  " << dummy2 << "     lossflag =  " << UserCommandFlag[CommNo]._FmapFlag_printloss << endl;
+			
+					 // FmapFlag = true;
+					strcpy(UserCommandFlag[CommNo].CommandStr,name);
+			}
       // FMA dp
       else if (strcmp("FmapdpFlag", name) == 0){         
     	  strcpy(dummy, "");
@@ -357,7 +355,7 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
 	   sscanf(line, "%*s %s %ld %ld %lf",  
 	          UserCommandFlag[CommNo]._EnergyTuneShift_nudp_file,
 	         &(UserCommandFlag[CommNo]._EnergyTuneShift_npoint),
-                 &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), 
+					 &(UserCommandFlag[CommNo]._EnergyTuneShift_nturn), 
 	         &(UserCommandFlag[CommNo]._EnergyTuneShift_deltamax));
         
 	strcpy(UserCommandFlag[CommNo].CommandStr,name);
@@ -386,13 +384,13 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
 		      UserCommandFlag[CommNo]._MomentumAccFlag_TrackDim,
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_istart),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_istop),
-        	      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminp),
+					&(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminp),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxp),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepp),
-        	      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminn),
+					&(UserCommandFlag[CommNo]._MomentumAccFlag_deltaminn),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_deltamaxn),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_nstepn),
-                      &(UserCommandFlag[CommNo]._MomentumAccFlag_nturn),
+					&(UserCommandFlag[CommNo]._MomentumAccFlag_nturn),
 		      &(UserCommandFlag[CommNo]._MomentumAccFlag_zmax));
 	 } 
 	 strcpy(UserCommandFlag[CommNo].CommandStr,name);   
@@ -406,15 +404,15 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
       }
       // 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",UserCommandFlag[CommNo].qf1,UserCommandFlag[CommNo].qf2,
-	UserCommandFlag[CommNo].qd1,UserCommandFlag[CommNo].qd2,
-	&(UserCommandFlag[CommNo].targetnux),&(UserCommandFlag[CommNo].targetnuz));
-        strcpy(UserCommandFlag[CommNo].CommandStr,name);
+         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",UserCommandFlag[CommNo].sxm1,UserCommandFlag[CommNo].sxm2,
-	&(UserCommandFlag[CommNo].targetksix),&(UserCommandFlag[CommNo].targetksiz));
-        strcpy(UserCommandFlag[CommNo].CommandStr,name);
+         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);
@@ -432,7 +430,7 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
 //         
 //       }
        else if (strcmp("InducedAmplitudeFlag", name) == 0){
-        strcpy(UserCommandFlag[CommNo].CommandStr,name); 
+         strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       }
 //        else if (strcmp("CodeComparaisonFlag", name) == 0){
 //        strcpy(UserCommand[CommandNo-1],name);
@@ -527,27 +525,27 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
       else if (strcmp("IDCorrFlag", name) == 0)
         strcpy(UserCommandFlag[CommNo].CommandStr,name); 
       else if (strcmp("scl_dbetax", name) == 0)//scaling of dbetax, for ID correction
-	sscanf(line, "%*s %lf", &scl_dbetax); 
+	      sscanf(line, "%*s %lf", &scl_dbetax); 
       else if (strcmp("scl_dbetay", name) == 0)//scaling of dbetay, for ID correction
-	sscanf(line, "%*s %lf", &scl_dbetay); 
+				sscanf(line, "%*s %lf", &scl_dbetay); 
       else if (strcmp("scl_dnux", name) == 0)//scaling of dnux, for ID correction
-	sscanf(line, "%*s %lf", &scl_dnux); 
+	      sscanf(line, "%*s %lf", &scl_dnux); 
       else if (strcmp("scl_dnuy", name) == 0)//scaling of dnuy, for ID correction
-	sscanf(line, "%*s %lf", &scl_dnuy); 
+	      sscanf(line, "%*s %lf", &scl_dnuy); 
       else if (strcmp("scl_nux", name) == 0) //scaling of nux, for ID correction
-	sscanf(line, "%*s %lf", &scl_nux); 
+	      sscanf(line, "%*s %lf", &scl_nux); 
       else if (strcmp("scl_nuy", name) == 0) //scaling of nuy, for ID correction
-	sscanf(line, "%*s %lf", &scl_nuy); 
+	      sscanf(line, "%*s %lf", &scl_nuy); 
       else if (strcmp("ID_step", name) == 0)//ID steps, for ID correction
-	sscanf(line, "%*s %lf", &ID_step); 
+	      sscanf(line, "%*s %lf", &ID_step); 
       else if (strcmp("N_calls", name) == 0) // ID correction parameters
-	sscanf(line, "%*s %d", &N_calls);
+	      sscanf(line, "%*s %d", &N_calls);
       else if (strcmp("N_steps", name) == 0)
-	sscanf(line, "%*s %d", &N_steps);
+	      sscanf(line, "%*s %d", &N_steps);
       else if (strcmp("N_Fam", name) == 0)
-	sscanf(line, "%*s %d", &N_Fam);
+	      sscanf(line, "%*s %d", &N_Fam);
       else if (strcmp("IDCquads", name) == 0) {
-	sscanf(line, "%*s %s %s %s %s %s %s %s %s %s %s %s",
+				sscanf(line, "%*s %s %s %s %s %s %s %s %s %s %s %s",
 	       IDCq_name[0], IDCq_name[1], IDCq_name[2], IDCq_name[3],
 	       IDCq_name[4], IDCq_name[5], IDCq_name[6], IDCq_name[7],
 	       IDCq_name[8], IDCq_name[9], IDCq_name[10]);
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index 1f6f106..705fd2f 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -223,7 +223,7 @@ Vector2 H;
        at position pos and for a energy offset dP
 
        For a givien delta
-       H = gamma xcod² + 2*alpha*xcod*xcod' + beta*xcod'*xcod'
+       H = gamma xcod² + 2*alpha*xcod*xcod' + beta*xcod'*xcod'
        
    Input:
        none
@@ -875,7 +875,7 @@ void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, do
  FILE *outf_loss; //file with the information of the lost particle
  long i = 0L, j = 0L;
  double Tab[DIM][NTURN], Tab0[DIM][NTURN];
- double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx, dfy;
+ double fx[NTERM2], fy[NTERM2], fx2[NTERM2], fy2[NTERM2], dfx = 0.0, dfy = 0.0;
  double x = 0.0, xp = 0.0, y = 0.0, zp = 0.0;
  double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
  const double ctau = 0.0;
@@ -892,12 +892,11 @@ void fmap(const char *FmapFile, long Nbx, long Nby, long Nbtour, double xmax, do
  long lastpos = 1;
  ss_vect<double> x1;
  
- 
- //if(loss){
- char FmapLossFile[S_SIZE+5]=" ";
- strcpy(FmapLossFile,FmapFile);
- strcat(FmapLossFile,".loss");
- //}
+ if(printloss){
+	 char FmapLossFile[S_SIZE+5]=" ";
+	 strcpy(FmapLossFile,FmapFile);
+	 strcat(FmapLossFile,".loss");
+ }
  
  /* Get time and date */
  time_t aclock;
@@ -1081,15 +1080,15 @@ void fmap_p(const char *FmapFile_p, long Nbx, long Nby, long Nbtour, double xmax
  char FmapFile[max_str];
  sprintf(FmapFile,"%d",myid);
  strcat(FmapFile,FmapFile_p);
- fprintf(stdout, "%s\n",FmapFile);
+ if (trace) fprintf(stdout, "%s\n",FmapFile);
 
 //variables of the returned the tracked particle
- long lastn = 0;
- long lastpos = 1;
+ long lastn = 0L;
+ long lastpos = 1L;
  ss_vect<double> x1;
 
  char FmapLossFile[S_SIZE+5]=" ";
- strcpy(FmapLossFile,FmapFile);
+ strcpy(FmapLossFile, FmapFile);
  strcat(FmapLossFile,".loss");
 
 
@@ -1484,7 +1483,7 @@ void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double
  long lastpos = 1;
  ss_vect<double> x1;
 
- char FmapdpLossFile[S_SIZE+5]=" ";
+ char FmapdpLossFile[max_str+5]=" ";
  strcpy(FmapdpLossFile,FmapdpFile);
  strcat(FmapdpLossFile,".loss");
 
-- 
GitLab