Changeset 23 in TRACY3 for branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc

Timestamp:

Dec 6, 2013, 5:12:43 PM (10 years ago)

Author:

zhangj

Message:

Clean version of Tracy: SoleilVersion at the end of 2011.Use this clean version to find the correct dipole fringe field to have the correct FMAP and FMAPDP of ThomX. Modified files: tpsa_lin.cc, soleillib.cc, prtmfile.cc, rdmfile.cc, read_script.cc, physlib.cc, tracy.cc, t2lat.cc, t2elem.cc, naffutils.cc in /tracy/tracy/src folder; naffutils.h, tracy_global.h, physlib.h, tracy.h, read_script.h, solielilib.h, t2elem.h in /tracy/tracy.inc folder; soltracy.cc in tracy/tools folder

File:

• branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc (modified) (87 diffs)

branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc

@@ -25 +25 @@
        none
 
-       
    Global variables:
        trace
@@ -38 +37 @@
 void Get_Disp_dp(void)
 {
-  long i=0L;
+  long i;
 //  long lastpos = 0;
   const char nomfic[] = "dispersion.out";
@@ -96 +95 @@
 {
   Vector        x1;     /* tracking coordinates */
-  long          i = 0L, k = 0L, imax = 50L;
+  long          i = 0L, k = 0L, imax = 50;
   FILE *        outf;
   double        dP = 0.0, dP20 = 0.0, dpmax = 0.06;
   Vector2       amp = {0.0, 0.0}, H = {0.0, 0.0};
   const char    nomfic[] = "amp_ind.out";
-  long          lastpos = 0L;
+  long          lastpos = 0;
   CellType      Celldebut, Cell;
   Vector        codvector[Cell_nLocMax];
@@ -201 +200 @@
 {
   CellType Cell;
-  long i=0L;
+  long i;
 Vector2 H;
 
@@ -252 +251 @@
 {
   CellType Cell;
-  long i=0L;
+  long i;
   long lastpos = 1L;
   Vector2 H;
@@ -462 +461 @@
        See also LoadApers in nsrl-ii.cc
        J.Zhang 07/10 soleil
-       
-      (1) 15/10/2013 Jianfeng Zhang @ LAL
-         Replace fgets() by getline() to fix the bug 
-         to read arbritrary line length.
 ****************************************************************************/
 void ReadCh(const char *AperFile)
 {
- char    Name1[max_str], Name2[max_str];
- char    *line = NULL;
- size_t  len = 0;
- ssize_t read;
+ char    in[max_str], Name1[max_str], Name2[max_str];
+ char    *line;
   int     Fnum1=0, Fnum2=0, Kidnum1=0, Kidnum2=0, k1=0, k2=0; 
   int     i=0, j=0,LineNum=0;
@@ -480 +473 @@
 
   fp = file_read(AperFile);
-  if(fp == NULL){
-    printf("ReadCh(): Error! Failure to read file %s\n", AperFile);
-    exit_(1);
-  }
-  
-  printf("\n Loading and setting vacuum apertures to lattice elements...\n");
-
-  while ((read=getline(&line, &len, fp)) != -1) {
-    
-      /* count the line number that has been read*/
-    LineNum++;
-    if(!prt){
-      printf("Line # %ld : \n",LineNum);
-      printf("Retrieved line of length %zu :\n",read);
-      printf("%s\n",line);      
-    }
-    
+
+  printf("\n");
+  printf("Loading and setting vacuum apertures to lattice elements...\n");
+
+  while (line=fgets(in, max_str, fp)) {
+  /* kill preceding whitespace generated by "table" key
+        or "space" key, but leave \n 
+        so we're guaranteed to have something*/
+     while(*line == ' ' || *line == '\t') {
+       line++;
+     }       
+    /* count the line number that has been read*/
+    LineNum++;  
     /* NOT read comment line or blank line with the end of line symbol '\n','\r' or '\r\n'*/
-    if (strstr(line, "#") == NULL && line[0] != '\n' &&
-        line[0] != '\r' &&strcmp(line,"\r\n") != 0)
+    if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0 &&
+        strcmp(line,"\r") != 0 &&strcmp(line,"\r\n") != 0)
     /* read the aperture setting */ 
     {
@@ -570 +559 @@
  //     printf("%s", line);
   } 
-  free(line);
   fclose(fp);
   // turn on the global flag for CheckAmpl()
@@ -617 +605 @@
   bool lostF = true; /* Lost particle Flag */
   Vector x1;            /* tracking coordinates */
-  long i=0L;
+  long i;
   Vector2  aperture;
   aperture[0] = 1e0; aperture[1] = 1e0;
@@ -811 +799 @@
 #undef nterm
 
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 double get_D(const double df_x, const double df_y)
 {
-  double  D=0.0;
+  double  D;
 
   const double D_min = -2.0, D_max = -10.0;
@@ -851 +834 @@
    Input:
        FmapFile file to save calculated frequency map analysis
-       Nbx                horizontal step number
-       Nby                vertical step number
-       xmax               horizontal maximum amplitude
-       zmax              vertical maximum amplitude
-       Nbtour            number of turn for tracking
-       energy            particle energy offset
-       diffusion        flag to calculate tune diffusion/ tune 
-                        difference between the first Nbtour and last Nbtour
+       Nbx    horizontal step number
+       Nby    vertical step number
+       xmax   horizontal maximum amplitude
+       zmax   vertical maximum amplitude
+       Nbtour number of turn for tracking
+       energy particle energy offset
        matlab  set file print format for matlab post-process; specific for nsrl-ii
        
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -895 +876 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2 = true;
+ bool status = true;
  struct tm *newtime;
 
@@ -946 +927 @@
  //    z  = z0 + sgn(j)*sqrt((double)abs(j))*zstep;
      // tracking around closed orbit
-     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
-     if (status2) { // if trajectory is stable
+     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
+     if (status) { // if trajectory is stable
        // gets frequency vectors
        Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
@@ -993 +974 @@
 #undef NTERM2
 
-/****************************************************************************/
-/* void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion, bool loss)
-
-   Purpose:
-       
-          
-       Compute a frequency map of Nbx x Nbz points
-       For each set of initial conditions the particle is tracked over
-       Nbtour for an energy offset dp
-
-       Frequency map is based on fixed beam energy, trace x versus z,
-       or, tracking transverse dynamic aperture for fixed momentum
-       (usually, on-momentum) particle.
-       
-       The stepsize follows a square root law
-
-       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
-       zmax              vertical maximum amplitude
-       Nbtour            number of turn for tracking
-       energy            particle energy offset
-       diffusion        flag to calculate tune diffusion/ tune 
-                        difference between the first Nbtour and last Nbtour
-       matlab  set file print format for matlab post-process; specific for nsrl-ii
-       
-   Output:
-       status2 true if stable
-              false otherwise
-
-   Return:
-       none
-
-   Global variables:
-       none
-
-   Specific functions:
-       Trac_Simple, Get_NAFF
-
-   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
-                
-       11/06/2013  Modified by Jianfeng Zhang @ LAL 
-       The same as famp(onst char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion); but with a flag to print/not print the final information of 
-          the particle; if the final turn is not the "Nbtour", then the particle is lost. 
-****************************************************************************/
-#define NTERM2  10
-void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion, bool printloss)          
-{
- FILE * outf;   //file with the tunes at the grid point
- 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], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
- double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
- double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
- const double ctau = 0.0;
- double xstep = 0.0, zstep = 0.0;
- double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
- int nb_freq[2] = {0, 0};
- long nturn = Nbtour;
- bool status2 = true;
- 
- struct tm *newtime;
- 
-//variables of the returned the tracked particle
- long lastn = 0;
- long lastpos = 1;
- ss_vect<double> x1;
- 
- 
- //if(loss){
- char FmapLossFile[S_SIZE+5]=" ";
- strcpy(FmapLossFile,FmapFile);
- strcat(FmapLossFile,".loss");
- //}
- 
- /* Get time and date */
- time_t aclock;
- time(&aclock);                 /* Get time in seconds */
- newtime = localtime(&aclock);  /* Convert time to struct */
-
- if (trace) printf("Entering fmap ... results in %s\n\n",FmapFile);
-
- /* Opening file */
- if ((outf = fopen(FmapFile, "w")) == NULL) {
-   fprintf(stdout, "fmap: error while opening file %s\n", FmapFile);
-   exit_(1);
- }
-
-
- 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");
-//
-         fprintf(outf,"#    x[m]          z[m]           fx             fz"
-         "            dfx            dfz\n");
-
-
-//file with lost particle information
-if(printloss){ 
-
-  if ((outf_loss = fopen(FmapLossFile, "w")) == NULL) {
-      fprintf(stdout, "fmap: error while opening file %s\n", FmapLossFile);
-      exit_(1);
-   }
-
- fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
- fprintf(outf_loss,"# information of the lost particle: "
-                   "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
-
-        fprintf(outf_loss,"#   x[m]    z[m]   Nturn   lostp    S[m]       x[m]"
-         "       xp[rad]    z[m]      zp[rad]     delta      ctau\n");
- }
-         
- if ((Nbx < 1) || (Nbz < 1))
-   fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
- 
- // steps in both planes
- xstep = xmax/sqrt((double)Nbx);
- zstep = zmax/sqrt((double)Nbz);
-
- // double number of turn if diffusion to compute
- if (diffusion) nturn = 2*Nbtour;
-
- // px and pz zeroed
- xp = xp0;
- zp = zp0;
-
-// Tracking part + NAFF 
- for (i = 0; i <= Nbx; i++) {
- x  = x0 + sqrt((double)i)*xstep;
-// for (i = -Nbx; i <= Nbx; i++) {
- //  x  = x0 + sgn(i)*sqrt((double)abs(i))*xstep;
-//   if (!matlab) fprintf(outf,"\n");
-   fprintf(stdout,"\n");
-   for (j = 0; j<= Nbz; j++) {
-   z  = z0 + sqrt((double)j)*zstep;
- //  for (j = -Nbz; j<= Nbz; j++) {
- //    z  = z0 + sgn(j)*sqrt((double)abs(j))*zstep;
-  
-   //print out the lost information
-   if(printloss)
-     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
-     else
-     // tracking around closed orbit
-     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
-
-   if (status2) { // if trajectory is stable
-       // gets frequency vectors
-       Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
-       Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
-       if (diffusion) { // diffusion
-         // shift data for second round NAFF
-         Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
-         // gets frequency vectors
-         Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq);
-         Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
-       }
-     } // unstable trajectory
-     else { //zeroing output
-      nux1 = 0.0; nuz1 = 0.0;
-      nux2 = 0.0; nuz2 = 0.0;
-     }
-     
-     // printout value
-     if (!diffusion){
-//       fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n",
-//             1e3*x, 1e3*z, nux1, nuz1);
-//       fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
-//             1e3*x, 1e3*z, nux1, nuz1);
-               fprintf(outf,"%10.6e %10.6e %10.6e %10.6e\n",
-               x, z, nux1, nuz1);
-       fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
-               x, z, nux1, nuz1);
-     }
-     else {
-       dfx = nux1 - nux2; dfz = nuz1 - nuz2;
-       fprintf(outf,"%10.6e %10.6e %10.6e %10.6e %10.6e %10.6e\n",
-               x, z, nux1, nuz1, dfx, dfz);
-       fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
-               x, z, nux1, nuz1, dfx, dfz);
-     }
-   
-     //print out the information of the lost particle
-     if(printloss)
-        fprintf(outf_loss, "%6.3e %6.3e %8ld %3d %+12.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n", 
-                            x, z, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1], 
-                            x1[2], x1[3], x1[4], x1[5]);
-
-   }
- }
-
- fclose(outf);
- if(printloss)
-   fclose(outf_loss);
-}
-#undef NTERM2
-
 
 /****************************************************************************/
@@ -1253 +1024 @@
 #define NTERM2  10
 void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax, 
-            double zmax, double energy, bool diffusion, bool printloss, int numprocs, int myid)   
+            double zmax, double energy, bool diffusion, int numprocs, int myid)   
 {
  FILE * outf;
- 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], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
+ double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
  double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
  double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
@@ -1267 +1037 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2 = true;
+ bool status = true;
  struct tm *newtime;
 
@@ -1274 +1044 @@
  strcat(FmapFile,FmapFile_p);
  printf("%s\n",FmapFile);
-
-//variables of the returned the tracked particle
- long lastn = 0;
- long lastpos = 1;
- ss_vect<double> x1;
-
- char FmapLossFile[S_SIZE+5]=" ";
- strcpy(FmapLossFile,FmapFile);
- strcat(FmapLossFile,".loss");
-
 
  /* Get time and date */
@@ -1299 +1059 @@
  }
 
-//file with lost particle information
- if(printloss){ 
-   if ((outf_loss = fopen(FmapLossFile, "w")) == NULL) {
-     fprintf(stdout, "fmap: error while opening file %s\n", FmapLossFile);
-     exit_(1);
-   }
- }
-
  if(myid==0)
  {
@@ -1312 +1064 @@
   fprintf(outf,"# nu = f(x) \n");
   fprintf(outf,"#    x[m]          z[m]           fx             fz            dfx            dfz\n");
- 
-  if(printloss){ 
-    fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
-    fprintf(outf_loss,"# information of the lost particle: "
-            "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
-    
-    fprintf(outf_loss,"#   x[m]    z[m]   Nturn   lostp    S[m]       x[m]"
-            "       xp[rad]    z[m]      zp[rad]     delta      ctau\n");
-  }
-}
+ }
  
  if ((Nbx < 1) || (Nbz < 1))
@@ -1370 +1113 @@
  
      // tracking around closed orbit
-     if(printloss)
-       Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
-     else
-       Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
-
-
-     if (status2) // if trajectory is stable
+     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
+     if (status) // if trajectory is stable
      { 
        // gets frequency vectors
@@ -1410 +1148 @@
      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", x, z, nux1, nuz1, dfx, dfz);
      }
-
- //print out the information of the lost particle
-     if(printloss)
-        fprintf(outf_loss, "%6.3e %6.3e %8ld %3d %+12.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n", 
-                            x, z, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1], 
-                            x1[2], x1[3], x1[4], x1[5]);
    }
  }
 
  fclose(outf);
-
- if(printloss)
-   fclose(outf_loss);
-
 }
 #undef NTERM2
@@ -1453 +1181 @@
        emax             maximum energy
        z                vertical amplitude
-       diffusion        flag to calculate tune diffusion/ tune 
-                        difference between the first Nbtour and last Nbtour
+       diffusion        flag to calculate tune diffusion
        matlab  set file print format for matlab post-process; specific for nsrl-ii
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -1474 +1201 @@
        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
-
-       18/06/2013   by Jianfeng Zhang @ LAL
-
-      Add bool flag to print out the last information      of the tracked particle
 ****************************************************************************/
 #define NTERM2  10
 void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
-              double z, bool diffusion, bool printloss)
+              double z, bool diffusion)
 {
  FILE * outf;
-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], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
+ double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
  double x = 0.0, xp = 0.0, zp = 0.0, dp = 0.0, ctau = 0.0;
  double x0 = 1e-6, xp0 = 0.0, zp0 = 0.0;
@@ -1495 +1217 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2=true;
+ bool status=true;
  struct tm *newtime;
-
- //variables of the returned the tracked particle
- long lastn = 0;
- long lastpos = 1;
- ss_vect<double> x1;
- 
- char FmapdpLossFile[S_SIZE+5]=" ";
- strcpy(FmapdpLossFile,FmapdpFile);
- strcat(FmapdpLossFile,".loss");
 
  /* Get time and date */
@@ -1514 +1227 @@
  if (diffusion && globval.Cavity_on == false) nturn = 2*Nbtour;
 
- if (trace) printf("Entering fmapdp ... results in %s\n\n",FmapdpFile);
+ if (trace) printf("Entering fmap ... results in %s\n\n",FmapdpFile);
 
  /* Opening file */
@@ -1527 +1240 @@
 fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx           dfz\n");
  
-
-//file with lost particle information
-if(printloss){ 
-
-  if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL) {
-      fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpLossFile);
-      exit_(1);
-   }
-
- fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
- fprintf(outf_loss,"# information of the lost particle: "
-                   "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
-
-        fprintf(outf_loss,"#  dp      x[m]    Nturn   lostp     S[m]       x[m]"
-         "      xp[rad]     z[m]     zp[rad]     delta      ctau\n");
- }
-
  if ((Nbx <= 1) || (Nbe <= 1))
    fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
@@ -1566 +1262 @@
         diffusion = false;
      }   
-     else{
+     else
       // x  = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
          x  = x0 + sqrt((double)j)*xstep;
-     }
-
-if(printloss)
-     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
-     else
-     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status2);
-
-     if (status2) {
+     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
+     if (status) {
        Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
        Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
@@ -1610 +1300 @@
         dp, x, nux1, nuz2, dfx, dfz);
      }
-   
-     if(printloss)
-       fprintf(outf_loss," %-+8.3f %-+8.3f %-8ld %2d %+12.3f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n", 
-               dp, x, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
-
    }
  }
 
  fclose(outf);
-
-if(printloss)
-   fclose(outf_loss);
 }
 #undef NTERM2
@@ -1656 +1338 @@
 
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -1671 +1353 @@
        14/11/2011  add features to parallel calculate fmapdp. 
                    Merged with the version written by Mao-sen Qiu at Taiwan light source.
-                   
-       18/06/2013   by Jianfeng Zhang @ LAL
-        Add bool flag to print out the last information of the tracked particle            
 ****************************************************************************/
 #define NTERM2  10
 void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double xmax, 
-              double emax, double z, bool diffusion, bool printloss, int numprocs, int myid)
+              double emax, double z, bool diffusion, int numprocs, int myid)
 {
  FILE * outf;
-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], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
+ double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
  double x = 0.0, xp = 0.0, zp = 0.0, dp = 0.0, ctau = 0.0;
  double x0 = 1e-6, xp0 = 0.0, zp0 = 0.0;
@@ -1691 +1369 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2=true;
+ bool status=true;
  struct tm *newtime;
 
@@ -1699 +1377 @@
  printf("%s\n",FmapdpFile);
 
-//variables of the returned the tracked particle
- long lastn = 0;
- long lastpos = 1;
- ss_vect<double> x1;
-
- char FmapdpLossFile[S_SIZE+5]=" ";
- strcpy(FmapdpLossFile,FmapdpFile);
- strcat(FmapdpLossFile,".loss");
-
  /* Get time and date */
  time_t aclock;
@@ -1721 +1390 @@
    fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpFile);
    exit_(1);
- }
-
-if(printloss){ 
-  if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL) {
-      fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpLossFile);
-      exit_(1);
-   }
  }
 
@@ -1736 +1398 @@
      // fprintf(outf,"#    dp[%%]         x[mm]          fx            fz           dfx           dfz\n");
      fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx           dfz\n");
- 
-
-     if(printloss){
-       fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
-       fprintf(outf_loss,"# information of the lost particle: "
-               "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
-
-       fprintf(outf_loss,"#  dp      x[m]    Nturn   lostp     S[m]       x[m]"
-               "      xp[rad]     z[m]     zp[rad]     delta      ctau\n");
-     }
-  }
+   }
  
 if ((Nbx <= 1) || (Nbe <= 1))
@@ -1796 +1448 @@
         diffusion = false;
      }   
-     else{
+     else
       // x  = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
          x  = x0 + sqrt((double)j)*xstep;
-     }
-
-     if(printloss)
-       Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
-     else
-       Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status2);
-
-     if (status2) {
+     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
+     if (status) {
        Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
        Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
@@ -1840 +1486 @@
         dp, x, nux1, nuz2, dfx, dfz);
      }
-
-     if(printloss)
-       fprintf(outf_loss," %-+8.3f %-+8.3f %-8ld %2d %+12.3f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n", 
-               dp, x, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
-
    }
  }
 
  fclose(outf);
-
- if(printloss)
-   fclose(outf_loss);
 }
 #undef NTERM2
@@ -1897 +1535 @@
   double nux1 = 0.0, nuz1 = 0.0;
   int nb_freq[2] = {0, 0};
-  bool status2 = true;
+  bool status = true;
   struct tm *newtime;
 
@@ -1929 +1567 @@
     ctau = ctau0;
     
-    Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,Nbtour,Tab,&status2); // tracking around closed orbit
-    if (status2) {
+    Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,Nbtour,Tab,&status); // tracking around closed orbit
+    if (status) {
        Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq); // get frequency vectors
        Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
@@ -1936 +1574 @@
     else {
        nux1 = 0.0; nuz1 = 0.0;
-       status2 = true;
+       status = true;
     }
     
@@ -1991 +1629 @@
   FILE *outf;
   const char *fic = phasefile; 
-  int i=0;
-  bool status2=false;
+  int i;
+  bool status;
   struct tm *newtime;
 
@@ -2025 +1663 @@
   }
   
-  Trac_Simple6DCOD(x,xp,y,yp,energy,ctau,Nbtour,Tab,&status2);
-
+  Trac_Simple6DCOD(x,xp,y,yp,energy,ctau,Nbtour,Tab,&status);
   for (i = 0; i < Nbtour; i++) {
     fprintf(outf,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
             Tab[0][i],Tab[1][i],Tab[2][i],Tab[3][i],Tab[4][i],Tab[5][i]);
   }
-  //  cout << "status is: " << status2 << endl;
   fclose(outf);
 }
@@ -2069 +1705 @@
   FILE *outf;
   const char  *fic="phasepoly.out";
-  long        lastpos = 0L,lastn = 0L;
-  int         i=0,j=0;
-  double      x=0.0, z=0.0, px=0.0, pz=0.0, delta=0.0, ctau=0.0;
+  long        lastpos = 0,lastn = 0;
+  int         i,j;
+  double      x, z, px, pz, delta, ctau;
   double      ex = 1368E-9, el = 1.78E-4;
   double      betax = 9.0, /*betaz = 8.2, */betal = 45.5;
@@ -2163 +1799 @@
   double start = 0.0, step = 0.0;
   double x = 0.0, px = 0.0, z = 0.0, pz = 0.0, delta = 0.0, ctau = 0.0;
-  bool status2 = true;
+  bool status = true;
   struct tm *newtime;
 
@@ -2222 +1858 @@
    fprintf(stdout,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
             x,px,z,pz,delta,ctau);
-    Trac_Simple4DCOD(x,px,z,pz,delta,ctau,Nbtour,Tab,&status2);
+    Trac_Simple4DCOD(x,px,z,pz,delta,ctau,Nbtour,Tab,&status);
    for (i = 0; i < Nbtour; i++) {
       fprintf(outf,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
@@ -2265 +1901 @@
   FILE *outf;
   const char fic[] = "check_ampl.out";
-  int i=0;
+  int i;
 
   if ((outf = fopen(fic, "w")) == NULL)
@@ -2320 +1956 @@
   FILE *outf;
   const char fic[] = "enveloppe.out";
-  int i=0,j=0 ;
+  int i,j ;
   CellType Cell;
 
@@ -3562 +3198 @@
   FILE *fi;
   const char fic_skew[] = "QT-solamor_2_3.dat";
-  int i=0;
+  int i;
   double theta[500]; /* array for skew quad tilt*/
-  double b2=0.0, mKL=0.0;
+  double b2, mKL;
   CellType Cell;
-  long mOrder=0L;
+  long mOrder;
 
   int nquad = 0;  /* Number of skew quadrupoles */
@@ -3767 +3403 @@
   struct tm     *newtime;  // for time
   Vector        codvector[Cell_nLocMax];
-  bool          cavityflag=false, radiationflag=false;
+  bool          cavityflag, radiationflag;
   bool          trace=true; 
   
@@ -4165 +3801 @@
   struct tm     *newtime;  // for time
   Vector        codvector[Cell_nLocMax];
-  bool          cavityflag=false, radiationflag=false;
+  bool          cavityflag, radiationflag;
   bool          trace=true; 
   
@@ -4650 +4286 @@
  const char xfic[] = "xspectrum.out";
  const char zfic[] = "zspectrum.out";
- long i=0L, j=0L, k=0L;
+ long i, j, k;
  #define nterm2  20
  double Tab[6][NTURN], fx[nterm2], fz[nterm2], fx2[nterm2], fz2[nterm2];
@@ -4658 +4294 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2=true;
+ bool status=true;
  struct tm *newtime;
 
@@ -4699 +4335 @@
    for (j = 0; j<= Nbz; j++) {
      z  = z0 + sqrt((double)j)*zstep;
-     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
-     if (status2) {
+     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+     if (status) {
       Get_NAFF(nterm2, Nbtour, Tab, fx, fz, nb_freq);
      }
@@ -4779 +4415 @@
             long nmax, long pos, long &lastn, long &lastpos, FILE *outf1)
 {
-  bool lostF=false; /* Lost particle Flag */
+  bool lostF; /* Lost particle Flag */
   Vector x1;     /* tracking coordinates */
   Vector2  aperture;
@@ -4792 +4428 @@
   getelem(pos-1,&Cell);
 
-  if (trace) printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
+  if (!trace) printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
              dp*1e2, Cell.BeamPos[0]*1e3, Cell.BeamPos[2]*1e3);
 
@@ -4863 +4499 @@
   CellType Cell;
   int qlist[320];
-  int nquad=0, i=0;
+  int nquad=0, i;
   double A = 0.0;
 
@@ -4897 +4533 @@
 /****************************************************************************/
 /* void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-   double energy, bool *status2)
+   double energy, bool *status)
 
    Purpose:
@@ -4917 +4553 @@
 
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -4939 +4575 @@
  FILE * outf;
  const char fic[] = "fmapfull.out";
- int i=0, j=0, k=0;
+ int i, j, k;
  double Tab[DIM][NTURN], Tab0[DIM][NTURN];
  double fx[NTERM], fz[NTERM], fx2[NTERM], fz2[NTERM];
@@ -4948 +4584 @@
  double nux1[NTERM], nuz1[NTERM],nux2[NTERM], nuz2[NTERM];
  long nturn = Nbtour;
- bool status2=true;
+ bool status=true;
  struct tm *newtime;
  char name[14];
@@ -5008 +4644 @@
    for (j = 0; j<= Nbz; j++) {
      z  = z0 + sqrt((double)j)*zstep;
-     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
-
-     if (status2) {
+     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+
+     if (status) {
        Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq);
 
@@ -5091 +4727 @@
 /****************************************************************************/
 /* void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-   double energy, bool *status2)
+   double energy, bool *status)
 
    Purpose:
@@ -5111 +4747 @@
 
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -5133 +4769 @@
   FILE * outf;
   const char fic[] = "dyna.out";
-  long i=0, j=0;
+  long i, j;
   double Tab[6][NTURN], fx[NTERM2], fz[NTERM2];
   double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
@@ -5140 +4776 @@
   int nb_freq[2] = {0, 0};
   long nturn = Nbtour;
-  bool status2=true;
+  bool status=true;
   struct tm *newtime;
 
@@ -5174 +4810 @@
     for (j = 0; j<= Nbz; j++) {
       z  = z0 + sqrt((double)j)*zstep;
-      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
-      if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+      if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
       else {
        fx[0] = 0.0; fz[0] = 0.0;
       }
-      fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
-      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
+      fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
+      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
       if (diffusion) {
-        if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
-        fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
-        fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
+        if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+        fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
+        fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
       }
     }
@@ -5196 +4832 @@
     for (j = 0; j<= Nbz; j++) {
       z  = z0 + sqrt((double)j)*zstep;
-      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
-      if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+      if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
       else {
        fx[0] = 0.0; fz[0] =0.0;
       }
-      fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
-      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
+      fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
+      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
       if (diffusion) {
-        if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
+        if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
         fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", x, z, fx[0], fz[0]);
         fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", x, z, fx[0], fz[0]);
@@ -5247 +4883 @@
   FILE *outf;
   const char fic[] = "phase2.out";
-  long lastpos = 0L,lastn = 0L;
+  long lastpos = 0,lastn = 0;
   struct tm *newtime;
 
@@ -5373 +5009 @@
   int i,j;
   bool chroma=true, trace=false;
-  bool radiationflag=false, cavityflag=false;
+  bool radiationflag, cavityflag;
   double dP      = 0.0;
   double diffcmu = 0.0,   /* cos(mu1) - cos(mu2)*/
@@ -5795 +5431 @@
 void PhaseLongitudinalHamiltonien(void)
 {
-  long i=0L,j=0L;
+  long i,j;
   const double t = T;        // To get a one turn map
-  double phi=0.0, delta=0.0, H0=0.0;
+  double phi, delta, H0;
   long imax = 1000L,         // turn number 
        jmax = 25L;          // starting condition number
@@ -5993 +5629 @@
 }
 
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 double EnergySmall(double *X, double irho)
 {
- double A=0.0, B=0.0;
+ double A, B;
  double h = irho;
 
@@ -6008 +5639 @@
  return (A+B);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 double EnergyDrift(double *X)
 {
- double A=0.0;
+ double A;
 
  A = (X[1]*X[1]+X[3]*X[3])/2.0/(1.0+X[4]);
@@ -6055 +5681 @@
   FILE *outf;
   const char fic[] = "enveloppe2.out";
-  int i=0,j=0 ;
+  int i,j ;
   CellType Cell;
   /* Array for Enveloppes */
@@ -6185 +5811 @@
 void set_RFVoltage(const int Fnum, const double V_RF){
 
-  int k=0, n = 0;
+  int k, n = 0;
   
   
@@ -6232 +5858 @@
        01/2011  Fix the bug for reading the end of line symbol "\n" , "\r",'\r\n' 
                 at different operation system
-       04/2011  Change the set of 'seed' for rms error in file, now it's mandatory.
-       10/2013 Jianfeng Zhang @ LAL
-               Fix the bug of fgets() to read lines of arbritrary length.
+       04/2011  Change the set of 'seed' for rms error in file, now it's mandatory.     
 *****************************************************************************************************/
 void ReadFieldErr(const char *FieldErrorFile) 
 {  
-  bool  rms=false, set_rnd = false;
-  char     name[max_str],keywrd[max_str], *prm;
-  char    *line = NULL;
-  size_t   len = 0;
-  ssize_t   read;
+  bool  rms, set_rnd = false;
+  char    in[max_str], name[max_str],keywrd[max_str], *prm;
+  char    *line;
   int     n = 0;    /* field error order*/
   int     LineNum = 0;
-  int     seed_val = 0; // random seed number for the rms error 
+  int     seed_val; // random seed number for the rms error 
   double  Bn = 0.0, An = 0.0, r0 = 0.0; /* field error components and radius when the field error is measured */
   /* conversion number from A to T.m for soleil*/
@@ -6254 +5876 @@
 
   inf = file_read(FieldErrorFile);
-  if(inf == NULL){
-    printf("ReadFieldErr(): Error! Failure to read file %s !\n",FieldErrorFile);
-    exit_(1);
-  }
-  
+
   printf("\n");
   /* read lines*/
-  while ((read=getline(&line, &len, inf)) != -1) {  
+  while (line=fgets(in, max_str, inf)) {
+  
+  /* kill preceding whitespace generated by "table" key
+        or "space" key, but leave \n 
+        so we're guaranteed to have something*/
+     while(*line == ' ' || *line == '\t') {
+       line++;
+     }    
     
     /* count line number for debug*/
     LineNum++;
-    if(prt){
-      printf("Line # %ld \n",LineNum);
-      printf("Retrieved line of length %zu : \n",read);
-      printf("%s\n",line);
-    }
     
     /* check the line is whether comment line or null line*/
@@ -6287 +5907 @@
         /*read and assign the key words and measure radius*/
           sscanf(line, " %*s %s %lf",keywrd, &r0);
-          if (!prt) printf("\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,r0);
+          if (prt) printf("\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,r0);
           
           rms = (strcmp("rms", keywrd) == 0)? true : false;
@@ -6311 +5931 @@
             sscanf(prm, "%lf", &An);
           
-            if (!prt)
+            if (prt)
               printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
           
@@ -6335 +5955 @@
      // printf("%s", line);
   }
-  free(line);
+
   fclose(inf);
 }
@@ -6393 +6013 @@
       AddFieldValues_fam(Fnum,keywrd, r0, n, Bn, An);
     else 
-      printf("AddFieldErrors(): Error! undefined lattice element %s !\n", name);
+      printf("SetFieldErrors: undefined element %s\n", name);
   }
 }
@@ -6674 +6294 @@
          which also functions as these correctors.  
                  
-       10/2010  Written by Jianfeng Zhang @ SOLEIL
+       10/2010  Written by Jianfeng Zhang
 **********************************************************************/
 void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const char *keywrd, const double r0,
@@ -6682 +6302 @@
   double  bnL = 0.0, anL = 0.0;
   double  brho = 0.0, conv_strength = 0.0;
-  double  corr = 0.0;   /* skew quadrupole horizontal or vertical corrector error, read from a file*/
+  double  corr;   /* skew quadrupole horizontal or vertical corrector error, read from a file*/
   int    corrlistThick[120];   /* index of associated sextupole*/
   
@@ -6798 +6418 @@
 void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy)
 {
-  long      i=0L;
+  long      i;
   iVector2  nq1 = {0,0},nq2 = {0,0}, nq={0,0};
   Vector2   nu = {0.0, 0.0};
@@ -6847 +6467 @@
       delta            initial delta relative to closed orbit
       ctau             initial c*tau relative to closed orbit
-      nmax             maximum number of tracking turns
+      nmax             maximum number of turns
       
       
@@ -6859 +6479 @@
 {
    
-    long int i=0L, pos = 1L;
-    long int lastn = 0L, lastpos = 0L;
+    long int i,pos = 1;
+    long int lastn = 0, lastpos = 0;
     Vector x0, x1, x2, xf,codvector[Cell_nLocMax];
     FILE *outf;
@@ -6905 +6525 @@
     do {
         pos = 1;//track from first element
-        (lastn)++; //number of turns
+        (lastn)++;
         for (i = 0; i < nv_; i++)  //nv_ = 6
             x1[i] = x2[i];
 
-        //tracking for one turn
         while( pos <= globval.Cell_nLoc){  
          
@@ -6933 +6552 @@
                      pos, lastpos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_], 
                      1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
-          }
-           fprintf(outf,"i      name     S    x(mm)      px(mrad)     y(mm)     py(mrad)    delta(%)    z(mm)    Nturn\n");
-          fprintf(outf,"%6d  %s   %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
+          }     
+          fprintf(outf,"%6d  %s %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
                      pos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_], 
                      1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);            
@@ -6942 +6560 @@
         }//finish of tracking and printing to file  
         
-    } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc)); //track for nmax turns
-
-//     if (globval.MatMeth)
-//         Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
-
-    fclose(outf);   
-}
-/**********************************************************************
-void PrintTrackElem(const char *TrackFile, double x, double px, double y,double py, 
-           double delta, double ctau, long int nelem1, long int nelem2)
-           
-  Purpose:         
-    Print the coordinates tracked around a lattice element by tracking around COD
-
-  Input:
-      TrackFile        file to be print
-      x                initial x relative to closed orbit
-      px               initial px relative to closed orbit
-      y                initial y relative to closed orbit
-      py               initial py relative to closed orbit
-      delta            initial delta relative to closed orbit
-      ctau             initial c*tau relative to closed orbit
-      nelem1           start lattice index of the tracked element
-      nelem2           end lattice index of the tracked element
-      
-      
-  Output:
-  
-  Comments:
-    Written by Jianfeng Zhang  @ LAL 04/2013    
-**********************************************************************/ 
-void PrintTrackElem(const char *TrackFile, double x, double px, double y,double py, 
-           double delta, double ctau, long int nelem1, long int nelem2)
-{
-   
-    long int i=0L, pos = 1L;
-    Vector x0, x1;
-    FILE *outf;
-    struct    tm *newtime;
-    
-    bool prt = true;
-           
-    outf = file_write(TrackFile);
-    /* Get time and date */
-    newtime = GetTime();
-
-    fprintf(outf, "# Element tracking using TRACY III-- %s -- %s\n",TrackFile, asctime2(newtime));
-    fprintf(outf, "#\n");  
-   // fprintf(outf, "# work tunes: %7.5f %7.5f\n",globval.TotalTune[0], globval.TotalTune[1]);                
-   // fprintf(outf, "#    i  ElemName            S           x            p_x           y          p_y");
-   // fprintf(outf, "         delta         cdt     NElem \n");
-   // fprintf(outf, "#                          [m]        [mm]         [mrad]        [mm]        [mrad]");      
-   // fprintf(outf, "        [%%]          [mm]\n");
-
-      fprintf(outf,"i      ElemName     S[m]    x(m)      px(rad)     y(m)     py(rad)    delta(-)    z(m)   \n");
-    //initial coordinates
-    x0[x_] = x;
-    x0[px_] = px;
-    x0[y_] = y;
-    x0[py_] = py;
-    x0[delta_] = delta;
-    x0[ct_] = ctau;          
-   
-     for (i = 0; i < nv_; i++)  //nv_ = 6
-            x1[i] = x0[i];
-
- //print the coordinates at the element   
-        pos = nelem1+1;//first element "1" of Tracy is a default "debut", it's diffirent from the real lattice index
-      
-        //tracking from the start element to the end element
-        while( pos <= nelem2+1){  
-         
-          Elem_Pass(pos,x1);            
-        
-           
-          if (prt) {
-            printf("%4ld %s %6.4f %16.8f %16.8f %16.8f %16.8f %16.8f %16.8f \n",
-                     pos-1,Cell[pos].Elem.PName,Cell[pos].S, x1[x_], x1[px_], 
-                     x1[y_], x1[py_],x1[delta_], x1[ct_]);
-          }
-         
-          fprintf(outf,"%6d  %s   %8.4e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
-                     pos-1,Cell[pos].Elem.PName,Cell[pos].S, x1[x_], x1[px_], 
-                     x1[y_], x1[py_],x1[delta_], x1[ct_]);          
-        
-        pos++;               
-        }//finish of tracking and printing to file  
-        
- 
+    } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc));
 
 //     if (globval.MatMeth)