diff --git a/tracy/tools/Input_checkcode.prm b/tracy/tools/Input_checkcode.prm
index 8196550ee44131de51fa1500826692fcfa18841f..5cb13ce61f3986a93b9cb0caeebd956f5054fa52 100644
--- a/tracy/tools/Input_checkcode.prm
+++ b/tracy/tools/Input_checkcode.prm
@@ -157,7 +157,7 @@ in_dir /Users/nadolski/codes/tracy/maille/soleil/
   EtaFlag  false
 
 # calculate phase space  
-  PhaseSpaceFlag true 0.001 0.0 0.002 0.0 0.0 0.0 100
+#  PhaseSpaceFlag 6D 1e-6 0.0 1e-6 0.0 0.012 0.0 1000 false
  
 #******to be obsoleted************************
 # Do not REMOVE
diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 2da6525d1263c189e9d8d0d87a8a1ed6e2cd3a28..af048af1e8a9c9c6c09d615075054602bf9ef133 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -258,22 +258,22 @@ int main(int argc, char *argv[]) {
   // computes TuneShift with amplitudes
   if (AmplitudeTuneShiftFlag == true) {
     if (ChamberFlag == true) {
-      NuDx(_AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint,
+      TunesShiftWithAmplitude(_AmplitudeTuneShift_nxpoint, _AmplitudeTuneShift_nypoint,
           _AmplitudeTuneShift_nturn, _AmplitudeTuneShift_xmax,
           _AmplitudeTuneShift_ymax, _AmplitudeTuneShift_delta);
       //NuDx(31L,21L,516L,0.025,0.005,dP);
     } else { // Utility ?
-      NuDx(50L, 30L, 516L, 0.035, 0.02, dP);
+      TunesShiftWithAmplitude(50L, 30L, 516L, 0.035, 0.02, dP);
     }
 
   }
   if (EnergyTuneShiftFlag == true) {
     if (ChamberFlag == true) {
-      NuDp(_EnergyTuneShift_npoint, _EnergyTuneShift_nturn,
+      TunesShiftWithEnergy(_EnergyTuneShift_npoint, _EnergyTuneShift_nturn,
           _EnergyTuneShift_deltamax);
       //NuDp(31L,1026L,0.06);
     } else { // utility ?
-      NuDp(31L, 1026L, 0.06);
+      TunesShiftWithEnergy(31L, 1026L, 0.06);
     }
 
   }
@@ -315,7 +315,7 @@ int main(int argc, char *argv[]) {
       exit_(1);
     };
 
-    /* calculate momentum accepetance*/
+    /* calculate momentum acceptance*/
     MomentumAcceptance(_MomentumAccFlag_istart, _MomentumAccFlag_istop,
         _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp,
         _MomentumAccFlag_nstepp, _MomentumAccFlag_deltaminn,
@@ -332,7 +332,7 @@ int main(int argc, char *argv[]) {
   }
 
   if (EtaFlag == true) {
-    // compute cod and twiss parameters for different energy offsets
+    // compute cod and Twiss parameters for different energy offsets
     for (int ii = 0; ii <= 40; ii++) {
       dP = -0.02 + 0.001 * ii;
       Ring_GetTwiss(chroma = false, dP); /* Compute and get Twiss parameters */
@@ -349,10 +349,35 @@ int main(int argc, char *argv[]) {
   }
 
   if (PhaseSpaceFlag == true) {
+    bool cavityflag, radiationflag;
+    /* record the initial values*/
+    cavityflag = globval.Cavity_on;
+    radiationflag = globval.radiation;
+
+    /* set the dimension for the momentum tracking*/
+    if (strncmp("6D", _Phase_Dim, 2) == 0) {
+      globval.Cavity_on = true;
+    } else if (strncmp("4D", _Phase_Dim, 2) == 0) {
+      globval.Cavity_on = false;
+    } else {
+      printf("MomentumAccFlag: Error!!! _Phase_Dim must be '4D' or '6D'\n");
+      exit_(1);
+    };
+    /* setting damping */
+    if (_Phase_Damping == true) {
+      globval.radiation = true;
+    } else  {
+      globval.radiation = false;
+    }
+
     start = stampstart();
     Phase(_Phase_X, _Phase_Px, _Phase_Y, _Phase_Py, _Phase_delta, _Phase_ctau, _Phase_nturn);
-    printf("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);
+
+    /* restore the initial values*/
+    globval.Cavity_on = cavityflag;
+    globval.radiation = radiationflag;
   }
 
   //  ????????????? NSRL version, Check with Soleil version "MomentumAcceptance"
diff --git a/tracy/tracy/inc/naffutils.h b/tracy/tracy/inc/naffutils.h
index f17811c341aeff33d08c4c6383c1599b58b4e71a..b7ade74289374fab46a56c9e85076ccc87ce9538 100644
--- a/tracy/tracy/inc/naffutils.h
+++ b/tracy/tracy/inc/naffutils.h
@@ -23,6 +23,8 @@ void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiz);
 
 void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz);
 
-void Trac_Simple(double x, double px, double y, double py, double dp,
+void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
 		 double ctau, long nmax, double Tx[][NTURN], bool *status);
 
+void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
+         double ctau, long nmax, double Tx[][NTURN], bool *status);
diff --git a/tracy/tracy/inc/read_script.h b/tracy/tracy/inc/read_script.h
index a3755e838f0682e24f02d0ec3db20ba607ae340b..2b5eed89e2ed065976a48f45b48aa319f0d27ae6 100644
--- a/tracy/tracy/inc/read_script.h
+++ b/tracy/tracy/inc/read_script.h
@@ -36,6 +36,8 @@ extern double _MomentumAccFlag_deltaminp, _MomentumAccFlag_deltamaxp;
 extern double _Phase_X, _Phase_Px, _Phase_Y, _Phase_Py,
     _Phase_delta, _Phase_ctau;
 extern long _Phase_nturn;
+extern char _Phase_Dim[3];
+extern bool _Phase_Damping;
 
 extern bool ErrorCouplingFlag ; extern long err_seed; extern double err_rms;
 extern bool CouplingFlag ;
diff --git a/tracy/tracy/inc/soleillib.h b/tracy/tracy/inc/soleillib.h
index f6214690dafb00de4d430a96592b020dee6b330a..e2a0f047fe63ffb07b0ed3d6ef29841f50981733 100644
--- a/tracy/tracy/inc/soleillib.h
+++ b/tracy/tracy/inc/soleillib.h
@@ -44,7 +44,7 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
                double energy, bool diffusion);
                             
 /* Frequency map analysis */
-void NuDp(long Nb, long Nbtour, double emax);
+void TunesShiftWithEnergy(long Nb, long Nbtour, double emax);
 //void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
 //                 double energy, bool diffusion, bool matlab);
 //void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
@@ -54,7 +54,7 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
 void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
               double z, bool diffusion);
 void Nu_Naff(void);
-void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double ymax,
+void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, double ymax,
                  double energy);
 
 /* Vacuum chamber */
diff --git a/tracy/tracy/src/naffutils.cc b/tracy/tracy/src/naffutils.cc
index 3a6fac050f5b5c158d2cb37bd5e99c137377dc15..d0d92979e0cf67cf2390e5619cfdf150d83a93ae 100644
--- a/tracy/tracy/src/naffutils.cc
+++ b/tracy/tracy/src/naffutils.cc
@@ -7,13 +7,17 @@
    J. Bengtsson  NSLS-II, BNL  2004 -        
 
 */
-
+/*
+ Current revision $Revision: 1.3 $
+ On branch $Name: not supported by cvs2svn $
+ Latest change $Date: 2010-12-01 15:29:26 $ by $Author: nadolski $
+*/
 #include "complexeheader.h"
 
 double  pi = M_PI;
 
 /****************************************************************************/
-/* void Trac_Simple(double x, double px, double y, double py, double dp, long nmax,
+/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, long nmax,
                  double Tx[][NTURN], bool *status)
 
    Purpose:
@@ -21,7 +25,7 @@ double  pi = M_PI;
        The 6D phase trajectory is saved in a array
 
    Input:
-       x, px, y, py 4 transverses coordinates
+       x, px, y, py 4 transverse coordinates
        dp           energy offset
        nmax         number of turns
        pos          starting position for tracking
@@ -46,7 +50,7 @@ double  pi = M_PI;
        useful for connection with NAFF
        19/01/03 tracking around the closed orbit
 ****************************************************************************/
-void Trac_Simple(double x, double px, double y, double py, double dp, 
+void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
                  double ctau, long nmax, double Tx[][NTURN], bool *status2)
 {
   bool             lostF = false; /* Lost particle Flag */
@@ -111,6 +115,104 @@ void Trac_Simple(double x, double px, double y, double py, double dp,
   }
 }
 
+/****************************************************************************/
+/* void Trac_Simple6DCOD(double x, double px, double y, double py, double dp, long nmax,
+                 double Tx[][NTURN], bool *status)
+
+   Purpose:
+       Single particle tracking around the 6D closed orbit for NTURN turns
+       The 6D phase trajectory is saved in a array
+
+   Input:
+       x, px, y, py 4 transverses 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
+****************************************************************************/
+void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
+                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
+{
+  bool             lostF = false; /* Lost particle Flag */
+  ss_vect<double>  x1;  /* Tracking coordinates */
+  long             lastpos = globval.Cell_nLoc;
+  long             lastn = 0;
+  Vector2          aperture = {1.0, 1.0};
+
+  *status2 = true; /* stable */
+
+  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 + globval.CODvect[4]; x1[5] = ctau + + globval.CODvect[5];
+
+  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 Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
diff --git a/tracy/tracy/src/physlib.cc b/tracy/tracy/src/physlib.cc
index e089630b0adb9063fd7b823632da723c69807985..14d06f6343c2d0c997cef41b4322ffbd24ed6e75 100644
--- a/tracy/tracy/src/physlib.cc
+++ b/tracy/tracy/src/physlib.cc
@@ -6,7 +6,7 @@
  L. Nadolski   SOLEIL        2002          Link to NAFF, Radia field maps
  J. Bengtsson  NSLS-II, BNL  2004 -
  */
-/* Current revision $Revision: 1.11 $
+/* Current revision $Revision: 1.12 $
  On branch $Name: not supported by cvs2svn $
  Latest change by $Author: nadolski $
 */
@@ -124,7 +124,7 @@ uint32_t stampstop(uint32_t start) {
     struct tm *tm;
     uint32_t stop;
     const bool timedebug = false;
-    bool prt = false;
+    bool prt = true;
     // get the time
     gettimeofday(&tv, &tz);
     tm = localtime(&tv.tv_sec);
@@ -2681,7 +2681,7 @@ void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiy) {
     y  = y0;
     yp = yp0;
 
-    Trac_Simple(x, xp, y, yp, emax, 0.0, Nbtour, Tab, &status);
+    Trac_Simple4DCOD(x, xp, y, yp, emax, 0.0, Nbtour, Tab, &status);
     if (status){
     Get_NAFF(nterm, Nbtour, Tab, fx, fy, nb_freq);
     nux1 = (fabs(fx[0]) > ZERO ? fx[0] : fx[1]);
@@ -2703,7 +2703,7 @@ void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiy) {
     y = y0;
     yp = yp0;
 
-    Trac_Simple(x, xp, y, yp, -emax, 0.0, Nbtour, Tab, &status);
+    Trac_Simple4DCOD(x, xp, y, yp, -emax, 0.0, Nbtour, Tab, &status);
     if (status){
         Get_NAFF(nterm, Nbtour, Tab, fx, fy, nb_freq);
     if (trace)
@@ -2796,7 +2796,7 @@ void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz) {
 
     double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0;
 
-    Trac_Simple(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status);
+    Trac_Simple4DCOD(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status);
     if (status){
        Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
         *nux = (fabs(fx[0]) > ZERO ? fx[0] : fx[1]);
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index b87083b7f3e1a02da489ff0c0a0b608cbcc6f37d..ec444de1479c5c56f3164460cdcab3696029faf1 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -66,6 +66,8 @@ double _MomentumAccFlag_deltaminp=0.01, _MomentumAccFlag_deltamaxp=0.05;
 double _Phase_X=0.0, _Phase_Px=0.0, _Phase_Y=0.0, _Phase_Py=0.0,
     _Phase_delta=0.0, _Phase_ctau=0.0;
 long _Phase_nturn=512L;
+char _Phase_Dim[3]="4D";
+bool _Phase_Damping = false;
 
 bool ReadMultipoleFlag = false;
 bool MultipoleFlag = false, ThinsextFlag = false;
@@ -511,24 +513,21 @@ void read_script(const char *param_file_name, bool rd_lat)
         else {
           printf("set boolean flag true or false for globval.Cavity_on \n" );
           exit_(1);
-        } 
-      }
-      else if (strcmp("PhaseSpaceFlag", name) == 0){
-        sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld", str,
-                        &_Phase_X, &_Phase_Px,
-                        &_Phase_Y, &_Phase_Py,
-                        &_Phase_delta, &_Phase_ctau,
-                        &_Phase_nturn);
-        if(strcmp(str, "true") == 0)
-          PhaseSpaceFlag = true;
-        else if(strcmp(str, "false") == 0)
-          PhaseSpaceFlag = false;
+        }
+      } else if (strcmp("PhaseSpaceFlag", name) == 0) {
+        sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld %s", _Phase_Dim,
+            &_Phase_X, &_Phase_Px, &_Phase_Y, &_Phase_Py, &_Phase_delta,
+            &_Phase_ctau, &_Phase_nturn, str);
+        if (strcmp(str, "true") == 0)
+          _Phase_Damping = true;
+        else if (strcmp(str, "false") == 0)
+          _Phase_Damping = false;
         else {
-          printf("set boolean flag true or false for PhaseSpaceFlag \n" );
+          printf("set boolean flag true or false for PhaseSpaceFlag \n");
           exit_(1);
-        } 
-      }
-       else if (strcmp("TouschekFlag", name) == 0){
+        }
+        PhaseSpaceFlag = true;
+      } else if (strcmp("TouschekFlag", name) == 0){
         sscanf(line, "%*s %s", str);
         if(strcmp(str, "true") == 0)
           TouschekFlag = true;
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index c9bfb1a1f0095c691e78d958f1726b3e9f71a2ea..e783b0cc859b07b51a2f4f0075c91425bef5fcd6 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -765,7 +765,7 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
 
 
 /****************************************************************************/
-/* void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+/* void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
                double energy)
 
    Purpose:
@@ -802,7 +802,7 @@ void Trac_Tab(double x, double px, double y, double py, double dp,
 
 ****************************************************************************/
 #define nterm  4
-void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
+void TunesShiftWithAmplitude(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
                double energy)
 {
   FILE * outf;
@@ -821,7 +821,7 @@ void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
   /* Get time and date */
   newtime = GetTime();
 
-    if (!trace) printf("\n Entering NuDx ... results in nudx.out\n\n");
+    if (!trace) printf("\n Entering TunesShiftWithAmplitude ... results in nudx.out\n\n");
 
     /////////////
     // H tuneshift
@@ -852,7 +852,7 @@ void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
       z  = z0  ;
       zp = zp0 ;
 
-      Trac_Simple(x,xp,z,zp,energy,0.0,Nbtour,Tab,&stable); // tracking around closed orbit
+      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,Nbtour,Tab,&stable); // tracking around closed orbit
       if (stable) {
         Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq); // gets frequency vectors
         Get_freq(fx,fz,&nux,&nuz);  // gets nux and nuz
@@ -893,7 +893,7 @@ void NuDx(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
       z  = z0 + i*zstep;
       zp = zp0;
 
-      Trac_Simple(x,xp,z,zp,energy,0.0,Nbtour,Tab,&stable);
+      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,Nbtour,Tab,&stable);
       if (stable) {
         Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
         Get_freq(fx,fz,&nux,&nuz);  // gets nux and nuz
@@ -971,8 +971,6 @@ double get_D(const double df_x, const double df_y)
        
 ****************************************************************************/
 #define NTERM2  10
-//void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
-//          double energy, bool diffusion, bool matlab)
 void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
           double energy, bool diffusion)	  
 {
@@ -1039,7 +1037,7 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
  //  for (j = -Nbz; j<= Nbz; j++) {
  //    z  = z0 + sgn(j)*sqrt((double)abs(j))*zstep;
      // tracking around closed orbit
-     Trac_Simple(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
+     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);
@@ -1131,8 +1129,6 @@ void fmap(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
 
 ****************************************************************************/
 #define NTERM2  10
-//void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
-//              double z, bool diffusion, bool matlab)
 void fmapdp(long Nbx, long Nbe, long Nbtour, double xmax, double emax,
               double z, bool diffusion)
 {
@@ -1196,7 +1192,7 @@ fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx
      else
       // x  = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
 	 x  = x0 + sqrt((double)j)*xstep;
-     Trac_Simple(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
+     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
@@ -1239,7 +1235,7 @@ fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx
 #undef NTERM2
 
 /****************************************************************************/
-/* void NuDp(long Nb, long Nbtour, double emax)
+/* void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
 
    Purpose:
        Computes tunes versus energy offset by tracking
@@ -1267,7 +1263,7 @@ fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx
 
 ****************************************************************************/
 #define NTERM  4
-void NuDp(long Nb, long Nbtour, double emax)
+void TunesShiftWithEnergy(long Nb, long Nbtour, double emax)
 {
   FILE * outf;
   const char fic[] = "nudp.out";
@@ -1285,7 +1281,7 @@ void NuDp(long Nb, long Nbtour, double emax)
   /* Get time and date */
   newtime = GetTime();
 
-  if (!trace) printf("\n Entering NuDp ...\n\n");
+  if (!trace) printf("\n Entering TunesShiftWithEnergy ...\n\n");
 
   /* Opening file */
   if ((outf = fopen(fic, "w")) == NULL) {
@@ -1311,7 +1307,7 @@ void NuDp(long Nb, long Nbtour, double emax)
     zp   = zp0 ;
     ctau = ctau0;
     
-    Trac_Simple(x,xp,z,zp,dp,ctau,Nbtour,Tab,&status); // tracking around closed orbit
+    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
@@ -1362,10 +1358,10 @@ void NuDp(long Nb, long Nbtour, double emax)
        trace
 
    Specific functions:
-       Trac_Simple, Get_NAFF
+       Trac_Simple6DCOD, Get_NAFF
 
    Comments:
-       none
+       1 December 2010, Call to a Tracking round around the 6D and not 4D closed orbit
 
 ****************************************************************************/
 void Phase(double x,double xp,double y, double yp,double energy, double ctau, long Nbtour)
@@ -1405,7 +1401,7 @@ void Phase(double x,double xp,double y, double yp,double energy, double ctau, lo
     Tab[5][i] = 0.0;
   }
   
-  Trac_Simple(x,xp,y,yp,energy,ctau,Nbtour,Tab,&status);
+  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]);
@@ -1579,7 +1575,7 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
       start = ctau0; break;
   }
 
-  /** Step between intila conditions **/
+  /** Step between initial conditions **/
   step = (end - start)/Nb;
 
   for (j = 0; j <= Nb; j++){
@@ -1600,7 +1596,7 @@ void PhasePortrait(double x0,double px0,double z0, double pz0, double delta0,
 
    fprintf(stdout,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
             x,px,z,pz,delta,ctau);
-    Trac_Simple(x,px,z,pz,delta,ctau,Nbtour,Tab,&status);
+    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",
             Tab[0][i],Tab[1][i],Tab[2][i],Tab[3][i],Tab[4][i],Tab[5][i]);
@@ -3611,7 +3607,7 @@ void spectrum(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    x  = x0 + sqrt((double)i)*xstep;
    for (j = 0; j<= Nbz; j++) {
      z  = z0 + sqrt((double)j)*zstep;
-     Trac_Simple(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
      if (status) {
       Get_NAFF(nterm2, Nbtour, Tab, fx, fz, nb_freq);
      }
@@ -3920,7 +3916,7 @@ void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    x  = x0 + sqrt((double)i)*xstep;
    for (j = 0; j<= Nbz; j++) {
      z  = z0 + sqrt((double)j)*zstep;
-     Trac_Simple(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
 
      if (status) {
        Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq);
@@ -4086,7 +4082,7 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
     x  = x0 + sqrt((double)i)*xstep;
     for (j = 0; j<= Nbz; j++) {
       z  = z0 + sqrt((double)j)*zstep;
-      Trac_Simple(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+      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;
@@ -4108,7 +4104,7 @@ void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
     x  = x0 - sqrt((double)i)*xstep;
     for (j = 0; j<= Nbz; j++) {
       z  = z0 + sqrt((double)j)*zstep;
-      Trac_Simple(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
+      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;