Index: /branches/tracy3-3.10.1b/tracy/tools/soltracy.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tools/soltracy.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tools/soltracy.cc	(revision 23)
@@ -1,22 +1,25 @@
-/************************************
+/*
  Current revision $Revision$
  On branch $Name$
  Latest change $Date$ by $Author$
-*************************************/
-#define ORDER 1    
-//#define ORDER 4   
+*/
+#define ORDER 1          
 
 int no_tps = ORDER; // arbitrary TPSA order is defined locally
 
+
+
 extern bool freq_map;
 
-// #if HAVE_CONFIG_H
-// #include <config.h>
-// #endif
-
-// #if MPI_EXEC
-//   //library of parallel computation,must be before "stdio.h"
-//   #include <mpi.h>
-// #endif
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/*
+#if MPI_EXEC
+  //library of parallel computation,must be before "stdio.h"
+  #include <mpi.h>
+#endif
+*/
 
   #include "tracy_lib.h"
@@ -41,10 +44,4 @@
   globval.H_exact = false;
 
-  
-  //output files
-  char fic_twiss[S_SIZE + 4] = "";   //twiss file
-  char fic_summary[S_SIZE + 4]="";   //summary file
-  
-  
   /* parameters to read the user input script .prm */
   long i=0L; //initialize the for loop to read command string
@@ -64,8 +61,11 @@
   UserCommand UserCommandFlag[NCOMMAND];
  
-// #if MPI_EXEC
-//   //Initialize parallel computation
-//     MPI_Init(&argc, &argv);
-// #endif
+  /*
+#if MPI_EXEC
+  //Initialize parallel computation
+    MPI_Init(&argc, &argv);
+#endif
+  */
+
   /************************************************************************
    read in files and flags
@@ -76,5 +76,5 @@
     read_script(argv[1], true,CommandNo, UserCommandFlag);
   } else {
-    fprintf(stdout, "Not enough parameters\n Syntax is program parameter file\n");
+    fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
     exit_(1);
   }
@@ -105,20 +105,9 @@
   else
     globval.ge = ElemIndex(ge_name);
- 
- 
+  
 //  globval.g = ElemIndex("g");  /* get family index of  girder*/
- strcpy(fic_twiss,lat_file);
-  strcpy(fic_summary,lat_file);
-  strcat(fic_twiss,".twi"); 
-  strcat(fic_summary,".sum");
-  
-    //generate the twiss and summary files with proper appendix  
-  /* print the summary of the element in lattice to the screen and an external file*/  
+   
+  /* print the summary of the element in lattice */  
   printglob();   
-  printglob2file(fic_summary);
- //print the twiss file
-   /* print out lattice functions, with all the optical information for the lattice with design values */
-  printlatt(fic_twiss);  
-  
   /************************************************************************
     print files, very important file for debug
@@ -263,15 +252,4 @@
    }
    
-   else if(strcmp(CommandStr,"PrintTrackElemFlag") == 0) {
-      cout << "\n";
-      cout << "print the tracked coordinates to file: "
-           << UserCommandFlag[i]._PrintTrackElem_track_file << "\n";
-      	
-      PrintTrackElem(UserCommandFlag[i]._PrintTrackElem_track_file,
-                 UserCommandFlag[i]._PrintTrackElem_x, UserCommandFlag[i]._PrintTrackElem_px, 
-                 UserCommandFlag[i]._PrintTrackElem_y, UserCommandFlag[i]._PrintTrackElem_py, 
-                 UserCommandFlag[i]._PrintTrackElem_delta,UserCommandFlag[i]._PrintTrackElem_ctau,
-		 UserCommandFlag[i]._PrintTrackElem_nelem1,UserCommandFlag[i]._PrintTrackElem_nelem2);  
-   }
     //print the girder
   // else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
@@ -468,6 +446,7 @@
   else if(strcmp(CommandStr,"FmapFlag") == 0) {
     printf("\n begin Fmap calculation for on momentum particles: \n");
-    
-  //  #if MPI_EXEC
+  
+  //   /*  
+ //   #if MPI_EXEC
 
  //    //initialization for parallel computing
@@ -490,5 +469,4 @@
  // 	   UserCommandFlag[i]._FmapFlag_delta, 
  // 	   UserCommandFlag[i]._FmapFlag_diffusion,
- //        UserCommandFlag[i]._FmapFlag_printloss,
  // 	   numprocs,myid);
 
@@ -543,7 +521,7 @@
 	   UserCommandFlag[i]._FmapFlag_ymax, 
 	   UserCommandFlag[i]._FmapFlag_delta, 
-	   UserCommandFlag[i]._FmapFlag_diffusion,
-	   UserCommandFlag[i]._FmapFlag_printloss);
+	   UserCommandFlag[i]._FmapFlag_diffusion);
       //#endif
+    
   }
 
@@ -551,5 +529,6 @@
   else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
     printf("\n begin Fmap calculation for off momentum particles: \n");
-    
+   
+  //   /* 
  //   #if MPI_EXEC
 
@@ -572,5 +551,4 @@
  // 	     UserCommandFlag[i]._FmapdpFlag_z,
  // 	     UserCommandFlag[i]._FmapdpFlag_diffusion,
- //            UserCommandFlag[i]._FmapdpFlag_printloss,
  // 	     numprocs,myid);
 
@@ -625,7 +603,7 @@
 	      UserCommandFlag[i]._FmapdpFlag_emax, 
 	      UserCommandFlag[i]._FmapdpFlag_z,
-	      UserCommandFlag[i]._FmapdpFlag_diffusion,
-	       UserCommandFlag[i]._FmapdpFlag_printloss);
-       //  #endif
+	      UserCommandFlag[i]._FmapdpFlag_diffusion);
+       // #endif
+    
   }
 
@@ -657,5 +635,6 @@
     printf("\n Calculate momentum acceptance: \n");
 
-//  #if MPI_EXEC
+   
+ // #if MPI_EXEC
     
 //     /* calculate momentum acceptance*/
@@ -789,4 +768,5 @@
 
 // #else
+
      MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
 			UserCommandFlag[i]._MomentumAccFlag_istart, 
@@ -800,5 +780,7 @@
 			UserCommandFlag[i]._MomentumAccFlag_nturn,
 			UserCommandFlag[i]._MomentumAccFlag_zmax);
-     //  #endif   
+     /*
+  #endif   
+     */
 
     /* restore the initial values*/
@@ -863,5 +845,5 @@
 	  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");
+    printf("the simulation time for phase space in tracy 3 is \n");
     stop = stampstop(start);
 
@@ -1021,8 +1003,10 @@
 
  
- // #if MPI_EXEC
- //  MPI_Finalize();   //finish parallel computation
- // #endif
-  
+  /*
+ #if MPI_EXEC
+  MPI_Finalize();   //finish parallel computation
+ #endif
+  */
+
   return 0;
 }//end of main() 
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/naffutils.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/naffutils.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/naffutils.h	(revision 23)
@@ -27,7 +27,4 @@
 		 double ctau, long nmax, double Tx[][NTURN], bool *status);
 
-void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
-                 double ctau, long nmax, double Tx[][NTURN], long& lastn, long& lastpos,ss_vect<double>& x1,bool *status2);
-
 void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
          double ctau, long nmax, double Tx[][NTURN], bool *status);
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/physlib.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/physlib.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/physlib.h	(revision 23)
@@ -55,6 +55,4 @@
 
 void printglob(void);
-
-void printglob2file(const char fic[]); 
 
 void printlatt(const char fic[]);
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/read_script.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/read_script.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/read_script.h	(revision 23)
@@ -37,11 +37,4 @@
    long  _PrintTrack_nmax; 
    
-   //start track particle coordinates; PrintTrackFlag
-   char _PrintTrackElem_track_file[max_str];
-   double _PrintTrackElem_x, _PrintTrackElem_px, _PrintTrackElem_y, 
-          _PrintTrackElem_py, _PrintTrackElem_delta, _PrintTrackElem_ctau;
-   long  _PrintTrackElem_nelem1, _PrintTrackElem_nelem2; 
-   
-   
    //twiss file 
    char _PrintTwiss_twiss_file[max_str];
@@ -65,7 +58,5 @@
    long _FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn;
    double _FmapFlag_xmax, _FmapFlag_ymax, _FmapFlag_delta;
-   bool _FmapFlag_diffusion;
-   bool _FmapFlag_printloss;
-   
+   bool _FmapFlag_diffusion; 
  
  //extern bool FmapdpFlag;
@@ -73,6 +64,5 @@
    long _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn;
    double _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z;
-   bool _FmapdpFlag_diffusion;
-   bool _FmapdpFlag_printloss;
+   bool _FmapdpFlag_diffusion;	
 
    
@@ -120,5 +110,5 @@
    _FmapFlag_nxpoint=31L, _FmapFlag_nypoint=21L, _FmapFlag_nturn=516L;
    _FmapFlag_xmax=0.025, _FmapFlag_ymax=0.005, _FmapFlag_delta=0.0;
-   _FmapFlag_diffusion = true, _FmapFlag_printloss = false;
+   _FmapFlag_diffusion = true;
 
 /*fmap for off momentum particle*/
@@ -126,6 +116,6 @@
  _FmapdpFlag_nxpoint=31L, _FmapdpFlag_nepoint=21L, _FmapdpFlag_nturn=516L;
  _FmapdpFlag_xmax=0.025, _FmapdpFlag_emax=0.005, _FmapdpFlag_z=0.0;
- _FmapdpFlag_diffusion = true, _FmapdpFlag_printloss = false;
- 
+ _FmapdpFlag_diffusion = true;			
+
 /* tune shift with amplitude*/
  strcpy(_AmplitudeTuneShift_nudx_file,"nudx.out");
@@ -172,7 +162,5 @@
 
 /***** file names************/
-extern char    lat_file[max_str];
-
-//files with multipole errors; for soleil lattice
+ //files with multipole errors; for soleil lattice
  extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
  extern char multipole_file[max_str];
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/soleillib.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/soleillib.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/soleillib.h	(revision 23)
@@ -58,13 +58,11 @@
 //              double z, bool diffusion, bool matlab);
 void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion);	
-void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion, bool printloss);	
+          double energy, bool diffusion);	  
 void fmap_p(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion, bool printloss, int numprocs, int myid);	  
+          double energy, bool diffusion, int numprocs, int myid);	  
 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);
 void fmapdp_p(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax, 
-            double z, bool diffusion, bool printloss, int numprocs, int myid);
+            double z, bool diffusion, int numprocs, int myid);
 void Nu_Naff(void);
 
@@ -130,10 +128,8 @@
 /* fit tunes for soleil lattice, in which each quadrupole is cut into two parts*/			
 void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy);			 
-//print the coordinates at lattice elements tracked for n turns
+//print the coordinates at lattice elements
 void PrintTrack(const char *TrackFile, double x, double px,double y,double py, 
                 double delta, double ctau, long int nmax);
-//print the tracked coordinates after a lattice element
-void PrintTrackElem(const char *TrackFile, double x, double px,double y,double py, 
-                double delta, double ctau, long int nelem1,long int nelem2);
+		
 
 
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/t2elem.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/t2elem.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/t2elem.h	(revision 23)
@@ -38,7 +38,4 @@
 
 template<typename T>
-void Drift(double L,double h_bend, ss_vect<T> &x);
-
-template<typename T>
 void Drift(double L, ss_vect<T> &x);
 
@@ -47,11 +44,5 @@
 
 template<typename T>
-void bend_fringe(double hb, double gap, ss_vect<T> &x);
-
-template<typename T>
-static void bendCurvature(double irho, double H, ss_vect<T> &x); 
-
-template<typename T>
-static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x, bool Entrance);
+static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x);
 
 template<typename T>
@@ -61,9 +52,4 @@
 template<typename T>
 void Drift_Pass(CellType &Cell, ss_vect<T> &x);
-
-template<typename T>
-void dipole_kick(double L, double h_ref, double h_bend, ss_vect<T> &x);
-template<typename T>
-void multipole_kick(int Order, double MB[], double L, double h_bend, ss_vect<T> &x);
 
 template<typename T>
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/tpsa_lin.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/tpsa_lin.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/tpsa_lin.h	(revision 23)
@@ -75,19 +75,11 @@
 void dacos(const tps_buf &x, tps_buf &z);
 
-void datan(const tps_buf &x, tps_buf &z);
-
-void dactan(const tps_buf &x, tps_buf &z);
-
 void dasinh(const tps_buf &x, tps_buf &z);
 
 void dacosh(const tps_buf &x, tps_buf &z);
 
-void daarcsin(const tps_buf &x, tps_buf &z);
-
-void daarccos(const tps_buf &x, tps_buf &z);
+void datan(const tps_buf &x, tps_buf &z);
 
 void daarctan(const tps_buf &x, tps_buf &z);
-
-void daarcctan(const tps_buf &x, tps_buf &z);
 
 void dafun_(const char *fun, const tps_buf &x, tps_buf &z);
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/tracy.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/tracy.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/tracy.h	(revision 23)
@@ -110,7 +110,6 @@
   bool tuneflag, chromflag, codflag, mapflag, passflag, overflag, chambre;
   int lossplane; /* lost in: horizontal    1
-                             vertical      2
-                             longitudinal  3 
-                              not lost in any of the above plane: 0*/
+		                     vertical      2
+			                 longitudinal  3 */
 } statusrec;
 
Index: /branches/tracy3-3.10.1b/tracy/tracy/inc/tracy_global.h
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/inc/tracy_global.h	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/inc/tracy_global.h	(revision 23)
@@ -79,6 +79,4 @@
   int         Pmethod;   // Integration Method
   int         PN;        // Number of integration steps
-  long        dipEdge_effect1; // dipole edge effect at the entrance
-  long        dipEdge_effect2; // dipole edge effect at the exit
   long        quadFF1;         /* Entrance quadrupole Fringe field flag */
   long        quadFF2;         /* Exit quadrupole Fringe field flag */
@@ -100,9 +98,9 @@
   double      PdTrnd;    // (normal)random scale factor to rotation error PdTrms
   // Multipole strengths
-  mpolArray   PBpar;     // design field gradient; bn, and an
-  mpolArray   PBsys;     // systematic multipole errors gradient, bn and an
-  mpolArray   PBrms;     // rms multipole field errors gradient, bn and an
-  mpolArray   PBrnd;     // random scale factor of rms error PBrms gradient, bn and an
-  mpolArray   PB;        // total field strength(design,sys,rms) gradient, bn and an
+  mpolArray   PBpar;     // design field strength
+  mpolArray   PBsys;     // systematic multipole errors
+  mpolArray   PBrms;     // rms multipole field errors
+  mpolArray   PBrnd;     // random scale factor of rms error PBrms
+  mpolArray   PB;        // total field strength(design,sys,rms)
   int         Porder;    // The highest order in PB
   int         n_design;  // multipole order (design, All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6)
@@ -111,6 +109,4 @@
   double PTx1;           // horizontal entrance angle [deg]
   double PTx2;           // horizontal exit angle [deg]
-  double PH1;            // bending curvature of the entrance pole face of dipole, see P116 SAC-75. 
-  double PH2;            // bending curvature of the exit pole face of dipole, see P116 SAC-75. 
   double Pgap;           // total magnet gap [m] 
   double Pirho;          // angle([radian], but in lattice definition, angle is with unit degree)/length of the dipole, 1/rho [1/m] 
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/naffutils.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/naffutils.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/naffutils.cc	(revision 23)
@@ -18,11 +18,10 @@
 
 /****************************************************************************/
-/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
-                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
+/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, long nmax,
+                 double Tx[][NTURN], bool *status)
 
    Purpose:
        Single particle tracking around the closed orbit for NTURN turns
-       The 6D phase trajectory is saved in a array, but the 
-       tracked dp and ctau is not about the COD.
+       The 6D phase trajectory is saved in a array
 
    Input:
@@ -70,5 +69,5 @@
 
   
-  if (!trace && status.codflag) 
+  if (trace && status.codflag) 
     printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
              dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
@@ -108,5 +107,4 @@
   } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
 
-  
   if (lastpos != globval.Cell_nLoc)
   { /* Particle lost: Error message section */
@@ -118,122 +116,7 @@
 }
 
-
-/****************************************************************************/
-/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
-                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
-
-   Purpose:
-       Single particle tracking around the closed orbit for NTURN turns
-       The 6D phase trajectory is saved in a array, but the 
-       tracked dp and ctau is not about the COD.
-
-   Input:
-       x, px, y, py 4 transverse coordinates
-       dp           energy offset
-       nmax         number of turns
-       pos          starting position for tracking
-       aperture     global physical aperture
-
-   Output:
-      lastn         last n (should be nmax if  not lost)
-      lastpos       last position in the ring
-      Tx            6xNTURN matrix of phase trajectory
-
-   Return:
-       none
-
-   Global variables:
-       NTURN number of turn for tracking
-       globval
-
-   Specific functions:
-       Cell_Pass
-
-   Comments:
-       useful for connection with NAFF
-       19/01/03 tracking around the closed orbit
-       
-       11/06/2013  Modified by Jianfeng Zhang @ LAL 
-       The same feature as function  
-       Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
-                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
-       but add the feature to extract the position of the location of 
-       the lost particle; 
-       called by function 
-       fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-          double energy, bool diffusion, bool loss)	  
-       in soleillib.cc.
-       
-****************************************************************************/
-void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
-                 double ctau, long nmax, double Tx[][NTURN], long& lastn, long& lastpos, ss_vect<double>& x1, bool *status2)
-{
-  bool             lostF = false; /* Lost particle Flag */
-  Vector2          aperture = {1.0, 1.0};
-
-   lastn = 0;
-   lastpos = globval.Cell_nLoc;
-  *status2 = true; /* stable */
-   x1[0]=0,x1[1]=0,x1[2]=0,x1[3]=0,x1[4]=0,x1[5]=0;
-  
-  if (globval.MatMeth) Cell_Concat(dp);
-
-  /* Get closed orbit */
-
-  getcod(dp, lastpos);
-
-  
-  if (!trace && status.codflag) 
-    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
-             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
-
-  /* Tracking coordinates around the closed orbit */
-  x1[0] =  x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
-  x1[2] =  y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
-  x1[4] = dp; x1[5] = ctau;
-
-  Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
-  Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
-  Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
-  lastn++;
-
-  do
-  { /* tracking through the ring */
-    if ((lastpos == globval.Cell_nLoc) &&
-        (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
-    {
-      if (globval.MatMeth)
-        Cell_fPass(x1, lastpos);
-      else
-        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
-      Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
-      Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
-      Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
-    }
-    else
-    {
-      printf("Trac_Simple: Particle lost \n");
-      fprintf(stdout, "%6ld plane: %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", 
-         lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
-      lostF = true;
-      *status2 = false;
-    }
-    lastn++;
-  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
-
-  
-  if (lastpos != globval.Cell_nLoc)
-  { /* Particle lost: Error message section */
-    *status2 = false;
-    printf("Trac_Simple: Particle lost \n");
-    fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
-             status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
-  }
-}
-
-
 /****************************************************************************/
 /* void Trac_Simple6DCOD(double x, double px, double y, double py, double dp, long nmax,
-                 double Tx[][NTURN], bool *status2)
+                 double Tx[][NTURN], bool *status)
 
    Purpose:
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/physlib.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/physlib.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/physlib.cc	(revision 23)
@@ -16,10 +16,5 @@
 /**************************/
 
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 /**** same as asctime in C without the \n at the end****/
 char *asctime2(const struct tm *timeptr) {
@@ -38,10 +33,5 @@
     return result;
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 /** Get time and date **/
 struct tm* GetTime() {
@@ -199,6 +189,4 @@
 
  Comments copied from Tracy 2.7(soleil),Written by L.Nadolski.
- 
- 03/06/2013 Add feature to print the summary on the sreen and an external file
 
  ****************************************************************************/
@@ -263,103 +251,5 @@
 
 /****************************************************************************/
-/* void printglob2file(void)
-
- Purpose:
- Print global variables on screen
- Print tunes and chromaticities
- Print Oneturn matrix
-
- Input:
- none
-
- Output:
- output on the screen
-
- Return:
- none
-
- Global variables:
- globval
-
- Specific functions:
- none
-
- Comments: 
- 03/06/2013  by Jianfeng Zhang @ LAL
- The same feature as the file printglob(void), but 
- added feature to print the lattice summary to an external file
-
- 
- ****************************************************************************/
-void printglob2file(const char fic[]) {
-  FILE *fout;
-  int i=0,j=0;
-  
-  fout = fopen(fic,"w");
-  
-    fprintf(fout, "\n***************************************************************"
-        "***************\n");
-    fprintf(fout,"\n");
-    fprintf(fout,"  dPcommon     =  %9.3e  dPparticle   =  %9.3e"
-        "  Energy [GeV] = %.3f\n", globval.dPcommon, globval.dPparticle,
-            globval.Energy);
-    fprintf(fout,"  MaxAmplx [m] = %9.3e   MaxAmply [m] = %9.3e"
-        "   RFAccept [%%] = +/- %4.2f\n", Cell[0].maxampl[X_][1],
-            Cell[0].maxampl[Y_][1], globval.delta_RF * 1e2);
-    fprintf(fout,"  MatMeth      =  %s     ", globval.MatMeth ? "TRUE " : "FALSE");
-    fprintf(fout," Cavity_On    =  %s    ", globval.Cavity_on ? "TRUE " : "FALSE");
-    fprintf(fout,"  Radiation_On = %s     \n", globval.radiation ? "TRUE " : "FALSE");
-    
-    if(globval.bpm == 0)
-      fprintf(fout,"  bpm          =  0");  
-    else
-      fprintf(fout,"  bpm          =  %3d", GetnKid(globval.bpm));
-    if(globval.qt == 0)
-      fprintf(fout,"                             qt           = 0           \n"); 
-    else
-      fprintf(fout,"                             qt           = %3d         \n", GetnKid(globval.qt));    
-    if(globval.hcorr == 0)
-      fprintf(fout,"  hcorr         =  0");  
-    else
-      fprintf(fout,"  hcorr         = %3d", GetnKid(globval.hcorr));
-    if(globval.vcorr == 0)
-      fprintf(fout,"                             vcorr        =  0        \n");  
-    else
-      fprintf(fout,"                             vcorr        = %3d      \n", GetnKid(globval.vcorr));
-	    	       	       	    	    
-    fprintf(fout," Chambre_On   = %s     \n", globval.Aperture_on ? "TRUE " : "FALSE");
-  
-    fprintf(fout,"  alphac       =   %8.4e\n", globval.Alphac);
-    fprintf(fout,"  nux          =   %8.6f      nuy  =  %8.6f",
-            globval.TotalTune[X_], globval.TotalTune[Y_]);
-    if (globval.Cavity_on)
-        fprintf(fout,"  omega  = %13.9f\n", globval.Omega);
-    else {
-        fprintf(fout,"\n");
-        fprintf(fout,"  ksix         =    %8.6f      ksiy  =   %8.6f\n",
-                globval.Chrom[X_], globval.Chrom[Y_]);
-    }
-    fprintf(fout,"\n");
-    
-    fprintf(fout,"  OneTurn matrix:\n");
-    for(i=0;i<ss_dim;i++){
-      fprintf(fout,"\n");
-      for(j=0;j<ss_dim;j++)
-	fprintf(fout,"%14.6e",globval.OneTurnMat[i][j]);
-    }
-
-    fprintf(fout,"\n\nTwiss parameters at entrance:\n");
-    fprintf(fout,
-            "   Betax [m] = % 9.3e  Alphax = % 9.3e Etax [m] = % 9.3e Etaxp = % 9.3e\n"
-                "   Betay [m] = % 9.3e  Alphay = % 9.3e Etay [m] = % 9.3e Etayp = % 9.3e\n\n",
-            Cell[1].Beta[X_], Cell[1].Alpha[X_], Cell[1].Eta[X_],
-            Cell[1].Etap[X_], Cell[1].Beta[Y_], Cell[1].Alpha[Y_],
-            Cell[1].Eta[Y_], Cell[1].Etap[Y_]);
-
-    fclose(fout);
-}
-
-/****************************************************************************/
-/* void printlatt(const char fic[]) 
+/* void printlatt(void)
 
  Purpose:
@@ -423,10 +313,5 @@
     fclose(outf);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 double int_curly_H1(long int n) {
     /* Integration with Simpson's Rule */
@@ -445,10 +330,5 @@
     return curly_H;
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void prt_sigma(void) {
     long int i;
@@ -597,20 +477,10 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void getabn(Vector2 &alpha, Vector2 &beta, Vector2 &nu) {
     Vector2 gamma;
     Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void TraceABN(long i0, long i1, const Vector2 &alpha, const Vector2 &beta,
         const Vector2 &eta, const Vector2 &etap, const double dP) {
@@ -743,12 +613,7 @@
     Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, ksidkpL, ksiimax);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void FitDisp(long q, long pos, double eta) {
-    long i=0L, nq=0L;
+    long i, nq;
     fitvect qbuf;
 
@@ -760,10 +625,5 @@
     Ring_FitDisp(pos, eta, dispeps, nq, qbuf, dispdkL, dispimax);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 #define nfloq           4
 
@@ -775,10 +635,4 @@
 #undef nfloq
 
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
 #define ntrack          4
 
@@ -797,6 +651,6 @@
 
      */
-    long int i=0L;
-    double twoJx=0.0, twoJy=0.0, phix=0.0, phiy=0.0, scl_1 = 1.0, scl_2 = 1.0;
+    long int i;
+    double twoJx, twoJy, phix, phiy, scl_1 = 1.0, scl_2 = 1.0;
     Vector x0, x1, x2, xf;
     FILE *outf;
@@ -956,15 +810,9 @@
 #undef ntrack
 
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
 #define step            0.1
 #define px              0.0
 #define py              0.0
 void track_(double r, struct LOC_getdynap *LINK) {
-    long i=0L, lastn=0L, lastpos=0L;
+    long i, lastn, lastpos;
     Vector x;
 
@@ -1041,6 +889,5 @@
 void Trac(double x, double px, double y, double py, double dp, double ctau,
         long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) {
-  
-    bool lostF = true; /* Lost particle Flag */
+    bool lostF; /* Lost particle Flag */
     Vector x1; /* tracking coordinates */
     Vector2 aperture;
@@ -1060,4 +907,5 @@
     lastn = 0L;
     (lastpos) = pos;
+    lostF = true;
 
     if (trace)
@@ -1128,6 +976,5 @@
 #define nfloq     4
 bool chk_if_lost(double x0, double y0, double delta, long int nturn, bool floqs) {
-  
-    long int i=0L, lastn=0L, lastpos=0L;
+    long int i, lastn, lastpos;
     Vector x;
 
@@ -1259,6 +1106,6 @@
  ****************************************************************************/
 void getcsAscr(void) {
-    long i=0L, j=0L;
-    double phi=0.0;
+    long i, j;
+    double phi;
     Matrix R;
 
@@ -1331,6 +1178,6 @@
      Assumes mid-plane symmetry.                    */
 
-    long int i=0L, lastpos=0L;
-    double phi=0.0, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0;
+    long int i, lastpos;
+    double phi, x_min, x_max, y_min, y_max;
 
     getcod(delta, lastpos);
@@ -1422,6 +1269,6 @@
  ****************************************************************************/
 double get_aper(int n, double x[], double y[]) {
-    int i=0;
-    double A=0.0;
+    int i;
+    double A;
 
     A = 0.0;
@@ -1436,14 +1283,7 @@
 }
 
-/**************************************************************************
- * void GetTrack(const char *file_name, long *n, double x[], double px[],
-        double y[], double py[])
- * 
- * 
- * 
- **************************************************************************/
 void GetTrack(const char *file_name, long *n, double x[], double px[],
         double y[], double py[]) {
-    int k=0;
+    int k;
     char line[200];
     FILE *inf;
@@ -1495,5 +1335,5 @@
  ****************************************************************************/
 void Getj(long n, double *x, double *px, double *y, double *py) {
-    long int i=0L;
+    long int i;
 
     for (i = 0; i < n; i++) {
@@ -1529,6 +1369,5 @@
  ****************************************************************************/
 double GetArg(double x, double px, double nu) {
-    
-  double phi=0.0, val=0.0;
+    double phi, val;
 
     phi = GetAngle(x, px);
@@ -1568,5 +1407,5 @@
 void GetPhi(long n, double *x, double *px, double *y, double *py) {
     /* Calculates the linear phase */
-    long i=0L;
+    long i;
 
     for (i = 1; i <= n; i++) {
@@ -1582,5 +1421,5 @@
 void Sinfft(int n, double xr[]) {
     /* DFT with sine window */
-    int i=0;
+    int i;
     double xi[n];
 
@@ -1611,13 +1450,8 @@
     free_dvector(xi, 1, 2 * n);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void sin_FFT(int n, double xr[], double xi[]) {
     /* DFT with sine window */
-    int i=0;
+    int i;
     double *xri;
 
@@ -1636,10 +1470,5 @@
     free_dvector(xri, 1, 2 * n);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void GetInd(int n, int k, int *ind1, int *ind3) {
     if (k == 1) {
@@ -1654,10 +1483,5 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void GetInd1(int n, int k, int *ind1, int *ind3) {
     if (k == 1) {
@@ -1672,15 +1496,11 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void GetPeak(int n, double *x, int *k) {
     /* Locate peak in DFT spectrum */
-    int ind1=0, ind2=0, ind3=0;
-    double peak=0.0;
-
+    int ind1, ind2, ind3;
+    double peak;
+
+    peak = 0.0;
     *k = 1;
     for (ind2 = 1; ind2 <= n / 2 + 1; ind2++) {
@@ -1693,15 +1513,11 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void GetPeak1(int n, double *x, int *k) {
     /* Locate peak in DFT spectrum */
-    int ind1=0, ind2=0, ind3=0;
-    double peak=0.0;
-
+    int ind1, ind2, ind3;
+    double peak;
+
+    peak = 0.0;
     *k = 1;
     for (ind2 = 1; ind2 <= n; ind2++) {
@@ -1714,10 +1530,5 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 double Int2snu(int n, double *x, int k) {
     /* Get frequency by nonlinear interpolation with two samples
@@ -1728,6 +1539,6 @@
      N           A(k-1) + A(k)   2
      */
-    int ind=0, ind1=0, ind3=0;
-    double ampl1=0.0, ampl2=0.0;
+    int ind, ind1, ind3;
+    double ampl1, ampl2;
 
     GetInd(n, k, &ind1, &ind3);
@@ -1823,5 +1634,5 @@
      2      pi ( k + 1/2 )          pi ( k - 1/2 )
      */
-    double corr=0.0;
+    double corr;
 
     corr = (Sinc(M_PI * (k - 1 + 0.5 - nu * n)) + Sinc(M_PI * (k - 1 - 0.5 - nu
@@ -1858,6 +1669,6 @@
     /* Get phase by linear interpolation for rectangular window
      with -pi <= phi <= pi */
-    int i=0;
-    double phi=0.0;
+    int i;
+    double phi;
     double xr[n], xi[n];
 
@@ -1901,6 +1712,6 @@
  ****************************************************************************/
 void FndRes(struct LOC_findres *LINK) {
-    int i=0, j=0, FORLIM=0, FORLIM1=0;
-    double delta=0.0;
+    int i, j, FORLIM, FORLIM1;
+    double delta;
 
     FORLIM = LINK->n;
@@ -1970,12 +1781,7 @@
     } while (!V.found);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void GetPeaks(int n, double *x, int nf, double *nu, double *A) {
-    int i=0, k=0, ind1=0, ind3=0;
+    int i, k, ind1, ind3;
 
     for (i = 0; i < nf; i++) {
@@ -1991,12 +1797,7 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void GetPeaks1(int n, double *x, int nf, double *nu, double *A) {
-    int i=0, k=0, ind1=0, ind3=0;
+    int i, k, ind1, ind3;
 
     for (i = 0; i < nf; i++) {
@@ -2018,6 +1819,6 @@
 
 void SetTol(int Fnum, double dxrms, double dyrms, double drrms) {
-    int i=0;
-    long k=0L;
+    int i;
+    long k;
 
     for (i = 1; i <= GetnKid(Fnum); i++) {
@@ -2033,13 +1834,8 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void Scale_Tol(int Fnum, double dxrms, double dyrms, double drrms) {
-    int Knum=0;
-    long int loc=0L;
+    int Knum;
+    long int loc;
 
     for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
@@ -2079,5 +1875,5 @@
  ****************************************************************************/
 void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) {
-    long int loc=0L;
+    long int loc;
 
     loc = Elem_GetPos(Fnum, Knum);
@@ -2091,13 +1887,8 @@
     Mpole_SetdT(Fnum, Knum);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void ini_aper(const double Dxmin, const double Dxmax, const double Dymin,
         const double Dymax) {
-    int k=0;
+    int k;
 
     for (k = 0; k <= globval.Cell_nLoc; k++) {
@@ -2108,14 +1899,9 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void set_aper(const int Fnum, const double Dxmin, const double Dxmax,
         const double Dymin, const double Dymax) {
-    int i=0;
-    long int loc=0L;
+    int i;
+    long int loc;
 
     for (i = 1; i <= GetnKid(Fnum); i++) {
@@ -2127,14 +1913,9 @@
     }
 }
-/*************************************************
- * void LoadApertures(const char *ChamberFileName)
- * 
- * 
- * 
- **************************************************/
+
 void LoadApertures(const char *ChamberFileName) {
     char line[128], FamName[32];
-    long Fnum=0L;
-    double Xmin=0.0, Xmax=0.0, Ymin=0.0, Ymax=0.0;
+    long Fnum;
+    double Xmin, Xmax, Ymin, Ymax;
     FILE *ChamberFile;
 
@@ -2154,22 +1935,13 @@
     fclose(ChamberFile);
 }
-/**********************************************************************************
- * void LoadTolerances(const char *TolFileName)
- * 
- * 
- * 
- **********************************************************************************/
+
 // Load tolerances from the file
 void LoadTolerances(const char *TolFileName) {
     char line[128], FamName[32];
-    int Fnum=0;
-    double dx=0.0, dy=0.0, dr=0.0;
+    int Fnum;
+    double dx, dy, dr;
     FILE *tolfile;
 
     tolfile = file_read(TolFileName);
-    if(tolfile == NULL){
-      printf("LoadTolerances(): Error! Failure to read file: %s \n",TolFileName);
-      exit_(1);
-    }
 
     do
@@ -2193,23 +1965,13 @@
 }
 
-/**************************************************************************************
- * void ScaleTolerances(const char *TolFileName, const double scl) 
- * 
- * 
- * 
- * ************************************************************************************/
 // Load tolerances from the file
 void ScaleTolerances(const char *TolFileName, const double scl) {
     char line[128], FamName[32];
-    int Fnum=0;
-    double dx=0.0, dy=0.0, dr=0.0;
+    int Fnum;
+    double dx, dy, dr;
     FILE *tolfile;
 
     tolfile = file_read(TolFileName);
-    if(tolfile == NULL){
-      printf("LoadTolerances(): Error! Failure to read file: %s \n",TolFileName);
-      exit_(1);
-    }
-    
+
     do
         fgets(line, 128, tolfile);
@@ -2230,10 +1992,5 @@
     fclose(tolfile);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetKpar(int Fnum, int Knum, int Order, double k) {
 
@@ -2241,32 +1998,17 @@
     Mpole_SetPB(Fnum, Knum, Order);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetL(int Fnum, int Knum, double L) {
 
     Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L;
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetL(int Fnum, double L) {
-    int i=0;
+    int i;
 
     for (i = 1; i <= GetnKid(Fnum); i++)
         Cell[Elem_GetPos(Fnum, i)].Elem.PL = L;
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetdKpar(int Fnum, int Knum, int Order, double dk) {
 
@@ -2274,12 +2016,7 @@
     Mpole_SetPB(Fnum, Knum, Order);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetKLpar(int Fnum, int Knum, int Order, double kL) {
-    long int loc=0L;
+    long int loc;
 
     if (abs(Order) > HOMmax) {
@@ -2296,12 +2033,7 @@
     Mpole_SetPB(Fnum, Knum, Order);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetdKLpar(int Fnum, int Knum, int Order, double dkL) {
-    long int loc=0L;
+    long int loc;
 
     loc = Elem_GetPos(Fnum, Knum);
@@ -2332,5 +2064,5 @@
 **************************************************************/
 void SetdKrpar(int Fnum, int Knum, int Order, double dkrel) {
-    long int loc=0L;
+    long int loc;
 
     loc = Elem_GetPos(Fnum, Knum);
@@ -2343,48 +2075,28 @@
     Mpole_SetPB(Fnum, Knum, Order);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void Setbn(int Fnum, int order, double bn) {
-    int i=0;
+    int i;
 
     for (i = 1; i <= GetnKid(Fnum); i++)
         SetKpar(Fnum, i, order, bn);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetbnL(int Fnum, int order, double bnL) {
-    int i=0;
+    int i;
 
     for (i = 1; i <= GetnKid(Fnum); i++)
         SetKLpar(Fnum, i, order, bnL);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void Setdbn(int Fnum, int order, double dbn) {
-    int i=0;
+    int i;
 
     for (i = 1; i <= GetnKid(Fnum); i++)
         SetdKpar(Fnum, i, order, dbn);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetdbnL(int Fnum, int order, double dbnL) {
-    int i=0;
+    int i;
 
     for (i = 1; i <= GetnKid(Fnum); i++) {
@@ -2415,18 +2127,13 @@
 //void Setbnr(int Fnum, long order, double bnr) {
 void Setbnr(int Fnum, int order, double bnr) {
-    int i=0;
+    int i;
 
     for (i = 1; i <= GetnKid(Fnum); i++)
         SetdKrpar(Fnum, i, order, bnr);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetbnL_sys(int Fnum, int Order, double bnL_sys) {
-    int Knum=0;
-    long int loc=0L;
+    int Knum;
+    long int loc;
 
     for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
@@ -2440,13 +2147,8 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void set_dbn_rel(const int type, const int n, const double dbn_rel) {
-    long int j=0L;
-    double dbn=0.0;
+    long int j;
+    double dbn;
 
     printf("\n");
@@ -2535,12 +2237,7 @@
         return (Cell[loc].Elem.M->PBpar[Order + HOMmax]);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetdKLrms(int Fnum, int Order, double dkLrms) {
-    long int Knum=0L, loc=0L;
+    long int Knum, loc;
 
     for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
@@ -2555,12 +2252,7 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void Setdkrrms(int Fnum, int Order, double dkrrms) {
-    long int Knum=0L, loc=0L;
+    long int Knum, loc;
 
     for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
@@ -2576,24 +2268,14 @@
     }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetKL(int Fnum, int Order) {
-    long int Knum=0L;
+    long int Knum;
 
     for (Knum = 1; Knum <= GetnKid(Fnum); Knum++)
         Mpole_SetPB(Fnum, Knum, Order);
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void set_dx(const int type, const double sigma_x, const double sigma_y) {
-    long int j=0L;
+    long int j;
 
     printf("\n");
@@ -2611,12 +2293,7 @@
         }
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void SetBpmdS(int Fnum, double dxrms, double dyrms) {
-    long int Knum, loc=0L;
+    long int Knum, loc;
 
     for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
@@ -2644,7 +2321,7 @@
 ****************************************************************************/
 void codstat(double *mean, double *sigma, double *xmax, long lastpos, bool all) {
-    long i=0L, j=0L, n=0L;
+    long i, j, n;
     Vector2 sum, sum2;
-    double TEMP=0.0;
+    double TEMP;
 
     n = 0;
@@ -2711,7 +2388,7 @@
 void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos,
         long bpmdis[mnp]) {
-    long i=0L, j=0L, m=0L, n=0L;
+    long i, j, m, n;
     Vector2 sum, sum2;
-    double TEMP=0.0;
+    double TEMP;
 
     m = n = 0;
@@ -2871,5 +2548,5 @@
  ****************************************************************************/
 void GetMean(long n, double *x) {
-    long i=0L;
+    long i;
     double mean = 0e0;
 
@@ -3280,7 +2957,7 @@
     const int ntrial = 40; // maximum number of trials for closed orbit
     const double tolx = 1e-8; // numerical precision
-    int k=0;
-    int dim=0; // 4D or 6D tracking
-    long lastpos=0L;
+    int k;
+    int dim; // 4D or 6D tracking
+    long lastpos;
 
     vcod = dvector(1, 6);
@@ -3362,6 +3039,6 @@
     const int ntrial = 40; // maximum number of trials for closed orbit
     const double tolx = 1e-10; // numerical precision
-    int k=0, dim = 0;
-    long lastpos=0L;
+    int k, dim = 0;
+    long lastpos;
 
     // initializations
@@ -3429,5 +3106,5 @@
 
 void computeFandJS(double *x, int n, double **fjac, double *fvect) {
-    int i=0, k=0;
+    int i, k;
     long lastpos = 0L;
     Vector x0, fx, fx1, fx2;
@@ -3501,5 +3178,5 @@
  ****************************************************************************/
 void computeFandJ(int n, Vector &x, Matrix &fjac, Vector &fvect) {
-    int i=0, k=0;
+    int i, k;
     long lastpos = 0;
     Vector x0, fx1, fx2;
@@ -3583,6 +3260,6 @@
 
 void Newton_RaphsonS(int ntrial, double x[], int n, double tolx) {
-    int k=0, i=0, *indx;
-    double errx=0.0, d=0.0, *bet, *fvect, **alpha;
+    int k, i, *indx;
+    double errx, d, *bet, *fvect, **alpha;
 
     errx = 0.0;
@@ -3681,6 +3358,6 @@
  ****************************************************************************/
 int Newton_Raphson(int n, Vector &x, int ntrial, double tolx) {
-    int k=0, i=0;
-    double errx=0.0;
+    int k, i;
+    double errx;
     Vector bet, fvect;
     Matrix alpha;
@@ -3724,14 +3401,10 @@
     return 0;
 }
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 void rm_mean(long int n, double x[]) {
-    long int i=0L;
-    double mean=0.0;
-
+    long int i;
+    double mean;
+
+    mean = 0.0;
     for (i = 0; i < n; i++)
         mean += x[i];
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/prtmfile.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/prtmfile.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/prtmfile.cc	(revision 23)
@@ -153,6 +153,5 @@
 
  Purpose:
- Print out flatfile;
- print the lattice parameters in an external file.
+ Print out flatfile
 
  Input:
@@ -182,7 +181,7 @@
             Cell[i].dS[Y_], Cell[i].Elem.M->PdTpar, Cell[i].Elem.M->PdTsys
                 + Cell[i].Elem.M->PdTrms * Cell[i].Elem.M->PdTrnd);
-        fprintf(mfile, " %23.16e %23.16e %23.16e %23.16e %23.16e  %23.16e %23.16e\n",
+        fprintf(mfile, " %23.16e %23.16e %23.16e %23.16e %23.16e\n",
             Cell[i].Elem.PL, Cell[i].Elem.M->Pirho, Cell[i].Elem.M->PTx1,
-            Cell[i].Elem.M->PTx2, Cell[i].Elem.M->PH1,Cell[i].Elem.M->PH2,Cell[i].Elem.M->Pgap);
+            Cell[i].Elem.M->PTx2, Cell[i].Elem.M->Pgap);
         prtHOM(mfile, Cell[i].Elem.M->n_design, Cell[i].Elem.M->PB,
             Cell[i].Elem.M->Porder);
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/rdmfile.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/rdmfile.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/rdmfile.cc	(revision 23)
@@ -217,7 +217,7 @@
                 if (prt)
                     printf("%s\n", line);
-                sscanf(line, "%lf %lf %lf %lf %lf  %lf %lf", &Cell[i].Elem.PL,
+                sscanf(line, "%lf %lf %lf %lf %lf", &Cell[i].Elem.PL,
                         &Cell[i].Elem.M->Pirho, &Cell[i].Elem.M->PTx1,
-                        &Cell[i].Elem.M->PTx2, &Cell[i].Elem.M->PH1, &Cell[i].Elem.M->PH2, &Cell[i].Elem.M->Pgap);
+                        &Cell[i].Elem.M->PTx2, &Cell[i].Elem.M->Pgap);
                 if (Cell[i].Elem.M->Pirho != 0.0)
                     Cell[i].Elem.M->Porder = 1;
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/read_script.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/read_script.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/read_script.cc	(revision 23)
@@ -25,14 +25,4 @@
    Comments:
        Written by Jianfeng Zhang 03/2011 soleil
-       
-       (1) 15/10/2013 by Jianfeng Zhang @ LAL
-           Fix the bugs of fgets():
-              to read lines with arbritracy length, 
-	      to read the empty space at the end of the line.
-	      
-	   Remove the rontine to move the address of "line" if the 
-	   contents of "line" is empty space; otherwise free() can't
-	   find the orignial address of "line". 
-           
 
 ****************************************************************************/
@@ -40,5 +30,4 @@
 /* files */ 
 
-char    lat_file[max_str]="voidlattice";
 // girder file
 char girder_file[max_str];
@@ -78,13 +67,11 @@
  
   
-  char    str[max_str]="voidstring", dummy[max_str]="voidstring",dummy2[max_str]="voidstring", nextpara[max_str]="voidpara";
+  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 = NULL; //line to store the command without preceding white space
-  size_t   len = 0;
-  ssize_t  read;
+  char    *line; //line to store the command without preceding white space
   char    name[max_str]="voidname";
- // char    lat_file[max_str]="voidlattice";
+  char    lat_file[max_str]="voidlattice";
   char    EndName[]="void";
-
+ 
   FILE    *inf;
   const bool  prt = false; // for debugging printout each line of input file
@@ -100,30 +87,34 @@
   // Manipulation of the parameter file
   strcpy(dummy, param_file_name);
-
-  sprintf(full_param_file_name,"%s",dummy);
-  if (!prt) printf("\n reading in script file: %s\n",full_param_file_name);
+  pch = strstr (dummy,".prm"); // search for extension .prm
+  /* remove additional .prm if exist */
+  if (pch != NULL) strncpy(pch,"\0",1);
+
+  sprintf(full_param_file_name,"%s.prm",dummy);
+  if (prt) printf("\n reading in script file: %s\n",full_param_file_name);
 
   //
   inf = file_read(full_param_file_name);
-   if(inf == NULL){
-     printf("read_script(): Error! Read File %s Failure!\n",full_param_file_name);
-     exit_(1);  
-  }
+
   
   // read parameter file, line by line 
-   while ((read = getline(&line, &len, inf)) != -1) {
-     
-     LineNum++;
-     if(prt){
-       printf("Line # %ld \n",LineNum);
-       printf("Retrieved line of length %zu : \n",read);
-       cout << line << endl;
-     }
-               
-     
+ //  while (fgets(line, max_str, inf) != NULL) {
+   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++;
+     }       
+    
+    if(prt)
+      printf("cfg: %s",line);
+      
+      LineNum++;
            
     /* read the end of line symbol '\n','\r' or '\r\n' at different operation system*/
-    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) {
       // get initial command token
       sscanf(line, "%s", name); 
@@ -260,21 +251,4 @@
 	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
       }  
-      //print the cooridinates using tracking at each element
-      else if (strcmp("PrintTrackElemFlag", name) == 0){
-	  sscanf(line, "%*s %s",nextpara);
-	 
-	  if(strcmp(nextpara,"voidpara")!=0)
-	    sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld %ld", 
-	                UserCommandFlag[CommNo]._PrintTrackElem_track_file,
-          		&(UserCommandFlag[CommNo]._PrintTrackElem_x), 
-			&(UserCommandFlag[CommNo]._PrintTrackElem_px),
-          		&(UserCommandFlag[CommNo]._PrintTrackElem_y), 
-			&(UserCommandFlag[CommNo]._PrintTrackElem_py), 
-			&(UserCommandFlag[CommNo]._PrintTrackElem_delta), 
-			&(UserCommandFlag[CommNo]._PrintTrackElem_ctau),
-			&(UserCommandFlag[CommNo]._PrintTrackElem_nelem1),
-	    		&(UserCommandFlag[CommNo]._PrintTrackElem_nelem2));
-	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
-      } 
       //print close orbit(COD) flag
       else if (strcmp("PrintGirderFlag", name) == 0){
@@ -293,10 +267,9 @@
       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", 
+            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s", 
               		  UserCommandFlag[CommNo]._FmapFlag_fmap_file,
 			  &(UserCommandFlag[CommNo]._FmapFlag_nxpoint), 
@@ -306,5 +279,5 @@
 			  &(UserCommandFlag[CommNo]._FmapFlag_ymax),
           		  &(UserCommandFlag[CommNo]._FmapFlag_delta), 
-			  dummy,dummy2);
+			  dummy);
          
         if(strcmp(dummy, "true") == 0)
@@ -313,15 +286,4 @@
           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);
@@ -330,9 +292,8 @@
       else if (strcmp("FmapdpFlag", name) == 0){         
     	  strcpy(dummy, "");
-	  strcpy(dummy2,"");
 	  sscanf(line, "%*s %s",nextpara);
-	  
+	 
 	  if(strcmp(nextpara,"voidpara")!=0)
-            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s", 
+            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s", 
           		UserCommandFlag[CommNo]._FmapdpFlag_fmapdp_file,
 			&(UserCommandFlag[CommNo]._FmapdpFlag_nxpoint),
@@ -342,5 +303,5 @@
 			&(UserCommandFlag[CommNo]._FmapdpFlag_emax),
           		&(UserCommandFlag[CommNo]._FmapdpFlag_z), 
-			dummy,dummy2);
+			dummy);
 	 
         if(strcmp(dummy, "true") == 0)
@@ -349,15 +310,4 @@
           UserCommandFlag[CommNo]._FmapdpFlag_diffusion = false;
 
-	if(strcmp(dummy2, "true") == 0)
-          UserCommandFlag[CommNo]._FmapdpFlag_printloss = true;
-        else if(strcmp(dummy2, "false") == 0)
-          UserCommandFlag[CommNo]._FmapdpFlag_printloss = false;
-	else
-	  UserCommandFlag[CommNo]._FmapdpFlag_printloss = false;
-	
-	
-	cout << "debug:      " << line << endl;
-	cout << "dummy2 =  " << dummy2 << "     lossflag =  " << UserCommandFlag[CommNo]._FmapdpFlag_printloss << endl;
-	
 	  //  FmapdpFlag = true;
          strcpy(UserCommandFlag[CommNo].CommandStr,name);
@@ -585,12 +535,9 @@
         exit_(1);
       }
-  }
+    }
     /* continue read in the line */ 
     else
-     continue;
-  }//end of while loop
-  
-  free(line);
-  
+      continue;
+  }
   fclose(inf);
 }  
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc	(revision 23)
@@ -25,5 +25,4 @@
        none
 
-       
    Global variables:
        trace
@@ -38,5 +37,5 @@
 void Get_Disp_dp(void)
 {
-  long i=0L;
+  long i;
 //  long lastpos = 0;
   const char nomfic[] = "dispersion.out";
@@ -96,10 +95,10 @@
 {
   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,5 +200,5 @@
 {
   CellType Cell;
-  long i=0L;
+  long i;
 Vector2 H;
 
@@ -252,5 +251,5 @@
 {
   CellType Cell;
-  long i=0L;
+  long i;
   long lastpos = 1L;
   Vector2 H;
@@ -462,15 +461,9 @@
        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,24 +473,20 @@
 
   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,5 +559,4 @@
  //     printf("%s", line);
   } 
-  free(line);
   fclose(fp);
   // turn on the global flag for CheckAmpl()
@@ -617,5 +605,5 @@
   bool lostF = true; /* Lost particle Flag */
   Vector x1;            /* tracking coordinates */
-  long i=0L;
+  long i;
   Vector2  aperture;
   aperture[0] = 1e0; aperture[1] = 1e0;
@@ -811,13 +799,8 @@
 #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,16 +834,14 @@
    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,5 +876,5 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2 = true;
+ bool status = true;
  struct tm *newtime;
 
@@ -946,6 +927,6 @@
  //    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,214 +974,4 @@
 #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,11 +1024,10 @@
 #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,5 +1037,5 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2 = true;
+ bool status = true;
  struct tm *newtime;
 
@@ -1274,14 +1044,4 @@
  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,12 +1059,4 @@
  }
 
-//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,14 +1064,5 @@
   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,11 +1113,6 @@
  
      // 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,18 +1148,8 @@
      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,9 +1181,8 @@
        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,18 +1201,13 @@
        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,15 +1217,6 @@
  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,5 +1227,5 @@
  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,21 +1240,4 @@
 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,15 +1262,9 @@
 	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,16 +1300,8 @@
         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,5 +1338,5 @@
 
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -1671,17 +1353,13 @@
        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,5 +1369,5 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2=true;
+ bool status=true;
  struct tm *newtime;
 
@@ -1699,13 +1377,4 @@
  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,11 +1390,4 @@
    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,15 +1398,5 @@
      // 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,15 +1448,9 @@
 	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,16 +1486,8 @@
         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,5 +1535,5 @@
   double nux1 = 0.0, nuz1 = 0.0;
   int nb_freq[2] = {0, 0};
-  bool status2 = true;
+  bool status = true;
   struct tm *newtime;
 
@@ -1929,6 +1567,6 @@
     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,5 +1574,5 @@
     else {
        nux1 = 0.0; nuz1 = 0.0;
-       status2 = true;
+       status = true;
     }
     
@@ -1991,6 +1629,6 @@
   FILE *outf;
   const char *fic = phasefile; 
-  int i=0;
-  bool status2=false;
+  int i;
+  bool status;
   struct tm *newtime;
 
@@ -2025,11 +1663,9 @@
   }
   
-  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,7 +1705,7 @@
   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,5 +1799,5 @@
   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,5 +1858,5 @@
    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,5 +1901,5 @@
   FILE *outf;
   const char fic[] = "check_ampl.out";
-  int i=0;
+  int i;
 
   if ((outf = fopen(fic, "w")) == NULL)
@@ -2320,5 +1956,5 @@
   FILE *outf;
   const char fic[] = "enveloppe.out";
-  int i=0,j=0 ;
+  int i,j ;
   CellType Cell;
 
@@ -3562,9 +3198,9 @@
   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,5 +3403,5 @@
   struct tm     *newtime;  // for time
   Vector        codvector[Cell_nLocMax];
-  bool          cavityflag=false, radiationflag=false;
+  bool          cavityflag, radiationflag;
   bool          trace=true; 
   
@@ -4165,5 +3801,5 @@
   struct tm     *newtime;  // for time
   Vector        codvector[Cell_nLocMax];
-  bool          cavityflag=false, radiationflag=false;
+  bool          cavityflag, radiationflag;
   bool          trace=true; 
   
@@ -4650,5 +4286,5 @@
  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,5 +4294,5 @@
  int nb_freq[2] = {0, 0};
  long nturn = Nbtour;
- bool status2=true;
+ bool status=true;
  struct tm *newtime;
 
@@ -4699,6 +4335,6 @@
    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,5 +4415,5 @@
 	    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,5 +4428,5 @@
   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,5 +4499,5 @@
   CellType Cell;
   int qlist[320];
-  int nquad=0, i=0;
+  int nquad=0, i;
   double A = 0.0;
 
@@ -4897,5 +4533,5 @@
 /****************************************************************************/
 /* void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-   double energy, bool *status2)
+   double energy, bool *status)
 
    Purpose:
@@ -4917,5 +4553,5 @@
 
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -4939,5 +4575,5 @@
  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,5 +4584,5 @@
  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,7 +4644,7 @@
    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,5 +4727,5 @@
 /****************************************************************************/
 /* void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-   double energy, bool *status2)
+   double energy, bool *status)
 
    Purpose:
@@ -5111,5 +4747,5 @@
 
    Output:
-       status2 true if stable
+       status true if stable
               false otherwise
 
@@ -5133,5 +4769,5 @@
   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,5 +4776,5 @@
   int nb_freq[2] = {0, 0};
   long nturn = Nbtour;
-  bool status2=true;
+  bool status=true;
   struct tm *newtime;
 
@@ -5174,15 +4810,15 @@
     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,13 +4832,13 @@
     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,5 +4883,5 @@
   FILE *outf;
   const char fic[] = "phase2.out";
-  long lastpos = 0L,lastn = 0L;
+  long lastpos = 0,lastn = 0;
   struct tm *newtime;
 
@@ -5373,5 +5009,5 @@
   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,7 +5431,7 @@
 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,13 +5629,8 @@
 }
 
-/*******************************************************************************
- * 
- * 
- * 
- * 
- ******************************************************************************/
+
 double EnergySmall(double *X, double irho)
 {
- double A=0.0, B=0.0;
+ double A, B;
  double h = irho;
 
@@ -6008,13 +5639,8 @@
  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,5 +5681,5 @@
   FILE *outf;
   const char fic[] = "enveloppe2.out";
-  int i=0,j=0 ;
+  int i,j ;
   CellType Cell;
   /* Array for Enveloppes */
@@ -6185,5 +5811,5 @@
 void set_RFVoltage(const int Fnum, const double V_RF){
 
-  int k=0, n = 0;
+  int k, n = 0;
   
   
@@ -6232,18 +5858,14 @@
        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,20 +5876,18 @@
 
   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,5 +5907,5 @@
 	/*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,5 +5931,5 @@
 	    sscanf(prm, "%lf", &An);
 	  
-	    if (!prt)
+	    if (prt)
 	      printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
 	  
@@ -6335,5 +5955,5 @@
      // printf("%s", line);
   }
-  free(line);
+
   fclose(inf);
 }
@@ -6393,5 +6013,5 @@
       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,5 +6294,5 @@
 	 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,5 +6302,5 @@
   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,5 +6418,5 @@
 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,5 +6467,5 @@
       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,6 +6479,6 @@
 {
    
-    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,9 +6525,8 @@
     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,7 +6552,6 @@
                      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,93 +6560,5 @@
 	}//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)
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/t2elem.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/t2elem.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/t2elem.cc	(revision 23)
@@ -1,3 +1,3 @@
-/*Tracy-3
+/* Tracy-3
 
  J. Bengtsson, CBP, LBL      1990 - 1994   Pascal version
@@ -9,5 +9,4 @@
  Element propagators.                                                      */
 
-
 double c_1, d_1, c_2, d_2, cl_rad, q_fluct;
 double I2, I4, I5, dcurly_H, dI4;
@@ -20,6 +19,4 @@
 
 // Dynamic model
-
-
 
 /****************************************************************************/
@@ -84,15 +81,12 @@
 
  Purpose:
- Get the longitudinal momentum ps; for the exact/approximation Hamiltonian
+ Get the longitudinal momentum ps
  
     for approximation Hamitonian:
     ps = 1+delta   
     
-    
     For exact hamitonian:
     
-    (1) ps = sqrt((1+delta)^2 - px_^2 - py_^2)
-    
-    (2) using the check of TEAPOT to check ps < 0 
+    ps 
  **********************************************************/
 template<typename T>
@@ -101,5 +95,5 @@
     T p_s, p_s2;
 
-   if (!globval.H_exact)
+    if (!globval.H_exact)
     p_s = 1.0+x[delta_];
     else {
@@ -109,10 +103,7 @@
         else {
             // avoid compile warning
-            //p_s = 0.0;
-	     p_s = 1.0e-20;
-	    cout << "        " << 1+x[delta_] << "   " << x[px_] << "   " << x[py_] << endl;
+            p_s = 0.0;
             printf("get_p_s: *** Speed of light exceeded!\n");
-           // exit_(1);
-	    
+            exit_(1);
         }
     }
@@ -124,9 +115,15 @@
 
  Purpose:
- Exact Drift pass in the curvilinear coordinate for the dipole
-
- See Eqn. (10.26) on page 307 in E. Forest's beam dynamics, a new attitude
- 
- H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta                                  
+ Drift pass
+
+ If H_exact = false, using approximation Hamiltonian:
+
+                                   px^2+py^2
+	      H(x,y,l,px,py,delta) = -------------
+	                               2(1+delta)     
+       
+
+ Otherwise, using exact Hamiltonian
+
 
  Input:
@@ -144,67 +141,4 @@
 
  Specific functions:
- 
- 
- 
- Comments:
- 
- Test version....rewritten by Jianfeng zhang @ LAL, 31/07/2013
-
- ****************************************************************************/
-
-
-template<typename T>
-void Drift(double L,double h_bend, ss_vect<T> &x) {
-    T u;
-    ss_vect<T> x0;
-
-    x0=x;
-    
-    if (globval.H_exact) {
-      if(get_p_s(x) == 0)
-	return;
-
-      x[x_] = x0[x_]/cos(L*h_bend)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
-      x[px_]= x0[px_]*cos(L*h_bend) + get_p_s(x)*sin(L*h_bend);
-      x[y_] = x0[y_]+x0[px_]*x0[x_]*tan(L*h_bend)/get_p_s(x)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
-      x[py_]= x0[py_];
-      x[ct_]= x0[ct_]+((1+x0[delta_])*x0[x_]*tan(L*h_bend))/get_p_s(x)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
- 
-    
-    }
-}
-
-/****************************************************************************/
-/* Drift(double L, ss_vect<T> &x)
-
- Purpose:
- Drift pass in the Cartesian coordinate
-
- If H_exact = false, using approximation Hamiltonian (only valid for big ring):
-
-                                       px^2+py^2
-	      H(x,y,cT,px,py,delta) = -------------
-	                               2(1+delta)     
-       
-
- Otherwise, using the exact Hamiltonian (valid for small ring) 
- (J. Bengtsson: Tracy-2 user's mannual):
- 
- H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta                                  
-
- Input:
- L:  drift length
- &x:  pointer to initial vector: x
-
- Output:
- none
-
- Return:
- none
-
- Global variables:
-
-
- Specific functions:
 
  ****************************************************************************/
@@ -215,15 +149,11 @@
     T u;
 
-    if (globval.H_exact) {
-      if(get_p_s(x) == 0)
-	return;
-      
-   u = L/get_p_s(x);
+    if (!globval.H_exact) {
+    u = L/(1.0+x[delta_]);
+    x[ct_] += u*(sqr(x[px_])+sqr(x[py_]))/(2.0*(1.0+x[delta_]));
+    } else {
+    u = L/get_p_s(x);
     x[ct_] += u*(1.0+x[delta_]) - L;
-    } else {
-     u = L/(1.0+x[delta_]);
-    x[ct_] += u*(sqr(x[px_])+sqr(x[py_]))/(2.0*(1.0+x[delta_]));
-    }
-    
+    }
     x[x_] += x[px_] * u;
     x[y_] += x[py_] * u;
@@ -671,8 +601,7 @@
     Correction for magnet gap (longitudinal fringe field)
 
-    Input: 
-     irho:         h = 1/rho [1/m]
-     phi:         dipole edge (entrance or exit) angle
-     gap:        full gap between poles
+     irho h = 1/rho [1/m]
+     phi  dipole edge angle
+     gap  full gap between poles
 
                              2
@@ -683,10 +612,6 @@
      K1 is usually 1/2
      K2 is zero here                                                  
-     For the type of Tanh-like fringe field
-     
-
-     18/06/2012  Jianfeng Zhang@LAL
-     
-     Fix the bug to calculate psi for the sector magnets which has phi=0
+
+
 
 *********************************************************************/
@@ -697,86 +622,11 @@
     const double k1 = 0.5, k2 = 0.0;
 
-//     if (phi == 0.0)  //NOT hard edge, but sector magnets!!!
-//         psi = 0.0;  
-//     else
+    if (phi == 0.0)  //hard edge
+        psi = 0.0;  
+    else
         psi = k1 * gap * irho * (1.0 + sqr(sin(dtor(phi)))) / cos(dtor(phi))
                 * (1.0 - k2 * gap * irho * tan(dtor(phi)));
 
     return psi;
-}
-
-/*****************************************************************
- * Exact map of the sector dipole
- * 
- * Forest's beam dynamics P.360 Eqn. (12.18)
- * Also see Laurent's thesis: P219-220.
- * 
- * Written by Jianfeng Zhang @ LAL, 01/08/2013
- * 
- * Test version........................................
- * 
- * ****************************************************************/
-template<typename T>
-void dipole_kick(double L, double h_ref, double h_bend, ss_vect<T> &x){
-
-  ss_vect<T> x0;
-  T ps, u, dpxf;
-    
-  x0 = x;
-  ps = get_p_s(x0);
-  
-  u = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
-  
-  x[py_] = x0[py_];
-  
-  x[delta_] =x0[delta_];
-  
-  x[px_]= x0[px_]*cos(L*h_ref) + (ps -h_bend*(1/h_ref+x0[x_]))*sin(L*h_ref);
-  
-  dpxf = h_ref*sqrt((1+x[delta_])*(1+x[delta_]) - x[px_]*x[px_] - x[py_]*x[py_])-h_bend*(1+h_ref*x[x_]);
-  
-  
-  x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
-  
-  x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x[delta_])-sqr(x[px_])-sqr(x[py_]))-dpxf - h_bend);
-  
-  x[ct_] += (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
- 
-}
-
-/*****************************************************************
- * pure multipole kick
- * 
- * Forest's beam dynamics P.360 Eqn. (12.18)
- * Also see Laurent's thesis: P219-220.
- * 
- * Written by Jianfeng Zhang @ LAL, 01/08/2013
- * 
- * Test version........................................
- * 
- * ****************************************************************/
-template<typename T>
-void multipole_kick(int Order, double MB[], double L, double h_bend, ss_vect<T> &x){
-
-  int j=0;
-  ss_vect<T> x0;
-  
-  T ByoBrho = 0, BxoBrho = 0, ByoBrho1=0 ;
-  
-  if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) {
-        x0 = x;
-        /* compute the total magnetic field Bx and By with Horner's rule */
-        ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e)
-        BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e)
-        for (j = Order - 1; j >= 1; j--) {
-            ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax];
-            BxoBrho = x0[y_] * ByoBrho + x0[x_] * BxoBrho + MB[HOMmax - j];
-            ByoBrho = ByoBrho1; 
-	}
-  }
-  
-  x[px_] -= L * (ByoBrho + h_bend*1); 
-	//vertical kick due to the magnets
-        x[py_] += L * BxoBrho;
 }
 
@@ -786,32 +636,24 @@
 
   Purpose:
-        Refer to Tracy 2 manual.
-        Calculate multipole kick. 
-  , 
-  (1)     The exact Hamiltonian is
+        Calculate multipole kick. The Hamiltonian is
 
             H = A + B where A and B are the kick part defined by
-                 
-                 (exact form of the H of the drift)
-                 A(x,y,-l,px,py,dP) = -sqrt( (1 + delta)^2 - px^2 - py^2 )  + delta 
-                                      
-                 (exact form of the H of the kick from the magnet)                                       
-                 B(x,y,-l,px,py,dP) = -h*x*sqrt( (1 + delta)^2 - px^2 - py^2 ) + h x + h^2*x^2/2 + int(Re(By+iBx)/Brho)   
-                 
-                 so this is the exact Hamitonian for small ring.
-                 
-
-	The kick is given by
-                                                                                                e                     
-         Dp_x = L* (h_ref*sqrt( (1 + delta)^2 - (px^2 + py^2) ) - h_bend - x*h_bend*h_ref -   ---- B_y ) 
-                                                                                               p_0                 
-         Dp_y = L* e/p_0 * B_x
-         
-         x    = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * p_x
-         
-         y    = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * p_y
-         
-         CT   = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * (1+delta)
-                                                                                           
+                                         2    2
+                                       px + py
+                 A(x,y,-l,px,py,dP) = ---------
+                                       2(1+dP)
+                                                   2 2
+                 B(x,y,-l,px,py,dP) = -h*x*dP + 0.5 h x + int(Re(By+iBx)/Brho)
+
+                 so this is the appproximation for large ring
+                 the chromatic term is missing hx*A
+
+
+		The kick is given by
+
+              e L       L delta    L x              e L
+     Dp_x = - --- B_y + ------- - ----- ,    Dp_y = --- B_x
+              p_0         rho     rho^2             p_0
+
     where
         e      1
@@ -823,51 +665,12 @@
                            /
                            ====
-  
-  (2)     The approximate Hamiltonian is
-
-            H = A + B where A and B are defined by
-                                         2    2
-                                       px + py
-                 A(x,y,-l,px,py,dP) = ---------                   (approximate form of the H of the drift)
-                                      2(1+delta)
-                                                        2 2
-                 B(x,y,-l,px,py,dP) = -h*x*delta + 0.5 h x + int(Re(By+iBx)/Brho)   (approximate form of the H of the kick from the magnet)
-
-                 so this is the appproximation for large ring
-                 the chromatic term is missing hx*A
-
-
-		The kick is given by
-
-              e L       L delta    L x              e L
-     Dp_x =  --- B_y + ------- - ----- ,    Dp_y = --- B_x
-              p_0         rho     rho^2             p_0
-
-     //second order of the Hamiltonian motion
-     x    =  h_ref * x_0 / (1+delta) * p_x
-     
-     y    =  h_ref * x_0 / (1+delta) * p_y
-     
-    CT    =  h_ref * x_0 / (1+delta) * (p_x^2 + p_y^2)/2/(1+delta)
-     
-    where
-        e      1
-      --- = -----
-      p_0   B rho   
-                           ====
-                           \ 
-      (B_y + iB_x) = B rho  >   (ia_n  + b_n ) (x + iy)^n-1
-                           /
-                           ====
-
-                    = B rho [(b_n*x - a_n*y) + i*(bn*y + an*x) + b_(n-1) + i*a_(n-1)] * (x+iy)^n-2
-                    =......
- 
+
  Input:
  Order  maximum non zero multipole component
  MB     array of an and bn, magnetic field components
  L      multipole length
- h_bend   B_dipole/Brho is the dipole field normalized by the reference momentum; one of the example is the combined dipole
- h_ref   1/rho [m^-1] is the curvature of the reference orbit in the dipoles which is used in curvilinear coordinate,
+ h_bend   1/rho in curvilinear coordinate,
+ 0    in cartisian cooridinate.
+ h_ref   1/rho [m^-1]
  x      initial coordinates vector
 
@@ -885,7 +688,5 @@
 
  Comments:
- 06/11/2012     Jianfeng Nadolski @ LAL
-                Fix the bug in the kick map of the exact Hamiltonian. 
-                Add the second order of the approximate Hamiltonian.
+ none
 
  **************************************************************************/
@@ -895,14 +696,12 @@
     int j;
     T BxoBrho, ByoBrho, ByoBrho1, B[3];
-    T sqrtpx, dpx, ps2new, psnew;
-    T n1, n2;
     ss_vect<T> x0, cod;
-    T u=0.0;
-    
+
     if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) {
         x0 = x;
-        /* compute the total magnetic field Bx and By with Horner's rule */
+        /* compute field with Horner's rule */
         ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e)
         BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e)
+
         for (j = Order - 1; j >= 1; j--) {
             ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax];
@@ -911,5 +710,4 @@
         }
 
-        
         if (globval.radiation || globval.emittance) {
             B[X_] = BxoBrho;
@@ -919,87 +717,12 @@
         }
 
-        //curvilinear cocodinates, from the dipole components
         if (h_ref != 0.0) {
-	  
-	  //exact Hamiltonian
-	  if (globval.H_exact) {
-	    
-	    if(get_p_s(x0)==0)
-	      return;
-	    
-	    //dipole kick map from the exact Hamiltonian; not symplectic
-// 	    u = L * h_ref * x0[x_] /get_p_s(x0);
-//             x[x_] += u * x0[px_];
-//             x[y_] += u * x0[py_];
-// 	    x[px_]-= L * (-h_ref * get_p_s(x0) + h_bend + h_ref * h_bend
-// 	             * x0[x_] + ByoBrho);
-//             x[ct_] += u * (1.0+x0[delta_]);
-	     
-	      //kick map due to the dipole vector potential
-	     // x[px_] -= L*(h_bend + h_ref*h_bend*x0[x_]+ByoBrho);
- 
-//          sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
-// // 	
-// 	x[px_] = x0[px_]*cos(L*h_ref) + (get_p_s(x0)-h_bend*(h_ref+x0[x_]))*sin(L*h_ref)-L*ByoBrho;
-// 	
-// 	dpx = -x0[px_]*h_ref*sin(L*h_ref)+(get_p_s(x0)-h_bend*(h_ref+x0[x_]))*h_ref*cos(L*h_ref);
-// 	
-// 	x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[py_]))-dpx - h_bend);
-// 	
-// 	x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(sqrtpx*x0[px_])-asin(sqrtpx*x[px_]));
-// 	
-// 	x[ct_]+= (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(sqrtpx*x0[px_])-asin(sqrtpx*x[px_]));
-	
-        //map of a wedge 
-//          x[px_] = x0[px_]*cos(L*h_bend) + (get_p_s(x0) - h_ref*x0[x_])*sin(L*h_bend) -L*ByoBrho;
-// 	 
-// 	 sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
-// 	 n1 = x0[px_]*sqrtpx;
-// 	 n2 = x[px_]*sqrtpx;
-// 	 
-//          x[x_] = x0[x_]*cos(L*h_bend) + (x0[x_]*x0[px_]*sin(2*L*h_bend)+sqr(sin(L*h_bend))*(2*x0[x_]*get_p_s(x0)-h_ref*sqr(x0[x_])))
-//                  /(sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[y_]))+get_p_s(x0)*cos(L*h_bend)-x0[px_]*sin(L*h_bend));
-// 
-//          x[y_] += x0[py_]*L*h_bend*h_ref + x0[py_]*h_ref*(asin(n1) - asin(n2));
-//          x[ct_]+= (1+x0[delta_])*L*h_bend*h_ref + (1+x0[delta_])*h_ref*(asin(n1) - asin(n2));
-//  
-          }//approximate Hamiltonian 
-          else {//first order map
             x[px_] -= L * (ByoBrho + (h_bend - h_ref) / 2.0 + h_ref * h_bend
                     * x0[x_] - h_ref * x0[delta_]);
             x[ct_] += L * h_ref * x0[x_];
-  
-	    //  second order, cubic term of H, H3; has problem when do tracking...works for THOMX!!!!!
-	    // not symplectic????
-         if(0){	
-	    u = L * h_ref * x0[x_] /(1.0+x0[delta_]);
-            x[x_] += u * x0[px_];
-            x[y_] += u * x0[py_];
-	 //   x[px_]+= 0.5*u/x0[x_]*(x0[px_]*x0[px_]+x0[py_]*x0[py_]);
-            x[ct_] += u*(sqr(x0[px_])+sqr(x0[py_]))/(2.0*(1.0+x0[delta_]));
-            
-	   // cout << "test......2.....  u = " << u<< endl;
-	  }	
-      }
-	  
-	}//from other straight magnets 
-	else 
-            x[px_] -= L * (ByoBrho + h_bend*1); //h_bend should be trigger on for combined dipoles
-
-	//vertical kick due to the magnets
+        } else
+            x[px_] -= L * (ByoBrho + h_bend);
+
         x[py_] += L * BxoBrho;
-      
-	/*
-	if(h_ref == 0 && h_bend !=0){
-	sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
-	ps2new = sqr(1+x0[delta_]) - sqr(x[px_]) -sqr(x0[py_]);
-	psnew = sqrt(ps2new);
-	
-	x[x_] += 1/h_bend*(psnew - get_p_s(x0));
-	x[y_] += x0[py_]/h_bend*(asin(x0[px_]*sqrtpx)-asin(x[px_]*sqrtpx));
-	x[ct_]+= (1+x0[delta_])/h_bend*(asin(x0[px_]*sqrtpx)-asin(x[px_]*sqrtpx));
-	}*/
-	
-	
     }
 }
@@ -1012,5 +735,5 @@
  Compute edge focusing for a dipole
  There is no radiation coming from the edge
- The standard formula used is : (SLAC-75, K. Brown's first order of fringe field but with the correction of 1/(1+delta) in R43)
+ The standard formula used is :
 
        px = px0 +  irho tan(phi) *x0
@@ -1029,7 +752,5 @@
  gap  = gap of the dipole for longitudinal fringing field (see psi)
  x    = input particle coordinates
-Entrance:   bool flag, true, then calculate edge focus and fringe field at the entrance of the dipole
-                       false,...........................................at the exit fo the dipole
- 
+
  Output:
  x    output particle coordinates
@@ -1053,121 +774,18 @@
 
  April 2011, No energy dependence for H-plane according to SSC-141 note (E. Forest)
- Agreement better with MAD8 (LNLS), SOLEIL small effect (.1). Modified by Laurent Nadolski @ soleil
- Comments: Jianfeng Zhang @ LAL, 05/06/2013
- This model only works for Soleil lattice, but not works for ThomX and SuperB DR lattice.
- 
- 18/06/2012  Jianfeng Zhang @ LAL
- Add the contribution of x' to the x[py_] for small ring like Thom-X, using Alexandre Loulergue's correction.
- This models works for ThomX, Soleil, and SuperB DR lattices; and 
- close to the Forest's model and K. Brown's models.
- 
+ Agreement better with MAD8 (LNLS), SOLEIL small effect (.1)
  ****************************************************************************/
 template<typename T>
-static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x, bool Entrance) {
-  
-
+static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x) {
     if (true) {
-      
-       // cout << "Dipole edge effect WITH Alex's  correction: " << endl;
-      
-      //with the contribution from the entrance angle of the particle at the edge
-// Original model, written by L. Nadolski, this is a reasonal model!!!! 
-//    warning: => diverging Taylor map (see SSC-141)
-//         x[px_] += irho * tan(dtor(phi)) * x[x_];
-//         x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]
-//                 / (1.0 + x[delta_]);
-
-
-// K.Brown 2rd order transfer matrix of the dipole fringe field for sector dipoles (phi=0)
-// compared the fringe field model used in Tracy
-// results:  K. Brown's model is wrong for small ring like Thom-X!!!!
-// The original model of tracy is correct for tracking!!!!!!!!!!
-  /*
-  double psi;
-    psi=gap*0.5*irho*(1+sin(dtor(phi))*sin(dtor(phi)))/cos(dtor(phi));
- 
-    if(Entrance){
-   x[x_] += irho/2*x[y_]*x[y_];
-   x[px_] += 0.5*x[x_]*x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5*x[y_]*x[y_];
-   x[py_] -= -irho*tan(-dtor(psi))*x[y_]+1*x[x_]*x[y_]+irho*x[px_]*x[y_]+ (irho*psi/cos(-dtor(psi))/cos(-dtor(psi))) * x[y_]
-                 / (1.0 + x[delta_]);
-    }
-    else
-    {
-       x[x_] -= irho/2*x[y_]*x[y_];
-   x[px_] += 0.5*x[x_]*x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5*x[y_]*x[y_];
-   x[py_] -= -irho*tan(-dtor(psi))*x[y_]+1*x[x_]*x[y_]-irho*x[px_]*x[y_]+ (irho*psi/cos(-dtor(psi))/cos(-dtor(psi))) * x[y_]
-                 / (1.0 + x[delta_]);
-    }
-    */
-
-    //add the contribution of the entrance momentum of the particle, from Alex Loulergue. 
-    if(Entrance){
-       x[px_] += irho * tan(dtor(phi) +x[py_]/(1+x[delta_]*1)) * x[x_]/ (1.0 + x[delta_]*0);
-       x[py_] -= irho * tan(dtor(phi) + x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)/(1+x[delta_]*0)) * x[y_]
-                / (1.0 + x[delta_]*0);
-    }else{
-       x[px_] += irho * tan(dtor(phi) - x[py_]/(1+x[delta_]*1)) * x[x_]/ (1.0 + x[delta_]*0);
-       x[py_] -= irho * tan(dtor(phi) - x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)/(1+x[delta_]*0)) * x[y_]
-                / (1.0 + x[delta_]*0);
-    }
-      
-    } else {//original one
-
-   //cout << "Dipole edge effect WITHOUT Alex's  correction: " << endl;
+        // warning: => diverging Taylor map (see SSC-141)
         x[px_] += irho * tan(dtor(phi)) * x[x_];
-        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]/ (1.0 + x[delta_]*1);
-    }
-}
-
-
-/******************************************************************
- template<typename T>
-static void BendCurvature(double irho, double H, ss_vect<T> &x)
-
-The leanding order term T211, T233, and T413 on page 117 & page 118
-of the dipole edge effect in SLAC-75, using the symplectic
-form in 
- E. Forest: "Fringe field in MAD, Part II: Bend curvature in MAD-X for the 
-             mudule PTC". 
-             P.10 Eq. (5) and (42).
- 
- For an rectangular dipole, the bending curvature term is zero, since 1/R = 0.
- 
- irho:   curvature of the central orbit inside the dipole
- H:      curvature of the entrance/exit pole face of the dipole
- 
- Comments:
- Written by Jianfeng Zhang @ LAL, 01/10/2012.
-******************************************************************/
-template<typename T>
-static void BendCurvature(double irho,  double H, ss_vect<T> &x) {
-  
-  if(trace)
-  cout << "Forest's bend curvature...."<<endl;
-  
-  T pm2,u, coeff1, coeff2,coeff3;
-  
-   ss_vect<T> x0;
-   
-   x0 = x;
-  
-   pm2 = sqr(1+x0[delta_]) - sqr(x0[px_]); 
-   u = irho*H/2/pm2;
-   
-   coeff1 = u * 2*x0[px_]*(1+x0[delta_])/pm2; //
-   coeff2 = u * (1+x0[delta_]); //
-   coeff3 = u * (pm2 - 2*(1+x0[delta_]))/pm2; 
-   
- 
-   x[x_]  += x0[x_]/(1-coeff1*x0[y_]*x0[y_]);
-   x[px_] -= coeff2*x0[y_]*x0[y_] - irho*H/2*sqr(x0[x_]);
-   x[py_] -= 2*coeff2*x[x_]*x0[y_];
-   x[ct_] -= coeff3*x[x_]*x0[y_]*x0[y_];
-}   
-   
-   
-      
- 
+        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]
+                / (1.0 + x[delta_]);
+    } else {
+        x[px_] += irho * tan(dtor(phi)) * x[x_];
+        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_];
+    }
+}
 
 /****************************************************************
@@ -1176,14 +794,8 @@
 
   Purpose:
-     Rotate from the beam coorinate to the dipole face.
-     This is a pure geometric operation, no physical meaning.
-     
-     See P.307 Eq. (10.26) in 
-     E. Forest "Beam dynamics: A new attitude and framework".
-     
+     the effect of edge focusing of the dipole
   Input:
      phi:  entrance or exit angle of the dipole edge
      x:    initial cooridinate of the particle
-     
   Output:
   
@@ -1207,5 +819,4 @@
 }
 
-
 /***************************************************************
 template<typename T>
@@ -1213,10 +824,6 @@
 
   Purpose:
-     the effect of longitudinal fringe field using exact Hamitonian.
-     This is a hard effect model, the fringe field and edge focus
-     happens at the entrance/exit point of the dipole pole face.
-
-      include only the first order of irho.
-     Input:
+     the effect of longitudinal fringe field using exact Hamitonian 
+  Input:
   
   Output:
@@ -1242,84 +849,8 @@
         x[ct_] = x1[ct_] - coeff * x1[px_] * sqr(x[y_]) * (1.0 + x1[delta_])
                 / ps3;
-    } else {//give the error message and set the unstable x value during the tracking
-             //this x value will be discard as an unstable value during the tracking
+    } else {
         printf("bend_fringe: *** Speed of light exceeded!\n");
-        x[y_] = 10;
-	x[x_] = 10;
-	x[py_] =10;
-	x[ct_] = 10;
-	//exit_(1);
-    }
-}
-
-/***************************************************************
-template<typename T>
-void bend_fringe(double hb, double gap, ss_vect<T> &x)
-
-  Purpose:
-     the effect of longitudinal fringe field using exact Hamitonian.
-     This is a hard effect model, the fringe field and edge focus
-     happens at the entrance/exit point of the dipole pole face.
-
-      include the second order of irho.
-      
-      See E. Forest and et al. 
-      "Fringe field in MAD Part I: Second order Fringe in MAD-X for the module PTC",
-      P. 8-9, Eq. (35) and (34).
- 
-   Input:
-  
-  Output:
-  
-  Return:
-  
-****************************************************************/
-template<typename T>
-void bend_fringe(double hb, double gap, ss_vect<T> &x) {
-  
-  bool track=false;
-  
-  if(track) cout << "bend_fringe(): Forest's model"<<endl;
-  
-    T coeff1,coeff2,coeff3,coeff4;
-    T ps, ps2, ps3, ps5;
-    
-    ss_vect<T> x1;
- 
-    gap = gap*0.5;
-    
-    x1 = x;
-    ps = get_p_s(x);
-    
-    if(ps==0)
-      return;
-    
-    ps2 = sqr(ps);
-    ps3 = ps * ps2;
-    ps5 = ps2*ps3; 
-    
-    coeff1 = 1.0 - 4.0*sqr(hb)*gap*x1[py_]*x1[y_]/ps3;
-    coeff2 = sqr(hb)*gap*2*x1[px_]*(2*sqr(x1[px_])-sqr(1.0+x1[delta_]))/ps5 +  
-            hb*(ps/(sqr(1+x1[delta_])-sqr(x1[px_])) + 2*sqr(x1[px_])*ps/sqr(sqr(1+x1[delta_])-sqr(x1[px_])));
-             
-    coeff3 = hb*(x1[px_]/ps)/(1+sqr(x1[py_])/ps2) - sqr(hb)*gap*((ps2+sqr(x1[px_]))/ps3+sqr(x1[px_])/ps2*(ps2+x1[py_])/ps3);
-    coeff4 = -sqr(hb)*gap*4*sqr(x1[px_])*(1+x1[delta_])/ps5;
-        
-    if (coeff1 >= 0.0) {
-        x[y_] = 2.0 * x1[y_] / (1.0 + sqrt(coeff1));
-    } else {
-      //give the error message and set the unstable x value during the tracking
-             //this x value will be discard as an unstable value during the tracking
-        printf("bend_fringe: *** Speed of light exceeded!\n");
-        x[y_] = 10;
-	x[x_] = 10;
-	x[py_] =10;
-	x[ct_] = 10;
-     // exit_(1);
-    }
-
-  x[x_] += 0.5*coeff2*sqr(x[y_]);
-  x[py_]-= coeff3*x[y_];
-  x[ct_] -= 0.5*coeff4*sqr(x[y_]);
+        exit_(1);
+    }
 }
 
@@ -1450,7 +981,4 @@
     elemtype *elemp;
     MpoleType *M;
-    
-    ss_vect<T> x1, x2;
-    
 
     elemp = &Cell.Elem;
@@ -1460,5 +988,4 @@
     GtoL(x, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1);
 
-    //entrance of the magnet
     if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) { /* fringe fields */
 
@@ -1468,23 +995,15 @@
         }
 
-        //fringe field and edge focusing of dipoles
-       if (M->Pirho != 0.0 && M->dipEdge_effect1 == 1){
-         if (!globval.H_exact) { //big ring */ modified linear term T43 on page 117 & 118 in SLAC-75.
-	               EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x,true);
-	   } else {//small and big rings; Forest's model
-	  if(M->PH1!=0){
-	    // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0.
-	   BendCurvature(M->Pirho, M->PH1, x); 
-	  }
-	   p_rot(M->PTx1,  x);  //rotate from cartesian cooridnate to curlinear curved beam coordinate; 
-             // since the map of a sector dipole is used in Tracy 3 
-	   bend_fringe(M->Pirho, M->Pgap, x); //fringe field
-         }
-       }
+        if (!globval.H_exact) {
+            if (M->Pirho != 0.0)
+                EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x);
+        } else {
+            p_rot(M->PTx1, x);
+            bend_fringe(M->Pirho, x);
+        }
     }
 
     if (M->Pthick == thick) {
-      if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates
-      // if (M->n_design == Dip || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates
+        if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates
             h_ref = M->Pirho;
             dL = elemp->PL / M->PN;
@@ -1494,5 +1013,5 @@
                 dL = elemp->PL / M->PN;
             else
-                dL = 1.0 / M->Pirho * (sin(dtor(M->PTx1)) + sin(dtor(M->PTx2)))  //straight length of the dipole
+                dL = 1.0 / M->Pirho * (sin(dtor(M->PTx1)) + sin(dtor(M->PTx2)))
                         / M->PN;
         }
@@ -1522,23 +1041,9 @@
 
     case Meth_Second:
-      
-      //  cout << "Mpole_Pass: Meth_Second not supported" << endl;
-       // exit_(0);
-       // break;
-
-    //specially for the test of the dipole map,see Forest's beam dynamics, p361. eqn.(12.19)
- 
-    dL1 = 0.5*dL;
-    dL2 = dL;
-           
-            for (seg = 1; seg <= M->PN; seg++) {
-	      dipole_kick(dL1,M->Pirho,h_ref,x);
-	      multipole_kick(M->Porder, M->PB, dkL1, M->Pirho, x);
-	      dipole_kick(dL1,M->Pirho,h_ref,x);
-            }
-
-            
-       
-    case Meth_Fourth:  //forth order symplectic integrator
+        cout << "Mpole_Pass: Meth_Second not supported" << endl;
+        exit_(0);
+        break;
+
+    case Meth_Fourth:
         if (M->Pthick == thick) {
             dL1 = c_1 * dL;
@@ -1556,8 +1061,8 @@
                 }
 
-		 Drift(dL1, x);
+                Drift(dL1, x);
                 thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
                 Drift(dL2, x);
-		thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
+                thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
 
                 if (globval.emittance && (!globval.Cavity_on) && (M->Pirho
@@ -1566,9 +1071,8 @@
                     dI4 += 4.0 * is_tps<tps>::get_dI4(x);
                 }
-                
-                
+
                 Drift(dL2, x);
-		thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
-		Drift(dL1, x);
+                thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
+                Drift(dL1, x);
 
                 if (globval.emittance && (!globval.Cavity_on) && (M->Pirho
@@ -1594,21 +1098,13 @@
     }
 
-    //exit of the magnets
     if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) {
-       
-      /* dipole fringe fields */
-             if (M->Pirho != 0.0 && M->dipEdge_effect2 == 1){
-	    if (!globval.H_exact) { //big ring, linear model correction
-	                       EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x,false);
-	       }else{
-             bend_fringe(-M->Pirho, M->Pgap,x); //Forest's fringe field of the dipole
-	      p_rot(M->PTx2, x); //rotate back to the cartesian cooridinate
-	     if(M->PH2!=0){
-            // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0.
-             BendCurvature(M->Pirho, M->PH2, x);
-	     }
-         }
-     }
-        //quadrupole fringe field
+        /* fringe fields */
+        if (!globval.H_exact) {
+            if (M->Pirho != 0.0)
+                EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x);
+        } else {
+            bend_fringe(-M->Pirho, x);
+            p_rot(M->PTx2, x);
+        }
         if (globval.quad_fringe && (M->PB[Quad + HOMmax] != 0.0) && (M->quadFF2
                 == 1))
@@ -2916,6 +2412,7 @@
     mergeto4by5(M, ah, av);
 }
-  /* **************************************************************  
-   * Compute transfert matrix for a quadrupole magnet
+
+void Mpole_Setmatrix(int Fnum1, int Knum1, double K) {
+    /*   Compute transfert matrix for a quadrupole magnet
      the transfert matrix A is plitted into two part
      A = AD55xAU55
@@ -2929,8 +2426,5 @@
      corresponding to a half magnet w/ an entrance
      angle and no exit angle.
-     The linear fringe field is taken into account                
-     **************************************************************/
-void Mpole_Setmatrix(int Fnum1, int Knum1, double K) {
-  
+     The linear frindge field is taken into account                    */
 
     CellType *cellp;
@@ -3426,5 +2920,5 @@
         printf("Elem_GetPos: there are no kids in family %d (%s)\n", Fnum1,
                 ElemFam[Fnum1 - 1].ElemF.PName);
-        exit_(1);
+        exit_(0);
     }
 
@@ -3457,8 +2951,4 @@
     d_2 = 1e0 - 2e0 * d_1;
 
-    
-    
-    
-    
     // classical radiation
     //  C_gamma = 8.846056192e-05;
@@ -3853,6 +3343,4 @@
     M->PTx1 = 0.0; /* Entrance angle */
     M->PTx2 = 0.0; /* Exit angle */
-    M->PH1 = 0.0; /* Entrance pole face curvature*/
-    M->PH2 = 0.0; /* Exit pole face curvature */
     M->Pgap = 0.0; /* Gap for fringe field ??? */
 
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/t2lat.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/t2lat.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/t2lat.cc	(revision 23)
@@ -44,5 +44,5 @@
   leq, pwrsym, lss, geo, gtr, period_, charcon, stringcon, ident, geq, lsym,
   bobrhosym, bobrhovsym, bobrhohsym, kxvsym, kxhsym, phisym, ksym,
-  tsym, t1sym, t2sym,h1sym,h2sym,
+  tsym, t1sym, t2sym,
   gapsym, thksym, invsym, thnsym,
   endsym, tsksym, bemsym, corsym, prnsym, tblsym, possym, prmsym,
@@ -51,10 +51,9 @@
   nbdsym, frgsym, latsym, mpsym, dbnsym, kssym, homsym, lmdsym, dtsym, xytsym,
   vrfsym, harnumsym, frqsym, gstsym, typsym, rollsym, idsym,
-  fnamesym1, fnamesym2, scalingsym1, scalingsym2, fmsym, harmsym, sprsym, recsym, solsym,edge_effect1sym,edge_effect2sym,
+  fnamesym1, fnamesym2, scalingsym1, scalingsym2, fmsym, harmsym, sprsym, recsym, solsym,
   ff1sym, ff2sym, ffscalingsym, tiltsym
 } Lat_symbol;   /*\*/
 // idsym fnamesym1 fnamesym2 scalingsym added for insertion
 // ring sym added
-  //edge_effect1sym, edge_effect2sym added, to turn on/off the edge effect of dipoles @ 05/10/2012 by Jianfeng Zhang
 // ff1sym ff2sym for quadrupole entrance and exit fringe field
 // ffscalingsym scaling factor for entrance and exit fringe field.  /*J.Zhang
@@ -658,8 +657,4 @@
       if (!strncmp(id, "t              ", sizeof(alfa_)))
 	*sym = tsym;
-      else if (!strncmp(id, "h1            ", sizeof(alfa_)))
-	*sym = h1sym;
-      else if (!strncmp(id, "h2            ", sizeof(alfa_)))
-	*sym = h2sym;
       else if (!strncmp(id, "gap            ", sizeof(alfa_)))
 	*sym = gapsym;
@@ -1059,8 +1054,4 @@
 	else if (!strncmp(fname, "t2             ", sizeof(partsName)))
 	  x = WITH1->PTx2;
-	else if (!strncmp(fname, "h1            ", sizeof(partsName)))
-	  x = WITH1->PH1;
-	else if (!strncmp(fname, "h2            ", sizeof(partsName)))
-	  x = WITH1->PH2;
 	else if (!strncmp(fname, "gap            ", sizeof(partsName)))
 	  x = WITH1->Pgap;
@@ -2026,15 +2017,11 @@
   struct LOC_Lat_DealElement V;
   bool           Result = false;
-  double         t=0.0, t1=0.0, t2=0.0, h1=0,h2=0,gap=0.0, QL = 0.0, QK=0.0;
-  double         QKick=0.0;  /* kick angle of the corrector [rad]*/
-  int            Kplane=0; /* kick plane of the corrector, 1 for H plane, -1 for V plane*/
-  double         QKV=0.0, QKH=0.0, QKxV=0.0, QKxH=0.0, QPhi=0.0, QKS=0.0;
-  double         dt=0.0, Frf=0.0, Vrf=0.0;
-  long            harnum=0L;
-  int            k1=1; //number of cut slice of the element
-  int            k2=4; //integration method; Meth_Linear=0, Meth_First=1,Meth_Second=2,Meth_Fourth=4
-  int            edge_effect1=0,edge_effect2=0; //flag to turn on/off dipole fringe field
-  int            FF1=0, FF2=0;
-  double         FFscaling=0.0;
+  double         t, t1, t2, gap, QL = 0.0, QK;
+  double         QKick;  /* kick angle of the corrector [rad]*/
+  int            Kplane; /* kick plane of the corrector, 1 for H plane, -1 for V plane*/
+  double         QKV, QKH, QKxV, QKxH, QPhi, QKS;
+  double         dt, Frf, Vrf;
+  long k1, k2, harnum, FF1, FF2;
+  double FFscaling;
   Lat_symbol     sym1;
   symset         mysys, SET;
@@ -2089,5 +2076,4 @@
     test__(P_expset(SET, 1 << ((long)semicolon)), "<;> expected", &V);
     GetSym__(&V);
-    QL = *V.rnum; 
     globval.Elem_nFam++;
     if (globval.Elem_nFam <= Elem_nFamMax) {
@@ -2095,6 +2081,5 @@
       WITH1 = &WITH->ElemF;
       memcpy(WITH1->PName, ElementName, sizeof(partsName));
-   //   WITH1->PL = *V.rnum;
-       WITH1->PL = QL;
+      WITH1->PL = *V.rnum;
       WITH1->Pkind = PartsKind(drift);
       Drift_Alloc(&WITH->ElemF);
@@ -2116,8 +2101,4 @@
             T2     = <exit angle>, ( [degree] )
             gap    = <total magnet gap>, ( [m] )
-   edge_effect1    = 1/0, turn on/ff the entrance edge effect, 
-                     and fringe field (if gap != 0)
-   edge_effect2    = 1/0, turn on/ff the exit edge effect, 
-                     and fringe field (if gap != 0)
             K      = <K-value>, ( [m^-2] )
                      ( K > 0 : focusing in horizontal plane )
@@ -2135,5 +2116,5 @@
     Example
 
-      B: bending, L=0.70, T=10.0, T1:=5.0, T2:=5.0, H1=0, H2=0, gap=0.03, edge_effect1=1, edge_effect2=1, K=-1.0, N=8, Method=2;
+      B: bending, L=0.70, T=10.0, T1:=5.0, T2:=5.0, K=-1.0, N=8, Method=2;
 
     *************************************************************************/
@@ -2143,14 +2124,10 @@
     QL = 0.0;   /* L */
     QK = 0.0;   /* K, quadrupole components*/
-    k1 = 1;   /* N */
+    k1 = 0;   /* N */
     t  = 0.0;   /* T */
     t1 = 0.0;   /* T1 */
     t2 = 0.0;   /* T2 */
-    h1 = 0.0;
-    h2 = 0.0;
     gap = 0.0;   /* gap */
-    edge_effect1 = 0; /* edge_effect effect at the entrance */ 
-    edge_effect2 = 0; /* edge effect at the exit */ 
-    k2 = Meth_Fourth;   /* method */
+    k2 = Meth_Linear;   /* method */
     dt = 0.0;
     ClearHOMandDBN(&V);
@@ -2163,9 +2140,5 @@
     P_addset(mysys, (long)t1sym);
     P_addset(mysys, (long)t2sym);
-    P_addset(mysys, (long)h1sym);
-    P_addset(mysys, (long)h2sym);
     P_addset(mysys, (long)gapsym);
-    P_addset(mysys, (long)edge_effect1sym);
-    P_addset(mysys, (long)edge_effect2sym);
     P_addset(mysys, (long)homsym);
     P_addset(mysys, (long)dbnsym);
@@ -2205,22 +2178,6 @@
 	break;
 
-	case h1sym:
-	h1 = EVAL_(&V);
-	break;
-
-      case h2sym:
-	h2 = EVAL_(&V);
-	break;
-	
       case gapsym:
 	gap = EVAL_(&V);
-	break;
-
-      case edge_effect1sym:
-	edge_effect1 = (long)EVAL_(&V);
-	break;
-
-      case edge_effect2sym:
-	edge_effect2 = (long)EVAL_(&V);
 	break;
 
@@ -2266,9 +2223,5 @@
       else
 	WITH2->Pirho = t * M_PI / 180.0;
-      WITH2->PTx1 = t1; WITH2->PTx2 = t2; //entrance and exit angle of the dipoles (relative to rectangular magnets)
-      WITH2->PH1 = h1; WITH2->PH2 = h2; //entrance and exit angle of the dipoles
-      WITH2->Pgap = gap;
-      WITH2->dipEdge_effect1 = edge_effect1;//edge_effect effect at the entrance
-      WITH2->dipEdge_effect2 = edge_effect2;//edge effect at the exit
+      WITH2->PTx1 = t1; WITH2->PTx2 = t2; WITH2->Pgap = gap;
       WITH2->n_design = Dip;
       AssignHOM(globval.Elem_nFam, &V);
@@ -2302,5 +2255,5 @@
 
       Q1: Quadrupole, L=0.5, K=2.213455, N=4, Method=4;
-      Q1 : Quadrupole, L=0.5, K=2.213455, N=4, FF1 =1, FF2=0, FFscaling=1, Method=4;
+      Q1 : Quadrupole, L=0.5, K=2.213455, N=4, FF1 =1, FF2, FFscaling=1, Method=4;
 
     **************************************************************/
@@ -2310,5 +2263,5 @@
     QL = 0.0;   /* L */
     QK = 0.0;   /* K */
-    k1 = 1;   /* N */
+    k1 = 0;   /* N */
     k2 = Meth_Fourth;   /* method */
     dt = 0.0;
@@ -2438,5 +2391,5 @@
     QL = 0.0;   /* L */
     QK = 0.0;   /* K */
-    k1 = 1;   /* N */
+    k1 = 0;   /* N */
     k2 = Meth_Fourth;   /* method */
     FF1 = 0;     /* Entrance Fringe field */
@@ -2653,5 +2606,5 @@
     QKick = 0.0; /* kick angle of the corrector [rad]*/
     Kplane = 0;   /* 1 is horizontal corrector, -1 is vertical corrector */
-    k1 = 1;     /* N */
+    k1 = 0;     /* N */
     k2 = Meth_Linear;   /* method */
     dt = 0.0;
@@ -2931,12 +2884,10 @@
     QL = 0.0;   /* L */
     QK = 0.0;   /* K */
-    k1 = 1;   /* N */
+    k1 = 0;   /* N */
     t = 0.0;   /* T */
     t1 = 0.0;   /* T1 */
     t2 = 0.0;   /* T2 */
-    h1 = 0.0;
-    h2 = 0.0;
     gap = 0.0;   /* gap */
-    k2 = Meth_Fourth;   /* method */
+    k2 = Meth_Linear;   /* method */
     dt = 0.0;
     ClearHOMandDBN(&V);
@@ -2947,6 +2898,4 @@
     P_addset(mysys, (long)t1sym);
     P_addset(mysys, (long)t2sym);
-    P_addset(mysys, (long)h1sym);
-    P_addset(mysys, (long)h2sym);
     P_addset(mysys, (long)gapsym);
     P_addset(mysys, (long)rollsym);
@@ -2983,12 +2932,4 @@
 	break;
 
-      case h1sym:
-	h1 = EVAL_(&V);
-	break;
-
-      case h2sym:
-	h2 = EVAL_(&V);
-	break;
-	
       case gapsym:
 	gap = EVAL_(&V);
@@ -3038,7 +2979,5 @@
       }
       WITH2->PN = k1; WITH2->Pmethod = k2;
-      WITH2->PTx1 = t1; WITH2->PTx2 = t2;
-      WITH2->PH1 = h1; WITH2->PH2 = h2;
-      WITH2->Pgap = gap;
+      WITH2->PTx1 = t1; WITH2->PTx2 = t2; WITH2->Pgap = gap;
       WITH2->PdTpar = dt;
       AssignHOM(globval.Elem_nFam, &V);
@@ -3083,15 +3022,6 @@
     getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
     GetSym__(&V);
-    QL = 0e0;
-    QK = 0e0; 
-    QKV = 0e0; 
-    QKH = 0e0; 
-    QKxV = 0e0; 
-    QKxH = 0e0;
-    QPhi = 0e0; 
-    QKS = 0e0; 
-    k1 = 1; 
-    k2 = Meth_Fourth; 
-    dt = 0e0;
+    QL = 0e0; QK = 0e0; QKV = 0e0; QKH = 0e0; QKxV = 0e0; QKxH = 0e0;
+    QPhi = 0e0; QKS = 0e0; k1 = 0; k2 = Meth_Linear; dt = 0e0;
     ClearHOMandDBN(&V);
     P_addset(P_expset(mysys, 0), (long)lsym);
@@ -3219,8 +3149,5 @@
     getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
     GetSym__(&V);
-    QL = 0.0; 
-    k1 = 1; 
-    strcpy(str1, ""); 
-    strcpy(str2, "");
+    QL = 0.0; k1 = 0; strcpy(str1, ""); strcpy(str2, "");
     P_addset(P_expset(mysys, 0), (long)lsym);
     P_addset(mysys, (long)nsym);
@@ -3572,5 +3499,5 @@
     QL = 0.0; /* L */
     QK = 0.0; /* K */
-    k1 = 1;   /* N */
+    k1 = 0;   /* N */
     P_addset(P_expset(mysys, 0), (long)lsym);
     P_addset(mysys, (long)bobrhosym);
@@ -3731,6 +3658,4 @@
   Reg("drift          ", drfsym, &V);
   Reg("dt             ", dtsym, &V);
-  Reg("edge_effect1   ", edge_effect1sym, &V);     /* Jianfeng Zhang*/
-  Reg("edge_effect2   ", edge_effect2sym, &V);     /* Jianfeng Zhang */
   Reg("end            ", endsym, &V);
   Reg("ff1            ", ff1sym, &V);     /* Laurent */
@@ -3746,6 +3671,4 @@
   Reg("gap            ", gapsym, &V);
   Reg("ghost          ", gstsym, &V);
-  Reg("h1             ", h1sym, &V);
-  Reg("h2             ", h2sym, &V);
   Reg("harm           ", harmsym, &V);
   Reg("harnum         ", harnumsym, &V);
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/tpsa_lin.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/tpsa_lin.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/tpsa_lin.cc	(revision 23)
@@ -414,7 +414,5 @@
 }
 
-/****************************
- * tan(x)
- * *************************/
+
 void datan(const tps_buf &x, tps_buf &z)
 {
@@ -424,88 +422,9 @@
 }
 
-/*************************
- *  cotan(x)
- * ***********************/
-void dactan(const tps_buf &x, tps_buf &z)
-{
-  tps_buf  c, s;
-
-  dacos(x, c); dasin(x, s); dadiv_(c, s, z);
-}
-
-/*************************************
- * 
- * Purpose:
- *     Get the derivatives of the arcsin(x)
- *    
- *    d(arcsin(u))/dx = (1/sqrt(1-u^2))*du/dx;
- * 
- *  Input:
- *         x: u
- *  Output: 
- *        z: the differential of x.
- *
- *       
- * Comments:
- *           13/11/2012
- *          Written by Jianfeng Nadolski @ LAL 
- *                      
- * ***********************************/
-void daarcsin(const tps_buf &x, tps_buf &z)
-{
-  int     i;
-  double  a;  //differential of x; a = sqrt(1/1-u^2)
-
-  a = x[0]; z[0] = asin(a); a = 1.0/sqrt(1.0-sqr(a));
-  for (i = 1; i <= ss_dim; i++)
-    z[i] = a*x[i];
-}
-
-/*************************************
- * 
- * Purpose:
- *     Get the derivatives of the arccos(x)
- *    
- *    d(arccos(u))/dx = -(1/sqrt(1-u^2))*du/dx;
- * 
- *  Input:
- *         x: u
- *  Output: 
- *        z: the differential of x.
- *
- *       
- * Comments:
- *           13/11/2012
- *          Written by Jianfeng Nadolski @ LAL 
- *                      
- * ***********************************/
-void daarccos(const tps_buf &x, tps_buf &z)
-{
-int     i;
-  double  a;  //differential of x; a = -sqrt(1/1-u^2)
-
-  a = x[0]; z[0] = asin(a); a = -1.0/sqrt(1.0-sqr(a));
-  for (i = 1; i <= ss_dim; i++)
-    z[i] = a*x[i];
-  
-}
-/*************************************
- * 
- * Purpose:
- *     Get the derivatives of the arctan(x)
- *    
- *    d(arctan(u))/dx = (1/1+u^2)*du/dx;
- * 
- *  Input:
- *         x: u
- *  Output: 
- *        z: the differential of x.
- *
- *         
- * ***********************************/
+
 void daarctan(const tps_buf &x, tps_buf &z)
 {
   int     i;
-  double  a;  //differential of x; a = (1/1+u^2)
+  double  a;
 
   a = x[0]; z[0] = atan(a); a = 1.0/(1.0+sqr(a));
@@ -514,32 +433,9 @@
 }
 
-/*************************************
- * 
- * Purpose:
- *     Get the derivatives of the arcctan(x)
- *    
- *    d(arcctan(u))/dx = -(1/1+u^2)*du/dx;
- * 
- *  Input:
- *         x: u
- *  Output: 
- *        z: the differential of x.
- *
- *         
- * ***********************************/
-void daarcctan(const tps_buf &x, tps_buf &z)
-{
-  int     i;
-  double  a;  //differential of x; a = -(1/1+u^2)
-
-  a = x[0]; z[0] = atan(a); a = -1.0/(1.0+sqr(a));
-  for (i = 1; i <= ss_dim; i++)
-    z[i] = a*x[i];
-}
 
 void dafun_(const char *fun, const tps_buf &x, tps_buf &z)
 {
   tps_buf  u;
-  
+
   if (!strncmp(fun, "INV ", dafunlen))
     dainv_(x, u);
@@ -554,20 +450,12 @@
   else if (!strncmp(fun, "COS ", dafunlen))
     dacos(x, u);
-  else if (!strncmp(fun, "TAN ", dafunlen))
-    datan(x, u);
-  else if (!strncmp(fun, "CTAN ", dafunlen))
-    dactan(x, u);
   else if (!strncmp(fun, "SINH ", dafunlen))
     dasinh(x, u);
   else if (!strncmp(fun, "COSH ", dafunlen))
     dacosh(x, u);
-  else if (!strncmp(fun, "ASIN", dafunlen))
-    daarcsin(x, u);
-  else if (!strncmp(fun, "ACOS", dafunlen))
-    daarccos(x, u);
+  else if (!strncmp(fun, "TAN ", dafunlen))
+    datan(x, u);
   else if (!strncmp(fun, "ATAN", dafunlen))
     daarctan(x, u);
-  else if (!strncmp(fun, "ACTAN", dafunlen))
-    daarcctan(x, u);
   else {
     printf("dafun: illegal function %s\n", fun);
Index: /branches/tracy3-3.10.1b/tracy/tracy/src/tracy.cc
===================================================================
--- /branches/tracy3-3.10.1b/tracy/tracy/src/tracy.cc	(revision 22)
+++ /branches/tracy3-3.10.1b/tracy/tracy/src/tracy.cc	(revision 23)
@@ -97,8 +97,4 @@
 template void radiate(ss_vect<tps> &, const double, const double, const tps[]);
 
-template void Drift(double, double, ss_vect<double> &);
-
-template void Drift(double, double, ss_vect<tps> &);
-
 template void Drift(double, ss_vect<double> &);
 
@@ -109,15 +105,7 @@
 template void bend_fringe(double, ss_vect<tps> &);
 
-template void bend_fringe(double, double, ss_vect<double> &);
-
-template void bend_fringe(double, double, ss_vect<tps> &);
-
-template static void BendCurvature(double, double, ss_vect<double> &);
-
-template static void BendCurvature(double, double, ss_vect<tps> &);
-
-template static void EdgeFocus(double, double, double, ss_vect<double> &, bool);
-
-template static void EdgeFocus(double, double, double, ss_vect<tps> &,bool);
+template static void EdgeFocus(double, double, double, ss_vect<double> &);
+
+template static void EdgeFocus(double, double, double, ss_vect<tps> &);
 
 template void quad_fringe(double, ss_vect<double> &);
@@ -128,10 +116,4 @@
 
 template void Drift_Pass(CellType &, ss_vect<tps> &);
-
-
-template void dipole_kick(double, double, double, ss_vect<double> &);
-template void dipole_kick(double, double, double, ss_vect<tps> &);
-template void multipole_kick(int, double[], double, double, ss_vect<double> &);
-template void multipole_kick(int, double[], double, double, ss_vect<tps> &);
 
 template void thin_kick(int, double[], double, double, double,
