source: TRACY3/branches/tracy3-3.10.1b/tracy/tools/soltracy.cc @ 23

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

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