diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index aded39b218c91de605796e487abed050467885de..22e7101fe73f1169d354a52415b7395c475bd067 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -5,8 +5,7 @@ int no_tps = ORDER; // arbitrary TPSA order is defined locally
 extern bool  freq_map;
 
 #include "tracy_lib.h"
-#include "time.h"
-#include <sys/time.h>
+
 
 //***************************************************************************************
 //
@@ -16,29 +15,30 @@ extern bool  freq_map;
  int main(int argc, char *argv[])
 {
   const long  seed = 1121;
-  const bool All = TRUE;
   iniranf(seed); setrancut(2.0);
 
   // turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
   // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
   globval.H_exact     = false; 
-  globval.quad_fringe = false;                // quadrupole fringe field
-  
-  globval.radiation   = false;                // synchrotron radiation
-  globval.emittance   = false;                 // emittance
+  globval.quad_fringe = false;              // quadrupole fringe field
+  globval.radiation   = false;             // synchrotron radiation
+  globval.emittance   = false;             // emittance
   globval.pathlength  = false; 
 //  globval.bpm         = 0;
   
-  const double  x_max_FMA = 10e-3, delta_FMA = 10e-2;
-  const int     n_x = 80, n_dp = 80, n_tr = 2048;
+//  const double  x_max_FMA = 10e-3,  delta_FMA = 10e-2;
+//  const int     n_x = 801, n_dp = 80, n_tr = 2048;
   
   double nux=0, nuz=0, ksix=0, ksiz=0;
-  
+
   bool chroma;
   double dP = 0.0;
   long lastpos = -1L;
   char str1[S_SIZE];
-  
+
+  // for time handling
+  uint32_t         start, stop;
+
   /************************************************************************
       start read in files and flags
   *************************************************************************/
@@ -74,7 +74,10 @@ extern bool  freq_map;
 
   // compute chromaticities by tracking (should be the same as by DA)
   if (ChromTracFlag == true){
-    GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
+	start = stampstart();
+	GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
+	stop = stampstop(start);
+
     fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
   }
 
@@ -191,35 +194,12 @@ extern bool  freq_map;
   
 
   if (PhaseSpaceFlag == true){
-     int   t_sec,t_usec,t_msec; 
-     struct timeval tv,tv1;
-     struct timezone tz,tz1 ;
-     struct tm *tm,*tm1;
-     gettimeofday(&tv, &tz);
-     tm=localtime(&tv.tv_sec);
-     
-     t_sec = - tm->tm_sec;
-     printf(" %d:%02d:%02d %d \n", tm->tm_hour, tm->tm_min,
-             tm->tm_sec, tv.tv_usec);
-
-      Phase(0.001,0,0.001,0,0.001,0, 1);
-      
-       gettimeofday(&tv1, &tz1);
-       tm1=localtime(&tv1.tv_sec);
-     
-       t_sec += tm->tm_sec;
-       t_usec = tv1.tv_usec -  tv.tv_usec; 
-       t_msec = t_sec*1000 + t_usec/1000;
-         
-       
-      printf(" %d:%02d:%02d %d \n", tm1->tm_hour, tm1->tm_min,
-             tm1->tm_sec, tv1.tv_usec);
-      
-     
-       printf("%d %d %d \n",t_sec*1000,t_usec/1000, t_msec);
-     printf("the simulation time for phase space in tracy 3 is %ld milliseconds\n",t_msec);  
-   
+		start = stampstart();
+		Phase(0.001,0,0.001,0,0.001,0, 1);
+	    printf("the simulation time for phase space in tracy 3 is \n");
+		stop = stampstop(start);
   }
+  return 0;
 }    
-  
-  
\ No newline at end of file
+
+
diff --git a/tracy/tracy/inc/field.h b/tracy/tracy/inc/field.h
index 78e5efdf2bd7b88c6ee62403f8341505d20c6c48..749b3a5dc7d968b1187c96a5ae82cd2453a5a2e6 100644
--- a/tracy/tracy/inc/field.h
+++ b/tracy/tracy/inc/field.h
@@ -11,7 +11,7 @@ const int  n_m2    = 21;  // no of 2nd moments
 const int  ss_dim  = 6;   // state space dimension
 #else
 const int  ss_dim  = 6+1; /* state space dimension:
-			     phase space and parameter dependance */
+			     phase space and parameter dependence */
 #endif
 
 // spatial components
diff --git a/tracy/tracy/inc/physlib.h b/tracy/tracy/inc/physlib.h
index c8f322c851bb83ac54327d41a4544d7f2a3f0c3e..e300f75ecfd94c860ae7d923f654b558f8d4b05f 100644
--- a/tracy/tracy/inc/physlib.h
+++ b/tracy/tracy/inc/physlib.h
@@ -8,6 +8,9 @@
 
 */
 
+#include "time.h"
+#include <sys/time.h>
+
 /* For tune fitting */
 #define nueps           1e-6      //precision
 #define nudkL           0.01   //step
diff --git a/tracy/tracy/inc/tracy.h b/tracy/tracy/inc/tracy.h
index ae8a5aa08e0650b9eb58bb24e9709e5799fce946..4a8b350ebc2efe6bafb6c9bfaa434caa03f593b6 100644
--- a/tracy/tracy/inc/tracy.h
+++ b/tracy/tracy/inc/tracy.h
@@ -8,10 +8,16 @@
 
 */
 
+// Intel compiler
+#pragma warning (disable:193)
+#pragma warning (disable:981)
+#pragma warning (disable:1572)
+// LSN
 
 #define DOF    (ss_dim/2)
- 
-#define M_PI   3.14159265358979323846  // pi
+
+// defined in math.h
+//#define M_PI   3.14159265358979323846  // pi
 
 #ifndef LONG_MAX
 # define LONG_MAX   ((long)(((unsigned long) -1) >> 1))
@@ -53,7 +59,7 @@
 
 
 const double  c0    = 2.99792458e8;             // speed of light in vacuum
-const double  q_e   = 1.602e-19;                // electron charge
+const double  q_e   = 1.60217733e-19;                // electron charge
 const double  m_e   = 0.51099906e6;             // electron rest mass [eV]
 const double  mu_0  = 4.0*M_PI*1e-7;            // permittivity of free space
 const double  eps_0 = 1.0/(sqr(c0)*mu_0);       // permeability of free space
@@ -101,8 +107,8 @@ typedef double mpolArray[HOMmax+HOMmax+1];
 typedef struct statusrec{
   bool tuneflag, chromflag, codflag, mapflag, passflag, overflag, chambre;
   int lossplane; /* lost in: horizontal    1
-		             vertical      2
-			     longitudinal  3 */
+		                     vertical      2
+			                 longitudinal  3 */
 } statusrec;
 
 
@@ -146,4 +152,6 @@ extern void lsoc(int niter, int bpm, int corr, int plane);
 
 /**** same as asctime in C without the \n at the end****/
 char *asctime2(const struct tm *timeptr);
-struct tm* GetTime();
+struct tm* GetTime(void);
+uint32_t stampstop(uint32_t start);
+uint32_t stampstart(void);
diff --git a/tracy/tracy/inc/tracy_global.h b/tracy/tracy/inc/tracy_global.h
index cfc6d5157dad79fe1d9f76bf189ee6fe4157f7a9..7f2b37f9a96c6ea328340bdf8260cbdcfd1e488b 100644
--- a/tracy/tracy/inc/tracy_global.h
+++ b/tracy/tracy/inc/tracy_global.h
@@ -73,6 +73,11 @@ struct DriftType {
 struct MpoleType {
   int         Pmethod;   // Integration Method
   int         PN;        // Number of integration steps
+  long quadFF1;         /* Entrance quadrupole Fringe field flag */
+  long quadFF2;         /* Exit quadrupole Fringe field flag */
+  double quadFFscaling;         /* quadrupole Fringe field scaling factor flag */
+  long sextFF1;         /* Entrance sextupole Fringe field flag */
+  long sextFF2;         /* Exit sextupole Fringe field flag */
   // Displacement Errors
   Vector2     PdSsys;    // systematic [m] 
   Vector2     PdSrms;    // rms [m]
diff --git a/tracy/tracy/inc/tracy_lib.h b/tracy/tracy/inc/tracy_lib.h
index 9db68d6139f2ecc52d8d613c40bf3bc1acc515b0..98fc5f09f0e8243f52621806d5e4fdef05866718 100644
--- a/tracy/tracy/inc/tracy_lib.h
+++ b/tracy/tracy/inc/tracy_lib.h
@@ -25,7 +25,9 @@
 #include <setjmp.h>
 #include <time.h>
 #include <memory.h>
+#if !defined(__APPLE__)
 #include <malloc.h>
+#endif
 //#include <execinfo.h>
 
 
diff --git a/tracy/tracy/src/physlib.cc b/tracy/tracy/src/physlib.cc
index 69bdc94edaf5e2554251c3e251439dfa36ca9bc9..f8feb3aba856371a84f253f78db326ed4d1bcd22 100644
--- a/tracy/tracy/src/physlib.cc
+++ b/tracy/tracy/src/physlib.cc
@@ -46,7 +46,122 @@ struct tm* GetTime()
   whattime = localtime(&aclock);  /* Convert time to struct */
   return whattime;
 }
+/****************************************************************************
+ * uint32_t stampstart()
 
+
+   Purpose: record time in millliseconds
+
+   Input:
+       none
+
+   Output:
+       non
+
+   Return:
+       time in milliseconds
+
+   Global variables:
+       none
+
+   specific functions:
+       none
+
+   Comments:
+       to be used with stampstop()
+
+****************************************************************************/
+
+uint32_t stampstart(void)
+{
+	struct timeval  tv;
+	struct timezone tz;
+	struct tm      *tm;
+	uint32_t         start;
+	const bool timedebug = false;
+
+	// get the time
+	gettimeofday(&tv, &tz);
+	tm = localtime(&tv.tv_sec);
+
+	// print detailed time in milliseconds
+    if (timedebug)
+	printf("TIMESTAMP-START\t  %d:%02d:%02d:%d (~%d ms)\n", tm->tm_hour,
+	       tm->tm_min, tm->tm_sec, tv.tv_usec,
+	       tm->tm_hour * 3600 * 1000 + tm->tm_min * 60 * 1000 +
+	       tm->tm_sec * 1000 + tv.tv_usec / 1000);
+
+
+	start = tm->tm_hour * 3600 * 1000 + tm->tm_min * 60 * 1000 +
+		tm->tm_sec * 1000 + tv.tv_usec / 1000;
+
+	return (start);
+
+}
+
+// compute time elapsed since start
+/****************************************************************************
+ * uint32_t stampstop(uint32_t start)
+
+   Purpose: compute time elapsed since start time
+
+   Input:
+       start starting time i millisecond
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   specific functions:
+       none
+
+   Comments:
+       to be used with stampstart
+
+****************************************************************************/
+uint32_t stampstop(uint32_t start)
+{
+
+	struct timeval  tv;
+	struct timezone tz;
+	struct tm      *tm;
+	uint32_t         stop;
+	const bool timedebug = false;
+
+	// get the time
+	gettimeofday(&tv, &tz);
+	tm = localtime(&tv.tv_sec);
+
+	stop = tm->tm_hour * 3600 * 1000 + tm->tm_min * 60 * 1000 +
+		tm->tm_sec * 1000 + tv.tv_usec / 1000;
+
+	// print detailed time in milliseconds
+	if (timedebug){
+	printf("TIMESTAMP-END\t  %d:%02d:%02d:%d (~%d ms) \n", tm->tm_hour,
+	       tm->tm_min, tm->tm_sec, tv.tv_usec,
+	       tm->tm_hour * 3600 * 1000 + tm->tm_min * 60 * 1000 +
+	       tm->tm_sec * 1000 + tv.tv_usec / 1000);
+	printf("ELAPSED\t  %d ms\n", stop - start);
+	}
+
+	uint32_t delta, hour, minute, second, millisecond;
+	delta = stop - start;
+	hour = delta/3600000;
+	minute = (delta-3600000*hour)/60000;
+	second = (delta-3600000*hour-minute*60000)/1000;
+	millisecond = delta-3600000*hour-minute*60000-second*1000;
+
+	printf("ELAPSED\t  %d h %d min %d s %d ms\n",
+			hour, minute, second, millisecond);
+
+	return (stop);
+
+}
 /****************************************************************************/
 /* void printglob(void)
 
diff --git a/tracy/tracy/src/radia2tracy.cc b/tracy/tracy/src/radia2tracy.cc
index cabd39a44bc0167820e3dc947f431be931d17a66..3f65ea812eaf0d22d34948da231b0e935ab5bcfb 100644
--- a/tracy/tracy/src/radia2tracy.cc
+++ b/tracy/tracy/src/radia2tracy.cc
@@ -8,10 +8,18 @@
 
 */
 
-
+/* Spline routine adapted from NR with template */
 template<typename T>
 void spline(const double x[], const T y[], int const n,
 	    double const yp1, const double ypn, T y2[])
+/* Function to be called only once
+Given arrays x[1..n] and y[1..n] containing a tabulated function, i.e., yi = f(xi), with
+x1 <x2 < :: : < xN, and given values yp1 and ypn for the first derivative of the interpolating
+function at points 1 and n, respectively, this routine returns an array y2[1..n] that contains
+the second derivatives of the interpolating function at the tabulated points xi. If yp1 and/or
+ypn are equal to 10e30 or larger, the routine is signaled to set the corresponding boundary
+condition for a natural spline, with zero second derivative on that boundary.
+*/
 {
   int     i, k;
   double  sig;
@@ -42,6 +50,11 @@ void spline(const double x[], const T y[], int const n,
 template<typename T, typename U>
 void splint(const double xa[], const U ya[], const U y2a[],
 	    const int n, const T &x, T &y)
+/* SPLine INTerpolation. Once spline are known, only evaluate polynome
+ * Given the arrays xa[1..n] and ya[1..n], which tabulate a function (with the xai's in order),
+and given the array y2a[1..n], which is the output from spline above, and given a value of
+x, this routine returns a cubic-spline interpolated value y.
+*/
 {
   int     klo, khi, k;
   double  h;
@@ -66,23 +79,34 @@ template<typename T>
 void splin2(const double x1a[], const double x2a[],
 	    double **ya, double **y2a, const int m, const int n,
 	    const T &x1, const T &x2, T &y)
-{
+/* Main function computing the bicubic spline
+Given x1a, x2a, ya, m, n as described in splie2 and y2a as produced by that routine; and
+given a desired interpolating point x1,x2; this routine returns an interpolated function value y
+by bicubic spline interpolation.
+ */{
   int  j;
-  T    ytmp[m+1], yytmp[m+1];
+  T    ytmp[m+1], yytmp[m+1]; // Perform m evaluations of the row splines constructed by
   
-  for (j = 1; j <= m; j++)
-    splint(x2a, ya[j], y2a[j], n, x2, yytmp[j]);
-  spline(x1a, yytmp, m, 1.0e30, 1.0e30, ytmp);
+  for (j = 1; j <= m; j++) //splie2, using the one-dimensional spline evaluator
+    splint(x2a, ya[j], y2a[j], n, x2, yytmp[j]); //splint.
+  spline(x1a, yytmp, m, 1.0e30, 1.0e30, ytmp); // Construct the one-dimensional column spline and evaluate it.
   splint(x1a, yytmp, ytmp, m, x1, y);
 }
 
 
 void splie2(double x1a[], double x2a[], double **ya,
 	    int m, int n, double **y2a)
+/* precompute the auxiliary second-derivative table
+Given an m by n tabulated function ya[1..m][1..n], and tabulated independent variables
+x2a[1..n], this routine constructs one-dimensional natural cubic splines of the rows of ya
+and returns the second-derivatives in the array y2a[1..m][1..n]. (The array x1a[1..m] is
+included in the argument list merely for consistency with routine splin2.)
+*/
 {
   int  j;
 
   for (j = 1; j <= m; j++)
+	  // Values 10e30 signal a natural spline
     spline(x2a, ya[j], n, 1.0e30, 1.0e30, y2a[j]);
 }
 
@@ -225,7 +249,6 @@ void Read_IDfile(char *fic_radia, double *L, int *pnx, int *pnz,
   
   /* Array creation for thetaz */
   for (i = 0; i < nz; i++) {
-//    fscanf(fi, "%*lf");
     fscanf(fi, "%*f");
     for (j = 0; j < nx; j++) {
       fscanf(fi, "%lf", &thetaz[i][j]);
@@ -238,7 +261,7 @@ void Read_IDfile(char *fic_radia, double *L, int *pnx, int *pnz,
     if (traceID) fprintf(stdout, "\n");
   }
   
-  /* For debbuging */
+  /* For debugging */
   if (traceID)
     for (j = 0; j < nx; j++) {
       fprintf(stdout, "tabx[%d] =% lf\n", j, tabx[j]);
@@ -282,7 +305,7 @@ void Read_IDfile(char *fic_radia, double *L, int *pnx, int *pnz,
    Specific functions:
  
    Comments:
-       none
+       Search for index could be speed up easily
  
 ****************************************************************************/
 template<typename T>
@@ -422,7 +445,7 @@ void LinearInterpolation2(T &X, T &Z, T &TX, T &TZ, CellType &Cell,
  
    Input:
        X, Z location of the interpolation
-       Cell elment containing ID device
+       Cell element containing ID device
  
    Output:
        TX, TZ thetax and thetaz interpolated at X and Z
@@ -453,7 +476,7 @@ void SplineInterpolation2(T &X, T &Z, T &thetax, T &thetaz,
     /* test whether X and Z within the transverse map area */
     if (X < WITH->tabx[0] || X > WITH->tabx[nx-1] ||
 	Z > WITH->tabz[0] || Z < WITH->tabz[nz-1]) {
-        printf("SplineInterpDeriv2: out of borders in element s= %4.2f %*s\n",
+        printf("SplineInterpolation2: out of borders in element s= %4.2f %*s\n",
 	       Cell.S, 5, Cell.Elem.PName);
         printf("X = % lf but tabx[0] = % lf and tabx[nx-1] = % lf\n",
 	       is_double<T>::cst(X), WITH->tabx[0], WITH->tabx[nx-1]);
@@ -464,12 +487,14 @@ void SplineInterpolation2(T &X, T &Z, T &thetax, T &thetaz,
     }
 
     out = false;
+    // H-plane
     splin2(WITH->tab2-1, WITH->tab1-1, WITH->tx, WITH->f2x, nz, nx,
 	   Z, X, thetax);
 /*    if (fabs(temp) > ZERO_RADIA)
       *thetax = (double) temp;
     else
       *thetax = 0.0;*/
+    // V-plane
     splin2(WITH->tab2-1, WITH->tab1-1, WITH->tz, WITH->f2z, nz, nx,
 	   Z, X, thetaz);
 /*    if (fabs(temp) > ZERO_RADIA)
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index b778bc6511e670c16389bdcdd82eb38e7b7464fc..9bc3ad5294eef550dc5b83ffa7e7cfcd5d28d990 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -122,6 +122,17 @@ void read_script(const char *param_file_name, bool rd_lat)
       } 
       
       /* read in bool flags */
+      else if (strcmp("globval.quad_fringe", name) == 0){
+        sscanf(line, "%*s %s", str);
+        if(strcmp(str, "true") == 0)
+        	globval.quad_fringe = true;
+        else if(strcmp(str, "false") == 0)
+        	globval.quad_fringe = false;
+        else {
+          printf("set boolean flag true or false for TuneTracFlag \n");
+          exit_(1);
+        }
+      }
       else if (strcmp("TuneTracFlag", name) == 0){
         sscanf(line, "%*s %s", str);
         if(strcmp(str, "true") == 0)
@@ -129,7 +140,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           TuneTracFlag = false;
         else {
-          printf("set bool flag true or false for TuneTracFlag \n");
+          printf("set boolean flag true or false for TuneTracFlag \n");
           exit_(1);
         } 
       }
@@ -140,7 +151,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           ChromTracFlag = false;
         else {
-          printf("set bool flag true or false for ChromTracFlag \n");
+          printf("set boolean flag true or false for ChromTracFlag \n");
           exit_(1);
         } 
       }
@@ -151,7 +162,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
          FmapFlag = false;
         else {
-          printf("set bool flag true or false for FmapFlag \n");
+          printf("set boolean flag true or false for FmapFlag \n");
           exit_(1);
         } 
       }
@@ -162,7 +173,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
            ExperimentFMAFlag= false;
         else {
-          printf("set bool flag true or false for ExperimentFMAFlag \n");
+          printf("set boolean flag true or false for ExperimentFMAFlag \n");
           exit_(1);
         } 
       }
@@ -173,7 +184,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           DetailedFMAFlag = false;
         else {
-          printf("set bool flag true or false for DetailedFMAFlag \n");
+          printf("set boolean flag true or false for DetailedFMAFlag \n");
           exit_(1);
         } 
       }
@@ -184,7 +195,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           TuneShiftFlag = false;
         else {
-          printf("set bool flag true or false for TuneShiftFlag \n");
+          printf("set boolean flag true or false for TuneShiftFlag \n");
           exit_(1);
         } 
       }
@@ -195,7 +206,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           ErrorCouplingFlag = false;
         else {
-          printf("set bool flag true or false for ErrorCouplingFlag \n");
+          printf("set boolean flag true or false for ErrorCouplingFlag \n");
           exit_(1);
         } 
       }
@@ -206,7 +217,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           CouplingFlag = false;
         else {
-          printf("set bool flag true or false for CouplingFlag \n");
+          printf("set boolean flag true or false for CouplingFlag \n");
           exit_(1);
         } 
       }
@@ -217,7 +228,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           MomentumAccFlag = false;
         else {
-          printf("set bool flag true or false for MomentumAccFlag \n");
+          printf("set boolean flag true or false for MomentumAccFlag \n");
           exit_(1);
         } 
       }
@@ -228,7 +239,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           MultipoleFlag = false;
         else {
-          printf("set bool flag true or false for MultipoleFlag \n");
+          printf("set boolean flag true or false for MultipoleFlag \n");
           exit_(1);
         } 
       }
@@ -239,7 +250,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           ThinsextFlag = false;
         else {
-          printf("set bool flag true or false for ThinsextFlag \n");
+          printf("set boolean flag true or false for ThinsextFlag \n");
           exit_(1);
         } 
       }
@@ -250,7 +261,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           FitTuneFlag = false;
         else {
-          printf("set bool flag true or false for FitTuneFlag \n");
+          printf("set boolean flag true or false for FitTuneFlag \n");
           exit_(1);
         } 
       }
@@ -261,7 +272,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           FitChromFlag = false;
         else {
-          printf("set bool flag true or false for FitChromFlag \n");
+          printf("set boolean flag true or false for FitChromFlag \n");
           exit_(1);
         } 
       }
@@ -272,7 +283,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           ChamberFlag = false;
         else {
-          printf("set bool flag true or false for ChamberFlag \n" );
+          printf("set boolean flag true or false for ChamberFlag \n" );
           exit_(1);
         } 
       }
@@ -283,7 +294,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           ChamberNoU20Flag = false;
         else {
-          printf("set bool flag true or false for ChamberNoU20Flag \n");
+          printf("set boolean flag true or false for ChamberNoU20Flag \n");
           exit_(1);
         } 
       }
@@ -294,7 +305,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           ReadChamberFlag = false;
         else {
-          printf("set bool flag true or false for ReadChamberFlag \n");
+          printf("set boolean flag true or false for ReadChamberFlag \n");
           exit_(1);
         } 
       }
@@ -305,7 +316,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           GirderErrorFlag = false;
         else {
-          printf("set bool flag true or false for GirderErrorFlag \n");
+          printf("set boolean flag true or false for GirderErrorFlag \n");
           exit_(1);
         } 
       }
@@ -316,7 +327,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
          SigmaFlag = false;
         else {
-          printf("set bool flag true or false for SigmaFlag \n");
+          printf("set boolean flag true or false for SigmaFlag \n");
           exit_(1);
         } 
       }
@@ -327,7 +338,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           PX2Flag = false;
         else {
-          printf("set bool flag true or false for PX2Flag \n");
+          printf("set boolean flag true or false for PX2Flag \n");
           exit_(1);
         } 
       }
@@ -338,7 +349,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           InducedAmplitudeFlag = false;
         else {
-          printf("set bool flag true or false for InducedAmplitudeFlag \n" );
+          printf("set boolean flag true or false for InducedAmplitudeFlag \n" );
           exit_(1);
         } 
       }
@@ -349,7 +360,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           CodeComparaisonFlag = false;
         else {
-          printf("set bool flag true or false for  CodeComparaisonFlag \n" );
+          printf("set boolean flag true or false for  CodeComparaisonFlag \n" );
           exit_(1);
         } 
       }
@@ -360,7 +371,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           EtaFlag = false;
         else {
-          printf("set bool flag true or false for EtaFlag \n" );
+          printf("set boolean flag true or false for EtaFlag \n" );
           exit_(1);
         } 
       }
@@ -373,7 +384,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           globval.Cavity_on = false;
         else {
-          printf("set bool flag true or false for globval.Cavity_on \n" );
+          printf("set boolean flag true or false for globval.Cavity_on \n" );
           exit_(1);
         } 
       }
@@ -384,7 +395,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           PhaseSpaceFlag = false;
         else {
-          printf("set bool flag true or false for PhaseSpaceFlag \n" );
+          printf("set boolean flag true or false for PhaseSpaceFlag \n" );
           exit_(1);
         } 
       }
@@ -395,7 +406,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           TouschekFlag = false;
         else {
-          printf("set bool flag true or false for TouschekFlag \n" );
+          printf("set boolean flag true or false for TouschekFlag \n" );
           exit_(1);
         } 
       }
@@ -406,7 +417,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           IBSFlag = false;
         else {
-          printf("set bool flag true or false for IBSFlag \n" );
+          printf("set boolean flag true or false for IBSFlag \n" );
           exit_(1);
         } 
       }
@@ -417,7 +428,7 @@ void read_script(const char *param_file_name, bool rd_lat)
         else if(strcmp(str, "false") == 0)
           TousTrackFlag = false;
         else {
-          printf("set bool flag true or false for TousTrackFlag \n" );
+          printf("set boolean flag true or false for TousTrackFlag \n" );
           exit_(1);
         } 
       }
diff --git a/tracy/tracy/src/soleilcommon.cc b/tracy/tracy/src/soleilcommon.cc
index 60b4d4328a8d0724729c27a3041012354a385778..31ccdf104ca8d795477180cb09b4fccdb0ad8939 100644
--- a/tracy/tracy/src/soleilcommon.cc
+++ b/tracy/tracy/src/soleilcommon.cc
@@ -1492,8 +1492,8 @@ long get_qt_number(void)
  
 // }
 
-/****************************************************************************/
-/* 
+/****************************************************************************
+ *
 
    Purpose: 
 
@@ -1516,3 +1516,4 @@ long get_qt_number(void)
        none
 
 ****************************************************************************/
+
diff --git a/tracy/tracy/src/t2elem.cc b/tracy/tracy/src/t2elem.cc
index b83be5af05d33e75b527a468a1b5e0cfd5dd36e4..fe2a04a69df2955c7b739be6a3c1f8e60d4bb6bc 100644
--- a/tracy/tracy/src/t2elem.cc
+++ b/tracy/tracy/src/t2elem.cc
@@ -547,7 +547,7 @@ void radiate(ss_vect<T> &x, const double L, const double h_ref, const T B[])
 
 static double get_psi(double irho, double phi, double gap)
 {
-  /* Correction for magnet gap (lonitudinal fringe field)
+  /* Correction for magnet gap (longitudinal fringe field)
 
        irho h = 1/rho [1/m]
        phi  dipole edge angle
@@ -767,8 +767,8 @@ void bend_fringe(double hb, ss_vect<T> &x)
   }
 }
 
-/****************************************************************************/
-/* template<typename T>
+/****************************************************************************
+ * template<typename T>
   void quad_fringe(double b2, ss_vect<T> &x)
 
    Purpose:
@@ -903,8 +903,9 @@ void Mpole_Pass(CellType &Cell, ss_vect<T> &x)
   if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) 
   {  /* fringe fields */ 
   
-    if (globval.quad_fringe && (M->PB[Quad+HOMmax] != 0.0))
-      quad_fringe(M->PB[Quad+HOMmax], x);
+    if (globval.quad_fringe && (M->PB[Quad+HOMmax] != 0.0) &&  (M->quadFF1 == 1)){
+      quad_fringe(M->quadFFscaling*M->PB[Quad+HOMmax], x);}
+
     
     if (!globval.H_exact) 
        {	 
@@ -1024,8 +1025,8 @@ void Mpole_Pass(CellType &Cell, ss_vect<T> &x)
     {
       bend_fringe(-M->Pirho, x); p_rot(M->PTx2, x);
     }
-    if (globval.quad_fringe && (M->PB[Quad+HOMmax] != 0.0))
-      quad_fringe(-M->PB[Quad+HOMmax], x);
+    if (globval.quad_fringe && (M->PB[Quad+HOMmax] != 0.0) &&  (M->quadFF2 == 1))
+      quad_fringe(-M->quadFFscaling*M->PB[Quad+HOMmax], x);
   } 
 
   /* Local -> Global */
@@ -1046,6 +1047,32 @@ void Marker_Pass(CellType &Cell, ss_vect<T> &X)
 }
 
 
+/****************************************************************************
+ * void Cav_Pass(CellType *Cell, double *X)
+
+   Purpose:
+       Tracking through a cavity
+
+   Input:
+       Cell cavity element to track through
+       X    input coordinates
+
+   Output:
+       X output coordinates
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 template<typename T>
 void Cav_Pass(CellType &Cell, ss_vect<T> &X)
 {
@@ -1550,7 +1577,7 @@ template<typename T>
 void Insertion_Pass(CellType &Cell, ss_vect<T> &x)
 {
   /* Purpose:
-       Track vector x through a insertion
+       Track vector x through an insertion
        If radiation or cavity on insertion is like a drift
 
    Input:
@@ -1602,8 +1629,9 @@ void Insertion_Pass(CellType &Cell, ss_vect<T> &x)
   for (i = 1; i <= Nslice; i++) {
     // First order kick
     if (elemp->ID->firstorder) {
-      if (!elemp->ID->linear)
-        SplineInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable);
+      if (!elemp->ID->linear){
+    	//cout << "InsertionPass: Spline\n";
+        SplineInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable);}
       else
         LinearInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable, 1);
       if (outoftable) {
@@ -1616,10 +1644,12 @@ void Insertion_Pass(CellType &Cell, ss_vect<T> &x)
       
     // Second order kick
     if (elemp->ID->secondorder){
-      if (!elemp->ID->linear)
-        SplineInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable);
-      else
-        LinearInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable, 2);
+      if (!elemp->ID->linear){
+    	//cout << "InsertionPass: Spline\n";
+        SplineInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable);}
+      else{
+      	//cout << "InsertionPass: Linear\n";
+        LinearInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable, 2);}
       if (outoftable) {
 	x[x_] = 1e30;
         return;
@@ -1933,6 +1963,53 @@ static void quadmat(Matrix &ahv, double L, double k)
 }
 
 
+/****************************************************************************/
+/* static void bendmat(vector *M, double L, double irho, double phi1,
+                    double phi2, double gap, double k)
+
+   Purpose:  called  by Mpole_Setmatrix
+
+       For a quadrupole  see quadmat routine for explanation
+
+       For a dipole
+
+                         (1            0 0)
+           Edge(theta) = (h*tan(theta) 1 0)
+                         (0            0 1)
+
+                         (1                 0 0)
+           Edge(theta) = (-h*tan(theta-psi) 1 0)
+                         (0                 0 1)
+
+                                    2
+                   K1*gap*h*(1 + sin phi)
+            psi = -----------------------, K1 = 1/2
+                        cos phi
+
+   Input:
+      L   : length [m]
+      irho: 1/rho [1/m]
+      phi1: entrance edge angle [degres]
+      phi2: exit edge angle [degres]
+      K   : gradient = n/Rho
+
+   Output:
+       M transfer matrix
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       quadmat, make3by3
+       UnitMat, MulRMat, psi
+
+   Comments:
+       none
+
+****************************************************************************/
 static void bendmat(Matrix &M, double L, double irho, double phi1,
                     double phi2, double gap, double k)
 {
@@ -2596,6 +2673,33 @@ static void Wiggler_Print(FILE *f, int Fnum1)
 }
 
 
+/****************************************************************************
+ * void Insertion_Print(FILE **f, long Fnum1)
+
+   Purpose: called by Elem_Print
+       Print out drift characteristics in a file
+           Name, kind, length
+
+   Input:
+       Fnum1 Family number
+       f     pointer on file id
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 static void Insertion_Print(FILE *f, int Fnum1)
 {
   elemtype *elemp;
@@ -2742,6 +2846,31 @@ void MulLsMat(Matrix &A, Matrix &B)
 #undef n
 
 
+/****************************************************************************
+ * void Drift_Alloc(elemtype *Elem)
+
+   Purpose:
+       Dynamic memory allocation for drift element
+
+   Input:
+       Pointer on element
+
+   Output:
+       memory  space for drift in Elem->UU.D
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 void Drift_Alloc(elemtype *Elem)
 {
   Elem->D = (DriftType *)malloc(sizeof(DriftType));
@@ -2777,9 +2906,40 @@ void Mpole_Alloc(elemtype *Elem)
   M->Pgap   = 0.0; /* Gap for fringe field ??? */
 
   M->Pc0 = 0.0; M->Pc1 = 0.0; M->Ps1 = 0.0;
+  M->quadFF1 = 0L;
+  M->quadFF2 = 0L;
+  M->sextFF1 = 0L;
+  M->sextFF2 = 0L;
+  M->quadFFscaling = 0.0;
+
 }
 
 
+/****************************************************************************
+ * void Cav_Alloc(elemtype *Elem)
+
+   Purpose:
+       Constructor for a cavity element
+
+   Input:
+       none
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 void Cav_Alloc(elemtype *Elem)
 {
   CavityType *C;
@@ -2844,6 +3004,34 @@ void FieldMap_Alloc(elemtype *Elem, const bool alloc_fm)
 //  free_matrix(Bz, 1, m_max_FM, 1, n_max_FM);
 }
 
+/****************************************************************************
+ * void Insertion_Alloc(elemtype *Elem, boolean firstflag, boolean secondflag)
+
+   Purpose: called by Insertion_Init and Lat_DealElement
+       Dynamic memory allocation for a Insertion element and various
+       initializations
+
+   Input:
+       Elem Element for memory allocation
+       firstflag true if first order kick map to be loaded
+       secondflag true if second order kick map to be loaded
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       10/01/05 First order kick part added
+
+****************************************************************************/
 
 void Insertion_Alloc(elemtype *Elem)
 {
@@ -2856,7 +3044,7 @@ void Insertion_Alloc(elemtype *Elem)
   ID->Pmethod = Meth_Linear; ID->PN = 0;
   ID->nx = 0; ID->nz = 0;
 
-  /* Initialisation thetax and thetaz to 0*/
+  /* Initialization thetax and thetaz to 0*/
   
   // first order kick map
   if (ID->firstorder){
@@ -2886,13 +3074,10 @@ void Insertion_Alloc(elemtype *Elem)
   // filenames
   strcpy(ID->fname1,""); strcpy(ID->fname2,"");
   
-//  ID->kx = 0e0;
   for (j = 0; j <= 1; j++) {
     ID->PdSsys[j] = 0.0; ID->PdSrnd[j] = 0.0;
   }
   ID->PdTpar = 0.0; ID->PdTsys = 0.0; ID->PdTrnd = 0.0;
-//  for (j = 0; j <= HOMmax; j++)
-//    ID->PBW[j+HOMmax] = 0.0;
   ID->Porder = 0;
 }
 
@@ -3170,6 +3355,32 @@ void FieldMap_Init(int Fnum1)
 }
 
 
+/****************************************************************************
+ * void Cav_Init(long Fnum1)
+
+   Purpose: called by Cell_Init()
+       Constructor for a cavity
+       set the RF voltage, frequency read from the lattice file
+
+   Input:
+       Fnum1 Family number
+
+   Output:
+       none
+
+   Return:
+       none
+
+   Global variables:
+       ElemFam, Cell
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+
+****************************************************************************/
 void Cav_Init(int Fnum1)
 {
   int          i;
diff --git a/tracy/tracy/src/t2lat.cc b/tracy/tracy/src/t2lat.cc
index dd6dc20ca3f72c5102371096409724cd66824d40..f290913da7452cf10e48dcc73586e42a5a49f9b8 100644
--- a/tracy/tracy/src/t2lat.cc
+++ b/tracy/tracy/src/t2lat.cc
@@ -46,10 +46,13 @@ typedef enum
   cctsym, usesym, andsym, dspsym, kicksym, wglsym, nsym, mrksym,
   nbdsym, frgsym, latsym, mpsym, dbnsym, kssym, homsym, lmdsym, dtsym, xytsym,
   vrfsym, harnumsym, frqsym, gstsym, typsym, rollsym, idsym,
-  fnamesym1, fnamesym2, scalingsym, fmsym, harmsym, sprsym, recsym, solsym
+  fnamesym1, fnamesym2, scalingsym, fmsym, harmsym, sprsym, recsym, solsym,
+  ff1sym, ff2sym, ffscalingsym
 } Lat_symbol;   /*\*/
 // idsym fnamesym1 fnamesym2 scalingsym added for insertion
 // ring sym added
+// ff1sym ff2sym for quadrupole entrance and exit fringe field
+// ffscalingsym scaling factor for entrance and exit fringe field.  /*J.Zhang
 /* p2c: t2lat.pas, line 52:
  * Note: Line breaker spent 0.0 seconds, 5000 tries on line 603 [251] */
 
@@ -2009,7 +2012,8 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
   double         t, t1, t2, gap, QL = 0.0, QK;
   double         QKV, QKH, QKxV, QKxH, QPhi, QKS;
   double         dt, Frf, Vrf;
-  long           k1, k2, harnum;
+  long k1, k2, harnum, FF1, FF2;
+  double FFscaling;
   Lat_symbol     sym1;
   symset         mysys, SET;
   long           SET1[(long)lsym / 32 + 2];
@@ -2022,8 +2026,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
   long           SET4[(long)dbnsym / 32 + 2];
   WigglerType    *WITH4;
   FieldMapType   *WITH6;
-  /* ID Laurent */
-  InsertionType  *WITH5;
+  InsertionType  *WITH5;   /* ID Laurent */
   SolenoidType   *WITH7;
   char str1[100] = "";
   char str2[100] = "";
@@ -2232,6 +2235,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
             K =<K-value>, ( [m-2] )
             N =<# of kicks>,
             Method=<method>,
+            FFscaling = Fringe field scaling factor, default 1.0,
             Roll=<roll angle>, ( [deg], design roll angle )
             HOM=(i, <Bi>, <Ai>, ( higher order component in USA notation )
                  j, <Bj>, <Aj>, ( Systematic error Only )
@@ -2241,6 +2245,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
     Example
 
       Q1: Quadrupole, L=0.5, K=2.213455, N=4, Method=4;
+      Q1 : Quadrupole, L=0.5, K=2.213455, N=4, FF =1, FFscaling=2, Method=4;
 
     **************************************************************/
   case qdsym:  /*4*/
@@ -2251,6 +2256,9 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
     k1 = 0;   /* N */
     k2 = Meth_Fourth;   /* method */
     dt = 0.0;
+    FF1 = 0;     /* Entrance Fringe field */
+    FF2 = 0;     /* Exit Fringe field */
+    FFscaling = 1.0; /*Fringe field scaling factor */
     ClearHOMandDBN(&V);
     P_addset(P_expset(mysys, 0), (long)lsym);
     P_addset(mysys, (long)ksym);
@@ -2259,6 +2267,9 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
     P_addset(mysys, (long)rollsym);
     P_addset(mysys, (long)homsym);
     P_addset(mysys, (long)dbnsym);
+    P_addset(mysys, (long)ff1sym);
+    P_addset(mysys, (long)ff2sym);
+    P_addset(mysys, (long)ffscalingsym);
     do {   /*5: read L, K, N, T, T1, T2 */
       test__(mysys, "illegal parameter", &V);
       sym1 = *V.sym;
@@ -2269,6 +2280,18 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
 	QL = EVAL_(&V);
 	break;
 
+      case ffscalingsym:     // read quad ff scaling factor, default 1.0
+      FFscaling = EVAL_(&V);
+      break;
+
+      case ff1sym:
+           FF1 = (long)EVAL_(&V);
+           break;
+
+      case ff2sym:
+           FF2 = (long)EVAL_(&V);
+           break;
+
       case ksym:
 	QK = EVAL_(&V);
 	break;
@@ -2316,6 +2339,12 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       Mpole_Alloc(&WITH->ElemF);
       WITH2 = WITH1->M;
       WITH2->Pmethod = k2; WITH2->PN = k1; WITH2->PdTpar = dt;
+      WITH2->quadFF1 = FF1; /* entrance fringe field flag */
+      WITH2->quadFF2 = FF2; /* exit fringe field flag */
+      WITH2->quadFFscaling = FFscaling; /*fringe field scaling factor flag */
+	  WITH2->sextFF1 = 0; /* entrance sextupole fringe field flag */
+       WITH2->sextFF2 = 0; /* exit sextupole fringe field flag */
+
       AssignHOM(globval.Elem_nFam, &V);
       SetDBN(&V);
       WITH2->n_design = Quad; WITH2->PBpar[HOMmax+2] = QK;
@@ -2349,6 +2378,8 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
     QK = 0.0;   /* K */
     k1 = 0;   /* N */
     k2 = Meth_Fourth;   /* method */
+    FF1 = 0;     /* Entrance Fringe field */
+    FF2 = 0;     /* Exit Fringe field */
     dt = 0.0;
     ClearHOMandDBN(&V);
     getest__(P_expset(SET, (1 << ((long)comma)) | (1 << ((long)semicolon))),
@@ -2362,6 +2393,8 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       P_addset(mysys, (long)rollsym);
       P_addset(mysys, (long)homsym);
       P_addset(mysys, (long)dbnsym);
+      P_addset(mysys, (long)ff1sym);
+      P_addset(mysys, (long)ff2sym);
       do {   /*5: read L, K, N, T, T1, T2 */
 	test__(mysys, "illegal parameter", &V);
 	sym1 = *V.sym;
@@ -2380,6 +2413,14 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
 	    k1 = (long)floor(EVAL_(&V) + 0.5);
 	    break;
 
+	  case ff1sym:
+        FF1 = (long)EVAL_(&V);
+        break;
+
+      case ff2sym:
+        FF2 = (long)EVAL_(&V);
+        break;
+
 	  case mthsym:
 	    k2 = (long)floor(EVAL_(&V) + 0.5);
 	    if (k2 != Meth_Linear) globval.MatMeth = false;
@@ -2426,6 +2467,10 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       WITH2 = WITH1->M;
       WITH2->Pmethod = k2;
       WITH2->PN = k1;
+      WITH2->quadFF1 = 0; /* entrance fringe field flag */
+      WITH2->quadFF2 = 0; /* exit fringe field flag */
+      WITH2->sextFF1 = FF1; /* entrance fringe field flag */
+      WITH2->sextFF2 = FF1; /* exit fringe field flag */
       if (WITH1->PL != 0.0)
 	WITH2->Pthick = pthicktype(thick);
       else
@@ -3562,6 +3607,9 @@ static void init_reserved_words(struct LOC_Lattice_Read *LINK)
   Reg("drift          ", drfsym, &V);
   Reg("dt             ", dtsym, &V);
   Reg("end            ", endsym, &V);
+  Reg("ff1            ", ff1sym, &V);     /* Laurent */
+  Reg("ff2            ", ff2sym, &V);     /* Laurent */
+  Reg("ffscaling      ", ffscalingsym, &V);/* quadff scaling  J.Zhang */
   Reg("fieldmap       ", fmsym, &V);
   Reg("file1          ", fnamesym1, &V);   /* ID Laurent */
   Reg("file2          ", fnamesym2, &V);   /* ID Laurent */