source: TRACY3/trunk/tracy/tools/max4.cc @ 3

Last change on this file since 3 was 3, checked in by zhangj, 12 years ago

Initiale import
  • Property svn:executable set to *
File size: 15.1 KB
Line 
1#define ORDER 1         
2
3int no_tps = ORDER; // arbitrary TPSA order is defined locally
4
5extern bool  freq_map;
6
7#include "tracy_lib.h"
8
9//***************************************************************************************
10//
11//  MAIN CODE
12//
13//****************************************************************************************
14 int main(int argc, char *argv[])
15{
16  const long  seed = 1121;
17  const bool All = TRUE;
18  iniranf(seed); setrancut(2.0);
19
20  // turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
21  // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
22  globval.H_exact     = false; 
23  globval.quad_fringe = false;                // quadrupole fringe field
24 
25  globval.radiation   = false;                // synchrotron radiation
26  globval.emittance   = false;                 // emittance
27  globval.pathlength  = false; 
28
29 
30 
31  // overview, on energy: 25-15
32  //const double  x_max_FMA = 20e-3, y_max_FMA = 10e-3; //const x_max_FMA = 25e-3, y_max_FMA = 15e-3;
33  //const int     n_x = 80, n_y = 80, n_tr = 2048;
34  // overview, off energy: 10-10
35  const double  x_max_FMA = 10e-3, delta_FMA = 10e-2;
36  const int     n_x = 80, n_dp = 80, n_tr = 2048;
37  //
38  // zoom, on energy: 8-2.5
39  //const double  x_max_FMA = 8e-3, y_max_FMA = 2.5e-3;
40  //const int     n_x = 64, n_y = 15, n_tr = 2048;
41  // zoom, off energy: 7-3
42  //const double  x_max_FMA = 3e-3, delta_FMA = 7e-2;
43  //const int     n_x = 28, n_dp = 56, n_tr = 2048;
44 
45  double nux=0.0, nuz=0.0, ksix=0.0, ksiz=0.0;
46 
47  bool chroma;
48  double dP = 0.0;
49  long lastpos = -1L;
50  char str1[S_SIZE];
51 
52  /************************************************************************
53      start read in files and flags
54  *************************************************************************/
55  read_script(argv[1], true);
56
57
58 /************************************************************************
59    end  read in files and flags
60  *************************************************************************/
61   
62 
63 
64 
65 
66 
67//  if (true)
68//  //  Read_Lattice("/home/nadolski/codes/tracy/maille/soleil/solamor2_tracy3");
69//    Read_Lattice(argv[1]); //sets some globval params
70//  else
71//    rdmfile("flat_file.dat"); //instead of reading lattice file, get data from flat file
72
73  //no_sxt(); //turns off sextupoles
74 // Ring_GetTwiss(true, 0e-2); //gettwiss computes one-turn matrix arg=(w or w/o chromat, dp/p)
75  //get_matching_params_scl();
76  //get_alphac2();
77  //GetEmittance(ElemIndex("cav"), true);
78 
79//prt_lat("linlat.out", globval.bpm, true);  /* print lattice file for nsrl-ii*/
80prtmfile("flat_file.dat"); // writes flat file   /* very important file for debug*/
81//prt_chrom_lat(); //writes chromatic functions into chromlat.out
82//  printlatt();  /* print out lattice functions */
83 
84
85/* print lattice file */
86//  prt_lat("linlatBNL.out", globval.bpm, All); // BNL print for all elements
87  printlatt();  /* SOLEIL print out lattice functions */
88  printglob();
89 
90 
91 
92    // Flag factory
93//  bool TuneTracFlag = true;
94//  bool ChromTracFlag = true;
95 
96 
97  //*************************************************************
98  //=============================================================
99
100   
101    // Chamber factory
102  if (ChamberFlag == false)
103     ChamberOff(); // turn off vacuum chamber setting, use the default one
104  else if (ChamberNoU20Flag == true)
105     DefineChNoU20();  // using vacuum chamber setting but without undulator U20
106  else if (ReadChamberFlag == true)
107     ReadCh("Apertures.dat"); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
108//LoadApers("Apertures.dat", 1.0, 1.0);  /* read vacuum chamber definition for bnl */
109  PrintCh();
110 
111 
112  // compute tunes by tracking (should be the same as by DA)
113  if (TuneTracFlag == true) {
114    GetTuneTrac(1026L, 0.0, &nux, &nuz);
115    fprintf(stdout,"From tracking: nux = % f nuz = % f \n",nux,nuz);
116  }
117
118  // compute chromaticities by tracking (should be the same as by DA)
119  if (ChromTracFlag == true){
120    GetChromTrac(2L, 1026L, 1e-5, &ksix, &ksiz);
121    fprintf(stdout,"From tracking: ksix= % f ksiz= % f \n",ksix,ksiz);
122  }
123
124
125  if (FitTuneFlag == true){
126    fprintf(stdout, "\n Fitting tunes\n");
127    FitTune(ElemIndex("qp7"),ElemIndex("qp9"), targetnux, targetnuz);
128    Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
129    printglob();                      /* print parameter list */
130  }
131
132  if (FitChromFlag == true){
133    fprintf(stdout, "\n Fitting chromaticities\n");
134    FitChrom(ElemIndex("sx9"),ElemIndex("sx10"), targetksix, targetksiz);
135    Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
136    printglob();                      /* print parameter list */
137  }
138
139    //SetKLpar(ElemIndex("QT"), 1, 2L,   0.001026770838382);
140 
141  // coupling calculation
142   if (CouplingFlag == true){
143     Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
144     printlatt();                      /* dump linear lattice functions into "linlat.dat" */
145  //   Coupling_Edwards_Teng();
146     printglob();   /* print parameter list */
147   }
148
149  // add coupling by random rotating of the quadrupoles
150  if (ErrorCouplingFlag == true){
151    SetErr();
152    Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
153    printlatt();                      /* dump linear lattice functions into "linlat.dat" */
154//    Coupling_Edwards_Teng();
155    printglob();   /* print parameter list */
156  }
157
158  // WARNING Fit tunes and chromaticities before applying errors !!!!
159  //set multipoles in all magnets
160  if (MultipoleFlag == true ){
161    if (ThinsextFlag ==true){
162      fprintf(stdout, "\n Setting Multipoles for lattice with thin sextupoles \n");
163      Multipole_thinsext();  /* for thin sextupoles */
164     
165      Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
166      printglob(); 
167     }
168    else{
169      fprintf(stdout, "\n Setting Multipoles for lattice with thick sextupoles \n");
170      Multipole_thicksext();  /* for thick sextupoles */
171     
172      Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
173      printglob(); 
174    }
175  }
176                       /* print parameter list */
177 
178
179  // PX2 chicane
180//  if (PX2Flag ==true){
181//  setPX2chicane();
182//  //get closed orbit   
183//  getcod (0.0, &lastpos);
184//  printcod();
185//  Ring_GetTwiss(chroma=true, 0.0);  /* Compute and get Twiss parameters */
186//  printglob();                      /* print parameter list */
187//  }
188 
189 // Computes FMA
190  if (FmapFlag == true){
191    if (ChamberFlag == true ){
192      if (ExperimentFMAFlag == true)
193         fmap(40,12,258,-20e-3,5e-3,0.0,true); // for experimental
194      if (DetailedFMAFlag == true)
195        fmap(100,50,1026,20e-3,5e-3,0.0,true);
196      }
197      else{
198        if (ExperimentFMAFlag == true)
199          fmap(40,12,258,-32e-3,5e-3,0.0,true);
200        if (DetailedFMAFlag == true)
201          fmap(200,100,1026,32e-3,7e-3,0.0,true);
202      }
203  }
204 
205  if (CodeComparaisonFlag){
206          // SOLEIL
207          fmap(100,50,1026,32e-3,7e-3,0.0,true);
208          //fmap(200,100,1026,-32e-3,7e-3,0.0,true);
209  }
210
211  if (MomentumAccFlag == true){
212    //MomentumAcceptance(10L, 28L, 0.01, 0.05, 4L, -0.01, -0.05, 4L);
213     MomentumAcceptance(1L, 28L, 0.01, 0.05, 40L, -0.01, -0.05, 40L);
214  //  MomentumAcceptance(1L, 108L, 0.01, 0.05, 100L, -0.01, -0.05, 100L);
215  }
216
217  // computes Tuneshift with amplitudes
218  if (TuneShiftFlag == true){
219    if (ChamberFlag == true ){
220      TunesShiftWithAmplitude(31L,21L,516L,0.025,0.005,dP);
221      //NuDp(31L,516L,0.06);
222      //NuDp(31L,1026L,0.06);
223      }
224      else{
225        TunesShiftWithAmplitude(50L,30L,516L,0.035,0.02,dP);
226        TunesShiftWithEnergy(31L,1026L,0.06);
227      }
228
229  }
230 
231//  if (SigmaFlag == true){printsigma();
232//  }
233// 
234
235  // induced amplitude
236  if (InducedAmplitudeFlag == true){
237      InducedAmplitude(193L); 
238  }
239 
240 if (EtaFlag == true){
241  // compute cod and twiss parameters for different energy offsets   
242    for (int ii=0; ii<=40; ii++) { 
243    dP = -0.02+ 0.001*ii;
244    Ring_GetTwiss(chroma=false, dP);   /* Compute and get Twiss parameters */
245    printlatt();                      /* dump linear lattice functions into "linlat.dat" */
246    getcod (dP, lastpos);
247//     printcod();
248    prt_cod("cod.out", globval.bpm, true);
249    //system("mv linlat.out linlat_ooo.out");
250    sprintf(str1, "mv cod.out cod_%02d.out", ii); 
251    system(str1);
252    sprintf(str1, "mv linlat.out linlat_%02d.out", ii); 
253    system(str1);
254    }
255}   
256 
257
258 
259  //*********************************************************************************
260  //---------------------------------------------------------------------------------------------------------------------------------
261 
262  // Delicated for max4 lattice. To load alignment error files and do correction
263  if (false) {
264    //prt_cod("cod.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
265   
266   
267    globval.bpm = ElemIndex("bpm_m");                   //definition for max4 lattice, bpm
268  //  globval.bpm = ElemIndex("bpm");                         
269    globval.hcorr = ElemIndex("corr_h"); globval.vcorr = ElemIndex("corr_v");    //definition for max4 lattice, correctors
270   // globval.hcorr = ElemIndex("ch"); globval.vcorr = ElemIndex("cv");
271    globval.gs = ElemIndex("GS"); globval.ge = ElemIndex("GE");   //definition for max4 lattice, girder maker
272   
273   
274    //compute response matrix (needed for OCO)
275    gcmat(globval.bpm, globval.hcorr, 1);  gcmat(globval.bpm, globval.vcorr, 2);
276     
277
278    //print response matrix (routine in lsoc.cc)
279    //prt_gcmat(globval.bpm, globval.hcorr, 1);  prt_gcmat(globval.bpm, globval.vcorr, 2);
280       
281    //gets response matrix, does svd, evaluates correction for N seed orbits
282    //get_cod_rms_scl(dx, dy, dr, n_seed)
283    //get_cod_rms_scl(100e-6, 100e-6, 1.0e-3, 100); //trim coils aren't reset when finished
284       
285
286    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
287    //CorrectCOD_N("/home/simon/projects/src/lattice/AlignErr.dat", 3, 1, 1);
288   
289    //for field errors check LoadFieldErr(in nsls_ii_lib.cc) and FieldErr.dat
290    //LoadFieldErr("/home/simon/projects/src/lattice/FieldErr.dat", true, 1.0, true);
291   
292    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
293    //LoadAlignTol("/home/simon/projects/src/lattice/AlignErr.dat", true, 1.0, true, 1);
294    //LoadAlignTol("/home/simon/projects/out/20091126/AlignErr.dat", true, 1.0, true, 1);
295    //prt_cod("cod_err.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
296   
297     
298    // delicated for max4 lattice
299    //load alignment errors and field errors, correct orbit, repeat N times, and get statistics
300    get_cod_rms_scl_new(1); //trim coils aren't reset when finished
301 
302   
303    //for aperture limitations use LoadApers (in nsls_ii_lib.cc) and Apertures.dat
304    //globval.Aperture_on = true;
305    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
306   
307  }
308
309
310
311//*******************************************************************************
312//------------------------------------------------------------------------------------------------------------------------------------- 
313
314  // Call nsls-ii_lib.cc
315  // tune shift with amplitude
316  double delta = 0;
317  if (false) {
318    cout << endl;
319    cout << "computing tune shifts" << endl;
320    dnu_dA(20e-3, 10e-3, 0.0);
321    get_ksi2(delta); // this gets the chromas and writes them into chrom2.out
322 // get_ksi2(5.0e-2); // this gets the chromas and writes them into chrom2.out
323  }
324 
325  if (false) {
326    //fmap(n_x, n_y, n_tr, x_max_FMA, y_max_FMA, 0.0, true, false);
327//    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true, false); // always use -delta_FMA (+delta_FMA appears broken)
328    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true); // always use -delta_FMA (+delta_FMA appears broken)
329  } else {
330    globval.Cavity_on = true; // this gives longitudinal motion
331    globval.radiation = false; // this adds ripple around long. ellipse (needs many turns to resolve damp.)
332    //globval.Aperture_on = true;
333    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
334    //get_dynap_scl(delta, 512);
335  }
336 
337 
338
339 
340 
341 
342 
343  //
344  // IBS & TOUSCHEK
345  //
346  int     k, n_turns;
347  double  sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
348  FILE    *outf;
349  const double  Qb   = 5e-9;
350 
351  if (false) {
352    double  sum_delta[globval.Cell_nLoc+1][2];
353    double  sum2_delta[globval.Cell_nLoc+1][2];
354   
355    GetEmittance(globval.cav, true);
356   
357    // initialize momentum aperture arrays
358    for(k = 0; k <= globval.Cell_nLoc; k++){
359      sum_delta[k][0] = 0.0; sum_delta[k][1] = 0.0;
360      sum2_delta[k][0] = 0.0; sum2_delta[k][1] = 0.0;
361    }
362   
363    //sigma_delta = 7.70e-04;      //410:7.70e-4,  411:9.57e-4,  412:9.12e-4
364    //globval.eps[X_] = 0.326e-9;  //410:3.26e-10, 411:2.63e-10, 412:2.01e-10
365    globval.eps[Y_] = 0.008e-9;
366    //sigma_s = 9.73e-3;           //410:9.73e-3,  411:12.39e-3, 412:12.50e-3/10.33e-3
367    //globval.eps[Z_] = sigma_delta*sigma_s;
368    //globval.delta_RF given by cav voltage in lattice file
369     //globval.delta_RF = 6.20e-2; //410:6.196e-2, 411:5.285e-2, 412:4.046e-2/5.786e-2
370    n_turns = 490;               //410:490(735), 411:503(755), 412:439(659)/529(794)
371   
372
373    alpha_z =
374      -globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
375      - globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
376    beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
377    gamma_z = (1+sqr(alpha_z))/beta_z;
378   
379    sigma_delta = sqrt(gamma_z*globval.eps[Z_]);
380    sigma_s = sqrt(beta_z*globval.eps[Z_]);//50e-3
381    beta_z = sqr(sigma_s)/globval.eps[Z_];
382    alpha_z = sqrt(beta_z*gamma_z-1);
383
384    // INCLUDE LC (LC changes sigma_s and eps_z, but has no influence on sigma_delta)
385    if (false) {
386      double  newLength, bunchLengthening;
387      newLength = 50e-3;
388      bunchLengthening = newLength/sigma_s;
389      sigma_s = newLength;
390      globval.eps[Z_] = globval.eps[Z_]*bunchLengthening;
391      beta_z = beta_z*bunchLengthening;
392      gamma_z = gamma_z/bunchLengthening;
393      alpha_z = sqrt(beta_z*gamma_z-1);  // this doesn't change
394    }
395   
396    //globval.eps[X_] = 0.362e-9;
397    //sigma_delta = 1.04e-3;
398    //sigma_s = 14.8e-3;
399    Touschek(Qb, globval.delta_RF, globval.eps[X_], globval.eps[Y_],
400             sigma_delta, sigma_s);
401   
402   
403    // IBS
404    if (false) {       
405      // initialize eps_IBS with eps_SR
406      for(k = 0; k < 3; k++)
407        eps[k] = globval.eps[k];
408      for(k = 0; k < 20; k++) //prototype (looping because IBS routine doesn't check convergence)
409        IBS(Qb, globval.eps, eps, alpha_z, beta_z);
410     }
411   
412
413    // TOUSCHEK TRACKING
414    if (false) {       
415      globval.Aperture_on = true;
416      LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
417      tau = Touschek(Qb, globval.delta_RF, false,
418                     globval.eps[X_], globval.eps[Y_],
419                     sigma_delta, sigma_s,
420                     n_turns, true, sum_delta, sum2_delta); //the TRUE flag requires apertures loaded
421     
422      printf("Touschek lifetime = %10.3e hrs\n", tau/3600.0);
423     
424      outf = file_write("mom_aper.out");
425      for(k = 0; k <= globval.Cell_nLoc; k++)
426        fprintf(outf, "%4d %7.2f %5.3f %6.3f\n",
427                k, Cell[k].S, 1e2*sum_delta[k][0], 1e2*sum_delta[k][1]);
428      fclose(outf);
429    }
430 
431  }
432}
Note: See TracBrowser for help on using the repository browser.