source: TRACY3/trunk/tracy/tracy/src/max4_lib.cc @ 32

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

Initiale import
  • Property svn:executable set to *
File size: 12.1 KB
Line 
1/*****************************************************
2*        MAX4 specific library                       * 
3*                                                    *   
4*           S.Leemann                                * 
5*                                                    *   
6 ****************************************************/   
7
8//copied here form nsls-ii_lib.cc; needed for LoadFieldErr_scl
9char* get_prm_scl(void)
10{
11  char  *prm;
12
13  prm = strtok(NULL, " \t");
14  if ((prm == NULL) || (strcmp(prm, "\r\n") == 0)) {
15    printf("get_prm: incorrect format\n");
16    exit_(1);
17  }
18
19  return prm;
20}
21
22
23
24
25//copied here form nsls-ii_lib.cc to add incrementing seed value
26void LoadFieldErr_scl(const char *FieldErrorFile, const bool Scale_it,
27                      const double Scale, const bool new_rnd, const int m) 
28{ 
29  bool    rms, set_rnd;
30  char    line[max_str], name[max_str], type[max_str], *prm;
31  int     k, n, seed_val;
32  double  Bn, An, r0;
33  FILE    *inf;
34
35  const bool  prt = true;
36
37  inf = file_read(FieldErrorFile);
38
39  set_rnd = false; 
40  printf("\n");
41  while (fgets(line, max_str, inf) != NULL) {
42    if (strstr(line, "#") == NULL) {
43      // check for whether to set new seed
44      sscanf(line, "%s", name); 
45      if (strcmp("seed", name) == 0) {
46        set_rnd = true;
47        sscanf(line, "%*s %d", &seed_val); 
48        seed_val += 2*m;
49        printf("setting random seed to %d\n", seed_val);
50        iniranf(seed_val); 
51      } else {
52        sscanf(line, "%*s %s %lf", type, &r0);
53        printf("%-4s %3s %7.1le", name, type, r0);
54        rms = (strcmp("rms", type) == 0)? true : false;
55        if (rms && !set_rnd) {
56          printf("LoadFieldErr: seed not defined\n");
57          exit_(1);
58        }
59        // skip first three parameters
60        strtok(line, " \t");
61        for (k = 1; k <= 2; k++)
62          strtok(NULL, " \t");
63        while (((prm = strtok(NULL, " \t")) != NULL) &&
64               (strcmp(prm, "\r\n") != 0)) {
65          sscanf(prm, "%d", &n);
66          prm = get_prm_scl();
67          sscanf(prm, "%lf", &Bn);
68          prm = get_prm_scl(); 
69          sscanf(prm, "%lf", &An);
70          if (Scale_it)
71            {Bn *= Scale; An *= Scale;}
72          if (prt)
73            printf(" %2d %9.1e %9.1e\n", n, Bn, An);
74          // convert to normalized multipole components
75          SetFieldErrors(name, rms, r0, n, Bn, An, true);
76        }
77      }
78    } else
79      printf("%s", line);
80  }
81 
82  fclose(inf);
83}
84
85
86
87
88void get_cod_rms_data(const int n_seed, const int nfam, const int fnums[], const double dx, const double dy, const double dr,
89                      double x_mean[][6], double x_sigma[][6], double theta_mean[][2], double theta_sigma[][2])
90{
91  const int  n_elem = globval.Cell_nLoc+1;
92
93  bool       cod;
94  long int   j;
95  int        i, k, ncorr[2];
96  double     x1[n_elem][6], x2[n_elem][6], theta1[n_elem][6], theta2[n_elem][6], a1L, b1L;;
97
98  for (j = 0; j <= globval.Cell_nLoc; j++){
99    for (k = 0; k < 6; k++) {
100      x1[j][k] = 0.0; x2[j][k] = 0.0;
101    }
102    for (k = 0; k < 2; k++){
103      theta1[j][k] = 0;
104      theta2[j][k] = 0;
105    }
106  } 
107   
108  for (i = 0; i < n_seed; i++) {
109    // reset orbit trims
110    set_bn_design_fam(globval.hcorr, Dip, 0.0, 0.0);
111    set_bn_design_fam(globval.vcorr, Dip, 0.0, 0.0);
112
113    for (j = 0; j < nfam; j++)
114      misalign_rms_fam(fnums[j], dx, dy, dr, true);
115   
116    cod = orb_corr(3);
117   
118    if (cod) {
119     
120      for (k = 0; k < 2; k++)
121        ncorr[k] = 0;
122     
123      for (j = 0; j <= globval.Cell_nLoc; j++){ // read back beam pos at bpm
124
125        if ( Cell[j].Elem.Pkind == Mpole ) {
126          for (k = 0; k < 6; k++) {
127            x1[j][k] += Cell[j].BeamPos[k];
128            x2[j][k] += sqr(Cell[j].BeamPos[k]);
129          }
130         
131          if ( (Cell[j].Fnum == globval.hcorr) || (Cell[j].Fnum == globval.vcorr) ){ // read back corr strength at corr
132            get_bnL_design_elem(Cell[j].Fnum, Cell[j].Knum, Dip, b1L, a1L);
133            if ( Cell[j].Fnum == globval.hcorr ){
134              ncorr[X_]++;
135              theta1[ncorr[X_]-1][X_] += b1L;
136              theta2[ncorr[X_]-1][X_] += sqr(b1L);
137            } else if ( Cell[j].Fnum == globval.vcorr ){
138              ncorr[Y_]++;
139              theta1[ncorr[Y_]-1][Y_] += a1L;
140              theta2[ncorr[Y_]-1][Y_] += sqr(a1L);
141            }
142          }
143        }
144      }
145     
146    } else
147      printf("orb_corr: failed\n");   
148  }
149 
150  for (k = 0; k < 2; k++)
151    ncorr[k] = 0;
152 
153  for (j = 0; j <= globval.Cell_nLoc; j++){
154       
155    if ( Cell[j].Elem.Pkind == Mpole ) {
156      for (k = 0; k < 6; k++) {
157        x_mean[j][k] = x1[j][k]/n_seed;
158        x_sigma[j][k] = sqrt((n_seed*x2[j][k]-sqr(x1[j][k]))
159                             /(n_seed*(n_seed-1.0)));
160      }
161     
162      if ( Cell[j].Fnum == globval.hcorr ){
163        ncorr[X_]++;
164        theta_mean[ncorr[X_]-1][X_] = theta1[ncorr[X_]-1][X_]/n_seed;
165        theta_sigma[ncorr[X_]-1][X_] = sqrt((n_seed*theta2[ncorr[X_]-1][X_]-sqr(theta1[ncorr[X_]-1][X_]))
166                                            /(n_seed*(n_seed-1.0)));
167       
168      } else if ( Cell[j].Fnum == globval.vcorr ){
169        ncorr[Y_]++;
170        theta_mean[ncorr[Y_]-1][Y_] = theta1[ncorr[Y_]-1][Y_]/n_seed;
171        theta_sigma[ncorr[Y_]-1][Y_] = sqrt((n_seed*theta2[ncorr[Y_]-1][Y_]-sqr(theta1[ncorr[Y_]-1][Y_]))
172                                            /(n_seed*(n_seed-1.0)));
173      }
174    }
175  }
176}
177
178
179
180
181void prt_cod_rms_data(const char name[], double x_mean[][6], double x_sigma[][6], double theta_mean[][2], double theta_sigma[][2])
182{
183  long     j;
184  int      k, ncorr[2];
185  FILE     *fp;
186  struct    tm *newtime;
187
188  fp = file_write(name);
189
190
191  /* Get time and date */
192  newtime = GetTime();
193
194  fprintf(fp,"# TRACY III v.3.5 -- %s -- %s \n",
195          name, asctime2(newtime));
196  fprintf(fp,"#         s   name              code  xcod_mean +/-  xcod_rms   ycod_mean +/-  ycod_rms"
197             "   dipx_mean +/-  dipx_sigma dipy_mean +/-  dipy_sigma\n");
198  fprintf(fp,"#        [m]                             [mm]           [mm]       [mm]           [mm]"
199             "      [mrad]         [mrad]     [mrad]         [mrad]\n");
200  fprintf(fp, "#\n");
201
202
203  for (k = 0; k < 2; k++)
204    ncorr[k] = 0;
205 
206  for (j = 0; j <= globval.Cell_nLoc; j++){
207   
208    fprintf(fp, "%4li %8.3f %s %6.2f %10.3e +/- %10.3e %10.3e +/- %10.3e",
209            j, Cell[j].S, Cell[j].Elem.PName, get_code(Cell[j]),
210            1e3*x_mean[j][x_], 1e3*x_sigma[j][x_],
211            1e3*x_mean[j][y_], 1e3*x_sigma[j][y_]);
212   
213    if ( Cell[j].Fnum == globval.hcorr ){
214      ncorr[X_]++;
215      fprintf(fp, " %10.3e +/- %10.3e %10.3e +/- %10.3e\n",
216              1e3*theta_mean[ncorr[X_]-1][X_], 1e3*theta_sigma[ncorr[X_]-1][X_], 0.0, 0.0);
217
218    } else if ( Cell[j].Fnum == globval.vcorr ){
219      ncorr[Y_]++;
220      fprintf(fp, " %10.3e +/- %10.3e %10.3e +/- %10.3e\n",
221              0.0, 0.0, 1e3*theta_mean[ncorr[Y_]-1][Y_], 1e3*theta_sigma[ncorr[Y_]-1][Y_]);
222    } else 
223      fprintf(fp, " %10.3e +/- %10.3e %10.3e +/- %10.3e\n",
224              0.0, 0.0, 0.0, 0.0);
225  }
226
227  fclose(fp);
228}
229
230
231
232
233void add_family( const char *name, int &nfam, int fnums[] )
234{
235  int fnum;
236  fnum = ElemIndex(name);
237  if (fnum != 0)
238    fnums[nfam++] = fnum;
239  else {
240    printf("Family not defined. %s %d\n", name, fnum); 
241    exit(1);
242  }
243}
244
245
246
247
248//copied here form nsls-ii_lib.cc to make changes
249void get_cod_rms_scl(const double dx, const double dy, const double dr,
250                     const int n_seed)
251{
252  const int  n_elem = globval.Cell_nLoc+1;
253
254  int      fnums[25], nfam; 
255  double   x_mean[n_elem][6], x_sigma[n_elem][6], theta_mean[n_elem][2], theta_sigma[n_elem][2];
256
257  nfam = 0;
258  add_family("qf", nfam, fnums);
259  add_family("qfm", nfam, fnums);
260  add_family("qfend", nfam, fnums);
261  add_family("qdend", nfam, fnums);
262  add_family("sfi", nfam, fnums);
263  add_family("sfo", nfam, fnums);
264  add_family("sfm", nfam, fnums);
265  add_family("sd", nfam, fnums);
266  add_family("sdend", nfam, fnums);
267
268  get_cod_rms_data(n_seed, nfam, fnums, dx, dy, dr, x_mean, x_sigma, theta_mean, theta_sigma);
269  prt_cod_rms_data("cod_rms.out", x_mean, x_sigma, theta_mean, theta_sigma);
270
271}
272
273
274
275
276/////////////////////////////////////////////////////////////////////////////////////////////////
277
278
279
280
281//this is a mashup of get_cod_rms_scl
282//instead of assigning errors to families we load and apply AlignErr.dat and FieldErr.dat
283//a correction is then applied and the results saved
284//this is repeated for n_seed seeds and overall statistics saved in cod_rms.out
285
286void get_cod_rms_scl_new(const int n_seed)
287{
288     
289  const int  n_elem = globval.Cell_nLoc+1;
290 
291  double     x_mean[n_elem][6], x_sigma[n_elem][6], theta_mean[n_elem][2], theta_sigma[n_elem][2];
292  bool       cod;
293  long int   j;
294  int        i, k, ncorr[2];
295  double     x1[n_elem][6], x2[n_elem][6], theta1[n_elem][6], theta2[n_elem][6], a1L, b1L;;
296 
297   
298 
299     
300  for (j = 0; j <= globval.Cell_nLoc; j++){
301    for (k = 0; k < 6; k++) {
302      x1[j][k] = 0.0; x2[j][k] = 0.0;
303    }
304    for (k = 0; k < 2; k++){
305      theta1[j][k] = 0;
306      theta2[j][k] = 0;
307    }
308//    printf("1  is one \n");
309  } 
310 
311 
312     
313  for (i = 0; i < n_seed; i++) {
314   
315    // reset orbit trims
316    set_bn_design_fam(globval.hcorr, Dip, 0.0, 0.0);
317    set_bn_design_fam(globval.vcorr, Dip, 0.0, 0.0);
318   
319    LoadFieldErr_scl("/home/zhang/codes/TracyIII/lattice/FieldErr.dat", false, 1.0, true, i);
320    LoadAlignTol("/home/zhang/codes/TracyIII/lattice/AlignErr.dat", false, 1.0, true, i);
321
322    cout << endl << "*** COD ITERATION " << i+1 << " ***" << endl;
323   
324    cod = orb_corr(3);
325   
326     
327    if (cod) {
328     
329      for (k = 0; k < 2; k++)
330        ncorr[k] = 0;
331     
332      for (j = 0; j <= globval.Cell_nLoc; j++){ // read back beam pos at bpm
333       
334        if ( Cell[j].Elem.Pkind == Mpole ) {
335          for (k = 0; k < 6; k++) {
336            x1[j][k] += Cell[j].BeamPos[k];
337            x2[j][k] += sqr(Cell[j].BeamPos[k]);
338          }
339         
340          if ( (Cell[j].Fnum == globval.hcorr) || (Cell[j].Fnum == globval.vcorr) ){ // read back corr strength at corr
341            get_bnL_design_elem(Cell[j].Fnum, Cell[j].Knum, Dip, b1L, a1L);
342            if ( Cell[j].Fnum == globval.hcorr ){
343              ncorr[X_]++;
344              theta1[ncorr[X_]-1][X_] += b1L;
345              theta2[ncorr[X_]-1][X_] += sqr(b1L);
346            } else if ( Cell[j].Fnum == globval.vcorr ){
347              ncorr[Y_]++;
348              theta1[ncorr[Y_]-1][Y_] += a1L;
349              theta2[ncorr[Y_]-1][Y_] += sqr(a1L);
350            }
351          }
352        }
353      }
354     
355    } else
356      printf("orb_corr: failed\n");   
357  }
358 
359  for (k = 0; k < 2; k++)
360    ncorr[k] = 0;
361 
362  for (j = 0; j <= globval.Cell_nLoc; j++){
363   
364    if ( Cell[j].Elem.Pkind == Mpole ) {
365      for (k = 0; k < 6; k++) {
366        x_mean[j][k] = x1[j][k]/n_seed;
367        x_sigma[j][k] = sqrt((n_seed*x2[j][k]-sqr(x1[j][k]))
368                             /(n_seed*(n_seed-1.0)));
369      }
370     
371      if ( Cell[j].Fnum == globval.hcorr ){
372        ncorr[X_]++;
373        theta_mean[ncorr[X_]-1][X_] = theta1[ncorr[X_]-1][X_]/n_seed;
374        theta_sigma[ncorr[X_]-1][X_] = sqrt((n_seed*theta2[ncorr[X_]-1][X_]-sqr(theta1[ncorr[X_]-1][X_]))
375                                            /(n_seed*(n_seed-1.0)));
376       
377      } else if ( Cell[j].Fnum == globval.vcorr ){
378        ncorr[Y_]++;
379        theta_mean[ncorr[Y_]-1][Y_] = theta1[ncorr[Y_]-1][Y_]/n_seed;
380        theta_sigma[ncorr[Y_]-1][Y_] = sqrt((n_seed*theta2[ncorr[Y_]-1][Y_]-sqr(theta1[ncorr[Y_]-1][Y_]))
381                                            /(n_seed*(n_seed-1.0)));
382      }
383    }
384  }
385 
386  prt_cod_rms_data("cod_rms.out", x_mean, x_sigma, theta_mean, theta_sigma);
387  // this writes the statictics over N seeds to file
388 
389}
390
391
392
393
394/////////////////////////////////////////////////////////////////////////////////////////////////
395
396
397
398// get dynamic aperture and output to file: dynap.out, dynap_dp.out
399double get_dynap_scl(const double delta, const int n_track2)
400{
401  char      str[max_str];
402  int       i;
403  double    x_aper[n_aper], y_aper[n_aper], DA;
404  FILE      *fp;
405
406  const int  prt = true;
407
408  fp = file_write("dynap.out");
409  dynap(fp, 5e-3, 0.0, 0.1e-3, n_aper, n_track2, x_aper, y_aper, false, prt);
410  fclose(fp);
411  DA = get_aper(n_aper, x_aper, y_aper);
412
413  if (true) {
414    sprintf(str, "dynap_dp%3.1f.out", 1e2*delta);
415    fp = file_write(str);
416    dynap(fp, 5e-3, delta, 0.1e-3, n_aper, n_track2,
417          x_aper, y_aper, false, prt);
418    fclose(fp);
419    DA += get_aper(n_aper, x_aper, y_aper);
420
421    for (i = 0; i < nv_; i++)
422      globval.CODvect[i] = 0.0;
423    sprintf(str, "dynap_dp%3.1f.out", -1e2*delta);
424    fp = file_write(str);
425    dynap(fp, 5e-3, -delta, 0.1e-3, n_aper,
426          n_track2, x_aper, y_aper, false, prt);
427    fclose(fp);
428    DA += get_aper(n_aper, x_aper, y_aper);
429  }
430
431  return DA/3.0;
432}
433
434
435
436
437/////////////////////////////////////////////////////////////////////////////////////////////////
438
439
440
441// get tunes and beta functions
442void get_matching_params_scl()
443{
444  double nux, nuy, betax, betay;
445 
446  nux = globval.TotalTune[X_];
447  nuy = globval.TotalTune[Y_];
448  betax = globval.OneTurnMat[0][1]/sin(2*pi*nux);
449  betay = globval.OneTurnMat[2][3]/sin(2*pi*nuy);
450 
451  printf("\n");
452  printf("beta_x* = %10.9e\n", betax);
453  printf("beta_y* = %10.9e\n", betay);
454  printf("\n");
455}
456
457
458
Note: See TracBrowser for help on using the repository browser.