diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 4ec0e6ab4559089a64072e265a3c66c1a0e31bf5..b6e6053f350e3ddce534fc8a2311bb0ea2341be4 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -282,6 +282,16 @@ int main(int argc, char *argv[]) {
       getcod(dP, lastpos);
       prt_cod(cod_file, globval.bpm, true);
    }
+       //print coordinates at lattice elements obtained by tracking
+   else if(strcmp(CommandStr,"PrintTrackFlag") == 0) {
+      cout << "\n";
+      cout << "print the coordinates to file: "<< track_file << "\n";
+      	
+      PrintTrack(track_file,UserCommandFlag[i]._x,UserCommandFlag[i]._px, 
+                 UserCommandFlag[i]._y,UserCommandFlag[i]._py, 
+                 UserCommandFlag[i]._delta,UserCommandFlag[i]._ctau,UserCommandFlag[i]._nmax);  
+   }
+   
     //print the girder
   // else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
    //   cout << "\n";
diff --git a/tracy/tracy/inc/read_script.h b/tracy/tracy/inc/read_script.h
index d33994ba311c26902ed4aea0d5ba7dc67a9832c8..8a2050f17ce0a594cb33cfe98a7efc5f4e883d21 100644
--- a/tracy/tracy/inc/read_script.h
+++ b/tracy/tracy/inc/read_script.h
@@ -5,6 +5,9 @@
    char  CommandStr[max_str];   // bool flag name
    
  /***** parameters *******/
+   //particle coordinates
+   double _x, _px, _y, _py, _delta, _ctau;
+   long  _nmax; 
     //RFVoltFlag
    double RFvolt;   // RF voltage
    
@@ -122,7 +125,7 @@
  extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
  extern char multipole_file[max_str];
 
-
+ extern char track_file[max_str];
  extern char twiss_file[max_str];
  extern char cod_file[max_str];
  extern char girder_file[max_str];
diff --git a/tracy/tracy/inc/soleillib.h b/tracy/tracy/inc/soleillib.h
index f3a240686c72265ad48e8e16ec1705470cdcc866..d721ef54a0cbc2cd5f9bc28fae3468886ab6a0fc 100644
--- a/tracy/tracy/inc/soleillib.h
+++ b/tracy/tracy/inc/soleillib.h
@@ -112,7 +112,10 @@ void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const
 			 
 /* 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
+void PrintTrack(const char *file_name, double x, double px,double y,double py, 
+                double delta, double ctau, long int nmax);
+		
 
 
 	      
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index ee0dc683e95d79ee1831c0726f423f72b613bb7d..f54966dab3e18a82720134e40e8c0f354b9e6bb4 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -28,6 +28,8 @@
 ****************************************************************************/
 
 /* files */ 
+ //file with coordinates obtained by tracking
+char track_file[max_str];
  //twiss file
 char twiss_file[max_str];
  // cod file
@@ -213,6 +215,19 @@ void read_script(const char *param_file_name, bool rd_lat, long& CommNo, UserCom
           sscanf(line, "%*s %s", cod_file);
 	  strcpy(UserCommandFlag[CommNo].CommandStr,name);  
       }
+      //print the cooridinates using tracking at each element
+      else if (strcmp("PrintTrackFlag", name) == 0){
+          sscanf(line, "%*s %s", track_file);
+	      sscanf(line, "%*s %*s %lf %lf %lf %lf %lf %lf %ld", 
+          		&(UserCommandFlag[CommNo]._x), 
+			&(UserCommandFlag[CommNo]._px),
+          		&(UserCommandFlag[CommNo]._y), 
+			&(UserCommandFlag[CommNo]._py), 
+			&(UserCommandFlag[CommNo]._delta), 
+			&(UserCommandFlag[CommNo]._ctau),
+			&(UserCommandFlag[CommNo]._nmax));
+	  strcpy(UserCommandFlag[CommNo].CommandStr,name);
+      }  
       //print close orbit(COD) flag
       else if (strcmp("PrintGirderFlag", name) == 0){
           sscanf(line, "%*s %s", girder_file);
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index 2cc605fb6cfec151bb07b9c3274bd31090e49de8..917d554b8dbed364243ba131f63c6b660bb520cb 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -5610,5 +5610,119 @@ void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy)
   Ring_Fittune(nu, nueps, nq, qfbuf, qdbuf, nudkL, nuimax);
 }
 
+/**********************************************************************
+void PrintTrack(const char *file_name, double x, double px, double y,double py, 
+           double delta, double ctau, long int nmax)
+	   
+  Purpose:	   
+    Print the coordinates at each lattice element using tracking
+
+  Input:
+      file_name        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
+      nmax             maximum number of turns
+      
+      
+  Output:
+  
+  Comments:
+    Written by Jianfeng Zhang  @ synchro. soleil 05/2011    
+**********************************************************************/	
+void PrintTrack(const char *file_name, double x, double px, double y,double py, 
+           double delta, double ctau, long int nmax)
+{
+   
+    long int i,pos = 1;
+    long int lastn = 0, lastpos = 0;
+    Vector x0, x1, x2, xf,codvector[Cell_nLocMax];
+    FILE *outf;
+    struct    tm *newtime;
+    
+    bool prt = false;
+           
+    outf = file_write(file_name);
+    /* Get time and date */
+    newtime = GetTime();
+
+    fprintf(outf, "# Tracking with TRACY III-- %s -- %s\n",file_name, 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     NTurns \n");
+    fprintf(outf, "#                          [m]        [mm]         [mrad]        [mm]        [mrad]");      
+    fprintf(outf, "        [%%]          [mm]\n");
+
+    
+    //initial coordinates
+    x0[x_] = x;
+    x0[px_] = px;
+    x0[y_] = y;
+    x0[py_] = py;
+    x0[delta_] = delta;
+    x0[ct_] = ctau;	     
+    //get the close orbit at each element around the ring 
+    set_vectorcod(codvector, delta);		     
+    //get the absolute initial coordinates		     
+    x2[x_] = x0[x_] + codvector[1][x_];
+    x2[px_] = x0[px_] + codvector[1][px_];
+    x2[y_] = x0[y_] + codvector[1][y_];
+    x2[py_] = x0[py_] + codvector[1][py_];
+    if (globval.Cavity_on) {
+        x2[delta_] = x0[delta_] + codvector[1][delta_];
+        x2[ct_] = x0[ct_] + codvector[1][ct_];
+    } else {
+        x2[delta_] = x0[delta_];
+        x2[ct_] = x0[ct_];
+    }
+     
+ //print the coordinates at each elements   
+    do {
+        pos = 1;//track from first element
+	(lastn)++;
+        for (i = 0; i < nv_; i++)  //nv_ = 6
+            x1[i] = x2[i];
+
+        while( pos <= globval.Cell_nLoc){  
+         
+          Cell_Pass(pos, pos, x2, lastpos); 	    	
+	  //check whether particle is lost
+	  if (!CheckAmpl(x2, pos)){
+	      fprintf(stderr,"Error!!! %d turn, Particle lost at element: %s!",
+	                      lastn, Cell[pos].Elem.PName);
+	       exit_(1);
+	  }
+	  //get the coordinates around the closed orbit   
+          for (i = x_; i <= py_; i++)   //x_=0,px_=1,y_=2,py_=3
+            xf[i] = x2[i] - codvector[pos][i];
+	    
+          for (i = delta_; i <= ct_; i++) //delta_=4,ct_=5
+            if (globval.Cavity_on && (i != ct_))
+              xf[i] = x2[i] - codvector[pos][i];
+            else
+              xf[i] = x2[i];
+           
+	  if (prt) {
+	    printf("%4ld %4ld %s %6.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %4ld \n",
+                     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,"%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);  	    
+	
+	pos++;	     	     
+	}//finish of tracking and printing to file  
+	
+    } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc));
 
+//     if (globval.MatMeth)
+//         Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
+
+    fclose(outf);   
+}