source: TRACY3/trunk/tracy/tools/soltracy.cc @ 32

Last change on this file since 32 was 32, checked in by zhangj, 10 years ago

active the transport of the twiss functions and orbits of the transfer line.
  • Property svn:executable set to *
File size: 38.3 KB
Line 
1/************************************
2 Current revision $Revision$
3 On branch $Name$
4 Latest change $Date$ by $Author$
5*************************************/
6#define ORDER 1   
7//#define ORDER 4   
8
9int no_tps = ORDER; // arbitrary TPSA order is defined locally
10
11extern bool freq_map;
12
13// #if HAVE_CONFIG_H
14// #include <config.h>
15// #endif
16
17// #if MPI_EXEC
18//   //library of parallel computation,must be before "stdio.h"
19//   #include <mpi.h>
20// #endif
21
22  #include "tracy_lib.h"
23
24
25//***************************************************************************************
26//
27//  MAIN CODE
28//
29//****************************************************************************************
30int main(int argc, char *argv[]) {
31 
32  /* for time handling */
33  uint32_t start, stop;
34 
35  /*set the random value for random generator*/
36  const long seed = 1121; //the default random seed number
37  iniranf(seed); //initialize the seed
38  setrancut(2.0); //default value of the normal cut for the normal distribution
39  // turn on globval.Cavity_on and globval.radiation to get proper synchro radiation damping
40  // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
41  globval.H_exact = false;
42
43 
44  //output files
45  char fic_twiss[S_SIZE + 4] = "";   //twiss file
46  char fic_summary[S_SIZE + 4]="";   //summary file
47 
48 
49  /* parameters to read the user input script .prm */
50  long i=0L; //initialize the for loop to read command string
51  long j=0L; //initialize the for loop to read the files from parallel computing
52  char CommandStr[max_str];
53  double nux = 0.0, nuy = 0.0, ksix = 0.0, ksiy = 0.0;
54  bool chroma=true;
55  double dP = 0.0;
56  long lastpos = -1L;
57  char str1[S_SIZE];
58 
59  // paramters to read the command line in user input script
60  long CommandNo = -1L;  //the number of commands, since this value is also
61                         // the index of array UserCommandFlag[], so this
62                         // value is always less than 1 in the really case.
63  const int NCOMMAND = 500;
64  UserCommand UserCommandFlag[NCOMMAND];
65 
66// #if MPI_EXEC
67//   //Initialize parallel computation
68//     MPI_Init(&argc, &argv);
69// #endif
70  /************************************************************************
71   read in files and flags
72   *************************************************************************/
73
74   
75  if (argc > 1) {
76    read_script(argv[1], true,CommandNo, UserCommandFlag);
77  } else {
78    fprintf(stdout, "Not enough parameters\n Syntax is program parameter file\n");
79    exit_(1);
80  }
81     
82  /* get the family index of the special elements, prepare for printglob()*/
83  if(strcmp(bpm_name,"")==0)
84    globval.bpm = 0;
85  else
86    globval.bpm = ElemIndex(bpm_name);
87  if(strcmp(skew_quad_name,"")==0)
88    globval.qt = 0; 
89  else
90    globval.qt = ElemIndex(skew_quad_name);
91  if(strcmp(hcorr_name,"")==0)
92    globval.hcorr = 0;
93  else
94    globval.hcorr = ElemIndex(hcorr_name);
95  if(strcmp(vcorr_name,"")==0)
96    globval.vcorr = 0;
97  else
98    globval.vcorr = ElemIndex(vcorr_name);
99  if(strcmp(gs_name,"")==0)
100    globval.gs = 0;
101  else
102    globval.gs = ElemIndex(gs_name);
103  if(strcmp(ge_name,"")==0)
104    globval.ge = 0;
105  else
106    globval.ge = ElemIndex(ge_name);
107 
108 
109//  globval.g = ElemIndex("g");  /* get family index of  girder*/
110 strcpy(fic_twiss,lat_file);
111  strcpy(fic_summary,lat_file);
112  strcat(fic_twiss,".twi"); 
113  strcat(fic_summary,".sum");
114 
115    //generate the twiss and summary files with proper appendix 
116  /* print the summary of the element in lattice to the screen and an external file*/ 
117  printglob();   
118  printglob2file(fic_summary);
119 //print the twiss file
120   /* print out lattice functions, with all the optical information for the lattice with design values */
121  printlatt(fic_twiss); 
122 
123  /************************************************************************
124    print files, very important file for debug
125   *************************************************************************/
126 
127  //print flat file with all the design values of the lattice,
128  prtmfile("flat_file.dat");
129 
130  //find the closed orbit of a ring
131  if(globval.RingType == 1){
132    getcod(dP, lastpos);
133   
134    // compute up to 3rd order momentum compact factor
135    get_alphac2(); 
136  }
137 
138 
139  // print location, twiss parameters and close orbit/orbit at the BPM positions to a file
140    if(globval.bpm !=0){
141      prt_beampos("6Dcod.out");
142    }
143 
144   
145   
146   
147 
148 
149   
150 
151  //get_matching_params_scl(); // get tunes and beta functions at entrance
152 
153  //cout << endl << "computing tune shifts" << endl;
154  //dnu_dA(10e-3, 5e-3, 0.002);
155  //get_ksi2(0.0); // this gets the chromas and writes them into chrom2.out
156 
157   
158 
159   
160 
161/*********************************************************************
162                            Flag factory
163*********************************************************************/
164   
165  for(i=0; i<=CommandNo; i++){//read user defined command by sequence
166   //assign user defined command
167    strcpy(CommandStr,UserCommandFlag[i].CommandStr); 
168   
169       //turn on flag for quadrupole fringe field
170   if(strcmp(CommandStr,"QuadFringeOnFlag") == 0) {
171      globval.quad_fringe = true;
172      cout << "\n";
173      cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
174   }
175   
176         //deactive quadrupole fringe field
177   else if(strcmp(CommandStr,"QuadFringeOffFlag") == 0) {
178      globval.quad_fringe = false;
179      cout << "\n";
180      cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
181   }
182   
183  //set RF voltage 
184  else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
185    printf("\nSetting RF voltage:\n");
186    printf("    Old RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
187        / 1e6);
188    set_RFVoltage(globval.cav, UserCommandFlag[i].RFvolt);
189    printf("    New RF voltage is: %lf [MV]\n", get_RFVoltage(globval.cav)
190        / 1e6);
191  }
192
193  // Chamber factory
194  else if(strcmp(CommandStr,"ReadChamberFlag") == 0) {
195    ReadCh(UserCommandFlag[i].chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
196    PrintCh(); // print chamber into chamber.out
197  }
198 
199  // read the misalignment errors to the elements
200  //  Based on the function error_and_correction() in nsls-ii_lib_templ.h
201  else if (strcmp(CommandStr,"ReadaefileFlag") == 0) {
202   // load misalignments
203    cout << "Read misalignment errors from file: " << UserCommandFlag[i].ae_file << endl;
204    LoadAlignTol(UserCommandFlag[i].ae_file, true, 1.0, true, 0);
205    Ring_GetTwiss(true, 0.0);
206  }
207    //do COD correction using SVD method.
208  else if(strcmp(CommandStr,"CODCorrectFlag") == 0) {
209    cout << "Begin Closed Orbit correction: " << endl;
210   
211    if(n_orbit < 1){
212    cout << "interation number n_obit should NOT small than 1!!!" << endl;
213    exit_(1);
214    }
215   
216    CODCorrect(hcorr_file,vcorr_file, n_orbit,nwh,nwv);
217    Ring_GetTwiss(true, 0.0);
218    prt_cod("summary_miserr_codcorr.out", globval.bpm, true);
219  }
220   // set the field error into the lattice
221  // The corresponding field error is replaced by the new value.
222  // This feature is generic, works for all lattices
223  else if (strcmp(CommandStr,"ReadfefileFlag") == 0) {
224    fprintf(stdout,"\n Read multipole field errors from file:      %s \n", UserCommandFlag[i].fe_file);
225    LoadFieldErrs(UserCommandFlag[i].fe_file, true, 1.0, true, 1);
226       
227    Ring_GetTwiss(true, 0.0);
228    prtmfile("flat_file_fefile.dat");
229  }
230
231  // read multipole errors from a file; specific for soleil lattice
232  else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
233    fprintf(stdout,"\n Read Multipoles field error for SOLEIL lattice with thick sextupoles, from file:\n");
234    fprintf(stdout,"        %s\n",multipole_file);
235    ReadFieldErr(multipole_file);
236    //first print the full lattice with error as a flat file
237    prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors  /* very important file for debug*/
238    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
239    printglob();
240  }
241   // read the sources of coupling from a file; for soleil lattice
242  else if(strcmp(CommandStr,"ReadVirtualSkewquadFlag") == 0) {
243    fprintf(stdout,"\n Read virtual skew quadrupole setting of SOLEIL lattice from file:\n");
244    fprintf(stdout,"    %s \n",virtualskewquad_file);
245   
246    ReadVirtualSkewQuad(virtualskewquad_file);
247    //first print the full lattice with sources of coupling as a flat file
248    prtmfile("flat_file_skewquad.dat");  /* very important file for debug*/
249    //get the coupling
250    Coupling_Edwards_Teng(); 
251  }
252
253     //print the twiss paramaters in a file defined by the name
254   else if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
255      cout << "\n";
256      cout << "print the twiss parameters to file: "
257           << UserCommandFlag[i]._PrintTwiss_twiss_file << "\n";
258           
259      printlatt(UserCommandFlag[i]._PrintTwiss_twiss_file); 
260   }
261   //print the close orbit
262   else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
263      cout << "\n";
264      cout << "print the close orbit to file: "
265           << UserCommandFlag[i]._PrintCOD_cod_file << "\n";
266      getcod(dP, lastpos);
267      prt_cod(UserCommandFlag[i]._PrintCOD_cod_file, globval.bpm, true);
268   }
269       //print coordinates at lattice elements obtained by tracking
270   else if(strcmp(CommandStr,"PrintTrackFlag") == 0) {
271      cout << "\n";
272      cout << "print the tracked coordinates to file: "
273           << UserCommandFlag[i]._PrintTrack_track_file << "\n";
274       
275      PrintTrack(UserCommandFlag[i]._PrintTrack_track_file,
276                 UserCommandFlag[i]._PrintTrack_x, UserCommandFlag[i]._PrintTrack_px, 
277                 UserCommandFlag[i]._PrintTrack_y, UserCommandFlag[i]._PrintTrack_py, 
278                 UserCommandFlag[i]._PrintTrack_delta,UserCommandFlag[i]._PrintTrack_ctau,
279                 UserCommandFlag[i]._PrintTrack_nmax); 
280   }
281   
282   else if(strcmp(CommandStr,"PrintTrackElemFlag") == 0) {
283      cout << "\n";
284      cout << "print the tracked coordinates to file: "
285           << UserCommandFlag[i]._PrintTrackElem_track_file << "\n";
286       
287      PrintTrackElem(UserCommandFlag[i]._PrintTrackElem_track_file,
288                 UserCommandFlag[i]._PrintTrackElem_x, UserCommandFlag[i]._PrintTrackElem_px, 
289                 UserCommandFlag[i]._PrintTrackElem_y, UserCommandFlag[i]._PrintTrackElem_py, 
290                 UserCommandFlag[i]._PrintTrackElem_delta,UserCommandFlag[i]._PrintTrackElem_ctau,
291                 UserCommandFlag[i]._PrintTrackElem_nelem1,UserCommandFlag[i]._PrintTrackElem_nelem2); 
292   }
293    //print the girder
294  // else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
295   //   cout << "\n";
296   //   cout << "print the information of girder to file: "<< girder_file << "\n";
297     
298     
299     // getcod(dP, lastpos);
300     // prt_cod(cod_file, globval.bpm, true);
301  // }
302   
303
304 
305  // compute tunes by tracking (should be the same as by tps)
306  else if (strcmp(CommandStr,"TuneTracFlag") == 0) {
307    GetTuneTrac(1026L, 0.0, &nux, &nuy);
308    fprintf(stdout, "From tracking: nux = % f nuz = % f \n", nux, nuy);
309    }
310
311  // compute chromaticities by tracking (should be the same as by DA)
312  else if(strcmp(CommandStr,"ChromTracFlag") == 0) {
313    start = stampstart();
314    GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiy);
315    stop = stampstop(start);
316    fprintf(stdout, "From tracking: ksix= % f ksiz= % f \n", ksix, ksiy);
317    }
318   
319
320  //generic function, to fit tunes using 1 family of quadrupoles
321 // if (FitTuneFlag == true) {
322  else if(strcmp(CommandStr,"FitTuneFlag") == 0) {
323    double qf_bn = 0.0, qd_bn = 0.0;
324    double qf_bn_fit = 0.0, qd_bn_fit = 0.0;
325
326    /* quadrupole field strength before fitting*/
327    qf_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
328    qd_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
329
330    /* fitting tunes*/
331    fprintf(stdout,
332        "\n Fitting tunes: %s %s, targetnux = %f, targetnuz = %f \n", UserCommandFlag[i].qf, 
333        UserCommandFlag[i].qd,UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
334    FitTune(ElemIndex(UserCommandFlag[i].qf), ElemIndex(UserCommandFlag[i].qd), 
335            UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
336
337    /* integrated field strength after fitting*/
338    qf_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf), 1)].Elem.M->PB[HOMmax + 2];
339    qd_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd), 1)].Elem.M->PB[HOMmax + 2];
340    /* print out the quadrupole strengths before and after the fitting*/
341    printf("Before the fitting, the quadrupole field strengths are: \n");
342    printf(" %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn, UserCommandFlag[i].qd, qd_bn);
343    printf("After the fitting, the quadrupole field strengths are: \n");
344    printf(" %s = %f,    %s = %f\n", UserCommandFlag[i].qf, qf_bn_fit, UserCommandFlag[i].qd, qd_bn_fit);
345
346    /* Compute and get Twiss parameters */
347    Ring_GetTwiss(chroma = true, 0.0);
348    printglob(); /* print parameter list */
349  }
350
351  /* specific for soleil ring in which the quadrupole is cut into 2 parts*/
352  //if (FitTune4Flag == true) {
353  else if(strcmp(CommandStr,"FitTune4Flag") == 0) {
354    double qf1_bn = 0.0, qf2_bn = 0.0, qd1_bn = 0.0, qd2_bn = 0.0;
355    double qf1_bn_fit = 0.0, qf2_bn_fit = 0.0, qd1_bn_fit = 0.0, qd2_bn_fit =
356        0.0;
357
358    /* quadrupole field strength before fitting*/
359    qf1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
360    qf2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
361    qd1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
362    qd2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
363
364    /* fitting tunes*/
365    fprintf(
366        stdout,
367        "\n Fitting tunes for Soleil ring: %s %s %s %s, targetnux = %f, targetnuz = %f \n",
368        UserCommandFlag[i].qf1, UserCommandFlag[i].qf2, UserCommandFlag[i].qd1, UserCommandFlag[i].qd2, 
369        UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
370    FitTune4(ElemIndex(UserCommandFlag[i].qf1), ElemIndex(UserCommandFlag[i].qf2),
371             ElemIndex(UserCommandFlag[i].qd1), ElemIndex(UserCommandFlag[i].qd2),
372             UserCommandFlag[i].targetnux, UserCommandFlag[i].targetnuz);
373
374    /* integrated field strength after fitting*/
375    qf1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf1), 1)].Elem.M->PB[HOMmax + 2];
376    qf2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qf2), 1)].Elem.M->PB[HOMmax + 2];
377    qd1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd1), 1)].Elem.M->PB[HOMmax + 2];
378    qd2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].qd2), 1)].Elem.M->PB[HOMmax + 2];
379    /* print out the quadrupole strengths before and after the fitting*/
380    printf("Before the fitting, the quadrupole field strengths are: \n");
381    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", UserCommandFlag[i].qf1, qf1_bn,
382        UserCommandFlag[i].qf2, qf2_bn, UserCommandFlag[i].qd1, qd1_bn, UserCommandFlag[i].qd2, qd2_bn);
383    printf("After the fitting, the quadrupole field strengths are: \n");
384    printf("    %s = %f,    %s = %f,    %s = %f,    %s = %f\n", UserCommandFlag[i].qf1,
385        qf1_bn_fit, UserCommandFlag[i].qf2, qf2_bn_fit, UserCommandFlag[i].qd1, qd1_bn_fit, 
386        UserCommandFlag[i].qd2, qd2_bn_fit);
387
388    /* Compute and get Twiss parameters */
389    Ring_GetTwiss(chroma = true, 0.0);
390    printglob(); /* print parameter list */
391  }
392
393  /* fit the chromaticities*/
394  else if(strcmp(CommandStr,"FitChromFlag") == 0) {
395 
396    double sxm1_bn = 0.0, sxm2_bn = 0.0;
397    double sxm1_bn_fit = 0.0, sxm2_bn_fit = 0.0;
398    sxm1_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
399    sxm2_bn = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
400
401    fprintf(stdout,"\n Fitting chromaticities: %s %s, targetksix = %f,  targetksiz = %f\n",
402        UserCommandFlag[i].sxm1, UserCommandFlag[i].sxm2, UserCommandFlag[i].targetksix,
403         UserCommandFlag[i].targetksiz);
404         
405    FitChrom(ElemIndex(UserCommandFlag[i].sxm1), ElemIndex(UserCommandFlag[i].sxm2), 
406             UserCommandFlag[i].targetksix, UserCommandFlag[i].targetksiz);
407
408    sxm1_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm1), 1)].Elem.M->PB[HOMmax + 3];
409    sxm2_bn_fit = Cell[Elem_GetPos(ElemIndex(UserCommandFlag[i].sxm2), 1)].Elem.M->PB[HOMmax + 3];
410    /* print out the sextupole strengths before and after the fitting*/
411    printf("Before the fitting, the sextupole field strengths are \n");
412    printf("    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn, UserCommandFlag[i].sxm2, sxm2_bn);
413    printf("After the fitting, the sextupole field strengths are: \n");
414    printf("    %s = %f,    %s = %f\n", UserCommandFlag[i].sxm1, sxm1_bn_fit, UserCommandFlag[i].sxm2, sxm2_bn_fit);
415
416    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
417    printglob(); /* print parameter list */
418  }
419
420  // coupling calculation
421  else if(strcmp(CommandStr,"CouplingFlag") == 0) {
422    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
423    printlatt("linlat_coupling.out"); /* dump linear lattice functions into "linlat.dat" */
424   // GetEmittance(ElemIndex("cav"), true);  //only for test
425    Coupling_Edwards_Teng();
426   // Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
427  //  printglob(); /* print parameter list */
428  }
429
430  // add coupling by random rotating of the full quadrupole magnets
431  //if (ErrorCouplingFlag == true) {
432  else if(strcmp(CommandStr,"ErrorCouplingFlag") == 0) {
433    prtmfile("flat_file.dat"); //print the elements without rotation errors
434    SetErr(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
435    prtmfile("flat_file_errcoupling_full.dat"); //print the elements with rotation errors
436    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
437    printlatt("linlat_errcoupling.out"); /* dump linear lattice functions into "linlat.dat" */
438    Coupling_Edwards_Teng();
439   // GetEmittance(ElemIndex("cav"), true); //only for test
440    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
441  //  printlatt();
442    printglob(); /* print parameter list */
443  }
444 
445//  add coupling by random rotating of the half quadrupole magnets, delicated for soleil
446  else if(strcmp(CommandStr,"ErrorCoupling2Flag") == 0) {
447    prtmfile("flat_file.dat"); //print the elements without rotation errors
448    SetErr2(UserCommandFlag[i].err_seed, UserCommandFlag[i].err_rms);
449    prtmfile("flat_file_errcoupling_half.dat"); //print the elements with rotation errors
450    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
451    printlatt("linlat_errcoupling2.out"); /* dump linear lattice functions into "linlat.dat" */
452    Coupling_Edwards_Teng();
453   // GetEmittance(ElemIndex("cav"), true);
454   Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
455    printglob(); /* print parameter list */
456  }
457
458
459 
460 
461
462  /******************************************************************************************/
463  // COMPUTATION PART after setting the model
464  /******************************************************************************************/
465  // computes TuneShift with amplitudes
466  else if(strcmp(CommandStr,"AmplitudeTuneShiftFlag") == 0) {
467      TunesShiftWithAmplitude(UserCommandFlag[i]._AmplitudeTuneShift_nudx_file,
468                              UserCommandFlag[i]._AmplitudeTuneShift_nudz_file,
469                              UserCommandFlag[i]._AmplitudeTuneShift_nxpoint,
470                              UserCommandFlag[i]._AmplitudeTuneShift_nypoint,
471                              UserCommandFlag[i]._AmplitudeTuneShift_nturn, 
472                              UserCommandFlag[i]._AmplitudeTuneShift_xmax,
473                              UserCommandFlag[i]._AmplitudeTuneShift_ymax, 
474                              UserCommandFlag[i]._AmplitudeTuneShift_delta);
475    } 
476  // compute tuneshift with energy
477 else if(strcmp(CommandStr,"EnergyTuneShiftFlag") == 0) {
478     TunesShiftWithEnergy(UserCommandFlag[i]._EnergyTuneShift_nudp_file,
479                          UserCommandFlag[i]._EnergyTuneShift_npoint, 
480                          UserCommandFlag[i]._EnergyTuneShift_nturn,
481                          UserCommandFlag[i]._EnergyTuneShift_deltamax);
482   }
483
484  // Computes FMA
485  else if(strcmp(CommandStr,"FmapFlag") == 0) {
486    printf("\n begin Fmap calculation for on momentum particles: \n");
487   
488  //  #if MPI_EXEC
489
490 //    //initialization for parallel computing
491 //    int myid=0, numprocs=0;
492 //    int namelen=0;
493 //    char processor_name[MPI_MAX_PROCESSOR_NAME]="";
494 //    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
495 //    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
496 //    MPI_Get_processor_name(processor_name, &namelen);
497 //    printf("\n Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); 
498   
499 //    printf("\n begin Fmap calculation for on momentum particles: \n");
500 //    //Each core or process, characterized by myid, will track different region of fmap
501 //    fmap_p(UserCommandFlag[i]._FmapFlag_fmap_file,
502 //        UserCommandFlag[i]._FmapFlag_nxpoint,
503 //        UserCommandFlag[i]._FmapFlag_nypoint,
504 //        UserCommandFlag[i]._FmapFlag_nturn,
505 //        UserCommandFlag[i]._FmapFlag_xmax,
506 //        UserCommandFlag[i]._FmapFlag_ymax,
507 //        UserCommandFlag[i]._FmapFlag_delta,
508 //        UserCommandFlag[i]._FmapFlag_diffusion,
509 //        UserCommandFlag[i]._FmapFlag_printloss,
510 //        numprocs,myid);
511
512 //    //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
513 //    MPI_Barrier(MPI_COMM_WORLD);
514   
515 //    //Collecting data from all files generated by cores, then write to fmap_p.out
516 //    if(myid==0)
517 //      {
518 //     //  system("cp 0fmap.out fmap_p.out");
519       
520 //     FILE *outf;
521 //     if ((outf = fopen(UserCommandFlag[i]._FmapFlag_fmap_file, "w")) == NULL)
522 //       {
523 //         fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapFlag_fmap_file);
524 //         exit_(1);
525 //       }
526
527 //     FILE *fp;
528 //     char buffer[81];
529 //     char FmapFile[50];
530 //     for(int j=0; j<numprocs; j++)
531 //       {
532 //         FmapFile[0]='\0';
533 //         sprintf(FmapFile,"%d",j);
534 //         strcat(FmapFile,UserCommandFlag[i]._FmapFlag_fmap_file);
535     
536 //         if((fp = fopen(FmapFile, "r")) == NULL) 
537 //           {
538 //             fprintf(stdout, "%s: error while opening file.\n", FmapFile);
539 //             exit_(1);
540 //           }
541
542 //         while(fgets(buffer, 80, fp) != NULL)
543 //           {
544 //             fputs(buffer, outf);
545 //             buffer[0]='\0';
546 //           }
547 //         fclose(fp);
548 //       }
549 //     fclose(outf);
550 //      }
551 //    MPI_Barrier(MPI_COMM_WORLD);
552 //    printf("Process %d finished for on momentum fmap.\n", myid);
553
554 // #else
555      fmap(UserCommandFlag[i]._FmapFlag_fmap_file,
556           UserCommandFlag[i]._FmapFlag_nxpoint, 
557           UserCommandFlag[i]._FmapFlag_nypoint, 
558           UserCommandFlag[i]._FmapFlag_nturn, 
559           UserCommandFlag[i]._FmapFlag_xmax,
560           UserCommandFlag[i]._FmapFlag_ymax, 
561           UserCommandFlag[i]._FmapFlag_delta, 
562           UserCommandFlag[i]._FmapFlag_diffusion,
563           UserCommandFlag[i]._FmapFlag_printloss);
564      //#endif
565  }
566
567  // Compute FMA dp
568  else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
569    printf("\n begin Fmap calculation for off momentum particles: \n");
570   
571 //   #if MPI_EXEC
572
573 //    //initialization for parallel computing
574 //    int myid=0, numprocs=0;
575 //    int namelen=0;
576 //    char processor_name[MPI_MAX_PROCESSOR_NAME]="";
577 //    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
578 //    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
579 //    MPI_Get_processor_name(processor_name, &namelen);
580 //    printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); 
581
582 //    printf("\n begin Fmap calculation for off momentum particles: \n");
583 //    fmapdp_p(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
584 //          UserCommandFlag[i]._FmapdpFlag_nxpoint,
585 //          UserCommandFlag[i]._FmapdpFlag_nepoint,
586 //          UserCommandFlag[i]._FmapdpFlag_nturn,
587 //          UserCommandFlag[i]._FmapdpFlag_xmax,
588 //          UserCommandFlag[i]._FmapdpFlag_emax,
589 //          UserCommandFlag[i]._FmapdpFlag_z,
590 //          UserCommandFlag[i]._FmapdpFlag_diffusion,
591 //            UserCommandFlag[i]._FmapdpFlag_printloss,
592 //          numprocs,myid);
593
594 //    //Synchronize all cores, till all cores finish fmap of different region,then all cores will proceed.
595 //    MPI_Barrier(MPI_COMM_WORLD);
596   
597 //    //Collecting data from all files generated by cores, then write to fmap_p.out
598 //    if(myid==0)
599 //     {
600 //       //  system("cp 0fmap.out fmap_p.out");
601
602 //       FILE *outf;
603 //       if ((outf = fopen(UserCommandFlag[i]._FmapdpFlag_fmapdp_file, "w")) == NULL)
604 //      {
605 //        fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
606 //        exit_(1);
607 //      }
608
609 //       FILE *fp;
610 //       char buffer[81];
611 //       char FmapFile[50];
612 //       for(int j=0; j<numprocs; j++)
613 //      {
614 //        FmapFile[0]='\0';
615 //        sprintf(FmapFile,"%d",j);
616 //        strcat(FmapFile,UserCommandFlag[i]._FmapdpFlag_fmapdp_file);
617     
618 //        if((fp = fopen(FmapFile, "r")) == NULL) 
619 //          {
620 //            fprintf(stdout, "%s: error while opening file.\n", FmapFile);
621 //            exit_(1);
622 //          }
623
624 //        while(fgets(buffer, 80, fp) != NULL)
625 //          {
626 //            fputs(buffer, outf);
627 //            buffer[0]='\0';
628 //          }
629 //        fclose(fp);
630 //      }
631 //       fclose(outf);
632 //     }
633 //    MPI_Barrier(MPI_COMM_WORLD);
634 //    printf("Process %d finished off momentum fmap.\n", myid);
635
636 // #else
637       fmapdp(UserCommandFlag[i]._FmapdpFlag_fmapdp_file,
638              UserCommandFlag[i]._FmapdpFlag_nxpoint, 
639              UserCommandFlag[i]._FmapdpFlag_nepoint, 
640              UserCommandFlag[i]._FmapdpFlag_nturn,
641              UserCommandFlag[i]._FmapdpFlag_xmax, 
642              UserCommandFlag[i]._FmapdpFlag_emax, 
643              UserCommandFlag[i]._FmapdpFlag_z,
644              UserCommandFlag[i]._FmapdpFlag_diffusion,
645               UserCommandFlag[i]._FmapdpFlag_printloss);
646       //  #endif
647  }
648
649//  // if (CodeComparaisonFlag) {
650//   else if(strcmp(CommandStr,"CodeComparaisonFlag") == 0) {
651//     fmap(200, 100, 1026, -32e-3, 7e-3, 0.0, true);
652//   }
653
654  // MOMENTUM ACCEPTANCE
655  else if(strcmp(CommandStr,"MomentumAccFlag") == 0) {
656    bool cavityflag, radiationflag;
657     
658    /* record the initial values*/
659    cavityflag = globval.Cavity_on;
660    radiationflag = globval.radiation; 
661   
662    if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"6D") == 0) {
663      globval.Cavity_on = true;
664      globval.radiation = true;
665    }else if (strcmp(UserCommandFlag[i]._MomentumAccFlag_TrackDim,"4D") == 0) {
666      globval.Cavity_on = false;
667      globval.radiation = false;
668    } else {
669      printf("MomentumAccFlag: Error!!! DimTrack must be '4D' or '6D'\n");
670      exit_(1);
671    };
672
673 /* calculate momentum acceptance*/
674    printf("\n Calculate momentum acceptance: \n");
675
676//  #if MPI_EXEC
677   
678//     /* calculate momentum acceptance*/
679//     //initialization for parallel computing
680//     int myid=0, numprocs=0;
681//     int namelen=0;
682//     char processor_name[MPI_MAX_PROCESSOR_NAME]="";
683//     MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
684//     MPI_Comm_rank(MPI_COMM_WORLD, &myid);
685//     MPI_Get_processor_name(processor_name, &namelen);
686//     printf("Beginning parallel computing: %d of %d is on %s\n", myid, numprocs, processor_name); 
687
688//     printf("\n Calculate momentum acceptance: \n");
689//     MomentumAcceptance_p(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
690//                        UserCommandFlag[i]._MomentumAccFlag_istart,
691//                        UserCommandFlag[i]._MomentumAccFlag_istop,
692//                        UserCommandFlag[i]._MomentumAccFlag_deltaminp,
693//                     UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
694//                        UserCommandFlag[i]._MomentumAccFlag_nstepp,
695//                     UserCommandFlag[i]._MomentumAccFlag_deltaminn,
696//                        UserCommandFlag[i]._MomentumAccFlag_deltamaxn,
697//                     UserCommandFlag[i]._MomentumAccFlag_nstepn,
698//                     UserCommandFlag[i]._MomentumAccFlag_nturn,
699//                     UserCommandFlag[i]._MomentumAccFlag_zmax,
700//                        numprocs,myid);
701//   //merge the files
702//     //Synchronize all cores, till finish cal of different region,
703//     //then files from the cores will be processed.
704//     MPI_Barrier(MPI_COMM_WORLD);
705
706//     //Collect data from all files generated by different cores, then write to user defined file
707//     if(myid==0)
708//       {
709//      //merge the phase.out from different cpus
710//      FILE *outf1;  //phase.out, for debugging
711//      FILE *fp1;
712//      char buffer1[81];
713//      char PhaseFile[50];
714//      if ((outf1 = fopen("phase_p.out", "w")) == NULL)
715//        {
716//          fprintf(stdout, "psoltracy: error while opening file %s\n", "phase_p.out");
717//          exit_(1);
718//        }
719
720//      for(int j=0; j<numprocs; j++)
721//        {
722//          PhaseFile[0]='\0';
723//          sprintf(PhaseFile,"%d",j);
724//          strcat(PhaseFile,"phase.out");
725           
726//          if((fp1 = fopen(PhaseFile, "r")) == NULL){
727//              fprintf(stdout, "%s: error while opening file.\n", PhaseFile);
728//              exit_(1);
729//            }
730
731//          while(fgets(buffer1, 80, fp1) != NULL)
732//            {
733//              fputs(buffer1, outf1);
734//              buffer1[0]='\0';
735//            }
736//          fclose(fp1);
737//        }
738
739//      //collect data for the momentum acceptance
740//      FILE *outf2; //momentum acceptance
741//      FILE *fp2;
742//      char buffer2[81];
743//      char MAFile[50];
744//      if ((outf2 = fopen(UserCommandFlag[i]._MomentumAccFlag_momacc_file, "w")) == NULL)
745//        {
746//          fprintf(stdout, "psoltracy: error while opening file %s\n", UserCommandFlag[i]._MomentumAccFlag_momacc_file);
747//          exit_(1);
748//        }
749
750//      //Collect data in Positive region
751//      for(int j=0; j<numprocs; j++)
752//        {
753//          MAFile[0]='\0';
754//          sprintf(MAFile,"%d",j);
755//          strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
756
757//          if((fp2 = fopen(MAFile, "r")) == NULL) 
758//            {
759//              fprintf(stdout, "%s: error while opening file.\n", MAFile);
760//              exit_(1);
761//            }
762         
763//          while(fgets(buffer2, 80, fp2) != NULL)
764//            {
765//              if(strstr(buffer2,"Negative")!=0) break;
766
767//              fputs(buffer2, outf2);
768//              buffer2[0]='\0'; //reset buffer to NULL
769//            }
770//          buffer2[0]='\0';
771//          fclose(fp2);
772//        }
773
774//      fprintf(outf2,"\n"); // A void line
775
776//      //Collect data in Nagative region
777//      for(int j=0; j<numprocs; j++)
778//        {
779//          MAFile[0]='\0';
780//          sprintf(MAFile,"%d",j);
781//          strcat(MAFile,UserCommandFlag[i]._MomentumAccFlag_momacc_file);
782           
783//          if((fp2 = fopen(MAFile, "r")) == NULL) 
784//            {
785//              fprintf(stdout, "%s: error while opening file.\n", MAFile);
786//              exit_(1);
787//            }
788
789//          bool found=false;
790//          while(fgets(buffer2, 80, fp2) != NULL)
791//            {
792//              if( (strstr(buffer2,"Negative")==0) && (!found)) {  buffer2[0]='\0'; continue;}
793//              if( (strstr(buffer2,"Negative")!=0) ) { found=true; buffer2[0]='\0'; continue;}
794               
795//              fputs(buffer2, outf2);
796//              buffer2[0]='\0';
797//            }
798//          buffer2[0]='\0';
799//          fclose(fp2);
800//        }
801//      fclose(outf2);
802       
803//       }
804//     MPI_Barrier(MPI_COMM_WORLD);
805//     printf("process %d finished momentum acceptance.\n", myid);
806
807// #else
808     MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
809                        UserCommandFlag[i]._MomentumAccFlag_istart, 
810                        UserCommandFlag[i]._MomentumAccFlag_istop,
811                        UserCommandFlag[i]._MomentumAccFlag_deltaminp, 
812                        UserCommandFlag[i]._MomentumAccFlag_deltamaxp,
813                        UserCommandFlag[i]._MomentumAccFlag_nstepp, 
814                        UserCommandFlag[i]._MomentumAccFlag_deltaminn,
815                        UserCommandFlag[i]._MomentumAccFlag_deltamaxn, 
816                        UserCommandFlag[i]._MomentumAccFlag_nstepn,
817                        UserCommandFlag[i]._MomentumAccFlag_nturn,
818                        UserCommandFlag[i]._MomentumAccFlag_zmax);
819     //  #endif   
820
821    /* restore the initial values*/
822    globval.Cavity_on = cavityflag;
823    globval.radiation = radiationflag;
824  }
825
826  // induced amplitude
827  else if(strcmp(CommandStr,"InducedAmplitudeFlag") == 0) {
828   printf("\n Calculate induced amplitude: \n");
829    InducedAmplitude(193L);
830  }
831
832
833  else if(strcmp(CommandStr,"EtaFlag") == 0) {
834    // compute cod and Twiss parameters for different energy offsets
835    for (int ii = 0; ii <= 40; ii++) {
836      dP = -0.02 + 0.001 * ii;
837      Ring_GetTwiss(chroma = false, dP); /* Compute and get Twiss parameters */
838      printlatt("linlat_eta.out"); /* dump linear lattice functions into "linlat.dat" */
839      getcod(dP, lastpos);
840      //     printcod();
841      prt_cod("cod.out", globval.bpm, true);
842      //system("mv linlat.out linlat_ooo.out");
843      sprintf(str1, "mv cod.out cod_%02d.out", ii);
844      system(str1);
845      sprintf(str1, "mv linlat.out linlat_%02d.out", ii);
846      system(str1);
847    }
848  }
849
850// track to get phase space
851  else if(strcmp(CommandStr,"PhaseSpaceFlag") == 0) {
852    bool cavityflag, radiationflag;
853    /* record the initial values*/
854    cavityflag = globval.Cavity_on;
855    radiationflag = globval.radiation;
856
857    /* set the dimension for the momentum tracking*/
858    if (strcmp("6D", UserCommandFlag[i]._Phase_Dim) == 0) {
859      globval.Cavity_on = true;
860    } else if (strcmp("4D", UserCommandFlag[i]._Phase_Dim) == 0) {
861      globval.Cavity_on = false;
862    } else {
863      printf("MomentumAccFlag: Error!!! _Phase_Dim must be '4D' or '6D'\n");
864      exit_(1);
865    };
866    /* setting damping */
867    if ( UserCommandFlag[i]._Phase_Damping == true) {
868      globval.radiation = true;
869    } else  {
870      globval.radiation = false;
871    }
872
873    start = stampstart();
874    Phase(UserCommandFlag[i]._Phase_phase_file,
875          UserCommandFlag[i]._Phase_X, 
876          UserCommandFlag[i]._Phase_Px, 
877          UserCommandFlag[i]._Phase_Y, 
878          UserCommandFlag[i]._Phase_Py, 
879          UserCommandFlag[i]._Phase_delta, 
880          UserCommandFlag[i]._Phase_ctau, 
881          UserCommandFlag[i]._Phase_nturn);
882    printf("6D phase space tracking: \n the simulation time for phase space in tracy 3 is \n");
883    stop = stampstop(start);
884
885    /* restore the initial values*/
886    globval.Cavity_on = cavityflag;
887    globval.radiation = radiationflag;
888  }
889 
890 
891 
892
893 else if (strcmp(CommandStr,"TouschekFlag") == 0 ||strcmp(CommandStr,"IBSFlag") == 0 ||
894          strcmp(CommandStr,"TousTrackFlag") == 0 ){ 
895  //  ????????????? NSRL version, Check with Soleil version "MomentumAcceptance"
896  // IBS & TOUSCHEK
897  int k, n_turns;
898  double sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
899  FILE *outf;
900  const double Qb = 5e-9; // e charge in one bunch
901
902 
903//  else if(strcmp(CommandStr,"TouschekFlag") == 0) {
904    double sum_delta[globval.Cell_nLoc + 1][2];
905    double sum2_delta[globval.Cell_nLoc + 1][2];
906
907    GetEmittance(globval.cav, true);
908
909    // initialize momentum aperture arrays
910    for (k = 0; k <= globval.Cell_nLoc; k++) {
911      sum_delta[k][0] = 0.0;
912      sum_delta[k][1] = 0.0;
913      sum2_delta[k][0] = 0.0;
914      sum2_delta[k][1] = 0.0;
915    }
916
917    globval.eps[Y_] = 0.008e-9;
918    n_turns = 40;
919
920    // get the twiss parameters
921    alpha_z = -globval.Ascr[ct_][ct_] * globval.Ascr[delta_][ct_]
922        - globval.Ascr[ct_][delta_] * globval.Ascr[delta_][delta_];
923    beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
924    gamma_z = (1 + sqr(alpha_z)) / beta_z;
925
926    sigma_delta = sqrt(gamma_z * globval.eps[Z_]);
927    sigma_s = sqrt(beta_z * globval.eps[Z_]);//50e-3
928    beta_z = sqr(sigma_s) / globval.eps[Z_];
929    alpha_z = sqrt(beta_z * gamma_z - 1);
930
931    // INCLUDE LC (LC changes sigma_s and eps_z, but has no influence on sigma_delta)
932    if (false) {
933      double newLength, bunchLengthening;
934      newLength = 50e-3;
935      bunchLengthening = newLength / sigma_s;
936      sigma_s = newLength;
937      globval.eps[Z_] = globval.eps[Z_] * bunchLengthening;
938      beta_z = beta_z * bunchLengthening;
939      gamma_z = gamma_z / bunchLengthening;
940      alpha_z = sqrt(beta_z * gamma_z - 1); // this doesn't change
941    }
942
943    Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
944        sigma_delta, sigma_s);
945
946 // Intra Beam Scattering(IBS)
947    if (strcmp(CommandStr,"IBSFlag") == 0) {
948      // initialize eps_IBS with eps_SR
949      for (k = 0; k < 3; k++)
950        eps[k] = globval.eps[k];
951      for (k = 0; k < 20; k++) //prototype (looping because IBS routine doesn't check convergence)
952        IBS(Qb, globval.eps, eps, alpha_z, beta_z);
953    }
954
955    // TOUSCHEK TRACKING
956    // Calculate Touschek lifetime
957    // with the momentum acceptance which is determined by
958    // the RF acceptance delta_RF and the momentum aperture
959    // at each element location which is tracked over n turns,
960    //the vacuum chamber is read from the file "chamber_file" 
961    // finally, write the momentum acceptance to file "mom_aper.out".
962    if (strcmp(CommandStr,"TousTrackFlag") == 0) {
963        // globval.Aperture_on = true;
964        // ReadCh(UserCommandFlag[i].chamber_file);
965      //      LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
966      tau = Touschek(Qb, globval.delta_RF, false, globval.eps[X_],
967          globval.eps[Y_], sigma_delta, sigma_s, n_turns, true, sum_delta,
968          sum2_delta); //the TRUE flag requires apertures loaded
969
970      printf("Touschek lifetime = %10.3e hrs\n", tau / 3600.0);
971
972      outf = file_write("mom_aper.out");
973      for (k = 0; k <= globval.Cell_nLoc; k++)
974        fprintf(outf, "%4d %7.2f %5.3f %6.3f\n", k, Cell[k].S, 1e2
975            * sum_delta[k][0], 1e2 * sum_delta[k][1]);
976      fclose(outf);
977    }
978
979  }
980 
981  /*
982To do ID correction.
983Based on parts of functions get_param( ) & ID_corr(), etc in nsls-ii_lib.cc
984*/
985  else if(strcmp(CommandStr,"IDCorrFlag") == 0) {
986   
987    int k = 0;
988 
989    cout << endl;
990    cout << "************  Begin ID matching:   **********" <<endl;
991 
992  // read the family index of quadrupoles used for ID correction
993    if (N_calls > 0) {
994      if (N_Fam > N_Fam_max) {
995        printf("get_param: N_Fam > N_Fam_max: %d (%d)\n",
996               N_Fam, N_Fam_max);
997        exit_(0);
998      }
999
1000    for (k = 0; k < N_Fam; k++)
1001      Q_Fam[k] = ElemIndex(IDCq_name[k]);
1002    }
1003
1004    //For debug
1005    if(!trace){
1006      cout << "scl_dbetax = " << scl_dbetax << "    scl_dbetay = " << scl_dbetay <<endl;
1007      cout << "scl_dnux = " << scl_dnux     << "    scl_dnuy = " << scl_dnuy << endl;
1008      cout << "scl_nux = " << scl_nux       << "    scl_nuy = " << scl_nuy << endl;
1009      cout << "ID_step = " << ID_step <<endl;
1010      cout << "N_calls = " << N_calls << "    N_steps = " << N_steps <<endl;
1011      cout << "Number of quadrupole families used for ID correction:   " << N_Fam << endl; 
1012      cout << "Quadrupoles used for ID correction: " << endl;
1013      for (int k = 0; k < N_Fam; k++)
1014        printf("%d\n",Q_Fam[k]);
1015    } 
1016   
1017    //initialization
1018    ini_ID_corr();
1019    printlatt("linlat00.out");
1020
1021    ID_corr0();
1022    printlatt("linlat01.out");
1023
1024  // exit(0);
1025
1026    ini_ID_corr();
1027    printlatt("linlat0.out");
1028
1029    ID_corr(N_calls,N_steps);
1030    printlatt("linlat1.out");
1031
1032  }
1033  else{
1034    printf("Command string %s is invalid!!!!!\n ",CommandStr);
1035    exit_(1);
1036   }//give an error message
1037 }//end of looking for user defined flag
1038
1039 
1040 // #if MPI_EXEC
1041 //  MPI_Finalize();   //finish parallel computation
1042 // #endif
1043 
1044  return 0;
1045}//end of main()
1046
Note: See TracBrowser for help on using the repository browser.