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

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