source: TRACY3/trunk/tracy/tracy/src/nsls-ii_lib.cc @ 3

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

Initiale import
  • Property svn:executable set to *
File size: 148.5 KB
Line 
1/* NSLS-II specific library
2
3   J. Bengtsson  NSLS-II, BNL  2004 -       
4
5   T. Shaftan, I. Pinayev, Y. Luo, C. Montag, B. Nash
6
7*/
8/* Current revision $Revision: 1.18 $
9 On branch $Name:  $
10 Latest change by $Author: zhang $
11*/
12
13// global params
14
15bool    DA_bare       = false,  // dynamic aperture
16        freq_map      = false, // include frequency map
17        tune_shift    = false; // to calculate tune shift with amplitude and energy
18int     n_aper        =  15, // no of dynamical aperture points
19        n_track       = 512, // no of turns for tracking
20        n_orbit       = 3,  //number of integration for orbit correction
21        n_scale       = 1; //number to scale the random orbit error
22
23bool    bba           = false;   // beam based alignment
24int     n_lin         =  3,
25        SQ_per_scell  =  2,
26        BPM_per_scell = 12,
27        HCM_per_scell = 12,
28        VCM_per_scell = 12;
29double  kick          = 0.01e-3; // 0.01 mrad kick for trims
30int     n_stat        = 1;       // number of statistics;when set the random error;
31
32double  VDweight      = 1e3,     // weight for vertical dispersion
33        HVweight      = 1e0,     // weight for coupling Htrim vertical BPM
34        VHweight      = 1e0;     // weight for coupling Vtrim horizontal BPM
35
36
37double  delta_DA_  = 5.0e-2; // Parameters for dynamic aperture
38
39// Parameters for frequency map
40// Note NTURN is set to 10000 (2*NTURN for diffusion)) in "naffutils.h".
41int     n_x = 50, n_y = 30, n_dp = 25, n_tr = 2064;
42double  x_max_FMA = 20e-3, y_max_FMA = 6e-3, delta_FMA = 3e-2;
43//double  x_max_FMA = 20e-3, y_max_FMA = 3e-3, delta_FMA = 3e-2;
44
45const int  max_corr = 100;
46const int  max_bpm  = 200;
47
48int     h_corr[max_corr], v_corr[max_corr], bpm_loc[max_bpm]; //non-global varible
49
50int     Fnum_Cart, n_iter_Cart;
51
52double  u_Touschek;  // argument for Touschek D(ksi)
53double  chi_m;       // argument for IBS D(ksi)
54
55char    in_dir[max_str]          = "";
56
57// Computation result files
58const char  beam_envelope_file[] = "beam_envelope.out";
59
60// Lattice error and correction files
61const char  CodCorLatFileName[]  = "codcorlat.out";
62
63const char  SkewMatFileName[]    = "skewmat.out";
64const char  eta_y_FileName[]     = "eta_y.out";
65const char  deta_y_FileName[]    = "deta_y.out";
66 
67
68const int  max_elem = Cell_nLocMax;
69
70
71
72char     ae_file[max_str], fe_file[max_str], ap_file[max_str];
73int      N_BPM, N_HCOR, N_VCOR, N_SKEW, N_COUPLE;
74
75/* ID correction**********/
76// maximum number of quadrupole families used for ID correction
77const int  N_Fam_max = 15;
78// ID correction params
79int N_calls = 1, N_steps = 1, N_Fam, Q_Fam[N_Fam_max]; 
80// Weights for ID correction
81double  scl_nux = 1e2, scl_nuy = 1e2, scl_dbetax = 1e0, scl_dbetay = 1e0;
82double  scl_dnux = 0.1, scl_dnuy = 0.1, ID_step = 0.7;
83
84int      n_sext, sexts[max_elem];
85double   beta0_[max_elem][2], nu0s[max_elem][2], nu0_[2], b2[N_Fam_max];
86double   **SkewRespMat, *VertCouple, *SkewStrengthCorr, *eta_y;
87double   *b, *w, **V, **U;
88double   disp_wave_y;
89Vector2  beta_ref;
90
91char lat_FileName[max_str];
92
93// ID_corr global variables
94
95const int  n_b2_max    = 500;  // max no of quad corrector families
96const int  n_b3_max    = 1000; // max no of sextupoles
97const int  max_ID_Fams = 25;   // max no of ID families
98
99int           Nsext, Nquad, Nconstr, NconstrO, quad_prms[n_b2_max], id_loc;
100int           n_ID_Fams, ID_Fams[max_ID_Fams];
101double        Ss[n_b3_max], Sq[n_b2_max], sb[2][n_b3_max], sNu[2][n_b3_max];
102double        qb[2][n_b2_max], qb0[2][n_b2_max], sNu0[2][n_b3_max];
103double        qNu0[2][n_b2_max], qNu[2][n_b2_max], IDb[2], IDNu[2];
104double        Nu_X, Nu_Y, Nu_X0, Nu_Y0;
105double        **A1, *Xsext, *Xsext0, *b2Ls_, *w1, **U1, **V1;
106double        *Xoct, *b4s, **Aoct;
107Vector2       dnu0, nu_0;
108
109ss_vect<tps>  map;
110MNF_struct    MNF;
111 
112
113// conversion
114
115void lwr_case(char str[])
116{
117  int  k;
118
119  for (k = 0; k < (int)strlen(str); k++)
120    str[k] = tolower(str[k]);
121}
122
123
124void upr_case(char str[])
125{
126  int  k;
127
128  for (k = 0; k < (int)strlen(str); k++)
129    str[k] = toupper(str[k]);
130}
131
132
133// only supported by Redhat
134#if 0
135// generate backtrace
136void prt_trace (void)
137{
138  const int  max_entries = 20;
139
140  void    *array[max_entries];
141  size_t  size;
142  char    **strings;
143  size_t  i;
144     
145  size = backtrace(array, max_entries);
146  strings = backtrace_symbols(array, size);
147     
148  printf("prt_trace: obtained %zd stack frames\n", size);
149     
150  for (i = 0; i < size; i++)
151    printf ("%s\n", strings[i]);
152     
153  free (strings);
154}
155#endif
156
157
158// file I/O
159
160// C++
161
162void file_rd(ifstream &inf, const char file_name[])
163{
164
165  inf.open(file_name, ios::in);
166  if (!inf.is_open()) {
167    printf("File not found: %s\n", file_name);
168    exit_(-1);
169  }
170}
171
172
173void file_wr(ofstream &outf, const char file_name[])
174{
175
176  outf.open(file_name, ios::out);
177  if (!outf.is_open()) {
178    printf("Could not create file: %s\n", file_name);
179    exit_(-1);
180  }
181}
182
183
184// C
185/***********************************************
186FILE* file_read(const char file_name[])
187
188   Purpose:
189      Open a file, and return the file ID
190************************************************/
191FILE* file_read(const char file_name[])
192{
193  FILE      *fp;
194 
195  fp = fopen(file_name, "r");
196  if (fp == NULL) {
197    printf("File not found: %s\n", file_name);
198    exit_(-1);
199  } else
200    return(fp);
201    // avoid compiler warning
202    return NULL;
203}
204
205
206FILE* file_write(const char file_name[])
207{
208  FILE      *fp;
209 
210  fp = fopen(file_name, "w");
211  if (fp == NULL) {
212    printf("Could not create file: %s\n", file_name);
213    exit_(-1);
214    // avoid compiler warning
215    return NULL;
216  } else
217    return(fp);
218}
219
220/* Check for closed orbit if not exit the program*/
221void chk_cod(const bool cod, const char *proc_name)
222{
223
224  if (!cod) {
225    printf("%s: closed orbit not found\n", proc_name);
226    exit_(1);
227  }
228}
229
230
231void no_sxt(void)
232{
233  int       k;
234
235  cout << endl;
236  cout << "zeroing sextupoles" << endl;
237  for (k = 0; k <= globval.Cell_nLoc; k++)
238    if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->Porder >= Sext))
239      SetKpar(Cell[k].Fnum, Cell[k].Knum, Sext, 0.0);
240}
241
242
243void get_map(void)
244{
245  long int  lastpos;
246
247
248  getcod(0.0, lastpos);
249
250  map.identity(); map += globval.CODvect;
251  Cell_Pass(0, globval.Cell_nLoc, map, lastpos);
252  map -= globval.CODvect;
253}
254
255
256#if ORDER > 1
257
258tps get_h(void)
259{
260  ss_vect<tps>  map1, R;
261
262  // Parallel transport nonlinear kick to start of lattice,
263  // assumes left to right evaluation.
264
265  if (true)
266    // Dragt-Finn factorization
267    return LieFact_DF(Inv(MNF.A0*MNF.A1)*map*MNF.A0*MNF.A1, R)*R;
268  else {
269    // single Lie exponent
270    danot_(1); map1 = map; danot_(no_tps);
271    return LieFact(Inv(MNF.A0*MNF.A1)*map*Inv(map1)*MNF.A0*MNF.A1);
272  }
273}
274
275#endif
276
277void get_m2(const ss_vect<tps> &ps, tps m2[])
278{
279  int  i, j, k;
280
281  k = 0;
282  for (i = 0; i < 2*nd_tps; i++)
283    for (j = i; j < 2*nd_tps; j++) {
284      k++; m2[k-1] = ps[i]*ps[j];
285    }
286}
287
288
289ss_vect<tps> get_S(const int n_DOF)
290{
291  int           j;
292  ss_vect<tps>  S;
293
294  S.zero();
295  for (j = 0; j < n_DOF; j++) {
296    S[2*j] = tps(0.0, 2*j+2); S[2*j+1] = -tps(0.0, 2*j+1);
297  }
298
299  return S;
300}
301
302
303ss_vect<tps> tp_S(const int n_DOF, const ss_vect<tps> &A)
304{
305  int           j, jj[ss_dim];
306  ss_vect<tps>  S;
307
308  for (j = 1; j <= ss_dim; j++)
309    jj[j-1] = (j <= 2*n_DOF)? 1 : 0;
310
311  S = get_S(n_DOF);
312
313  return -S*PInv(A, jj)*S;
314}
315
316
317void get_dnu(const ss_vect<tps> &A, double dnu[])
318{
319  int  k;
320
321  for (k = 0; k <= 1; k++) {
322    dnu[k] = atan2(A[2*k][2*k+1], A[2*k][2*k])/(2.0*M_PI);
323    if (dnu[k] < 0.0) dnu[k] += 1.0;
324  }
325}
326
327
328ss_vect<tps> get_A_CS(const ss_vect<tps> &A, double dnu[])
329{
330  int           k;
331  double        c, s;
332  ss_vect<tps>  Id, R;
333
334  Id.identity(); R.identity(); get_dnu(A, dnu);
335  for (k = 0; k <= 1; k++) {
336
337    c = cos(2.0*M_PI*dnu[k]); s = sin(2.0*M_PI*dnu[k]);
338    R[2*k] = c*Id[2*k] - s*Id[2*k+1]; R[2*k+1] = s*Id[2*k] + c*Id[2*k+1];
339  }
340
341  return A*R;
342}
343
344
345void get_twoJ(const int n_DOF, const ss_vect<double> &ps,
346              const ss_vect<tps> &A, double twoJ[])
347{
348  int              j;
349  iVector          jj;
350  ss_vect<double>  z;
351
352  for (j = 0; j < nv_tps; j++)
353    jj[j] = (j < 2*n_DOF)? 1 : 0;
354
355  z = (PInv(A, jj)*ps).cst();
356
357  for (j = 0; j < n_DOF; j++)
358    twoJ[j] = sqr(z[2*j]) + sqr(z[2*j+1]);
359}
360
361
362double get_curly_H(const double alpha_x, const double beta_x,
363                   const double eta_x, const double etap_x)
364{
365  double  curly_H, gamma_x;
366
367  gamma_x = (1.0+pow(alpha_x, 2))/beta_x;
368
369  curly_H = gamma_x*sqr(eta_x) + 2.0*alpha_x*eta_x*etap_x + beta_x*sqr(etap_x);
370
371  return curly_H;
372}
373
374
375double get_eps_x(void)
376{
377  bool             cav, emit;
378  long int         lastpos;
379  double           eps_x;
380  ss_vect<tps>     A;
381
382  /* Note:
383
384        T
385       M  J M = J,
386
387        -1       T           |  0  I |        T   | beta   -alpha |
388       A   = -J A  J,    J = |       |,    A A  = |               |
389                             | -I  0 |            | -alpha  gamma |
390
391     Transform to Floquet Space:
392
393        -1           T
394       A   eta = -J A  J eta,
395
396               -1      T  -1                T    T
397       H~ = ( A   eta )  A   eta = ( J eta )  A A  ( J eta )
398
399  */
400
401  cav = globval.Cavity_on; emit = globval.emittance;
402
403  globval.Cavity_on = false; globval.emittance = false;
404
405  Ring_GetTwiss(false, 0.0);
406
407  putlinmat(6, globval.Ascr, A); A += globval.CODvect;
408
409  globval.emittance = true;
410
411  Cell_Pass(0, globval.Cell_nLoc, A, lastpos);
412
413  eps_x = 1470.0*pow(globval.Energy, 2)*I5/(I2-I4);
414
415  printf("\n");
416  printf("eps_x = %5.3f nm.rad\n", eps_x);
417  printf("J_x   = %5.3f, J_z = %5.3f\n", 1.0-I4/I2, 2.0+I4/I2);
418
419  globval.Cavity_on = cav; globval.emittance = emit;
420
421  return eps_x;
422}
423
424/****************************************************************************/
425/* void GetEmittance(const int Fnum, const bool prt)
426
427   Purpose:
428       get the emittance
429
430   Input:
431      Fnum  family index
432      prt   bool flag, true = print the calcuate information
433   Output:
434       none
435
436   Return:
437       emittance
438       
439   Global variables:
440       
441
442   specific functions:
443       
444
445   Comments:
446       
447
448****************************************************************************/
449void GetEmittance(const int Fnum, const bool prt)
450{
451  bool          emit, rad, cav, path;
452  int           i, j, h_RF;
453  long int      lastpos, loc;
454  double        C, theta, V_RF, phi0, alpha_z, beta_z, gamma_z;
455  double        sigma_s, sigma_delta;
456  Vector3       nu;
457  Matrix        Ascr;
458  ss_vect<tps>  Ascr_map;
459
460  // save state
461  rad = globval.radiation; emit = globval.emittance;
462  cav = globval.Cavity_on; path = globval.pathlength;
463
464  C = Cell[globval.Cell_nLoc].S; /* Circumference of the ring*/
465
466  // damped system
467  globval.radiation = true; globval.emittance  = true;
468  globval.Cavity_on = true; globval.pathlength = false;
469
470  Ring_GetTwiss(false, 0.0);
471
472  // radiation loss is computed in Cav_Pass
473
474  // Sum_k( alpha_k) = 2 * U_0 / E
475//  printf("\n");
476//  printf("%24.16e %24.16e\n",
477//       globval.dE,
478//       (globval.alpha_rad[X_]+globval.alpha_rad[Y_]
479//       +globval.alpha_rad[Z_])/2.0);
480
481  globval.U0 = 1e9*globval.dE*globval.Energy;
482  V_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Pvolt; //RF voltage
483  h_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Ph;  // RF cavity harmonic number
484  phi0 = fabs(asin(globval.U0/V_RF));
485  globval.delta_RF =
486    sqrt(-V_RF*cos(M_PI-phi0)*(2.0-(M_PI-2.0*(M_PI-phi0))
487    *tan(M_PI-phi0))/(abs(globval.Alphac)*M_PI*h_RF*1e9*globval.Energy));
488
489  // compute diffusion coeffs. for eigenvectors [sigma_xx, sigma_yy, sigma_zz]
490  putlinmat(6, globval.Ascr, Ascr_map); Ascr_map += globval.CODvect;
491
492  Cell_Pass(0, globval.Cell_nLoc, Ascr_map, lastpos);
493
494  for (i = 0; i < DOF; i++) {
495    // partition numbers
496    globval.J[i] = 2.0*(1.0+globval.CODvect[delta_])*globval.alpha_rad[i]
497                   /globval.dE;
498    // damping times
499    globval.tau[i] = -C/(c0*globval.alpha_rad[i]);
500    // diffusion coeff. and emittance
501    globval.eps[i] = -globval.D_rad[i]/(2.0*globval.alpha_rad[i]);
502    // fractional tunes
503    nu[i]  = atan2(globval.wi[i*2], globval.wr[i*2])/(2.0*M_PI);
504    if (nu[i] < 0.0) nu[i] = 1.0 + nu[i];
505  }
506
507  // Note, J_x + J_y + J_z not exactly 4 (1st order perturbations)
508//  printf("\n");
509//  printf("%24.16e\n", globval.J[X_]+globval.J[Y_]+globval.J[Z_]);
510 
511  // undamped system
512  globval.radiation = false; globval.emittance = false;
513
514  Ring_GetTwiss(false, 0.0);
515
516  /* compute the sigmas arround the lattice:
517
518       sigma = A diag[J_1, J_1, J_2, J_2, J_3, J_3] A^T
519
520  */
521  for (i = 0; i < 6; i++) {
522    Ascr_map[i] = tps(globval.CODvect[i]);
523    for (j = 0; j < 6; j++)
524      Ascr_map[i] += globval.Ascr[i][j]*sqrt(globval.eps[j/2])*tps(0.0, j+1);
525  }
526  for (loc = 0; loc <= globval.Cell_nLoc; loc++) {
527    Elem_Pass(loc, Ascr_map);
528    // sigma = A x A^tp
529    getlinmat(6, Ascr_map, Cell[loc].sigma); TpMat(6, Cell[loc].sigma);
530    getlinmat(6, Ascr_map, Ascr); MulLMat(6, Ascr, Cell[loc].sigma);
531  }
532
533  theta = atan2(2e0*Cell[0].sigma[x_][y_],
534          (Cell[0].sigma[x_][x_]-Cell[0].sigma[y_][y_]))/2e0;
535
536  // longitudinal alpha and beta
537  alpha_z =
538    -globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
539    - globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
540  beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
541  gamma_z = (1.0+sqr(alpha_z))/beta_z;
542
543  // bunch size
544  sigma_s = sqrt(beta_z*globval.eps[Z_]);
545  sigma_delta = sqrt(gamma_z*globval.eps[Z_]);
546
547  if (prt) {
548    printf("\n");
549    printf("Emittance:\n");
550    printf("\n");
551    printf("Energy loss per turn [keV]:     "
552           "U0          = %3.1f\n",
553           1e-3*globval.U0);
554    printf("Synchronous phase [deg]:        "
555           "phi0        = 180 - %4.2f\n",
556           phi0*180.0/M_PI);
557    printf("RF bucket height [%%]:           "
558           "delta_RF    = %4.2f\n", 1e2*globval.delta_RF);
559    printf("\n");
560    printf("Equilibrium emittance [m.rad]:  "
561           "eps_x       = %9.3e, eps_y  = %9.3e, eps_z  = %9.3e\n",
562            globval.eps[X_], globval.eps[Y_], globval.eps[Z_]);
563    printf("Bunch length [mm]:              "
564           "sigma_s     = %5.3f\n", 1e3*sigma_s);
565    printf("Momentum spread:                "
566           "sigma_delta = %9.3e\n", sigma_delta);
567    printf("Partition numbers:              "
568           "J_x         = %5.3f,     J_y    = %5.3f,     J_z    = %5.3f\n",
569            globval.J[X_], globval.J[Y_], globval.J[Z_]);
570    printf("Damping times [msec]:           "
571           "tau_x       = %3.1f,      tau_y  = %3.1f,      tau_z  = %3.1f\n",
572           1e3*globval.tau[X_], 1e3*globval.tau[Y_], 1e3*globval.tau[Z_]);
573    printf("\n");
574    printf("alphac:                         "
575           "alphac      = %8.4e\n", globval.Alphac); 
576    printf("\n");
577    printf("Fractional tunes:               "
578           "nu_x        = %7.5f, nu_y   = %7.5f, nu_z   = %7.5f\n",
579           nu[X_], nu[Y_], nu[Z_]);
580    printf("                                "
581           "1-nu_x      = %7.5f, 1-nu_y = %7.5f, 1-nu_z = %7.5f\n",
582           1e0-nu[X_], 1e0-nu[Y_], 1e0-nu[Z_]);
583    printf("\n");
584    printf("sigmas:                         "
585           "sigma_x     = %5.1f microns, sigma_px    = %5.1f urad\n",
586           1e6*sqrt(Cell[0].sigma[x_][x_]), 1e6*sqrt(Cell[0].sigma[px_][px_]));
587    printf("                                "
588           "sigma_y     = %5.1f microns, sigma_py    = %5.1f urad\n",
589           1e6*sqrt(Cell[0].sigma[y_][y_]), 1e6*sqrt(Cell[0].sigma[py_][py_]));
590    printf("                                "
591           "sigma_s     = %5.2f mm,      sigma_delta = %8.2e\n",
592           1e3*sqrt(Cell[0].sigma[ct_][ct_]),
593           sqrt(Cell[0].sigma[delta_][delta_]));
594
595    printf("\n");
596    printf("Beam ellipse twist [rad]:       tw = %5.3f\n", theta);
597    printf("                   [deg]:       tw = %5.3f\n", theta*180.0/M_PI);
598  }
599
600  // restore state
601  globval.radiation = rad; globval.emittance  = emit;
602  globval.Cavity_on = cav; globval.pathlength = path;
603}
604
605
606// output
607
608void prt_lat(const char *fname, const int Fnum, const bool all)
609{
610  long int      i = 0;
611  double        I2, I5, code = 0.0;
612  FILE          *outf;
613
614  outf = file_write(fname);
615  fprintf(outf, "#        name           s   code"
616                "  alphax  betax   nux   etax   etapx");
617  fprintf(outf, "  alphay  betay   nuy   etay   etapy    I5\n");
618  fprintf(outf, "#                      [m]"
619                "                 [m]           [m]");
620  fprintf(outf, "                   [m]           [m]\n");
621  fprintf(outf, "#\n");
622
623  I2 = 0.0; I5 = 0.0;
624  for (i = 0; i <= globval.Cell_nLoc; i++) {
625    if (all || (Cell[i].Fnum == Fnum)) {
626      switch (Cell[i].Elem.Pkind) {
627      case drift:
628        code = 0.0;
629        break;
630      case Mpole:
631        if (Cell[i].Elem.M->Pirho != 0.0)
632          code = 0.5;
633        else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0)
634          code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]);
635        else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0)
636          code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]);
637        else if (Cell[i].Fnum == globval.bpm)
638          code = 2.0;
639        else
640          code = 0.0;
641        break;
642      default:
643        code = 0.0;
644        break;
645      }
646      fprintf(outf, "%4ld %15s %6.2f %4.1f"
647              " %7.3f %6.3f %6.3f %6.3f %6.3f"
648              " %7.3f %6.3f %6.3f %6.3f %6.3f  %8.2e\n",
649              i, Cell[i].Elem.PName, Cell[i].S, code,
650              Cell[i].Alpha[X_], Cell[i].Beta[X_], Cell[i].Nu[X_],
651              Cell[i].Eta[X_], Cell[i].Etap[X_],
652              Cell[i].Alpha[Y_], Cell[i].Beta[Y_], Cell[i].Nu[Y_],
653              Cell[i].Eta[Y_], Cell[i].Etap[Y_], I5);
654    }
655  }
656
657//  fprintf(outf, "\n");
658//  fprintf(outf, "# emittance: %5.3f nm.rad\n", get_eps_x());
659
660  fclose(outf);
661}
662
663
664void prt_chrom_lat(void)
665{
666  long int  i;
667  double    dbeta_ddelta[Cell_nLocMax][2], detax_ddelta[Cell_nLocMax];
668  double    ksi[Cell_nLocMax][2];
669  double    code = 0.0;
670  FILE      *outf;
671
672  printf("\n");
673  printf("prt_chrom_lat: calling Ring_GetTwiss with delta != 0\n");
674  Ring_GetTwiss(true, globval.dPcommon);
675  for (i = 0; i <= globval.Cell_nLoc; i++) {
676    dbeta_ddelta[i][X_] = Cell[i].Beta[X_];
677    dbeta_ddelta[i][Y_] = Cell[i].Beta[Y_];
678    detax_ddelta[i] = Cell[i].Eta[X_];
679  }
680  printf("prt_chrom_lat: calling Ring_GetTwiss with delta != 0\n");
681  Ring_GetTwiss(true, -globval.dPcommon);
682  ksi[0][X_] = 0.0; ksi[0][Y_] = 0.0;
683  for (i = 0; i <= globval.Cell_nLoc; i++) {
684    dbeta_ddelta[i][X_] -= Cell[i].Beta[X_];
685    dbeta_ddelta[i][Y_] -= Cell[i].Beta[Y_];
686    detax_ddelta[i] -= Cell[i].Eta[X_];
687    dbeta_ddelta[i][X_] /= 2.0*globval.dPcommon;
688    dbeta_ddelta[i][Y_] /= 2.0*globval.dPcommon;
689    detax_ddelta[i] /= 2.0*globval.dPcommon;
690    if (i != 0) {
691      ksi[i][X_] = ksi[i-1][X_]; ksi[i][Y_] = ksi[i-1][Y_];
692    }
693    if (Cell[i].Elem.Pkind == Mpole) {
694        ksi[i][X_] -= Cell[i].Elem.M->PBpar[Quad+HOMmax]
695                     *Cell[i].Elem.PL*Cell[i].Beta[X_]/(4.0*M_PI);
696        ksi[i][Y_] += Cell[i].Elem.M->PBpar[Quad+HOMmax]
697                     *Cell[i].Elem.PL*Cell[i].Beta[Y_]/(4.0*M_PI);
698    }
699  }
700
701  outf = file_write("chromlat.out");
702  fprintf(outf, "#     name              s    code"
703                "  bx*ex  sqrt(bx*by)  dbx/dd*ex  bx*dex/dd"
704                "  by*ex  dby/dd*ex by*dex/dd  ksix  ksiy"
705                "  dbx/dd  bx/dd dex/dd\n");
706  fprintf(outf, "#                      [m]          [m]"
707                "      [m]          [m]       [m]");
708  fprintf(outf, "       [m]      [m]       [m]\n");
709  fprintf(outf, "#\n");
710  for (i = 0; i <= globval.Cell_nLoc; i++) {
711    switch (Cell[i].Elem.Pkind) {
712    case drift:
713      code = 0.0;
714      break;
715    case Mpole:
716      if (Cell[i].Elem.M->Pirho != 0)
717        code = 0.5;
718      else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0)
719        code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]);
720      else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0)
721        code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]);
722      else if (Cell[i].Fnum == globval.bpm)
723        code = 2.0;
724      else
725        code = 0.0;
726      break;
727    default:
728      code = 0.0;
729      break;
730    }
731    fprintf(outf, "%4ld %15s %6.2f %4.1f"
732                  "  %6.3f  %8.3f    %8.3f   %8.3f"
733                  "   %6.3f %8.3f   %8.3f  %5.2f  %5.2f"
734                  "  %6.3f  %6.3f  %6.3f\n",
735            i, Cell[i].Elem.PName, Cell[i].S, code,
736            Cell[i].Beta[X_]*Cell[i].Eta[X_],
737            sqrt(Cell[i].Beta[X_]*Cell[i].Beta[Y_]),
738            dbeta_ddelta[i][X_]*Cell[i].Eta[X_],
739            detax_ddelta[i]*Cell[i].Beta[X_],
740            Cell[i].Beta[Y_]*Cell[i].Eta[X_],
741            dbeta_ddelta[i][Y_]*Cell[i].Eta[X_],
742            detax_ddelta[i]*Cell[i].Beta[Y_],
743            ksi[i][X_], ksi[i][Y_],
744            dbeta_ddelta[i][X_], dbeta_ddelta[i][Y_], detax_ddelta[i]);
745  }
746  fclose(outf);
747}
748
749/****************************************************************************/
750/* void prt_cod(const char *file_name, const int Fnum, const bool all)
751
752   Purpose:
753       print the location, twiss parameters, displacement errors and close orbit
754       at the the location of Fnum or all lattice elements
755
756   Input:
757      file_name    file to save the data
758                    code in the file:
759                      0   drfit
760                      0.5  dipole
761                      -1    defocusing quadrupole
762                      1     focusing quadrupole
763                      1.5  sextupole
764                      2.0  bpm
765      Fnum        family index
766      all         true, print the cod at all elements
767                  false, print the cod at the family elements
768   Output:
769       none
770
771   Return:
772       none
773       
774   Global variables:
775       
776
777   specific functions:
778       
779
780   Comments:
781       
782
783****************************************************************************/
784void prt_cod(const char *file_name, const int Fnum, const bool all)
785{
786  long      i;
787  double    code = 0.0;
788  FILE      *outf;
789  long      FORLIM;
790  struct    tm *newtime;
791
792  outf = file_write(file_name);
793
794  /* Get time and date */
795  newtime = GetTime();
796
797  fprintf(outf,"# TRACY III -- %s -- %s \n",
798          file_name, asctime2(newtime));
799
800  fprintf(outf, "#       name             s  code  betax   nux   betay   nuy"
801          "   xcod   ycod    dSx    dSy   dipx   dipy\n");
802  fprintf(outf, "#                       [m]        [m]           [m]       "
803          "   [mm]   [mm]    [mm]   [mm] [mrad]  [mrad]\n");
804  fprintf(outf, "#\n");
805
806  FORLIM = globval.Cell_nLoc;
807  for (i = 0L; i <= FORLIM; i++) {
808    if (all || (Cell[i].Fnum == Fnum)) {
809      switch (Cell[i].Elem.Pkind) {
810      case drift:
811        code = 0.0;
812        break;
813      case Mpole:
814        if (Cell[i].Elem.M->Pirho != 0)
815          code = 0.5;
816        else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0)
817          code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]);
818        else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0)
819          code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]);
820        else if (Cell[i].Fnum == globval.bpm)
821          code = 2.0;
822        else
823          code = 0.0;
824        break;
825      default:
826        code = 0.0;
827        break;
828      }
829      /* COD is in local coordinates */
830      fprintf(outf, "%4ld %.*s %6.2f %4.1f %6.3f %6.3f %6.3f %6.3f"
831              " %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n",
832              i, SymbolLength, Cell[i].Elem.PName, Cell[i].S, code,
833              Cell[i].Beta[X_], Cell[i].Nu[X_],
834              Cell[i].Beta[Y_], Cell[i].Nu[Y_], 
835              1e3*Cell[i].BeamPos[x_], 1e3*Cell[i].BeamPos[y_],
836              1e3*Cell[i].dS[X_], 1e3*Cell[i].dS[Y_], //displacement error
837              -1e3*Elem_GetKval(Cell[i].Fnum, Cell[i].Knum, Dip),  //H kick angle, for correctors
838              1e3*Elem_GetKval(Cell[i].Fnum, Cell[i].Knum, -Dip)); //kick angle, for correctors
839    }
840  }
841  fclose(outf);
842}
843
844
845void prt_beampos(const char *file_name)
846{
847  int       k;
848  ofstream  outf;
849
850  file_wr(outf, file_name);
851
852  outf << "# k  s  name x   px    y   py   delta ct" << endl;
853  outf << "#" << endl;
854
855  for (k = 0; k <= globval.Cell_nLoc; k++)
856    outf << scientific << setprecision(5)
857         << setw(5) << k << setw(11) << Cell[k].Elem.PName
858         << setw(13) << Cell[k].BeamPos << endl;
859
860  outf.close();
861}
862
863
864/****************************************************************
865void CheckAlignTol(const char *OutputFile)
866
867  Purpose:
868    check aligment errors of individual magnets on giders
869    the dT and roll angle are all printed out
870****************************************************************/
871// misalignments
872void CheckAlignTol(const char *OutputFile)
873  // check aligment errors of individual magnets on giders
874  // the dT and roll angle are all printed out
875{
876  int i,j;
877  int  n_girders;
878  int  gs_Fnum, ge_Fnum;
879  int  gs_nKid, ge_nKid;
880  int  dip_Fnum,dip_nKid;
881  int  loc, loc_gs, loc_ge;
882  char * name;
883  double s;
884  double PdSsys[2], PdSrms[2], PdSrnd[2], dS[2], dT[2];
885  fstream fout;
886 
887  gs_Fnum = globval.gs;   gs_nKid = GetnKid(gs_Fnum);
888  ge_Fnum = globval.ge;   ge_nKid = GetnKid(ge_Fnum);
889  if (gs_nKid == ge_nKid)
890    n_girders= gs_nKid;
891  else {
892    cout << " The numbers of GS and GE not same. " << endl;
893    exit (1);
894  }
895
896  fout.open(OutputFile,ios::out);
897  if(!fout) {
898    cout << "error in opening the file  " << endl;
899    exit_(0);
900  }
901 
902  fout << "Girders, Quads, Sexts:  " << endl; 
903  for (i = 1; i <= n_girders; i++){
904    fout << i << ":" << endl;
905    loc_gs = Elem_GetPos(gs_Fnum, i); loc_ge = Elem_GetPos(ge_Fnum, i);
906   
907    loc = loc_gs;
908    PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
909    PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
910    PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
911    PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
912    PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
913    PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
914    dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
915    dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
916    s = Cell[loc].S; name = Cell[loc].Elem.PName;
917    fout << "  " << name << "  " << loc << "   " << s
918         << "  " <<  PdSsys[X_] << "  " <<  PdSsys[Y_]
919         << "   " << PdSrms[X_] << "  " <<  PdSrms[Y_]
920         << "   " << PdSrnd[X_] << "  " <<  PdSrnd[Y_]
921         << "   " << Cell[loc].Elem.M->PdTrms << "  "
922         << Cell[loc].Elem.M->PdTrnd << "   " << dS[X_]     << "  " <<  dS[Y_]
923         << "   " << atan2( dT[1], dT[0] )  << endl;
924   
925    for (j = loc_gs+1; j < loc_ge; j++) {
926      if ((Cell[j].Elem.Pkind == Mpole) &&
927          (Cell[j].Elem.M->n_design >= Quad || 
928           Cell[j].Elem.M->n_design >= Sext)) {
929        loc = j;
930        PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
931        PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
932        PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
933        PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
934        PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
935        PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
936        dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
937        dT[0] = Cell[loc].dT[0];   dT[1] = Cell[loc].dT[1];
938        s = Cell[loc].S; name=Cell[loc].Elem.PName;
939        fout << "  " << name << "  " << loc << "   " << s
940             << "  " <<  PdSsys[X_] << "  " <<  PdSsys[Y_]
941             << "   " << PdSrms[X_] << "  " <<  PdSrms[Y_]
942             << "   " << PdSrnd[X_] << "  " <<  PdSrnd[Y_]
943             << "   " << Cell[loc].Elem.M->PdTrms << "  "
944             << Cell[loc].Elem.M->PdTrnd
945             << "   " << dS[X_] << "  " <<  dS[Y_]
946             << "   " << atan2( dT[1], dT[0] )  << endl;
947      }
948    }
949   
950    loc = loc_ge;
951    PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
952    PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
953    PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
954    PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
955    PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
956    PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
957    dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
958    dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
959    s=Cell[loc].S; name=Cell[loc].Elem.PName;
960    fout << "  " << name << "  " << loc << "   " << s
961         << "  " <<  PdSsys[X_] << "  " <<  PdSsys[Y_]
962         << "   " << PdSrms[X_] << "  " <<  PdSrms[Y_]
963         << "   " << PdSrnd[X_] << "  " <<  PdSrnd[Y_]
964         << "   " << Cell[loc].Elem.M->PdTrms
965         << "  " << Cell[loc].Elem.M->PdTrnd
966         << "   " << dS[X_]     << "  " <<  dS[Y_]
967         << "   " << atan2( dT[1], dT[0] )  << endl;
968
969  }
970
971  fout << "  " << endl; 
972  fout << "Dipoles:  " << endl;
973  dip_Fnum = ElemIndex("B1"); dip_nKid = GetnKid(dip_Fnum);
974  for (i = 1; i <= dip_nKid; i++){
975    loc = Elem_GetPos(dip_Fnum, i);
976    PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
977    PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
978    PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
979    PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
980    PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
981    PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
982    dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
983    dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
984    s = Cell[loc].S; name = Cell[loc].Elem.PName;
985    fout << "  " << name << "  " << loc << "   " << s
986         << "  " <<  PdSsys[X_] << "  " <<  PdSsys[Y_]
987         << "   " << PdSrms[X_] << "  " <<  PdSrms[Y_]
988         << "   " << PdSrnd[X_] << "  " <<  PdSrnd[Y_]
989         << "   " << Cell[loc].Elem.M->PdTrms
990         << "  " << Cell[loc].Elem.M->PdTrnd
991         << "   " << dS[X_]     << "  " <<  dS[Y_]
992         << "   " << atan2( dT[1], dT[0] )  << endl;
993  }
994 
995  fout.close();
996} 
997
998/*********************************************************************
999void misalign_rms_elem(const int Fnum, const int Knum,
1000                       const double dx_rms, const double dy_rms,
1001                       const double dr_rms, const bool new_rnd)
1002                       
1003  Purpose:
1004     Set random misalignment to the element with Fnum and Knum
1005 
1006  Input:   
1007     Fnum          family number
1008     Knum          kid number
1009     dx_rms         rms value of the error in x
1010     dy_rms         rms value of the error in y
1011     dr_rms         rms value of the error in rotation around s
1012     new_rnd        flag to turn on/off using the new random seed
1013     
1014**********************************************************************/
1015void misalign_rms_elem(const int Fnum, const int Knum,
1016                       const double dx_rms, const double dy_rms,
1017                       const double dr_rms, const bool new_rnd)
1018{
1019  long int   loc;
1020  MpoleType  *mp;
1021
1022  const bool  normal = true;
1023
1024  loc = Elem_GetPos(Fnum, Knum); mp = Cell[loc].Elem.M;
1025
1026  mp->PdSrms[X_] = dx_rms; mp->PdSrms[Y_] = dy_rms; mp->PdTrms = dr_rms; 
1027  if (new_rnd) {
1028    if (normal) {
1029      mp->PdSrnd[X_] = normranf(); mp->PdSrnd[Y_] = normranf();
1030      mp->PdTrnd = normranf();
1031    } else {
1032      mp->PdSrnd[X_] = ranf(); mp->PdSrnd[Y_] = ranf();
1033      mp->PdTrnd = ranf();
1034    }
1035  }
1036
1037  Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum);
1038}
1039
1040
1041/*********************************************************************
1042void misalign_sys_elem(const int Fnum, const int Knum,
1043                       const double dx_sys, const double dy_sys,
1044                       const double dr_sys)
1045                       
1046  Purpose:
1047     Set systematic misalignment to the elements in "type"
1048 
1049  Input:   
1050     Fnum          family number
1051     Knum          kid number         
1052     dx_sys         sys value of the error in x
1053     dy_sys         sys value of the error in y
1054     dr_sys         sys value of the error in rotation around s
1055     
1056**********************************************************************/
1057void misalign_sys_elem(const int Fnum, const int Knum,
1058                       const double dx_sys, const double dy_sys,
1059                       const double dr_sys)
1060{
1061  long int   loc;
1062  MpoleType  *mp;
1063
1064  loc = Elem_GetPos(Fnum, Knum); mp = Cell[loc].Elem.M;
1065
1066  mp->PdSsys[X_] = dx_sys; mp->PdSsys[Y_] = dy_sys; mp->PdTsys = dr_sys; 
1067
1068  Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum);
1069}
1070
1071/*********************************************************************
1072void misalign_rms_fam(const int Fnum,
1073                      const double dx_rms, const double dy_rms,
1074                      const double dr_rms, const bool new_rnd)
1075                       
1076  Purpose:
1077     Set systematic misalignment to the elements on the girders
1078 
1079  Input:   
1080     Fnum           family number                     
1081     dx_rms         rms value of the error in x
1082     dy_rms         rms value of the error in y
1083     dr_rms         rms value of the error in rotation around s
1084     new_rnd        flag to generate random number
1085**********************************************************************/
1086void misalign_rms_fam(const int Fnum,
1087                      const double dx_rms, const double dy_rms,
1088                      const double dr_rms, const bool new_rnd)
1089{
1090  int  i;
1091
1092  for (i = 1; i <= GetnKid(Fnum); i++)
1093    misalign_rms_elem(Fnum, i, dx_rms, dy_rms, dr_rms, new_rnd);
1094}
1095
1096void misalign_sys_fam(const int Fnum,
1097                      const double dx_sys, const double dy_sys,
1098                      const double dr_sys)
1099{
1100  int  i;
1101
1102  for (i = 1; i <= GetnKid(Fnum); i++)
1103    misalign_sys_elem(Fnum, i, dx_sys, dy_sys, dr_sys);
1104}
1105
1106/*********************************************************************
1107void misalign_rms_type(const int type,
1108                       const double dx_rms, const double dy_rms,
1109                       const double dr_rms, const bool new_rnd)
1110                       
1111  Purpose:
1112     Set random misalignment to the elements in "type"
1113 
1114  Input:   
1115     type           element type,has the following value
1116                       All
1117                       dipole
1118                       quad
1119                       sext
1120                       bpm
1121                       girder
1122                       
1123     dx_rms         rms value of the error in x
1124     dy_rms         rms value of the error in y
1125     dr_rms         rms value of the error in rotation around s
1126     new_rnd        flag to turn on/off using the new random seed
1127     
1128**********************************************************************/
1129void misalign_rms_type(const int type,
1130                       const double dx_rms, const double dy_rms,
1131                       const double dr_rms, const bool new_rnd)
1132{
1133  long int   k;
1134
1135  if ((type >= All) && (type <= HOMmax)) {
1136    for (k = 1; k <= globval.Cell_nLoc; k++) {
1137      if ((Cell[k].Elem.Pkind == Mpole) &&
1138          ((type == Cell[k].Elem.M->n_design) ||
1139          ((type == All) &&
1140           ((Cell[k].Fnum != globval.gs) && (Cell[k].Fnum != globval.ge))))) {
1141        // if all: skip girders
1142        misalign_rms_elem(Cell[k].Fnum, Cell[k].Knum,
1143                          dx_rms, dy_rms, dr_rms, new_rnd);
1144      }
1145    }
1146  } else {
1147    printf("misalign_rms_type: incorrect type %d\n", type); exit_(1);
1148  }
1149}
1150
1151/*********************************************************************
1152void misalign_sys_type(const int type,
1153                       const double dx_sys, const double dy_sys,
1154                       const double dr_sys)
1155                       
1156  Purpose:
1157     Set systematic misalignment to the elements in "type"
1158 
1159  Input:   
1160     type           element type,has the following value
1161                       All
1162                       dipole
1163                       quad
1164                       sext
1165                       bpm
1166                       girder
1167                       
1168     dx_sys         sys value of the error in x
1169     dy_sys         sys value of the error in y
1170     dr_sys         sys value of the error in rotation around s
1171     
1172**********************************************************************/
1173void misalign_sys_type(const int type,
1174                       const double dx_sys, const double dy_sys,
1175                       const double dr_sys)
1176{
1177  long int   k;
1178
1179  if ((type >= All) && (type <= HOMmax)) {
1180    for (k = 1; k <= globval.Cell_nLoc; k++) {
1181      if ((Cell[k].Elem.Pkind == Mpole) &&
1182          ((type == Cell[k].Elem.M->n_design) ||
1183          ((type == All) &&
1184           ((Cell[k].Fnum != globval.gs) && (Cell[k].Fnum != globval.ge))))) {
1185        // if all: skip girders
1186        misalign_sys_elem(Cell[k].Fnum, Cell[k].Knum,
1187                          dx_sys, dy_sys, dr_sys);
1188      }
1189    }
1190  } else {
1191    printf("misalign_sys_type: incorrect type %d\n", type); exit_(1);
1192  }
1193}
1194
1195/*********************************************************************
1196void misalign_rms_girders(const int gs, const int ge,
1197                          const double dx_rms, const double dy_rms,
1198                          const double dr_rms, const bool new_rnd)
1199                       
1200  Purpose:
1201     Set systematic misalignment to the elements on the girders
1202 
1203  Input:   
1204     gs             girder start marker
1205     ge             girder end marker                 
1206     dx_rms         rms value of the error in x
1207     dy_rms         rms value of the error in y
1208     dr_rms         rms value of the error in rotation around s
1209     
1210**********************************************************************/
1211void misalign_rms_girders(const int gs, const int ge,
1212                          const double dx_rms, const double dy_rms,
1213                          const double dr_rms, const bool new_rnd)
1214{
1215  int       i, k, n_girders, n_ge, n_gs;
1216  long int  loc_gs, loc_ge, j;
1217  double    s_gs, s_ge, dx_gs[2], dx_ge[2], s;
1218
1219  n_gs = GetnKid(gs); n_ge = GetnKid(ge);
1220
1221  if (n_gs == n_ge)
1222    n_girders = n_gs;
1223  else {
1224    cout << "set_girders: no of GS != no of GE" << endl;
1225    exit (1);
1226  }
1227 
1228  misalign_rms_fam(gs, dx_rms, dy_rms, dr_rms, new_rnd);
1229  misalign_rms_fam(ge, dx_rms, dy_rms, dr_rms, new_rnd);
1230
1231  for (i = 1; i <= n_girders; i++) {
1232    loc_gs = Elem_GetPos(gs, i); loc_ge = Elem_GetPos(ge, i);
1233    s_gs = Cell[loc_gs].S; s_ge = Cell[loc_ge].S;
1234
1235    // roll for a rigid boby
1236    // Note, girders needs to be introduced as gs->ge pairs
1237    Cell[loc_ge].Elem.M->PdTrnd = Cell[loc_gs].Elem.M->PdTrnd;
1238    Mpole_SetdT(ge, i);
1239
1240    for (k = 0; k <= 1; k++) {
1241      dx_gs[k] = Cell[loc_gs].dS[k]; dx_ge[k] = Cell[loc_ge].dS[k];
1242    }
1243
1244    // move elements onto mis-aligned girder
1245    for (j = loc_gs+1; j < loc_ge; j++) {
1246      if ((Cell[j].Elem.Pkind == Mpole) || (Cell[j].Fnum == globval.bpm)) {
1247        s = Cell[j].S;
1248        for (k = 0; k <= 1; k++)
1249          Cell[j].Elem.M->PdSsys[k]
1250            = dx_gs[k] + (dx_ge[k]-dx_gs[k])*(s-s_gs)/(s_ge-s_gs);
1251        Cell[j].Elem.M->PdTsys = 
1252          Cell[loc_gs].Elem.M->PdTrms*Cell[loc_gs].Elem.M->PdTrnd;
1253      }
1254    }
1255  }
1256}
1257
1258/*********************************************************************
1259void misalign_sys_girders(const int gs, const int ge,
1260                          const double dx_sys, const double dy_sys,
1261                          const double dr_sys)
1262                       
1263  Purpose:
1264     Set systematic misalignment to the elements on the girders
1265 
1266  Input:   
1267     gs             girder start marker
1268     ge             girder end marker                 
1269     dx_sys         sys value of the error in x
1270     dy_sys         sys value of the error in y
1271     dr_sys         sys value of the error in rotation around s
1272     
1273**********************************************************************/
1274void misalign_sys_girders(const int gs, const int ge,
1275                          const double dx_sys, const double dy_sys,
1276                          const double dr_sys)
1277{
1278  int       i, k, n_girders, n_ge, n_gs;
1279  long int  loc_gs, loc_ge, j;
1280  double    s_gs, s_ge, dx_gs[2], dx_ge[2], s;
1281
1282  n_gs = GetnKid(gs); n_ge = GetnKid(ge);
1283
1284  if (n_gs == n_ge)
1285    n_girders = n_gs;
1286  else {
1287    cout << "set_girders: no of GS != no of GE" << endl;
1288    exit (1);
1289  }
1290 
1291  misalign_sys_fam(gs, dx_sys, dy_sys, dr_sys);
1292  misalign_sys_fam(ge, dx_sys, dy_sys, dr_sys);
1293
1294  for (i = 1; i <= n_girders; i++) {
1295    loc_gs = Elem_GetPos(gs, i); loc_ge = Elem_GetPos(ge, i);
1296    s_gs = Cell[loc_gs].S; s_ge = Cell[loc_ge].S;
1297
1298    // roll for a rigid boby
1299    // Note, girders needs to be introduced as gs->ge pairs
1300    Cell[loc_ge].Elem.M->PdTrnd = Cell[loc_gs].Elem.M->PdTrnd;
1301    Mpole_SetdT(ge, i);
1302
1303    for (k = 0; k <= 1; k++) {
1304      dx_gs[k] = Cell[loc_gs].dS[k]; dx_ge[k] = Cell[loc_ge].dS[k];
1305    }
1306
1307    // move elements onto mis-aligned girder
1308    for (j = loc_gs+1; j < loc_ge; j++) {
1309      if ((Cell[j].Elem.Pkind == Mpole) || (Cell[j].Fnum == globval.bpm)) {
1310        s = Cell[j].S;
1311        for (k = 0; k <= 1; k++)
1312          Cell[j].Elem.M->PdSsys[k]
1313            = dx_gs[k] + (dx_ge[k]-dx_gs[k])*(s-s_gs)/(s_ge-s_gs);
1314        Cell[j].Elem.M->PdTsys = 
1315          Cell[loc_gs].Elem.M->PdTrms*Cell[loc_gs].Elem.M->PdTrnd;
1316      }
1317    }
1318  }
1319}
1320
1321/***************************************************************************
1322void LoadAlignTol(const char *AlignFile, const bool Scale_it,
1323                  const double Scale, const bool new_rnd, const int k)
1324
1325  Purpose:
1326        Load misalignment error from an external file, then set the error
1327        to the corresponding elements.
1328 
1329  Input:
1330      AlignFile    file from which the misalignment error is readed                       
1331      Scale_it     flag to turn on/off the scale of the error
1332      Scale        value to scale the error
1333      new_rnd      use new random number or not
1334      k   
1335****************************************************************************/
1336void LoadAlignTol(const char *AlignFile, const bool Scale_it,
1337                  const double Scale, const bool new_rnd, const int k)
1338{
1339  char    line[max_str], Name[max_str],  type[max_str];
1340  int     Fnum, seed_val;
1341  double  dx, dy, dr;  // x and y misalignments [m] and roll error [rad]
1342  double  dr_deg;
1343  bool    rms = false, set_rnd;
1344  FILE    *fp;
1345
1346  fp = file_read(AlignFile);
1347
1348  printf("\n");
1349  if (new_rnd)
1350    printf("set alignment errors\n");
1351  else
1352    printf("scale alignment errors: %4.2f\n", Scale);
1353
1354  set_rnd = false;
1355  while (fgets(line, max_str, fp) != NULL) {
1356    //check for whether to set seed
1357    if ((strstr(line, "#") == NULL) && (strcmp(line, "\r\n") != 0)) {
1358      sscanf(line, "%s", Name);
1359      if (strcmp("seed", Name) == 0) {
1360        set_rnd = true;
1361        sscanf(line, "%*s %d", &seed_val);
1362        seed_val += 2*k;
1363        printf("setting random seed to %d\n", seed_val);
1364        iniranf(seed_val); 
1365      } else {
1366        sscanf(line,"%*s %s %lf %lf %lf", type, &dx, &dy, &dr);
1367        dr_deg = dr*180.0/M_PI;
1368
1369        if (strcmp(type, "rms") == 0){
1370          rms = true;
1371          printf("<rms> ");
1372        }
1373        else if (strcmp(type, "sys") == 0){
1374          rms = false;
1375          printf("<sys> ");
1376        }
1377        else {
1378          printf("LoadAlignTol: element %s:  need to specify rms or sys\n",
1379                 Name);
1380          exit_(1);
1381        }
1382
1383        if (rms && !set_rnd) {
1384          printf("LoadAlignTol: seed not defined\n");
1385          exit_(1);
1386        }
1387
1388        if (Scale_it) {
1389          dx *= Scale; dy *= Scale; dr *= Scale;
1390        } 
1391
1392        if (strcmp("all", Name) == 0) {
1393          printf("misaligning all:         dx = %e, dy = %e, dr = %e\n",
1394                 dx, dy, dr);
1395          if(rms)
1396            misalign_rms_type(All, dx, dy, dr_deg, new_rnd);
1397          else
1398            misalign_sys_type(All, dx, dy, dr_deg);
1399        } else if (strcmp("girder", Name) == 0) {
1400          printf("misaligning girders:     dx = %e, dy = %e, dr = %e\n",
1401                 dx, dy, dr);
1402          if (rms)
1403            misalign_rms_girders(globval.gs, globval.ge, dx, dy, dr_deg,
1404                                 new_rnd);
1405          else
1406            misalign_sys_girders(globval.gs, globval.ge, dx, dy, dr_deg);
1407        } else if (strcmp("dipole", Name) == 0) {
1408          printf("misaligning dipoles:     dx = %e, dy = %e, dr = %e\n",
1409                 dx, dy, dr);
1410          if (rms)
1411            misalign_rms_type(Dip, dx, dy, dr_deg, new_rnd);
1412          else
1413            misalign_sys_type(Dip, dx, dy, dr_deg);
1414        } else if (strcmp("quad", Name) == 0) {
1415          printf("misaligning quadrupoles: dx = %e, dy = %e, dr = %e\n",
1416                 dx, dy, dr);
1417          if (rms)
1418            misalign_rms_type(Quad, dx, dy, dr_deg, new_rnd);
1419          else
1420            misalign_sys_type(Quad, dx, dy, dr_deg);
1421        } else if (strcmp("sext", Name) == 0) {
1422          printf("misaligning sextupoles:  dx = %e, dy = %e, dr = %e\n",
1423                 dx, dy, dr);
1424          if (rms)
1425            misalign_rms_type(Sext, dx, dy, dr_deg, new_rnd);
1426          else
1427            misalign_sys_type(Sext, dx, dy, dr_deg);
1428        } else if (strcmp("bpm", Name) == 0) {
1429          printf("misaligning bpms:        dx = %e, dy = %e, dr = %e\n",
1430                 dx, dy, dr);
1431          if (rms)
1432            misalign_rms_fam(globval.bpm, dx, dy, dr_deg, new_rnd);
1433          else
1434            misalign_sys_fam(globval.bpm, dx, dy, dr_deg);
1435        } else {
1436          Fnum = ElemIndex(Name);
1437          if(Fnum > 0) {
1438            printf("misaligning all %s:  dx = %e, dy = %e, dr = %e\n",
1439                   Name, dx, dy, dr);
1440            if (rms)
1441              misalign_rms_fam(Fnum, dx, dy, dr_deg, new_rnd);
1442            else
1443              misalign_sys_fam(Fnum, dx, dy, dr_deg);
1444          } else 
1445            printf("LoadAlignTol: undefined element %s\n", Name);
1446        }
1447      }
1448    } else
1449      printf("%s", line);
1450  }
1451
1452  fclose(fp);
1453}
1454
1455/****************************************************************************/
1456/* void set_aper_elem(const int Fnum, const int Knum,
1457                   const double Dxmin, const double Dxmax,
1458                   const double Dymin, const double Dymax)
1459
1460   Purpose:
1461      Define vacuum chamber to the element 
1462     
1463
1464   Input:
1465      none
1466
1467   Output:
1468       none
1469
1470   Return:
1471       none
1472
1473   Global variables:
1474       none
1475
1476   Specific functions:
1477       none
1478
1479   Comments:
1480      Add k=0 for the default start point "start" in tracy.
1481
1482****************************************************************************/
1483// apertures
1484
1485void set_aper_elem(const int Fnum, const int Knum, 
1486                   const double Dxmin, const double Dxmax, 
1487                   const double Dymin, const double Dymax) 
1488{ 
1489  int  k; 
1490 
1491  if(Fnum==0 && Knum==0)
1492    k=0;
1493  else
1494    k = Elem_GetPos(Fnum, Knum);
1495     
1496    Cell[k].maxampl[X_][0] = Dxmin; 
1497    Cell[k].maxampl[X_][1] = Dxmax; 
1498    Cell[k].maxampl[Y_][0] = Dymin; 
1499    Cell[k].maxampl[Y_][1] = Dymax; 
1500 } 
1501
1502 /****************************************************************************/
1503/* void set_aper_fam(const int Fnum,
1504                  const double Dxmin, const double Dxmax,
1505                  const double Dymin, const double Dymax)
1506   Purpose:
1507      Define vacuum chamber to the family "Fnum" 
1508     
1509
1510   Input:
1511      none
1512     
1513   Output:
1514       none
1515
1516   Return:
1517       none
1518
1519   Global variables:
1520       none
1521
1522   Specific functions:
1523       none
1524
1525   Comments:
1526       Enumeration of special variables:
1527                 enum { All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6 }
1528
1529****************************************************************************/
1530void set_aper_fam(const int Fnum,
1531                  const double Dxmin, const double Dxmax, 
1532                  const double Dymin, const double Dymax)
1533{
1534  int k;
1535
1536  for (k = 1; k <= GetnKid(Fnum); k++)
1537    set_aper_elem(Fnum, k, Dxmin, Dxmax, Dymin, Dymax);
1538}
1539
1540/****************************************************************************/
1541/* void set_aper_type(const int type, const double Dxmin, const double Dxmax,
1542                   const double Dymin, const double Dymax)
1543
1544   Purpose:
1545      Define vacuum chamber to the elements with the same "type", either for "all" element, or "quad","sext", etc. 
1546     
1547
1548   Input:
1549       none
1550       
1551   Output:
1552       none
1553
1554   Return:
1555       none
1556
1557   Global variables:
1558       none
1559
1560   Specific functions:
1561       none
1562
1563   Comments:
1564      1) Enumeration of special variables:
1565                 enum { All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6 }
1566      2) change the start point from k=1 to k=0, in order to set the aperture for
1567         the default start point "begin" which is not a lattice element.
1568
1569****************************************************************************/
1570void set_aper_type(const int type, const double Dxmin, const double Dxmax, 
1571                   const double Dymin, const double Dymax)
1572{
1573  long int   k;
1574
1575  if (type >= All && type <= HOMmax) 
1576    {
1577 //   for(k = 1; k <= globval.Cell_nLoc; k++)
1578      for(k = 0; k <= globval.Cell_nLoc; k++)
1579        if (((Cell[k].Elem.Pkind == Mpole) &&
1580           (Cell[k].Elem.M->n_design == type)) || (type == All))
1581        set_aper_elem(Cell[k].Fnum, Cell[k].Knum, Dxmin, Dxmax, Dymin, Dymax);
1582     } 
1583  else
1584    printf("set_aper_type: bad design type %d\n", type);
1585}
1586
1587
1588/****************************************************************************/
1589/* void LoadApers(const char *AperFile, const double scl_x, const double scl_y)
1590
1591   Purpose:
1592      Read vacuum chamber from file "AperFile", and scale the values with scl_x and scl_y.
1593      Assign vacuum chamber to all elements, if there are definition of
1594      chamber at quad and sext, then assign the specific values at these elements.
1595     
1596
1597   Input:
1598       AperFile:  file with chamber definition
1599       scl_x   :  scale factor of x limitation
1600       scl_y   : scale factor of y limitation
1601
1602   Output:
1603       none
1604
1605   Return:
1606       none
1607
1608   Global variables:
1609       none
1610
1611   Specific functions:
1612       none
1613
1614   Comments:
1615       Enumeration of special variables:
1616                 enum { All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6 }
1617                 
1618    10/2010   Comments by Jianfeng Zhang
1619
1620****************************************************************************/
1621void LoadApers(const char *AperFile, const double scl_x, const double scl_y) 
1622{
1623  char    line[max_str], Name[max_str];
1624  int     Fnum; 
1625  double  dxmin, dxmax, dymin, dymax;  // min and max x and apertures
1626  FILE    *fp;
1627
1628  bool  prt = true;
1629
1630  fp = file_read(AperFile);
1631
1632  printf("\n");
1633  printf("...Load and Set Apertures.\n");
1634
1635  while (fgets(line, max_str, fp) != NULL) {
1636    if (strstr(line, "#") == NULL) 
1637    {
1638      sscanf(line,"%s %lf %lf %lf %lf",
1639             Name, &dxmin, &dxmax, &dymin, &dymax);
1640      dxmin *= scl_x; 
1641      dxmax *= scl_x; 
1642      dymin *= scl_y; 
1643      dymax *= scl_y;
1644     
1645      if (strcmp("all", Name)==0) {
1646        if(prt)
1647          printf("setting all apertures to"
1648                 " dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1649                 dxmin, dxmax, dymin, dymax);
1650        set_aper_type(All, dxmin, dxmax, dymin, dymax);
1651        //      ini_aper(dxmin, dxmax, dymin, dymax);
1652       }
1653       
1654      else if (strcmp("quad", Name)==0) {
1655        if(prt)
1656          printf("setting apertures at all quads to"
1657                 " dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1658                 dxmin, dxmax, dymin, dymax); 
1659        set_aper_type(Quad, dxmin, dxmax, dymin, dymax);
1660       }
1661       
1662      else if (strcmp("sext", Name) == 0) {
1663        if(prt)
1664          printf("setting apertures at all sextupoles to"
1665                 " dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1666                 dxmin, dxmax, dymin, dymax);
1667        set_aper_type(Sext, dxmin, dxmax, dymin, dymax);
1668       }
1669       
1670      else {
1671        Fnum = ElemIndex(Name);
1672        if(Fnum > 0) {
1673          if(prt)
1674            printf("setting apertures at all %s to"
1675                   " dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1676                   Name, dxmin, dxmax, dymin, dymax);
1677          set_aper_fam(Fnum, dxmin, dxmax, dymin, dymax);
1678          } 
1679         else 
1680          printf("LoadApers: lattice does not contain element %s\n", Name);
1681      }
1682    } 
1683   else
1684      printf("%s", line);
1685  }
1686   
1687  fclose(fp);
1688}
1689
1690
1691double get_L(const int Fnum, const int Knum)
1692{
1693  return Cell[Elem_GetPos(Fnum, Knum)].Elem.PL;
1694}
1695
1696
1697void set_L(const int Fnum, const int Knum, const double L)
1698{
1699
1700  Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L;
1701}
1702
1703
1704void set_L(const int Fnum, const double L)
1705{
1706  int  k;
1707
1708  for (k = 1; k <= GetnKid(Fnum); k++)
1709    set_L(Fnum, k, L);
1710}
1711
1712
1713void set_dL(const int Fnum, const int Knum, const double dL)
1714{
1715
1716  Cell[Elem_GetPos(Fnum, Knum)].Elem.PL += dL;
1717}
1718
1719
1720// multipole components
1721/*************************************************************
1722void get_bn_design_elem(const int Fnum, const int Knum,
1723                        const int n, double &bn, double &an)
1724                       
1725   Purpose:
1726       Get the n-th order design value (bn, an) of the element
1727       with family index "Fnum" and the kid number "Knum".
1728       
1729  Input:
1730      Fnum:   family index
1731      Knum:   kid index
1732         n:   n-th design order of the element
1733        bn:   bn component of the element
1734        an:   an component of the element
1735
1736 Output:
1737    None
1738   
1739 Return:
1740          bn, an
1741         
1742 Comments:
1743 
1744            11/2010 comments by jianfeng Zhang
1745**************************************************************/
1746void get_bn_design_elem(const int Fnum, const int Knum,
1747                        const int n, double &bn, double &an)
1748{
1749  elemtype  elem;
1750
1751  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1752
1753  bn = elem.M->PBpar[HOMmax+n]; an = elem.M->PBpar[HOMmax-n];
1754}
1755
1756
1757void get_bnL_design_elem(const int Fnum, const int Knum,
1758                         const int n, double &bnL, double &anL)
1759{
1760  elemtype  elem;
1761
1762  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1763
1764  bnL = elem.M->PBpar[HOMmax+n]; anL = elem.M->PBpar[HOMmax-n];
1765
1766  if (elem.PL != 0.0) {
1767    bnL *= elem.PL; anL *= elem.PL;
1768  }
1769}
1770
1771
1772void set_bn_design_elem(const int Fnum, const int Knum,
1773                        const int n, const double bn, const double an)
1774{
1775  elemtype  elem;
1776
1777  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1778
1779  elem.M->PBpar[HOMmax+n] = bn; elem.M->PBpar[HOMmax-n] = an;
1780
1781  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1782}
1783
1784
1785void set_dbn_design_elem(const int Fnum, const int Knum,
1786                         const int n, const double dbn, const double dan)
1787{
1788  elemtype  elem;
1789
1790  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1791
1792  elem.M->PBpar[HOMmax+n] += dbn; elem.M->PBpar[HOMmax-n] += dan;
1793
1794  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1795}
1796
1797
1798void set_bn_design_fam(const int Fnum,
1799                       const int n, const double bn, const double an)
1800{
1801  int k;
1802
1803  for (k = 1; k <= GetnKid(Fnum); k++)
1804    set_bn_design_elem(Fnum, k, n, bn, an);
1805}
1806
1807
1808void set_dbn_design_fam(const int Fnum,
1809                        const int n, const double dbn, const double dan)
1810{
1811  int k;
1812
1813  for (k = 1; k <= GetnKid(Fnum); k++)
1814    set_dbn_design_elem(Fnum, k, n, dbn, dan);
1815}
1816
1817
1818void set_bnL_design_elem(const int Fnum, const int Knum,
1819                         const int n, const double bnL, const double anL)
1820{
1821  elemtype  elem;
1822
1823  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1824
1825  if (elem.PL != 0.0) {
1826    elem.M->PBpar[HOMmax+n] = bnL/elem.PL;
1827    elem.M->PBpar[HOMmax-n] = anL/elem.PL;
1828  } else {
1829    // thin kick
1830    elem.M->PBpar[HOMmax+n] = bnL; elem.M->PBpar[HOMmax-n] = anL;
1831  }
1832
1833  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1834}
1835
1836
1837void set_dbnL_design_elem(const int Fnum, const int Knum,
1838                          const int n, const double dbnL, const double danL)
1839{
1840  elemtype  elem;
1841
1842  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1843
1844  if (elem.PL != 0.0) {
1845    elem.M->PBpar[HOMmax+n] += dbnL/elem.PL;
1846    elem.M->PBpar[HOMmax-n] += danL/elem.PL;
1847  } else {
1848    // thin kick
1849    elem.M->PBpar[HOMmax+n] += dbnL; elem.M->PBpar[HOMmax-n] += danL;
1850  }
1851
1852  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1853}
1854
1855
1856void set_dbnL_design_fam(const int Fnum,
1857                         const int n, const double dbnL, const double danL)
1858{
1859  int k;
1860
1861  for (k = 1; k <= GetnKid(Fnum); k++)
1862    set_bnL_design_elem(Fnum, k, n, dbnL, danL);
1863}
1864
1865
1866void set_bnL_design_fam(const int Fnum,
1867                        const int n, const double bnL, const double anL)
1868{
1869  int k;
1870
1871  for (k = 1; k <= GetnKid(Fnum); k++)
1872    set_bnL_design_elem(Fnum, k, n, bnL, anL);
1873}
1874
1875
1876void set_bnL_design_type(const int type,
1877                         const int n, const double bnL, const double anL)
1878{
1879  long int  k;
1880
1881  if ((type >= Dip) && (type <= HOMmax)) {
1882    for (k = 1; k <= globval.Cell_nLoc; k++)
1883      if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
1884        set_bnL_design_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL, anL);
1885  } else 
1886    printf("Bad type argument to set_bnL_design_type()\n");
1887}
1888
1889/***********************************************************************
1890void set_bnL_sys_elem(const int Fnum, const int Knum,
1891                      const int n, const double bnL, const double anL)
1892                         
1893   Purpose:
1894       Set field error to the kid of a family, errors are absolute
1895        integrated field strength, replace the previous value!!!
1896       
1897   
1898   Input:
1899      Fnum           family index
1900      Knum           kids index   
1901      n              order of the error
1902      bnL            absolute integrated B component for the n-th error
1903      anL            absolute integrated A component for the n-th error
1904           
1905
1906     
1907   Output:
1908      None
1909     
1910  Return:
1911      None
1912     
1913  Global variables
1914      None
1915   
1916  Specific functions:
1917     None   
1918     
1919 Comments:
1920     None
1921**********************************************************************/
1922void set_bnL_sys_elem(const int Fnum, const int Knum,
1923                      const int n, const double bnL, const double anL)
1924{
1925  elemtype  elem;
1926
1927  const bool  prt = false;
1928
1929  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1930
1931  if (elem.PL != 0.0) {
1932    elem.M->PBsys[HOMmax+n] = bnL/elem.PL;
1933    elem.M->PBsys[HOMmax-n] = anL/elem.PL;
1934  } else {
1935    // thin kick
1936    elem.M->PBsys[HOMmax+n] = bnL; elem.M->PBsys[HOMmax-n] = anL;
1937  }
1938 
1939  Mpole_SetPB(Fnum, Knum, n);    //set for Bn component
1940  Mpole_SetPB(Fnum, Knum, -n);   //set for An component
1941
1942  if (prt)
1943    printf("set the n=%d component of %s to %e %e\n",
1944           n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
1945           bnL, elem.M->PBsys[HOMmax+n]);
1946}
1947
1948/***********************************************************************
1949void set_bnL_sys_fam(const int Fnum,
1950                     const int n, const double bnL, const double anL)
1951                         
1952   Purpose:
1953       Set field error to all the kids of a family,errors are absolute
1954        integrated field strength
1955       
1956       
1957   
1958   Input:
1959      Fnum           family index
1960      n              order of the error
1961      bnL            absolute integrated B component for the n-th error
1962      anL            absolute integrated A component for the n-th error
1963           
1964
1965     
1966   Output:
1967      None
1968     
1969  Return:
1970      None
1971     
1972  Global variables
1973      None
1974   
1975  Specific functions:
1976     None   
1977     
1978 Comments:
1979     None
1980**********************************************************************/
1981void set_bnL_sys_fam(const int Fnum,
1982                     const int n, const double bnL, const double anL)
1983{
1984  int k;
1985
1986  for (k = 1; k <= GetnKid(Fnum); k++)
1987    set_bnL_sys_elem(Fnum, k, n, bnL, anL);
1988}
1989
1990/***********************************************************************
1991void set_bnL_sys_type(const int type,
1992                      const int n, const double bnL, const double anL)
1993                         
1994   Purpose:
1995       Set field error to the type, errors are absolute
1996        integrated field strength
1997   
1998       
1999       
2000   Input:
2001      type          type name
2002      n             order of the error
2003      bnL           absolute integrated B component for the n-th error
2004      anL           absolute integrated A component for the n-th error
2005           
2006
2007     
2008   Output:
2009      None
2010     
2011  Return:
2012      None
2013     
2014  Global variables
2015      None
2016   
2017  Specific functions:
2018     None   
2019     
2020 Comments:
2021      Attension:
2022          This code has a bug, in the for loop, when looking for the
2023          quadrupole family, it will also find the skew quadrupole family,
2024          since their Pkind == Mpole & n_design =2.
2025          It means the multipole error is additionally set N times,
2026          and N is the number of skew quadrupole in the lattice.
2027         
2028                 10/2010   Comments by Jianfeng Zhang
2029     
2030**********************************************************************/
2031void set_bnL_sys_type(const int type,
2032                      const int n, const double bnL, const double anL)
2033{
2034  long int   k;
2035
2036  if (type >= Dip && type <= HOMmax) {
2037    for(k = 1; k <= globval.Cell_nLoc; k++)
2038      if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2039        set_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL, anL);
2040  } else
2041    printf("Bad type argument to set_bnL_sys_type()\n");
2042}
2043
2044
2045void set_bnL_rms_elem(const int Fnum, const int Knum,
2046                      const int n, const double bnL, const double anL,
2047                      const bool new_rnd)
2048{
2049  elemtype  elem;
2050
2051  bool  normal = true, prt = false;
2052 
2053  elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
2054
2055  if (elem.PL != 0.0) {
2056    elem.M->PBrms[HOMmax+n] = bnL/elem.PL;
2057    elem.M->PBrms[HOMmax-n] = anL/elem.PL;
2058  } else {
2059    // thin kick
2060    elem.M->PBrms[HOMmax+n] = bnL; elem.M->PBrms[HOMmax-n] = anL;
2061  }
2062
2063  if(new_rnd){
2064    if (normal) {
2065      elem.M->PBrnd[HOMmax+n] = normranf();
2066      elem.M->PBrnd[HOMmax-n] = normranf();
2067    } else {
2068      elem.M->PBrnd[HOMmax+n] = ranf(); elem.M->PBrnd[HOMmax-n] = ranf();
2069    }
2070  }
2071 
2072  if (prt)
2073    printf("set_bnL_rms_elem:  Fnum = %d, Knum = %d"
2074           ", bnL = %e, anL = %e %e %e\n",
2075           Fnum, Knum, bnL, anL,
2076           elem.M->PBrms[HOMmax+n], elem.M->PBrms[HOMmax-n]);
2077
2078  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
2079}
2080
2081
2082void set_bnL_rms_fam(const int Fnum,
2083                     const int n, const double bnL, const double anL,
2084                     const bool new_rnd)
2085{
2086  int k;
2087
2088  for (k = 1; k <= GetnKid(Fnum); k++)
2089    set_bnL_rms_elem(Fnum, k, n, bnL, anL, new_rnd);
2090}
2091
2092
2093void set_bnL_rms_type(const int type,
2094                      const int n, const double bnL, const double anL,
2095                      const bool new_rnd)
2096{
2097  long int   k;
2098 
2099  if (type >= Dip && type <= HOMmax) {
2100    for(k = 1; k <= globval.Cell_nLoc; k++)
2101      if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2102        set_bnL_rms_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL, anL, new_rnd);
2103  } else
2104    printf("Bad type argument to set_bnL_rms_type()\n");
2105}
2106
2107/***********************************************************************
2108void set_bnr_sys_elem(const int Fnum, const int Knum,
2109                      const int n, const double bnr, const double anr)
2110                         
2111   Purpose:
2112       Set field error to the kid of a family, errors are relative to
2113        design values for (Dip, Quad, Sext, ...)
2114   
2115        If r0 != 0, call this function
2116       
2117   Input:
2118      Fnum           family index
2119      Knum           kids index   
2120      n              order of the error
2121      bnL            relative integrated B component for the n-th error
2122      anL            relative integrated A component for the n-th error
2123           
2124
2125     
2126   Output:
2127      None
2128     
2129  Return:
2130      None
2131     
2132  Global variables
2133      None
2134   
2135  Specific functions:
2136     None   
2137     
2138 Comments:
2139      !!!!!!!!!
2140      This function desn't work for dipole, since   
2141      mp->PBpar[HOMmax+1] is not assigned a value when dipole is
2142      read in t2lat.cc. 
2143     
2144      comments by Jianfeng Zhang  10/2010  @ soleil   
2145**********************************************************************/
2146void set_bnr_sys_elem(const int Fnum, const int Knum,
2147                      const int n, const double bnr, const double anr)
2148{
2149  int        nd;
2150  MpoleType  *mp;
2151  bool prt = false;
2152
2153  mp = Cell[Elem_GetPos(Fnum, Knum)].Elem.M; nd = mp->n_design;
2154  // errors are relative to design values for (Dip, Quad, Sext, ...)
2155  mp->PBsys[HOMmax+n] = bnr*mp->PBpar[HOMmax+nd];
2156  mp->PBsys[HOMmax-n] = anr*mp->PBpar[HOMmax+nd];
2157  //update the total field strength PB and maximum order p_order
2158  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
2159
2160  if (prt)
2161    printf("set the n=%d component of %s to %e %e %e\n",
2162           n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
2163           bnr, mp->PBpar[HOMmax+nd], mp->PBsys[HOMmax+n]);
2164}
2165
2166/***********************************************************************
2167void set_bnr_sys_fam(const int Fnum,
2168                     const int n, const double bnr, const double anr)
2169                         
2170   Purpose:
2171       Set field error to all the kids of a family, errors are relative
2172        to design values for (Dip, Quad, Sext, ...).
2173       
2174        If r0 != 0, call this function
2175   
2176   Input:
2177      Fnum           family index
2178      n              order of the error
2179      bnL            relative integrated B component for the n-th error
2180      anL            relative integrated A component for the n-th error
2181           
2182
2183     
2184   Output:
2185      None
2186     
2187  Return:
2188      None
2189     
2190  Global variables
2191      None
2192   
2193  Specific functions:
2194     None   
2195     
2196 Comments:
2197     None
2198**********************************************************************/
2199void set_bnr_sys_fam(const int Fnum,
2200                     const int n, const double bnr, const double anr)
2201{
2202  int k;
2203
2204  for (k = 1; k <= GetnKid(Fnum); k++)
2205    set_bnr_sys_elem(Fnum, k, n, bnr, anr);
2206}
2207
2208/***********************************************************************
2209void set_bnr_sys_type(const int type,
2210                      const int n, const double bnr, const double anr)
2211                         
2212   Purpose:
2213       Set field error to the type, errors are relative
2214        to design values for (Dip, Quad, Sext, ...)
2215   
2216        If r0 != 0, call this function
2217       
2218   Input:
2219      type          type name
2220      n             order of the error
2221      bnL           relative integrated B component for the n-th error
2222      anL           relative integrated A component for the n-th error
2223           
2224
2225     
2226   Output:
2227      None
2228     
2229  Return:
2230      None
2231     
2232  Global variables
2233      None
2234   
2235  Specific functions:
2236     None   
2237     
2238 Comments:
2239     None
2240**********************************************************************/
2241void set_bnr_sys_type(const int type,
2242                      const int n, const double bnr, const double anr)
2243{
2244  long int   k;
2245
2246  if (type >= Dip && type <= HOMmax) {
2247    for(k = 1; k <= globval.Cell_nLoc; k++)
2248      if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2249        set_bnr_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, bnr, anr);
2250  } else
2251    printf("Bad type argument to set_bnr_sys_type()\n");
2252}
2253
2254/***********************************************************************
2255void set_bnr_rms_elem(const int Fnum, const int Knum,
2256                      const int n, const double bnr, const double anr,
2257                      const bool new_rnd)
2258                         
2259   Purpose:
2260       Update rms field error PBrms and random value of the rms field error PBrnd
2261       to a elment, errors are relative to design values PBpar for (Dip, Quad,
2262       Sext, ...).
2263       
2264        If r0 != 0, call this function
2265   
2266   Input:
2267      Fnum           family index
2268      Knum           kid index
2269      n              order of the error
2270      bnL            relative integrated B component for the n-th error
2271      anL            relative integrated A component for the n-th error
2272      new_rnd        true, generate random value for the rms field error
2273                     false, NOT generate random value for the rms field error       
2274                     
2275   Output:
2276      None
2277     
2278  Return:
2279      None
2280     
2281  Global variables
2282      None
2283   
2284  Specific functions:
2285     None   
2286     
2287 Comments:
2288     None
2289**********************************************************************/
2290void set_bnr_rms_elem(const int Fnum, const int Knum,
2291                      const int n, const double bnr, const double anr,
2292                      const bool new_rnd)
2293{
2294  int        nd;
2295  MpoleType  *mp;
2296
2297  bool  normal = true, prt = false;
2298 
2299  mp = Cell[Elem_GetPos(Fnum, Knum)].Elem.M; nd = mp->n_design;
2300  // errors are relative to design values for (Dip, Quad, Sext, ...)
2301  mp->PBrms[HOMmax+n] = bnr*mp->PBpar[HOMmax+nd];
2302  mp->PBrms[HOMmax-n] = anr*mp->PBpar[HOMmax+nd];
2303  if(new_rnd){
2304    if (normal) {
2305      mp->PBrnd[HOMmax+n] = normranf(); mp->PBrnd[HOMmax-n] = normranf();
2306    } else {
2307      mp->PBrnd[HOMmax+n] = ranf(); mp->PBrnd[HOMmax-n] = ranf();
2308    }
2309  }
2310 
2311  if (prt)
2312    printf("set_bnr_rms_elem:  Fnum = %d, Knum = %d"
2313           ", bnr = %e, anr = %e %e %e\n",
2314           Fnum, Knum, bnr, anr, mp->PBrms[HOMmax+n], mp->PBrms[HOMmax-n]);
2315
2316  Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
2317}
2318
2319/***********************************************************************
2320void set_bnr_rms_fam(const int Fnum,
2321                     const int n, const double bnr, const double anr,
2322                     const bool new_rnd)
2323                         
2324   Purpose:
2325       Update rms field error PBrms and random value of the rms field error PBrnd
2326       to all the kid elments in a family, errors are relative to design values PBpar
2327       for (Dip, Quad, Sext, ...).
2328       
2329        If r0 != 0, call this function
2330   
2331   Input:
2332      Fnum           family index
2333      n              order of the error
2334      bnr            relative integrated B component for the n-th error
2335      anr            relative integrated A component for the n-th error
2336      new_rnd        true, generate random value for the rms field error
2337                     false, NOT generate random value for the rms field error       
2338                     
2339   Output:
2340      None
2341     
2342  Return:
2343      None
2344     
2345  Global variables
2346      None
2347   
2348  Specific functions:
2349     None   
2350     
2351 Comments:
2352     None
2353**********************************************************************/
2354void set_bnr_rms_fam(const int Fnum,
2355                     const int n, const double bnr, const double anr,
2356                     const bool new_rnd)
2357{
2358  int k;
2359
2360  for (k = 1; k <= GetnKid(Fnum); k++)
2361    set_bnr_rms_elem(Fnum, k, n, bnr, anr, new_rnd);
2362}
2363
2364/***********************************************************************
2365void set_bnr_rms_type(const int type,
2366                      const int n, const double bnr, const double anr,
2367                      const bool new_rnd)
2368                         
2369   Purpose:
2370       Update rms field error PBrms and random value of the rms field error PBrnd
2371       to all the kid elments in a family, errors are relative to design values PBpar
2372       for (Dip, Quad, Sext, ...).
2373       
2374        If r0 != 0, call this function
2375   
2376   Input:
2377      type           type of the elment,maybe 'dipole', 'quadrupole', 'sextupole', 'all'
2378      n              order of the error
2379      bnr            relative integrated B component for the n-th error
2380      anr            relative integrated A component for the n-th error
2381      new_rnd        true, generate random value for the rms field error
2382                     false, NOT generate random value for the rms field error       
2383                     
2384   Output:
2385      None
2386     
2387  Return:
2388      None
2389     
2390  Global variables
2391      None
2392   
2393  Specific functions:
2394     None   
2395     
2396 Comments:
2397     None
2398**********************************************************************/
2399void set_bnr_rms_type(const int type,
2400                      const int n, const double bnr, const double anr,
2401                      const bool new_rnd)
2402{
2403  long int   k;
2404 
2405  if (type >= Dip && type <= HOMmax) {
2406    for(k = 1; k <= globval.Cell_nLoc; k++)
2407      if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2408        set_bnr_rms_elem(Cell[k].Fnum, Cell[k].Knum, n, bnr, anr, new_rnd);
2409  } else
2410    printf("Bad type argument to set_bnr_rms_type()\n");
2411}
2412
2413
2414double get_Wiggler_BoBrho(const int Fnum, const int Knum)
2415{
2416  return Cell[Elem_GetPos(Fnum, Knum)].Elem.W->BoBrhoV[0];
2417}
2418
2419/****************************************************************************/
2420/* void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
2421
2422   Purpose:
2423       Set the integrated field strength for the specific element in the family
2424       of wiggers.       
2425       
2426   Input:
2427       Fnum     family name
2428       BoBrhoV  integrated field strength
2429
2430   Output:
2431      none
2432
2433   Return:
2434       none
2435
2436   Global variables:
2437       none
2438
2439   Specific functions:
2440       none
2441
2442   Comments:
2443       Called by set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
2444       
2445****************************************************************************/
2446void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
2447{
2448  Cell[Elem_GetPos(Fnum, Knum)].Elem.W->BoBrhoV[0] = BoBrhoV;
2449  Cell[Elem_GetPos(Fnum, Knum)].Elem.W->PBW[HOMmax+Quad] = -sqr(BoBrhoV)/2.0;
2450  Wiggler_SetPB(Fnum, Knum, Quad);
2451}
2452
2453/****************************************************************************/
2454/* void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
2455
2456   Purpose:
2457       Set the integrated field strength for the family of wiggers       
2458
2459   Input:
2460       Fnum     family name
2461       BoBrhoV  integrated field strength
2462
2463   Output:
2464      none
2465
2466   Return:
2467       none
2468
2469   Global variables:
2470       none
2471
2472   Specific functions:
2473       none
2474
2475   Comments:
2476       none
2477       
2478****************************************************************************/
2479void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
2480{
2481  int  k;
2482
2483  for (k = 1; k <= GetnKid(Fnum); k++)
2484    set_Wiggler_BoBrho(Fnum, k, BoBrhoV);
2485}
2486
2487/****************************************************************************/
2488/* void set_ID_scl(const int Fnum, const int Knum, const double scl)
2489
2490   Purpose:
2491       Set the scale factor for the spicified element of wigger, insertion devices,
2492       or fieldMap.
2493       Called by set_ID_scl(const int Fnum, const double scl).       
2494
2495   Input:
2496       Fnum  family name
2497       Knum  kid index of Fname
2498       scl   scaling factor
2499
2500   Output:
2501      none
2502
2503   Return:
2504       none
2505
2506   Global variables:
2507       none
2508
2509   Specific functions:
2510       none
2511
2512   Comments:
2513       none
2514       
2515****************************************************************************/
2516void set_ID_scl(const int Fnum, const int Knum, const double scl)
2517{
2518  int           k;
2519  WigglerType*  W;
2520
2521  switch (Cell[Elem_GetPos(Fnum, Knum)].Elem.Pkind) {
2522  case Wigl:
2523    // scale the ID field
2524    W = Cell[Elem_GetPos(Fnum, Knum)].Elem.W;
2525    for (k = 0; k < W->n_harm; k++) {
2526      W->BoBrhoH[k] = scl*ElemFam[Fnum-1].ElemF.W->BoBrhoH[k];
2527      W->BoBrhoV[k] = scl*ElemFam[Fnum-1].ElemF.W->BoBrhoV[k];
2528    }
2529    break;
2530  case Insertion:
2531    Cell[Elem_GetPos(Fnum, Knum)].Elem.ID->scaling1 = scl;
2532    Cell[Elem_GetPos(Fnum, Knum)].Elem.ID->scaling2 = scl;
2533    break;
2534  case FieldMap:
2535    Cell[Elem_GetPos(Fnum, Knum)].Elem.FM->scl = scl;
2536    break;
2537  default:
2538    cout << "set_ID_scl: unknown element type" << endl;
2539    exit_(1);
2540    break;
2541  }
2542}
2543
2544/****************************************************************************/
2545/* void set_ID_scl(const int Fnum, const double scl)
2546
2547   Purpose:
2548       Set the scale factor for the specified wigger, insertion devices, or fieldMap.       
2549
2550   Input:
2551       Fnum  family name
2552       scl   scaling factor
2553
2554   Output:
2555      none
2556
2557   Return:
2558       none
2559
2560   Global variables:
2561       none
2562
2563   Specific functions:
2564       none
2565
2566   Comments:
2567       none
2568       
2569****************************************************************************/
2570void set_ID_scl(const int Fnum, const double scl)
2571{
2572  int  k;
2573
2574  for (k = 1; k <= GetnKid(Fnum); k++)
2575    set_ID_scl(Fnum, k, scl);
2576}
2577
2578/***********************************************************************
2579void SetFieldValues_fam(const int Fnum, const bool rms, const double r0,
2580                        const int n, const double Bn, const double An,
2581                        const bool new_rnd)
2582                         
2583   Purpose:
2584       Set field error to all the kids in a family
2585   
2586   Input:
2587      Fnum            family name
2588      rms             flag to set rms field error
2589      r0              radius at which error is measured
2590      n               order of the error
2591      Bn              integrated B component for the n-th error
2592      An              integrated A component for the n-th error
2593      new_rnd         bool flag to set new random number           
2594
2595     
2596   Output:
2597      None
2598     
2599  Return:
2600      None
2601     
2602  Global variables
2603      None
2604   
2605  Specific functions:
2606     None   
2607     
2608 Comments:
2609     None
2610**********************************************************************/
2611void SetFieldValues_fam(const int Fnum, const bool rms, const double r0,
2612                        const int n, const double Bn, const double An,
2613                        const bool new_rnd)
2614{
2615  int     N;
2616  double  bnr, anr;
2617
2618  N = Cell[Elem_GetPos(Fnum, 1)].Elem.M->n_design;
2619  if (r0 == 0.0) {
2620    // input is: (b_n L), (a_n L)
2621    if(rms)
2622      set_bnL_rms_fam(Fnum, n, Bn, An, new_rnd);
2623    else
2624      set_bnL_sys_fam(Fnum, n, Bn, An);
2625  } else {
2626    bnr = Bn/pow(r0, n-N); anr = An/pow(r0, n-N);
2627    if(rms)
2628      set_bnr_rms_fam(Fnum, n, bnr, anr, new_rnd);
2629    else
2630      set_bnr_sys_fam(Fnum, n, bnr, anr);
2631  }
2632}
2633
2634
2635/***********************************************************************
2636void SetFieldValues_type(const int N, const bool rms, const double r0,
2637                         const int n, const double Bn, const double An,
2638                         const bool new_rnd)
2639                         
2640   Purpose:
2641       Set field error to the type
2642   
2643   Input:
2644      N            type name
2645      rms          flag to set rms field error
2646      r0           radius at which error is measured
2647      n            order of the error
2648      Bn           integrated B component for the n-th error
2649      An           integrated A component for the n-th error
2650      new_rnd      bool flag to set new random number       
2651
2652     
2653   Output:
2654      None
2655     
2656  Return:
2657      None
2658     
2659  Global variables
2660      None
2661   
2662  Specific functions:
2663     None   
2664     
2665 Comments:
2666     None
2667**********************************************************************/
2668void SetFieldValues_type(const int N, const bool rms, const double r0,
2669                         const int n, const double Bn, const double An,
2670                         const bool new_rnd)
2671{
2672  double  bnr, anr;
2673
2674  if (r0 == 0.0) {
2675    // input is: (b_n L), (a_n L)
2676    if(rms)
2677      set_bnL_rms_type(N, n, Bn, An, new_rnd);
2678    else
2679      set_bnL_sys_type(N, n, Bn, An);
2680  } else {
2681    bnr = Bn/pow(r0, n-N); anr = An/pow(r0, n-N);
2682    if(rms)
2683      set_bnr_rms_type(N, n, bnr, anr, new_rnd);
2684    else
2685      set_bnr_sys_type(N, n, bnr, anr);
2686  }
2687}
2688
2689
2690/***********************************************************************
2691void SetFieldErrors(const char *name, const bool rms, const double r0,
2692                    const int n, const double Bn, const double An,
2693                    const bool new_rnd)
2694                   
2695   Purpose:
2696       Set field error to the type or element
2697   
2698   Input:
2699      name         type name or element name
2700      rms          flag to set rms field error
2701      r0           radius at which error is measured
2702      n            order of the error
2703      Bn           integrated B component for the n-th error
2704      An           integrated A component for the n-th error
2705      new_rnd      flag to set new random number for the error     
2706
2707     
2708   Output:
2709      None
2710     
2711  Return:
2712      None
2713     
2714  Global variables
2715      None
2716   
2717  Specific functions:
2718     None   
2719     
2720 Comments:
2721       10/2010   Comments by Jianfeng Zhang
2722**********************************************************************/
2723void SetFieldErrors(const char *name, const bool rms, const double r0,
2724                    const int n, const double Bn, const double An,
2725                    const bool new_rnd) 
2726{
2727  int     Fnum;
2728
2729  if (strcmp("all", name) == 0) {
2730    printf("all: not yet implemented\n");
2731  } else if (strcmp("dip", name) == 0) {
2732    SetFieldValues_type(Dip, rms, r0, n, Bn, An, new_rnd);
2733  } else if (strcmp("quad", name) == 0) {
2734    SetFieldValues_type(Quad, rms, r0, n, Bn, An, new_rnd);
2735  } else if (strcmp("sext", name) == 0) {
2736    SetFieldValues_type(Sext, rms, r0, n, Bn, An, new_rnd);
2737  } else {
2738    Fnum = ElemIndex(name);
2739    if(Fnum > 0)
2740      SetFieldValues_fam(Fnum, rms, r0, n, Bn, An, new_rnd);
2741    else 
2742      printf("SetFieldErrors: undefined element %s\n", name);
2743  }
2744}
2745
2746
2747
2748/********************************************************************
2749char* get_prm(void)
2750
2751    Purpose:
2752        Find the stoke talbe symbol "\t" in the string, and
2753        then saved and return the position of the pointer.
2754     
2755    symbol for end of the line:
2756        \n   
2757    Comments:
2758       10/2010  By Jianfeng Zhang
2759        Fix the bug for reading the end of line symbol "\n" , "\r",'\r\n'
2760        at different operation system
2761*********************************************************************/
2762char* get_prm(void)
2763{
2764  char  *prm;
2765
2766  prm = strtok(NULL, " \t");
2767
2768  //if ((prm == NULL) || (strcmp(prm, "\r\n") == 0)) {
2769  if ((prm == NULL) || (strcmp(prm, "\n") == 0)|| (strcmp(prm, "\r") == 0)|| (strcmp(prm, "\r\n") == 0)) {
2770    printf("get_prm: incorrect format\n");
2771    exit_(1);
2772  }
2773
2774  return prm;
2775}
2776
2777
2778void LoadFieldErr(const char *FieldErrorFile, const bool Scale_it,
2779                  const double Scale, const bool new_rnd) 
2780{ 
2781  bool    rms, set_rnd;
2782  char    line[max_str], name[max_str], type[max_str], *prm;
2783  int     k, n, seed_val;
2784  double  Bn, An, r0;
2785  FILE    *inf;
2786
2787  const bool  prt = true;
2788
2789  inf = file_read(FieldErrorFile);
2790
2791  set_rnd = false; 
2792  printf("\n");
2793  while (fgets(line, max_str, inf) != NULL) {
2794    if (strstr(line, "#") == NULL) {
2795      // check for whether to set new seed
2796      sscanf(line, "%s", name); 
2797      if (strcmp("seed", name) == 0) {
2798        set_rnd = true;
2799        sscanf(line, "%*s %d", &seed_val); 
2800//      printf("setting random seed to %d\n", seed_val);
2801//      iniranf(seed_val);
2802        printf("\n");
2803        printf("setting new random numbers (%d)\n", seed_val);
2804      } else {
2805        sscanf(line, "%*s %s %lf", type, &r0);
2806        printf("%-4s %3s %7.1le", name, type, r0);
2807        rms = (strcmp("rms", type) == 0)? true : false;
2808        if (rms && !set_rnd) {
2809          printf("LoadFieldErr: seed not defined\n");
2810          exit_(1);
2811        }
2812        // skip first three parameters
2813        strtok(line, " \t");
2814        for (k = 1; k <= 2; k++)
2815          strtok(NULL, " \t");
2816        while (((prm = strtok(NULL, " \t")) != NULL) &&
2817               (strcmp(prm, "\r\n") != 0)) {
2818          sscanf(prm, "%d", &n);
2819          prm = get_prm();
2820          sscanf(prm, "%lf", &Bn);
2821          prm = get_prm(); 
2822          sscanf(prm, "%lf", &An);
2823          if (Scale_it)
2824            {Bn *= Scale; An *= Scale;}
2825          if (prt)
2826            printf(" %2d %9.1e %9.1e\n", n, Bn, An);
2827          // convert to normalized multipole components
2828          SetFieldErrors(name, rms, r0, n, Bn, An, true);
2829        }
2830      }
2831    } else
2832      printf("%s", line);
2833  }
2834
2835  fclose(inf);
2836}
2837
2838
2839// correction algorithms
2840
2841// control of orbit
2842
2843// closed orbit correction by n_orbit iterations
2844bool CorrectCOD(int n_orbit)
2845{
2846  bool      cod;
2847  int       i;
2848  long int  lastpos;
2849  Vector2   mean, sigma, max;
2850
2851  cod = getcod(0.0, lastpos);
2852  if (cod) {
2853    codstat(mean, sigma, max, globval.Cell_nLoc, false); // get orbit stats
2854    printf("\n");
2855    printf("Initial RMS orbit (BPMs):   x = %7.1e mm, y = %7.1e mm\n",
2856           1e3*sigma[X_], 1e3*sigma[Y_]);
2857    codstat(mean, sigma, max, globval.Cell_nLoc, true);
2858    printf("Initial RMS orbit (all):    x = %7.1e mm, y = %7.1e mm\n",
2859           1e3*sigma[X_], 1e3*sigma[Y_]);
2860
2861    for (i = 0; i < n_orbit; i++){
2862      lsoc(1, globval.bpm, globval.hcorr, 1); // correct horizontal orbit
2863      lsoc(1, globval.bpm, globval.vcorr, 2); // correct vertical orbit
2864      cod = getcod(0.0, lastpos);             // find closed orbit
2865      if (!cod) break;
2866    }
2867
2868    if (cod) {
2869      codstat(mean, sigma, max, globval.Cell_nLoc, false);
2870      printf("Corrected RMS orbit (BPMs): x = %7.1e mm, y = %7.1e mm\n",
2871             1e3*sigma[X_], 1e3*sigma[Y_]);
2872      codstat(mean, sigma, max, globval.Cell_nLoc, true);
2873      printf("Corrected RMS orbit (all):  x = %7.1e mm, y = %7.1e mm\n",
2874             1e3*sigma[X_], 1e3*sigma[Y_]);
2875    }
2876  }
2877
2878  return cod;
2879}
2880
2881
2882void Align_BPMs(const int n)
2883{
2884  // Align BPMs to adjacent multipoles.
2885
2886  bool      aligned;
2887  int       i, j, k;
2888  long int  loc;
2889
2890  const int  n_step = 5;
2891
2892  printf("\n");
2893  for (i = 1; i <= GetnKid(globval.bpm); i++) {
2894    loc = Elem_GetPos(globval.bpm, i);
2895
2896    if ((loc == 1) || (loc == globval.Cell_nLoc)) {
2897      printf("Align_BPMs: BPM at entrance or exit of lattice: %ld\n", loc);
2898      exit_(1);
2899    }
2900
2901    j = 1; aligned = false;
2902    do {
2903      if ((Cell[loc-j].Elem.Pkind == Mpole) &&
2904          (Cell[loc-j].Elem.M->n_design == n)) {
2905        for (k = 0; k <= 1; k++)
2906          Cell[loc].Elem.M->PdSsys[k] = Cell[loc-j].dS[k];
2907        printf("aligned BPM no %1d to %s\n", i, Cell[loc-j].Elem.PName);
2908        aligned = true; break;
2909      } else if ((Cell[loc+j].Elem.Pkind == Mpole) &&
2910                 (Cell[loc+j].Elem.M->n_design == n)) {
2911        for (k = 0; k <= 1; k++)
2912          Cell[loc].Elem.M->PdSsys[k] = Cell[loc+j].dS[k];
2913        printf("aligned BPM no %1d to %s\n", i, Cell[loc+j].Elem.PName);
2914        aligned = true; break;
2915      }
2916
2917      j++;
2918    } while (j <= n_step);
2919     
2920    if (aligned)
2921      Mpole_SetdS(globval.bpm, i);
2922    else
2923      printf("Align_BPMs: no multipole adjacent to BPM no %d\n", i);
2924  }
2925}
2926
2927/* *****************************************************
2928 void get_bare()
2929 
2930  Purpose:
2931           store values of the optics function at the sextupoles
2932*/
2933void get_bare()
2934{
2935  long int  j, k;
2936
2937  n_sext = 0;
2938  for (j = 0; j <= globval.Cell_nLoc; j++) {
2939    if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design >= Sext)) {
2940      n_sext++; sexts[n_sext-1] = j;
2941      for (k = 0; k <= 1; k++) {
2942        beta0_[n_sext-1][k] = Cell[j].Beta[k];
2943        nu0s[n_sext-1][k] = Cell[j].Nu[k];
2944      }
2945    }
2946  }
2947
2948  // beta-functions for normalization
2949  beta_ref[X_] = Cell[globval.Cell_nLoc].Beta[X_];
2950  beta_ref[Y_] = Cell[globval.Cell_nLoc].Beta[Y_];
2951
2952  nu0_[X_] = globval.TotalTune[X_]; nu0_[Y_] = globval.TotalTune[Y_];
2953}
2954
2955
2956void get_dbeta_dnu(double m_dbeta[], double s_dbeta[],
2957                   double m_dnu[], double s_dnu[])
2958{
2959  int       k;
2960  long int  j, ind;
2961  double    dbeta, dnu;
2962
2963  Ring_GetTwiss(false, 0.0);
2964 
2965  for (k = 0; k <= 1; k++) {
2966    m_dbeta[k] = 0.0; s_dbeta[k] = 0.0; m_dnu[k] = 0.0; s_dnu[k] = 0.0;
2967  }
2968 
2969  for (j = 0; j < n_sext; j++) {
2970    ind = sexts[j];
2971    for (k = 0; k <= 1; k++) {
2972      dbeta = (Cell[ind].Beta[k]-beta0_[j][k])/beta0_[j][k];
2973      m_dbeta[k] += dbeta; s_dbeta[k] += sqr(dbeta);
2974      dnu = Cell[ind].Nu[k] - nu0s[j][k];
2975      m_dnu[k] += dnu; s_dnu[k] += sqr(dnu);
2976    }
2977  }
2978
2979  for (k = 0; k <= 1; k++) {
2980    m_dbeta[k] /= n_sext; m_dnu[k] /= n_sext;
2981    s_dbeta[k] = sqrt((s_dbeta[k]-n_sext*sqr(m_dbeta[k]))/(n_sext-1));
2982    s_dnu[k] = sqrt((s_dnu[k]-n_sext*sqr(m_dnu[k]))/(n_sext-1));
2983  }
2984}
2985
2986
2987bool CorrectCOD_N(const char *ae_file, const int n_orbit,
2988                  const int n, const int k)
2989{
2990  bool    cod = false;
2991  int     i;
2992  double  m_dbeta[2], s_dbeta[2], m_dnu[2], s_dnu[2];
2993
2994  // Clear trim setpoints
2995  set_bnL_design_fam(globval.hcorr, Dip, 0.0, 0.0);
2996  set_bnL_design_fam(globval.vcorr, Dip, 0.0, 0.0);
2997
2998  // load misalignments
2999  LoadAlignTol(ae_file, true, 1.0, true, k);
3000  for (i = 1; i <= n; i++) {
3001    // Scale the rms values
3002    LoadAlignTol(ae_file, true, (double)i/(double)n, false, k);
3003   
3004    if (bba) {
3005      // Beam based alignment
3006      Align_BPMs(Quad);
3007    }
3008
3009    cod = CorrectCOD(n_orbit); 
3010   
3011    if (!cod) break;
3012   
3013    get_dbeta_dnu(m_dbeta, s_dbeta, m_dnu, s_dnu);
3014    printf("\n");
3015    printf("RMS dbeta_x/beta_x = %4.2f%%,   dbeta_y/beta_y = %4.2f%%\n",
3016           1e2*s_dbeta[X_], 1e2*s_dbeta[Y_]);
3017    printf("RMS dnu_x          = %7.5f, dnu_y          = %7.5f\n",
3018           s_dnu[X_], s_dnu[Y_]);
3019  }
3020
3021  return cod;
3022}
3023
3024
3025// Control of vertical beam size
3026
3027void FindSQ_SVDmat(double **SkewRespMat, double **U, 
3028                   double **V, double *w, int N_COUPLE, int N_SKEW)
3029{
3030  int i, j;
3031
3032  const double  cut = 1e-10; // cut value for SVD singular values
3033
3034  for (i = 1; i <= N_COUPLE; i++)
3035    for (j = 1; j <= N_SKEW; j++)
3036      U[i][j] = SkewRespMat[i][j];
3037
3038  // prepare matrices for SVD
3039  dsvdcmp(U, N_COUPLE, N_SKEW, w, V);
3040
3041  // zero singular values
3042  printf("\n");
3043  printf("singular values:\n");
3044  printf("\n");
3045  for (i = 1; i <= N_SKEW; i++) {
3046    printf("%11.3e", w[i]);
3047    if (w[i] < cut ) {
3048      w[i] = 0.0;
3049      printf(" (zeroed)");
3050      if (i % 8 == 0) printf("\n");
3051    }
3052  }
3053  if (i % 8 != 0) printf("\n");
3054}
3055
3056
3057void FindMatrix(double **SkewRespMat, const double deta_y_max)
3058{
3059  //  Ring_GetTwiss(true, 0.0) should be called in advance
3060  int       i, j, k;
3061  long int  loc;
3062  double    nuX, nuY, alpha, eta_y_max;
3063  double    *etaSQ;
3064  double    **betaSQ, **nuSQ, **betaBPM, **nuBPM;
3065  double    **betaHC, **nuHC, **betaVC, **nuVC;
3066  FILE      *SkewMatFile, *fp;
3067
3068  const int     Xi = 1, Yi = 2;
3069  const double  pi = M_PI, twopi = 2.0*M_PI;
3070
3071
3072  etaSQ = dvector(1, N_SKEW); betaSQ = dmatrix(1, N_SKEW, 1, 2);
3073  nuSQ = dmatrix(1, N_SKEW, 1, 2);
3074  betaBPM = dmatrix(1, N_BPM, 1, 2); nuBPM = dmatrix(1, N_BPM, 1, 2);
3075  betaHC = dmatrix(1, N_HCOR, 1, 2); nuHC = dmatrix(1, N_HCOR, 1, 2);
3076  betaVC = dmatrix(1, N_VCOR, 1, 2); nuVC = dmatrix(1, N_VCOR, 1, 2);
3077
3078  nuX = globval.TotalTune[X_]; nuY = globval.TotalTune[Y_];
3079
3080  for (i = 1; i <= N_SKEW; i++) {
3081    loc = Elem_GetPos(globval.qt, i);
3082    etaSQ[i] = Cell[loc].Eta[X_];
3083    betaSQ[i][Xi] = Cell[loc].Beta[X_]; betaSQ[i][Yi] = Cell[loc].Beta[Y_];
3084    nuSQ[i][Xi] = Cell[loc].Nu[X_]; nuSQ[i][Yi] = Cell[loc].Nu[Y_];
3085  } // for i=1..N_SKEW
3086
3087  for (i = 1; i <= N_BPM; i++) {
3088    betaBPM[i][Xi] = Cell[bpm_loc[i-1]].Beta[X_];
3089    betaBPM[i][Yi] = Cell[bpm_loc[i-1]].Beta[Y_];
3090    nuBPM[i][Xi] = Cell[bpm_loc[i-1]].Nu[X_];
3091    nuBPM[i][Yi] = Cell[bpm_loc[i-1]].Nu[Y_];
3092  } // for i=1..N_BPM
3093
3094  for (i = 1; i <= N_HCOR; i++) {
3095    betaHC[i][Xi] = Cell[h_corr[i-1]].Beta[X_];
3096    betaHC[i][Yi] = Cell[h_corr[i-1]].Beta[Y_];
3097    nuHC[i][Xi] = Cell[h_corr[i-1]].Nu[X_];
3098    nuHC[i][Yi] = Cell[h_corr[i-1]].Nu[Y_];
3099  } // for i=1..N_HCOR
3100
3101  for (i = 1; i <= N_VCOR; i++) {
3102    betaVC[i][Xi] = Cell[v_corr[i-1]].Beta[X_];
3103    betaVC[i][Yi] = Cell[v_corr[i-1]].Beta[Y_];
3104    nuVC[i][Xi] = Cell[v_corr[i-1]].Nu[X_];
3105    nuVC[i][Yi] = Cell[v_corr[i-1]].Nu[Y_];
3106  } // for i=1..N_VCOR
3107
3108
3109  for (i = 1; i <= N_SKEW; i++) {
3110    // looking for term for vertical dispersion
3111    alpha = etaSQ[i];
3112    // printf("For skew quad %3d kick is %9.2e\n",i,alpha);
3113    for (j = 1; j <= N_BPM; j++) {
3114      SkewRespMat[j][i] = VDweight*0.5*alpha*sqrt(betaSQ[i][Yi]*betaBPM[j][Yi])
3115        *cos(twopi*fabs(nuSQ[i][Yi]-nuBPM[j][Yi])-pi*nuY)/sin(pi*nuY);
3116    } // for (j=1; j<=N_BPM; j++)
3117
3118    // loking for coupling of horizontal trim to vertical BPM
3119    for (k = 1; k <= N_HCOR; k++) {
3120      // find v-kick by i-th skew quad due to the k-th h-trim
3121      alpha = 0.5*sqrt(betaSQ[i][Xi]*betaHC[k][Xi])*
3122        cos(twopi*fabs(nuSQ[i][Xi]-nuHC[k][Xi])-pi*nuX)/sin(pi*nuX);
3123      // find vertical orbit due to the kick
3124      for (j = 1; j <= N_BPM; j++) 
3125        SkewRespMat[N_BPM+(k-1)*N_HCOR+j][i] = 
3126          HVweight*0.5*alpha*sqrt(betaSQ[i][Yi]*betaBPM[j][Yi])*
3127          cos(twopi*fabs(nuSQ[i][Yi]-nuBPM[j][Yi])-pi*nuY)/sin(pi*nuY);
3128    } //for (k=1; k<=N_HCOR; k++)
3129
3130   //loking for coupling of vertical trim to horizontal BPM
3131    for (k = 1; k <= N_VCOR; k++) {
3132      // find h-kick by i-th skew quad due to the k-th v-trim
3133      alpha = 0.5*sqrt(betaSQ[i][Yi]*betaVC[k][Yi])*
3134        cos(twopi*fabs(nuSQ[i][Yi]-nuVC[k][Yi])-pi*nuY)/sin(pi*nuY);
3135      // find horizontal orbit due to the kick
3136      for (j = 1; j <= N_BPM; j++) 
3137        SkewRespMat[N_BPM+N_BPM*N_HCOR+(k-1)*N_VCOR+j][i] = 
3138          VHweight*0.5*alpha*sqrt(betaSQ[i][Xi]*betaBPM[j][Xi])*
3139          cos(twopi*fabs(nuSQ[i][Xi]-nuBPM[j][Xi])-pi*nuX)/sin(pi*nuX);
3140    } //for (k=1; k<=N_VCOR; k++)
3141  } // for i=1..N_SKEW
3142
3143
3144  SkewMatFile = file_write(SkewMatFileName);
3145  for (i = 1; i <= N_SKEW; i++) {
3146    for (j = 1; j <= N_COUPLE; j++) 
3147      fprintf(SkewMatFile, "%9.2e ", SkewRespMat[j][i]);
3148    fprintf(SkewMatFile, "\n");
3149  }
3150  fclose(SkewMatFile);
3151
3152  fp = file_write(deta_y_FileName);
3153  eta_y_max = 0.0;
3154  for (j = 1; j <= N_BPM; j++) {
3155    eta_y[j] = 0.0;
3156    for (i = 1; i <= N_SKEW; i++)
3157      if (i % SQ_per_scell == 0) {
3158        eta_y[j] += 0.5*etaSQ[i]*sqrt(betaSQ[i][Yi]*betaBPM[j][Yi])
3159            *cos(twopi*fabs(nuSQ[i][Yi]-nuBPM[j][Yi])-pi*nuY)/sin(pi*nuY);
3160      }
3161    eta_y_max = max(fabs(eta_y[j]), eta_y_max);
3162  }
3163  for (j = 1; j <= N_BPM; j++)
3164    eta_y[j] /= eta_y_max;
3165
3166  for (j = 1; j <= N_BPM; j++) {
3167    eta_y[j] *= deta_y_max;
3168    fprintf(fp, "%6.3f %10.3e\n", Cell[bpm_loc[j-1]].S, 1e3*eta_y[j]);
3169  }
3170  fclose(fp);
3171
3172  free_dvector(etaSQ, 1, N_SKEW); free_dmatrix(betaSQ, 1, N_SKEW, 1, 2);
3173  free_dmatrix(nuSQ, 1, N_SKEW, 1, 2);
3174  free_dmatrix(betaBPM, 1, N_BPM, 1, 2); free_dmatrix(nuBPM, 1, N_BPM, 1, 2);
3175  free_dmatrix(betaHC, 1, N_HCOR, 1, 2); free_dmatrix(nuHC, 1, N_HCOR, 1, 2);
3176  free_dmatrix(betaVC, 1, N_VCOR, 1, 2); free_dmatrix(nuVC, 1, N_VCOR, 1, 2);
3177} // FindMatrix
3178
3179
3180void ini_skew_cor(const double deta_y_max)
3181{
3182  int  k;
3183
3184  // No of skew quads, BPMs, and correctors
3185  N_SKEW = GetnKid(globval.qt);
3186
3187  N_BPM = 0;
3188  for (k = 1; k <= GetnKid(globval.bpm); k++) {
3189    N_BPM++;
3190
3191    if (N_BPM > max_bpm) {
3192      printf("ini_skew_cor: max no of BPMs exceeded %d (%d)\n",
3193             N_BPM, max_bpm);
3194      exit_(1);
3195    }
3196
3197    bpm_loc[N_BPM-1] = Elem_GetPos(globval.bpm, k);
3198  }
3199
3200  N_HCOR = 0;
3201  h_corr[N_HCOR++] = Elem_GetPos(globval.hcorr, 1);
3202  h_corr[N_HCOR++] = Elem_GetPos(globval.hcorr, GetnKid(globval.hcorr)/3);
3203  h_corr[N_HCOR++] = Elem_GetPos(globval.hcorr, 2*GetnKid(globval.hcorr)/3);
3204
3205  N_VCOR = 0;
3206  v_corr[N_VCOR++] = Elem_GetPos(globval.vcorr, 1);
3207  v_corr[N_VCOR++] = Elem_GetPos(globval.vcorr, GetnKid(globval.vcorr)/3);
3208  v_corr[N_VCOR++] = Elem_GetPos(globval.vcorr, 2*GetnKid(globval.vcorr)/3);
3209
3210  N_COUPLE = N_BPM*(1+N_HCOR+N_VCOR);
3211
3212  SkewRespMat = dmatrix(1, N_COUPLE, 1, N_SKEW);
3213  VertCouple = dvector(1, N_COUPLE);
3214  SkewStrengthCorr = dvector(1, N_SKEW);
3215  b = dvector(1, N_COUPLE); w = dvector(1, N_SKEW);
3216  V = dmatrix(1, N_SKEW, 1, N_SKEW); U = dmatrix(1, N_COUPLE, 1, N_SKEW);
3217  eta_y = dvector(1, N_BPM);
3218
3219  printf("\n");
3220  printf("Number of trims:                   horizontal = %d, vertical = %d\n",
3221         N_HCOR, N_VCOR);
3222  printf("Number of BPMs:                    %6d\n", N_BPM);
3223  printf("Number of skew quads:              %6d\n", N_SKEW);
3224  printf("Number of elements in skew vector: %6d\n", N_COUPLE);
3225
3226  // find matrix
3227  Ring_GetTwiss(true, 0.0);
3228
3229  printf("\n");
3230  printf("Looking for response matrix\n");   
3231  FindMatrix(SkewRespMat, deta_y_max);
3232
3233  printf("Looking for SVD matrices\n"); 
3234  FindSQ_SVDmat(SkewRespMat, U, V, w, N_COUPLE, N_SKEW);
3235}
3236
3237
3238void FindCoupVector(double *VertCouple)
3239{
3240  bool      cod;
3241  long      i,j;
3242  long      lastpos;
3243  double    *orbitP, *orbitN;
3244
3245  orbitP = dvector(1, N_BPM); orbitN = dvector(1, N_BPM);
3246
3247  // Find vertical dispersion
3248  Cell_Geteta(0, globval.Cell_nLoc, true, 0e0);
3249
3250  for (i = 1; i <= N_BPM; i++)
3251    VertCouple[i] = VDweight*Cell[bpm_loc[i-1]].Eta[Y_];
3252  // Finished finding vertical dispersion
3253
3254  // Find off diagonal terms for horizontal trims
3255  for (j = 1; j <= N_HCOR; j++) {
3256    // positive kick: "+Dip" for horizontal
3257    SetdKLpar(Cell[h_corr[j-1]].Fnum, Cell[h_corr[j-1]].Knum, +Dip, kick);
3258    cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3259    for (i = 1; i <= N_BPM; i++)
3260      orbitP[i] = Cell[bpm_loc[i-1]].BeamPos[y_];
3261
3262    //negative kick: "+Dip" for horizontal
3263    SetdKLpar(Cell[h_corr[j-1]].Fnum, Cell[h_corr[j-1]].Knum, +Dip, -2*kick);
3264    cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3265    for (i = 1; i <= N_BPM; i++)
3266      orbitN[i] = Cell[bpm_loc[i-1]].BeamPos[y_];
3267
3268    // restore trim valueL: "+Dip" for horizontal
3269    SetdKLpar(Cell[h_corr[j-1]].Fnum, Cell[h_corr[j-1]].Knum, +Dip, kick);
3270
3271    for (i = 1; i <= N_BPM; i++)
3272      VertCouple[N_BPM+(j-1)*N_HCOR+i] = 
3273        HVweight*(orbitN[i]-orbitP[i])*0.5/kick; // sign reversal
3274  } // hcorr cycle
3275
3276
3277  // Find off diagonal terms for vertical trims
3278  for (j = 1; j <= N_VCOR; j++){
3279    // positive kick: "-Dip" for vertical
3280    SetdKLpar(Cell[v_corr[j-1]].Fnum, Cell[v_corr[j-1]].Knum, -Dip, kick);
3281    cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3282    for (i = 1;  i <= N_BPM; i++)
3283      orbitP[i] = Cell[bpm_loc[i-1]].BeamPos[x_];
3284
3285    // negative kick: "-Dip" for vertical
3286    SetdKLpar(Cell[v_corr[j-1]].Fnum, Cell[v_corr[j-1]].Knum, -Dip, -2*kick);
3287    cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3288    for (i = 1; i <= N_BPM; i++)
3289      orbitN[i] = Cell[bpm_loc[i-1]].BeamPos[x_];
3290
3291    // restore corrector: "-Dip" for vertical
3292    SetdKLpar(Cell[v_corr[j-1]].Fnum, Cell[v_corr[j-1]].Knum, -Dip, kick);
3293
3294    for (i = 1; i <= N_BPM; i++) 
3295      VertCouple[N_BPM+N_BPM*N_HCOR+(j-1)*N_VCOR+i] = 
3296        VHweight*(orbitP[i]-orbitN[i])*0.5/kick;
3297  } // vcorr cycle
3298
3299  free_dvector(orbitP, 1, N_BPM); free_dvector(orbitN, 1, N_BPM);
3300} // FindCoupVector
3301
3302
3303void SkewStat(double VertCouple[])
3304{
3305  int     i;
3306  double  max, rms, sk;
3307
3308  // statistics for skew quadrupoles
3309  max = 0.0; rms = 0.0;
3310  for(i = 1; i <= N_SKEW; i++) {
3311    sk = GetKLpar(globval.qt, i, -Quad);
3312    if (fabs(sk) > max) max = fabs(sk);
3313    rms += sqr(sk);
3314  }
3315  rms = sqrt(rms/N_SKEW);
3316  printf("Rms skew strength:       %8.2e+/-%8.2e\n", max, rms);
3317
3318  // statistics for vertical dispersion function
3319  max = 0.0; rms = 0.0;
3320  for(i = 1; i <= N_BPM; i++) {
3321    if (fabs(VertCouple[i]) > max) max = fabs(VertCouple[i]/VDweight);
3322    rms += sqr(VertCouple[i]/VDweight);
3323  }
3324  rms = sqrt(rms/N_BPM);
3325  printf("Max vertical dispersion: %8.2e+/-%8.2e mm\n", 1e3*max, 1e3*rms);
3326
3327  // statistics for off diagonal terms of response matrix (trims->bpms)
3328  max = 0.0; rms = 0.0;
3329  for(i = N_BPM+1; i <= N_BPM*(1+N_HCOR); i++) {
3330    if (fabs(VertCouple[i]) > max) max = fabs(VertCouple[i]/HVweight);
3331    rms += sqr(VertCouple[i]/HVweight);
3332  }
3333  rms = sqrt(rms/(N_HCOR*N_BPM));
3334  printf("Max horizontal coupling: %8.2e+/-%8.2e mm/mrad\n", max, rms);
3335
3336  max = 0.0; rms = 0.0;
3337  for(i = N_BPM*(1+N_HCOR)+1; i <= N_COUPLE; i++) {
3338    if (fabs(VertCouple[i]) > max) max = fabs(VertCouple[i]/VHweight);
3339    rms += sqr(VertCouple[i]/VHweight);
3340  }
3341  rms = sqrt(rms/(N_VCOR*N_BPM));
3342  printf("Max vertical coupling:   %8.2e+/-%8.2e mm/mrad\n", max, rms);
3343}
3344
3345
3346void corr_eps_y(void)
3347{
3348  int   i, j;
3349  FILE  *outf;
3350
3351  // Clear skew quad setpoints
3352  set_bnL_design_fam(globval.qt, Quad, 0.0, 0.0);
3353
3354  // Find coupling vector
3355  printf("\n");
3356  printf("Looking for coupling error\n");
3357  FindCoupVector(VertCouple);
3358
3359  //Find and print coupling statistics
3360  printf("\n");
3361  printf("Before correction\n");
3362  SkewStat(VertCouple);
3363
3364  // Coupling Correction
3365  printf("\n");
3366  for (i = 1; i <= n_lin; i++) {
3367    printf("Looking for correction\n");
3368
3369    //Find Correcting Settings to skew quadrupoles
3370    for (j = 1; j <= N_BPM; j++)
3371      b[j] = VDweight*eta_y[j] - VertCouple[j];
3372
3373    for (j = N_BPM+1; j <= N_COUPLE; j++)
3374      b[j] = -VertCouple[j];
3375
3376    dsvbksb(U, w, V, N_COUPLE, N_SKEW, b, SkewStrengthCorr);
3377
3378    printf("Applying correction\n");
3379    // Add correction
3380    for (j = 1; j <= N_SKEW; j++) 
3381      SetdKLpar(globval.qt, j, -Quad, SkewStrengthCorr[j]);
3382
3383    printf("\n");
3384    printf("Looking for coupling error\n");
3385    // Find coupling vector
3386    FindCoupVector(VertCouple);
3387
3388    printf("\n");
3389    printf("After run %d of correction\n", i);
3390    // Find and print coupling statistics
3391    SkewStat(VertCouple);
3392  } // End of coupling Correction
3393
3394  outf = file_write(eta_y_FileName);
3395  for (i = 0; i <= globval.Cell_nLoc; i++)
3396    fprintf(outf, "%4d %7.3f %s %6.3f %10.3e %10.3e\n",
3397            i, Cell[i].S, Cell[i].Elem.PName,
3398            Cell[i].Nu[Y_], 1e3*Cell[i].Eta[Y_], 1e3*Cell[i].Etap[Y_]);
3399  fclose(outf);
3400
3401  FindCoupVector(VertCouple);
3402}
3403
3404
3405/*******************************************************************
3406void get_IDs(void)
3407
3408  Purpose:
3409    Get the ID.
3410 
3411  Input:
3412     none
3413   
3414  Output:
3415    none   
3416 
3417  Comments:
3418    Print out the scaling 2 of Insertion.
3419   
3420********************************************************************/
3421// Control of IDs
3422
3423void get_IDs(void)
3424{
3425  int  k;
3426
3427  printf("\n");
3428  n_ID_Fams = 0;
3429  for (k = 0; k < globval.Elem_nFam; k++)
3430    switch (ElemFam[k].ElemF.Pkind) {
3431    case Wigl:
3432      printf("found ID family:   %s %12.5e\n",
3433             ElemFam[k].ElemF.PName, ElemFam[k].ElemF.W->BoBrhoV[0]);
3434      n_ID_Fams++; ID_Fams[n_ID_Fams-1] = k + 1;
3435      break;
3436    case Insertion:
3437      printf("found ID family:   %s scaling1 = %12.5e,  scaling2 = %12.5e\n",
3438             ElemFam[k].ElemF.PName, ElemFam[k].ElemF.ID->scaling1, ElemFam[k].ElemF.ID->scaling2);
3439      n_ID_Fams++; ID_Fams[n_ID_Fams-1] = k + 1;
3440      break;
3441    case FieldMap:
3442      printf("found ID family:   %s %12.5e\n",
3443             ElemFam[k].ElemF.PName, ElemFam[k].ElemF.FM->scl);
3444      n_ID_Fams++; ID_Fams[n_ID_Fams-1] = k + 1;
3445      break;
3446    default:
3447      break;
3448    }
3449}
3450
3451/****************************************************************************/
3452/* void set_IDs(const double scl)
3453
3454   Purpose:
3455       Set the scale factor for all the wigger, insertion devices, or fieldMap.       
3456
3457   Input:
3458       scl   scaling factor
3459
3460   Output:
3461      none
3462
3463   Return:
3464       none
3465
3466   Global variables:
3467       none
3468
3469   Specific functions:
3470       none
3471
3472   Comments:
3473       See also:
3474                 set_ID_scl(const int Fnum, const double scl)
3475                 set_ID_scl(const int Fnum, const int knum, const double scl)
3476       
3477****************************************************************************/
3478void set_IDs(const double scl)
3479{
3480  int  k;
3481
3482  printf("\n");
3483  for (k = 0; k < n_ID_Fams; k++) {
3484    switch (ElemFam[ID_Fams[k]-1].ElemF.Pkind) {
3485    case Wigl:
3486      printf("setting ID family: %s %12.5e\n",
3487             ElemFam[ID_Fams[k]-1].ElemF.PName,
3488             scl*ElemFam[ID_Fams[k]-1].ElemF.W->BoBrhoV[0]);
3489
3490      set_Wiggler_BoBrho(ID_Fams[k],
3491                         scl*ElemFam[ID_Fams[k]-1].ElemF.W->BoBrhoV[0]);
3492      break;
3493    case Insertion:
3494      printf("setting ID family: %s %12.5e\n",
3495             ElemFam[ID_Fams[k]-1].ElemF.PName, scl);
3496
3497      set_ID_scl(ID_Fams[k], scl);
3498      break;
3499    case FieldMap:
3500      printf("setting ID family: %s %12.5e\n",
3501             ElemFam[ID_Fams[k]-1].ElemF.PName, scl);
3502
3503      set_ID_scl(ID_Fams[k], scl);
3504      break;
3505    default:
3506      cout << "set_IDs: unknown element type" << endl;
3507      exit_(1);
3508      break;
3509    }
3510  }
3511}
3512
3513
3514void reset_quads(void)
3515{
3516  int       k;
3517
3518  if (N_Fam > N_Fam_max) {
3519    printf("reset_quads: N_Fam > N_Fam_max: %d (%d)\n", N_Fam, N_Fam_max);
3520    exit_(0);
3521  }
3522
3523  for (k = 0; k < N_Fam; k++) {
3524    // Note, actual values can differ from the original values
3525/*    printf("setting quad family: %s %12.5e\n",
3526           ElemFam[Q_Fam[k]-1].ElemF.PName,
3527           ElemFam[Q_Fam[k]-1].ElemF.M->PBpar[HOMmax+Quad]);
3528
3529    set_bn_design_fam(Q_Fam[k], Quad,
3530                       ElemFam[Q_Fam[k]-1].ElemF.M->PBpar[HOMmax+Quad], 0.0);*/
3531
3532    printf("setting quad family: %s %12.5e\n",
3533           ElemFam[Q_Fam[k]-1].ElemF.PName, b2[k]);
3534
3535    set_bn_design_fam(Q_Fam[k], Quad, b2[k], 0.0);
3536  }
3537}
3538
3539
3540void SVD(const int m, const int n, double **M, double beta_nu[],
3541         double b2Ls_[], const double s_cut, const bool first)
3542{
3543  int     i, j;
3544
3545  const bool    prt   = true;
3546
3547  if (first) {
3548    for (i = 1; i <= m; i++)
3549      for (j = 1; j <= n; j++)
3550        U1[i][j] = M[i][j];
3551
3552    dsvdcmp(U1, m, n, w1, V1);
3553
3554    if (prt) { 
3555      printf("\n");
3556      printf("singular values:\n");
3557      printf("\n");
3558    }
3559
3560    for (i = 1; i <= n; i++) {
3561      if (prt) printf("%11.3e", w1[i]);
3562      if (w1[i] < s_cut) {
3563        w1[i] = 0.0;
3564        if (prt) printf(" (zeroed)");
3565      }
3566      if (prt) if (i % 8 == 0) printf("\n");
3567    }
3568    if (prt) if (n % 8 != 0) printf("\n");
3569  }
3570 
3571  dsvbksb(U1, w1, V1, m, n, beta_nu, b2Ls_);
3572}
3573
3574
3575void quad_config()
3576{
3577  int     i, j;
3578  double  an;
3579
3580  if (N_Fam > N_Fam_max) {
3581    printf("quad_config: N_Fam > N_Fam_max: %d (%d)\n", N_Fam, N_Fam_max);
3582    exit_(0);
3583  }
3584
3585  Nquad = 0;
3586  for (i = 0; i < N_Fam; i++) {
3587    for (j = 1; j <= GetnKid(Q_Fam[i]); j++) {
3588      Nquad++;
3589
3590      if (Nquad > n_b2_max) {
3591        printf("quad_config: max no of quadrupoles exceeded %d (%d)\n",
3592               Nquad, n_b2_max);
3593        exit_(1);
3594      }
3595
3596      quad_prms[Nquad-1] = Elem_GetPos(Q_Fam[i], j);
3597
3598      if (j == 1) get_bn_design_elem(Q_Fam[i], j, Quad, b2[i], an);
3599    }
3600  }
3601
3602  printf("\n");
3603  printf("quad_config: Nquad = %d\n", Nquad);
3604}
3605
3606
3607bool get_SQ(void)
3608{
3609  int      j, k;
3610  Vector2  alpha3[3], beta3[3], nu3[3], eta3[3], etap3[3];
3611  FILE     *outf;
3612
3613  const bool  prt = false;
3614 
3615  /* Note, IDs are split for evaluation of the driving terms at the center:
3616       id1  1, 2
3617       id2  1, 2
3618       ...                                                                  */
3619
3620  // Get Twiss params, no dispersion
3621  Ring_GetTwiss(false, 0.0);
3622
3623  if (!status.codflag || !globval.stable) return false;
3624
3625  // Get global tunes
3626  Nu_X = globval.TotalTune[X_]; Nu_Y = globval.TotalTune[Y_];
3627
3628  if (prt) {
3629    printf("\n");
3630    printf("nu_x = %8.12f, nu_y = %8.12f\n", Nu_X, Nu_Y);
3631
3632    // Get Twiss params in sext
3633    printf("\n");
3634    printf("listing lattice functions at sextupoles\n");
3635    printf("\n");
3636
3637    outf = file_write("latfunS.out");
3638
3639    fprintf(outf, "s betax nux betay nuy\n");
3640  }
3641
3642  printf("\n");
3643  Nsext = 0;
3644  for (k = 0; k < globval.Cell_nLoc; k++) {
3645    if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == Sext)) {
3646      Nsext++;
3647
3648      if (Nsext > n_b3_max) {
3649        printf("get_SQ: max no of sextupoles exceeded %d (%d)\n",
3650               Nsext, n_b3_max);
3651        exit_(1);
3652      }
3653
3654      Ss[Nsext-1] = Cell[k].S;
3655
3656      for (j = 0; j <= 1; j++) {
3657        sb[j][Nsext-1] = Cell[k].Beta[j];
3658        sNu[j][Nsext-1] = Cell[k].Nu[j] - nu_0[j];
3659      }
3660
3661      if (prt) {
3662        printf("%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3663               Cell[k].Elem.PName, Ss[Nsext-1],
3664               sb[X_][Nsext-1], sNu[X_][Nsext-1]-nu_0[X_],
3665               sb[Y_][Nsext-1], sNu[Y_][Nsext-1]-nu_0[Y_]);
3666        fprintf(outf, "%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3667                Cell[k].Elem.PName, Ss[Nsext-1],
3668                sb[X_][Nsext-1], sNu[X_][Nsext-1]-nu_0[X_],
3669                sb[Y_][Nsext-1], sNu[Y_][Nsext-1]-nu_0[Y_]);
3670      }
3671    }
3672  }
3673
3674  if (prt) fclose(outf);
3675
3676  // Number of sexts in the ring
3677  printf("No of sextupoles = %d\n", Nsext);
3678
3679  if (prt) {
3680    // Get Twiss params in quads
3681    printf("\n");
3682    printf("listing lattice functions at quadrupoles\n");
3683    printf("\n");
3684
3685    outf = file_write("latfunQ.out");
3686
3687    fprintf(outf, "s name betax nux betay nuy\n");
3688  }
3689
3690  for (k = 0; k < Nquad; k++) {
3691    Sq[k] = Cell[quad_prms[k]].S;
3692    for (j = 0; j <= 1; j++) {
3693      if (Cell[quad_prms[k]].Elem.M->Pthick == thick) {
3694        get_twiss3(quad_prms[k], alpha3, beta3, nu3, eta3, etap3);
3695        qb[j][k] = beta3[Y_][j]; qNu[j][k] = nu3[Y_][j] - nu_0[j];
3696      } else {
3697        qb[j][k] = Cell[quad_prms[k]].Beta[j];
3698        qNu[j][k] = Cell[quad_prms[k]].Nu[j] - nu_0[j];
3699      }
3700    }
3701
3702    if (prt) {
3703      printf("%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3704             Cell[quad_prms[k]].Elem.PName, Sq[k], qb[X_][k],
3705             qNu[X_][k], qb[Y_][k], qNu[Y_][k]);
3706
3707      fprintf(outf, "%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3708              Cell[quad_prms[k]].Elem.PName, Sq[k], qb[X_][k],
3709              qNu[X_][k], qb[Y_][k], qNu[Y_][k]);
3710    }
3711  }
3712
3713  if (prt) fclose(outf);
3714
3715  // Number of quads in the ring
3716  printf("No of quads      = %d\n", Nquad);
3717
3718  return true;
3719}
3720
3721
3722double Bet(double bq, double nus, double nuq, double NuQ)
3723{
3724  return bq*cos(2.0*M_PI*(2.0*fabs(nus-nuq)-NuQ))/(2.0*sin(2.0*M_PI*NuQ));
3725}
3726
3727
3728double Nus(double bq, double nus, double nuq, double NuQ)
3729{
3730  double  Nu, sgn;
3731
3732  sgn = ((nus-nuq) <= 0)? -1: 1;
3733
3734  Nu = -bq*sgn*(sin(2.0*M_PI*NuQ)+sin(2.0*M_PI*(2.0*fabs(nus-nuq)-NuQ)))
3735       /(8.0*M_PI*sin(2.0*M_PI*NuQ));
3736
3737  return Nu;
3738}
3739
3740/*****************************************************************
3741void A_matrix(void)
3742
3743Purpose:
3744    solove A in the matrix  b = A.x
3745    using SVD.
3746    Used for ID correction.
3747
3748    Comments:
3749      change scl_dbeta  to  scl_dbetax  and scl_dbetay
3750      change scl_dnu  to  scl_dnux  and scl_dnuy
3751      change scl_nu  to  scl_nux  and scl_nuy
3752 
3753     
3754      See also X_vector ()   
3755******************************************************************/     
3756void A_matrix(void)
3757{
3758  int     k, j;
3759  double  BtX, BtY, NuX, NuY;
3760
3761  const bool  prt = false;
3762
3763  // Defining Twiss in undisturbed quads
3764  for (k = 0; k < Nquad; k++)
3765    for (j = 0; j <= 1; j++) {
3766      qb0[j][k] = qb[j][k]; 
3767      qNu0[j][k] = qNu[j][k];
3768    }
3769 
3770  // Defining Twiss in undisturbed sexts
3771  for (k = 0; k < Nsext; k++)
3772    for (j = 0; j <= 1; j++)
3773      sNu0[j][k] = sNu[j][k];
3774 
3775  // Now creating matrix A in X=A*B2L
3776  for (k = 1; k <= Nsext; k++) {
3777    for (j = 1; j <= Nquad; j++) {
3778      BtX = Bet(qb0[X_][j-1], sNu0[X_][k-1], qNu0[X_][j-1], Nu_X0);
3779      NuX = -Nus(qb0[X_][j-1], sNu0[X_][k-1], qNu0[X_][j-1], Nu_X0);
3780      BtY = -Bet(qb0[Y_][j-1], sNu0[Y_][k-1], qNu0[Y_][j-1], Nu_Y0);
3781      NuY = Nus(qb0[Y_][j-1], sNu0[Y_][k-1], qNu0[Y_][j-1], Nu_Y0);
3782      A1[k][j] = scl_dbetax*BtX;
3783      A1[k+Nsext][j] = scl_dbetay*BtY;
3784      A1[k+2*Nsext][j] = scl_dnux*NuX;
3785      A1[k+3*Nsext][j] = scl_dnuy*NuY;
3786    }
3787  }
3788  // Now adding 2 more constraints for global tunes
3789  for (j = 1; j <= Nquad; j++) {
3790    A1[4*Nsext+1][j] = -scl_nux*qb0[X_][j-1]/(4.0*M_PI);
3791    A1[4*Nsext+2][j] =  scl_nuy*qb0[Y_][j-1]/(4.0*M_PI);
3792  }
3793
3794  if (prt) {
3795    printf("\n");
3796    printf("AA:\n");
3797    printf("\n");
3798    for (k = 1; k <= Nconstr; k++) {
3799      for (j = 1; j <= Nquad; j++)
3800        printf(" %10.3e", A1[k][j]);
3801      printf("\n");
3802    }
3803  }
3804}
3805
3806/****************************************************
3807void X_vector(const bool first)
3808
3809  Purpose:
3810    solove x in the matrix  b = A.x
3811    using SVD.
3812    Used for ID correction.
3813
3814    Comments:
3815      change scl_dbeta  to  scl_dbetax  and scl_dbetay
3816      change scl_dnu  to  scl_dnux  and scl_dnuy
3817      change scl_nu  to  scl_nux  and scl_nuy
3818   
3819      See also: 
3820      void A_matrix(void)
3821*****************************************************/
3822void X_vector(const bool first)
3823{
3824  int  k;
3825
3826  const bool  prt = false;
3827
3828  dnu0[X_] = globval.TotalTune[X_] - Nu_X0;
3829  dnu0[Y_] = globval.TotalTune[Y_] - Nu_Y0;
3830
3831  if (first) {
3832    // Initial fill of X
3833    for (k = 1; k <= Nsext; k++) {
3834      Xsext0[k]         = sb[X_][k-1];  Xsext0[k+Nsext]   = sb[Y_][k-1];
3835      Xsext0[k+2*Nsext] = sNu[X_][k-1]; Xsext0[k+3*Nsext] = sNu[Y_][k-1];
3836    }
3837    Xsext0[4*Nsext+1] = 0.0; Xsext0[4*Nsext+2] = 0.0;
3838  } else { 
3839    // Now substracting from X in X=A*B2L
3840    for (k = 1; k <= Nsext; k++) {
3841      Xsext[k]         = scl_dbetax*(Xsext0[k]-sb[X_][k-1])/sb[X_][k-1];
3842      Xsext[k+Nsext]   = scl_dbetay*(Xsext0[k+Nsext]-sb[Y_][k-1])/sb[Y_][k-1];
3843      Xsext[k+2*Nsext] = scl_dnux*(Xsext0[k+2*Nsext]-sNu[X_][k-1]+dnu0[X_]/2.0);
3844      Xsext[k+3*Nsext] = scl_dnuy*(Xsext0[k+3*Nsext]-sNu[Y_][k-1]+dnu0[Y_]/2.0);
3845    }
3846    Xsext[4*Nsext+1] = scl_nux*(Nu_X0-globval.TotalTune[X_]);
3847    Xsext[4*Nsext+2] = scl_nuy*(Nu_Y0-globval.TotalTune[Y_]);
3848  }
3849
3850  if (prt) {
3851    printf("\n");
3852    printf("X:\n");
3853    printf("\n");
3854    if (first) {
3855      for (k = 1; k <= Nconstr; k++) {
3856        printf(" %10.3e", Xsext0[k]);
3857        if (k % 10 == 0)  printf("\n");
3858      }
3859      if (Nconstr % 10 != 0) printf("\n");
3860    } else {
3861      for (k = 1; k <= Nconstr; k++) {
3862        printf(" %10.3e", Xsext[k]);
3863        if (k % 10 == 0)  printf("\n");
3864      }
3865      if (Nconstr % 10 != 0) printf("\n");
3866    }
3867  }
3868}
3869
3870
3871/*
3872void ini_ID_corr(void)
3873
3874  Purpose:
3875    Initial ID, used for insertion Device correction
3876
3877   
3878*/
3879void ini_ID_corr(void)
3880{
3881
3882  // store ID families
3883  get_IDs();
3884
3885  // zero ID's
3886  set_IDs(0.0);
3887
3888  // Configuring quads (1 --> C means thin quads located in the middle of 1s)
3889  quad_config();
3890
3891  // Configuring quads (1 --> C means thin quads located in the middle of 1s)
3892  // Read Betas and Nus
3893  get_SQ(); 
3894  Nconstr = 4*Nsext + 2;
3895
3896  // Note, allocated vectors and matrices are deallocated in ID_corr
3897  Xsext = dvector(1, Nconstr); Xsext0 = dvector(1, Nconstr);
3898  b2Ls_ = dvector(1, Nquad); A1 = dmatrix(1, Nconstr, 1, Nquad);
3899  U1 = dmatrix(1, Nconstr, 1, Nquad); w1 = dvector(1, Nquad);
3900  V1 = dmatrix(1, Nquad, 1, Nquad);
3901
3902  // shift zero point to center of ID
3903//  nu_0[X_] = Cell[id_loc].Nu[X_]; nu_0[Y_] = Cell[id_loc].Nu[Y_];
3904  nu_0[X_] = 0.0; nu_0[Y_] = 0.0;
3905
3906  // Defining undisturbed tunes
3907  Nu_X0 = globval.TotalTune[X_]; Nu_Y0 = globval.TotalTune[Y_];
3908
3909  // Set-up matrix A in X=A*b2Ls_
3910  A_matrix();
3911
3912  // Now fill the X in X=A*b2Ls_
3913  X_vector(true);
3914}
3915
3916
3917void W_diag(void)
3918{
3919  double    bxf, byf, nxf, nyf, b2Lsum;
3920  int       k;
3921
3922  bxf = 0.0; byf = 0.0; nxf = 0.0; nyf = 0.0;
3923  for (k = 1; k <= Nsext; k++) {
3924    bxf += sqr(Xsext[k]);
3925    byf += sqr(Xsext[k+Nsext]);
3926    nxf += sqr(Xsext[k+2*Nsext]);
3927    nyf += sqr(Xsext[k+3*Nsext]);
3928  }
3929
3930  dnu0[X_] = globval.TotalTune[X_] - Nu_X0;
3931  dnu0[Y_] = globval.TotalTune[Y_] - Nu_Y0;
3932
3933  b2Lsum = 0.0;
3934  for (k = 1; k <= Nquad; k++) 
3935    b2Lsum += sqr(b2Ls_[k]);
3936
3937  printf("\n");
3938  printf("Residuals: beta [%%], dnu : \n");
3939  printf("dbeta_x: %6.2f dbeta_y: %6.2f nu_x: %12.6e nu_y: %12.6e\n",
3940         sqrt(bxf)/Nsext*1e2, sqrt(byf)/Nsext*1e2,
3941         sqrt(nxf)/Nsext, sqrt(nyf)/Nsext);
3942  printf("tune shift: dnu_x = %8.5f, dnu_y = %8.5f\n", dnu0[X_], dnu0[Y_]);
3943  printf("Sum b2Ls_: %12.6e\n", sqrt(b2Lsum)/Nquad);
3944}
3945
3946
3947
3948
3949
3950
3951/*
3952   For ID correction, from Mao-Sen Qiu at Taiwan light source.
3953*/
3954
3955bool ID_corr0(void)
3956{
3957  int     i, j, k;
3958  double  b2, a2;
3959  FILE    *outf;
3960
3961  printf("\n");
3962  printf("Before correction!\n");
3963
3964  set_IDs(1.0);
3965  //set_IDs((double) i);
3966
3967  get_SQ();                               // Read Betas and Nus
3968  X_vector(false);                        // Fill in dX in dX=A*db2Ls_
3969  W_diag();                               // Get statistics
3970
3971  outf = file_write("ID_before_corr.out");
3972  fprintf(outf, "# dbeta_x/beta_x  dbeta_y/beta_y  dnu_x  dnu_y\n");
3973  fprintf(outf, "#      [%%]             [%%]\n");
3974  fprintf(outf, "#\n");
3975  for (k = 1; k <= Nsext; k++)
3976    fprintf(outf, "%6.1f %6.2f %6.2f %10.3e %10.3e\n", Ss[k-1], 1e2*Xsext[k], 1e2*Xsext[k+Nsext], Xsext[k+2*Nsext], Xsext[k+3*Nsext]);
3977  fclose(outf);
3978
3979  // Allow for repeated calls to ID_corr, allocation is done in ini_ID_corr.
3980  if (true) {
3981    free_dvector(Xsext, 1, Nconstr);
3982    free_dvector(Xsext0, 1, Nconstr);
3983    free_dvector(b2Ls_, 1, Nquad); 
3984    free_dmatrix(A1, 1, Nconstr, 1, Nquad);
3985    free_dmatrix(U1, 1, Nconstr, 1, Nquad); 
3986    free_dvector(w1, 1, Nquad);
3987    free_dmatrix(V1, 1, Nquad, 1, Nquad);
3988  }
3989
3990  printf("\n");
3991  printf("End of before correction!\n");
3992
3993  return true;
3994}
3995
3996
3997/******************************************************************
3998bool ID_corr(const int N_calls, const int N_steps)
3999
4000  Purpose:
4001    Compensate the field error introduced by insertion device.
4002 
4003******************************************************************/
4004bool ID_corr(const int N_calls, const int N_steps)
4005{
4006  int     i, j, k;
4007  double  b2, a2;
4008  FILE    *outf;
4009
4010  printf("\n");
4011  printf("ID matching begins!\n");
4012
4013  outf = file_write("ID_corr.out");
4014  for (i = 1; i <= N_steps; i++) { //This brings ID strength in steps
4015    set_IDs((double)i/(double)N_steps);
4016
4017    get_SQ();                               // Read Betas and Nus
4018    X_vector(false);                        // Fill in dX in dX=A*db2Ls_
4019    W_diag();                               // Get statistics
4020    for (j = 1; j <= N_calls; j++) {
4021      SVD(Nconstr, Nquad, A1, Xsext, b2Ls_, 1e0, j == 1);
4022
4023      if ((i == N_steps) && (j == N_calls)) {
4024        fprintf(outf, "b_2:\n");
4025        fprintf(outf, "\n");
4026      }
4027
4028      // add quad strengths (db2Ls_)
4029      for (k = 1; k <= Nquad; k++) {
4030        set_dbnL_design_elem(Cell[quad_prms[k-1]].Fnum,
4031                             Cell[quad_prms[k-1]].Knum, Quad,
4032                             -ID_step*b2Ls_[k], 0.0);
4033
4034        if ((i == N_steps) && (j == N_calls)) {
4035          get_bn_design_elem(Cell[quad_prms[k-1]].Fnum,
4036                             Cell[quad_prms[k-1]].Knum, Quad, b2, a2);
4037          fprintf(outf, "%10s %2d %8.5f\n",
4038                  Cell[quad_prms[k-1]].Elem.PName, k, b2);
4039        }
4040      }
4041
4042      printf("\n");
4043      printf("Iteration: %2d\n", j);
4044      if (get_SQ()) {
4045        X_vector(false);                    // Fill in dX in dX=A*db2Ls_
4046        W_diag();                           // Get statistics
4047
4048        printglob();
4049      } else {
4050        printf("ID_corr: correction failed\n");
4051        // restore lattice
4052        set_IDs(0.0); reset_quads();
4053        return false;
4054      }
4055    }
4056  }
4057  fclose(outf);
4058
4059  outf = file_write("ID_corr_res.out");
4060  fprintf(outf, "# dbeta_x/beta_x  dbeta_y/beta_y  dnu_x  dnu_y\n");
4061  fprintf(outf, "#      [%%]             [%%]\n");
4062  fprintf(outf, "#\n");
4063  for (k = 1; k <= Nsext; k++)
4064    fprintf(outf, "%6.1f %6.2f %6.2f %10.3e %10.3e\n",
4065            Ss[k-1], 1e2*Xsext[k], 1e2*Xsext[k+Nsext],
4066            Xsext[k+2*Nsext], Xsext[k+3*Nsext]);
4067  fclose(outf);
4068
4069  // Allow for repeated calls to ID_corr, allocation is done in ini_ID_corr.
4070  if (false) {
4071    free_dvector(Xsext, 1, Nconstr); free_dvector(Xsext0, 1, Nconstr);
4072    free_dvector(b2Ls_, 1, Nquad); free_dmatrix(A1, 1, Nconstr, 1, Nquad);
4073    free_dmatrix(U1, 1, Nconstr, 1, Nquad); free_dvector(w1, 1, Nquad);
4074    free_dmatrix(V1, 1, Nquad, 1, Nquad);
4075  }
4076
4077  printf("\n");
4078  printf("ID matching ends!\n");
4079
4080  return true;
4081}
4082
4083/*  End ID correction functions *****/
4084
4085
4086void get_param(const char *param_file)
4087/*
4088 * Read parameter file in input file xxxx.prm
4089 *
4090 */
4091{
4092  char    line[max_str], name[max_str], str[max_str], s_prm[max_str];
4093  char    lat_file[max_str], flat_file[max_str], IDCq_name[max_str][8];
4094  int     k;
4095  double  f_prm;
4096  FILE    *inf;
4097
4098  const bool  prt = true;
4099
4100  if (prt) {
4101    printf("\n");
4102    printf("reading in %s\n", param_file);
4103  }
4104
4105  inf = file_read(param_file);
4106
4107  // read parameter file
4108  // initialization
4109  strcpy(ae_file, ""); strcpy(fe_file, ""); strcpy(ap_file, "");
4110
4111  if (prt) printf("\n");
4112  // loop over all lines of the file
4113  while (fgets(line, max_str, inf) != NULL) {
4114    if (prt) printf("%s", line);
4115    if (strstr(line, "#") == NULL) {
4116      // get initial command token
4117      sscanf(line, "%s", name);
4118        // input files *************************************
4119      if (strcmp("in_dir", name) == 0)
4120        sscanf(line, "%*s %s", in_dir);
4121      else if (strcmp("ae_file", name) == 0){ 
4122        sscanf(line, "%*s %s", str);
4123        sprintf(ae_file,"%s/%s", in_dir, str);
4124      } else if (strcmp("fe_file", name) == 0) { 
4125        sscanf(line, "%*s %s", str);
4126        sprintf(fe_file, "%s/%s", in_dir, str);
4127      } else if (strcmp("ap_file", name) == 0) {
4128        sscanf(line, "%*s %s", str);
4129        sprintf(ap_file, "%s/%s", in_dir, str);
4130      } else if (strcmp("lat_file", name) == 0) {
4131        sscanf(line, "%*s %s", lat_file);
4132        sprintf(lat_FileName, "%s/%s", in_dir, lat_file);
4133        Read_Lattice(lat_FileName);
4134      } else if (strcmp("flat_file", name) == 0) {
4135        sscanf(line, "%*s %s", flat_file);
4136        strcat(flat_file, "flat_file.dat"); rdmfile(flat_file);
4137      // **********< parameters >****************************
4138      } else if (strcmp("s_cut", name) == 0) {
4139        sscanf(line, "%*s %lf", &f_prm);
4140        setrancut(f_prm);
4141      } else if (strcmp("n_stat", name) == 0)
4142        sscanf(line, "%*s %d", &n_stat);
4143      else if (strcmp("n_scale", name) == 0)
4144        sscanf(line, "%*s %d", &n_scale);
4145      else if (strcmp("n_orbit", name) == 0)
4146        sscanf(line, "%*s %d", &n_orbit);
4147      else if (strcmp("bpm_name", name) == 0) {
4148        sscanf(line, "%*s %s", s_prm);
4149        globval.bpm = ElemIndex(s_prm);
4150      } else if (strcmp("h_corr", name) == 0) {
4151        sscanf(line, "%*s %s", s_prm);
4152        globval.hcorr = ElemIndex(s_prm);
4153      } else if (strcmp("v_corr", name) == 0) {
4154        sscanf(line, "%*s %s", s_prm);
4155        globval.vcorr = ElemIndex(s_prm);
4156      } else if (strcmp("gs", name) == 0) {
4157        sscanf(line, "%*s %s", s_prm);
4158        globval.gs = ElemIndex(s_prm);
4159      } else if (strcmp("ge", name) == 0) {
4160        sscanf(line, "%*s %s", s_prm);
4161        globval.ge = ElemIndex(s_prm);
4162      } else if (strcmp("DA_bare", name) == 0) {
4163        sscanf(line, "%*s %s", s_prm);
4164        DA_bare = (strcmp(s_prm, "true") == 0)? true : false;
4165      } else if (strcmp("delta", name) == 0) { 
4166        sscanf(line, "%*s %lf", &delta_DA_);
4167      } else if (strcmp("freq_map", name) == 0) {
4168        sscanf(line, "%*s %s %d %d %d %d %lf %lf %lf",
4169               s_prm, &n_x, &n_y, &n_dp, &n_tr,
4170               &x_max_FMA, &y_max_FMA, &delta_FMA);
4171        freq_map = (strcmp(s_prm, "true") == 0)? true : false;
4172      } else if (strcmp("tune_shift", name) == 0) {
4173        sscanf(line, "%*s %s", s_prm);
4174        tune_shift = (strcmp(s_prm, "true") == 0)? true : false;
4175      } else if (strcmp("qt", name) == 0) {
4176        sscanf(line, "%*s %s", s_prm);
4177        globval.qt = ElemIndex(s_prm);
4178      } else if (strcmp("disp_wave_y", name) == 0) 
4179        sscanf(line, "%*s %lf", &disp_wave_y);
4180      else if (strcmp("n_lin", name) == 0)
4181        sscanf(line, "%*s %d", &n_lin);
4182      else if (strcmp("VDweight", name) == 0) 
4183        sscanf(line, "%*s %lf", &VDweight);
4184      else if (strcmp("HVweight", name) == 0) 
4185        sscanf(line, "%*s %lf", &HVweight);
4186      else if (strcmp("VHweight", name) == 0)
4187        sscanf(line, "%*s %lf", &VHweight);
4188      else if (strcmp("N_calls", name) == 0) // ID correction parameters
4189        sscanf(line, "%*s %d", &N_calls);
4190      else if (strcmp("N_steps", name) == 0)
4191        sscanf(line, "%*s %d", &N_steps);
4192      else if (strcmp("N_Fam", name) == 0)
4193        sscanf(line, "%*s %d", &N_Fam);
4194      else if (strcmp("IDCquads", name) == 0) {
4195        sscanf(line, "%*s %s %s %s %s %s %s %s %s",
4196               IDCq_name[0], IDCq_name[1], IDCq_name[2], IDCq_name[3],
4197               IDCq_name[4], IDCq_name[5], IDCq_name[6], IDCq_name[7]);
4198//      } else if (strcmp("scl_nu", name) == 0)
4199//      sscanf(line, "%*s %lf", &scl_nu);
4200      } else {
4201        printf("bad line in %s: %s\n", param_file, line);
4202        exit_(1);
4203      }
4204    } else
4205      printf("%s", line);
4206  }
4207
4208  fclose(inf);
4209
4210  if (N_calls > 0) {
4211    if (N_Fam > N_Fam_max) {
4212      printf("get_param: N_Fam > N_Fam_max: %d (%d)\n",
4213             N_Fam, N_Fam_max);
4214      exit_(0);
4215    }
4216
4217    for (k = 0; k < N_Fam; k++)
4218      Q_Fam[k] = ElemIndex(IDCq_name[k]);
4219  }
4220}
4221
4222
4223// control of vertical beam size
4224
4225// output procedures
4226
4227// Finding statistics on orbit in the sextupoles and maximal trim settings
4228void Orb_and_Trim_Stat(void)
4229{
4230  int     i, j, N;
4231  int     SextCounter = 0;
4232  int     bins[5] = { 0, 0, 0, 0, 0 };
4233  double  bin = 40.0e-6; // bin size for stat
4234  double  tr; // trim strength
4235
4236  Vector2   Sext_max, Sext_sigma, TrimMax, orb;
4237 
4238  Sext_max[X_] = 0.0; Sext_max[Y_] = 0.0; 
4239  Sext_sigma[X_] = 0.0; Sext_sigma[Y_] = 0.0;
4240  TrimMax[X_] = 0.0; TrimMax[Y_] = 0.0;
4241  N = globval.Cell_nLoc; SextCounter = 0;
4242  for (i = 0; i <= N; i++) {
4243    if ((Cell[i].Elem.Pkind == Mpole) && (Cell[i].Elem.M->n_design == Sext)) {
4244      SextCounter++;
4245      orb[X_] = Cell[i].BeamPos[x_]; orb[Y_] = Cell[i].BeamPos[y_];
4246      Sext_sigma[X_] += orb[X_]*orb[X_]; Sext_sigma[Y_] += orb[Y_]*orb[Y_];
4247      if (fabs(orb[X_]) > Sext_max[X_]) Sext_max[X_] = fabs(orb[X_]);
4248      if (fabs(orb[Y_]) > Sext_max[Y_]) Sext_max[Y_] = fabs(orb[Y_]);
4249      j = (int) (sqrt(orb[X_]*orb[X_]+orb[Y_]*orb[Y_])/bin);
4250      if (j > 4) j = 4;
4251      bins[j]++;
4252    } // sextupole handling
4253
4254    if (Cell[i].Fnum == globval.hcorr) {
4255      tr = Cell[i].Elem.M->PBpar[HOMmax+Dip];
4256      if (fabs(tr) > TrimMax[X_]) TrimMax[X_] = fabs(tr);
4257    } // horizontal trim handling
4258    if (Cell[i].Fnum == globval.vcorr) {
4259      tr = Cell[i].Elem.M->PBpar[HOMmax-Dip];
4260      if (fabs(tr) > TrimMax[Y_]) TrimMax[Y_] = fabs(tr);
4261    } // vertical trim handling
4262  } // looking throught the cells
4263
4264  Sext_sigma[X_] = sqrt(Sext_sigma[X_]/SextCounter);
4265  Sext_sigma[Y_] = sqrt(Sext_sigma[Y_]/SextCounter);
4266   
4267  printf("In sextupoles maximal horizontal orbit is:"
4268         " %5.3f mm with sigma %5.3f mm\n",
4269          1e3*Sext_max[X_], 1e3*Sext_sigma[X_]);
4270  printf("and maximal vertical orbit is:            "
4271         " %5.3f mm with sigma %5.3f mm.\n",
4272         1e3*Sext_max[Y_], 1e3*Sext_sigma[Y_]);
4273
4274  for (i = 0; i < 4;  i++) {
4275    printf("There are %3d sextupoles with offset between ", bins[i]);
4276    printf(" %5.3f mm and %5.3f mm\n", i*bin*1e3, (i+1)*bin*1e3);
4277  }
4278  printf("There are %3d sextupoles with offset ", bins[4]);
4279  printf("more than %5.3f mm \n", 4e3*bin);
4280  printf("Maximal hcorr is %6.3f mrad, maximal vcorr is %6.3f mrad\n",
4281         1e3*TrimMax[X_], 1e3*TrimMax[Y_]);
4282}
4283
4284
4285void prt_codcor_lat(void)
4286{
4287  int   i;
4288  FILE  *CodCorLatFile;
4289
4290
4291  CodCorLatFile = file_write(CodCorLatFileName);
4292
4293  fprintf(CodCorLatFile, "#    name     s   sqrt(BxBy) betaX    nuX"
4294          "    betaY    nuY     etaX etaX*betaY nuX-nuY \n");
4295  fprintf(CodCorLatFile, "#            [m]     [m]      [m]             [m]"
4296          "              [m]     [m*m] \n");
4297
4298  for (i = 0; i <= globval.Cell_nLoc; i++){
4299    fprintf(CodCorLatFile, "%4d:", i);
4300
4301    if (i == 0)
4302      fprintf(CodCorLatFile, "%.*s", 6, "begin ");
4303    else
4304      fprintf(CodCorLatFile, "%.*s", 6, Cell[i].Elem.PName);
4305
4306    fprintf(CodCorLatFile, "%7.3f  %5.2f    %5.2f  %7.4f  %5.2f  %7.4f"
4307            "  %6.3f  %6.3f  %6.3f\n",
4308            Cell[i].S, sqrt(Cell[i].Beta[X_]*Cell[i].Beta[Y_]), 
4309            Cell[i].Beta[X_], Cell[i].Nu[X_], Cell[i].Beta[Y_], Cell[i].Nu[Y_],
4310            Cell[i].Eta[X_], Cell[i].Eta[X_]*Cell[i].Beta[Y_],
4311            Cell[i].Nu[X_]-Cell[i].Nu[Y_]);
4312  }
4313  fclose(CodCorLatFile);
4314}
4315
4316void prt_beamsizes()
4317{
4318  int   k;
4319  FILE  *fp;
4320
4321  fp = file_write(beam_envelope_file);
4322
4323  fprintf(fp,"# k  s  name s_xx s_pxpx s_xpx s_yy s_pypy s_ypy theta_xy\n");
4324  for(k = 0; k <= globval.Cell_nLoc; k++){
4325    fprintf(fp,"%4d %10s %e %e %e %e %e %e %e %e\n",
4326            k, Cell[k].Elem.PName, Cell[k].S,
4327            Cell[k].sigma[x_][x_], Cell[k].sigma[px_][px_],
4328            Cell[k].sigma[x_][px_],
4329            Cell[k].sigma[y_][y_], Cell[k].sigma[py_][py_],
4330            Cell[k].sigma[y_][py_],
4331            atan2(2e0*Cell[k].sigma[x_][y_],
4332                  Cell[k].sigma[x_][x_]-Cell[k].sigma[y_][y_])/2e0*180.0/M_PI);
4333  }
4334
4335  fclose(fp);
4336}
4337
4338
4339float f_int_Touschek(const float u)
4340{
4341  double  f;
4342
4343  if (u > 0.0)
4344    f = (1.0/u-log(1.0/u)/2.0-1.0)*exp(-u_Touschek/u); 
4345  else
4346    f = 0.0;
4347
4348  return f;
4349} 
4350
4351/****************************************************************************/
4352/* double Touschek(const double Qb, const double delta_RF,
4353                const double eps_x, const double eps_y,
4354                const double sigma_delta, const double sigma_s)
4355
4356   Purpose:
4357       Calculate Touschek lifetime
4358       using the normal theoretical formula
4359       the momentum acceptance is RF acceptance delta_RF
4360     
4361   Parameters:   
4362             Qb  electric charge per bunch
4363       delta_RF  RF momentum acceptance
4364          eps_x  emittance x
4365          eps_y  emittance y
4366    sigma_delta  longitudinal beam size
4367        sigma_s   
4368   
4369   Input:
4370       none
4371               
4372   Output:
4373       
4374
4375   Return:
4376      Touschek lifetime [second]
4377
4378   Global variables:
4379       
4380
4381   Specific functions:
4382     
4383
4384   Comments:
4385       none
4386
4387****************************************************************************/
4388double Touschek(const double Qb, const double delta_RF,
4389                const double eps_x, const double eps_y,
4390                const double sigma_delta, const double sigma_s)
4391{
4392  long int  k;
4393  double    tau, sigma_x, sigma_y, sigma_xp, L, curly_H;
4394
4395  const double  gamma = 1e9*globval.Energy/m_e, N_e = Qb/q_e;
4396
4397  printf("\n");
4398  printf("Qb = %4.2f nC, delta_RF = %4.2f%%"
4399         ", eps_x = %9.3e nm.rad, eps_y = %9.3e nm.rad\n",
4400         1e9*Qb, 1e2*delta_RF, 1e9*eps_x, 1e9*eps_y);
4401  printf("sigma_delta = %8.2e, sigma_s = %4.2f mm\n",
4402         sigma_delta, 1e3*sigma_s);
4403
4404  tau = 0.0;
4405  for(k = 1; k <= globval.Cell_nLoc; k++) {
4406    L = Cell[k].S - Cell[k-1].S;
4407
4408    curly_H = get_curly_H(Cell[k].Alpha[X_], Cell[k].Beta[X_],
4409                          Cell[k].Eta[X_], Cell[k].Etap[X_]);
4410
4411    // compute estimated beam sizes for given
4412    // hor.,  ver. emittance, sigma_s, and sigma_delta (x ~ 0)
4413    sigma_x = sqrt(Cell[k].Beta[X_]*eps_x+sqr(sigma_delta*Cell[k].Eta[X_])); 
4414    sigma_y = sqrt(Cell[k].Beta[Y_]*eps_y); 
4415    sigma_xp = (eps_x/sigma_x)*sqrt(1.0+curly_H*sqr(sigma_delta)/eps_x); 
4416
4417    u_Touschek = sqr(delta_RF/(gamma*sigma_xp));
4418
4419    tau += qromb(f_int_Touschek, 0.0, 1.0)/(sigma_x*sigma_y*sigma_xp)*L; 
4420    fflush(stdout);
4421  }
4422
4423  tau *= N_e*sqr(r_e)*c0/(8.0*M_PI*cube(gamma)*sigma_s)
4424         /(sqr(delta_RF)*Cell[globval.Cell_nLoc].S);
4425
4426  printf("\n"); 
4427  printf("Touschek lifetime [hrs]: %10.3e\n", 1.0/(3600.0*tau)); 
4428  return(1.0/(tau));
4429}
4430
4431/****************************************************************************/
4432/* void mom_aper(double &delta, double delta_RF, const long int k,
4433              const int n_turn, const bool positive)
4434
4435
4436   Purpose:
4437      Using Binary search to determine momentum aperture at location k.
4438      If the particle is stable at the last step of tracking, so the
4439      momentum acceptance is RF momentum accpetance.
4440     
4441   Input:
4442           delta       particle energy deviation
4443           delta_RF    RF momentum acceptance
4444           k           element index
4445           n_turn      number of tracking turns 
4446
4447   Output:
4448           none
4449   Return:
4450           delta       searched momentum acceptance
4451
4452   Global variables:
4453         none
4454
4455   specific functions:
4456       Cell_Pass
4457
4458   Comments:
4459       nsls-ii version.
4460       See also:
4461       soleil version MomentumAcceptance( ) in soleillib.cc
4462       
4463****************************************************************************/
4464void mom_aper(double &delta, double delta_RF, const long int k,
4465              const int n_turn, const bool positive)
4466{
4467  // Binary search to determine momentum aperture at location k.
4468  int       j;
4469  long int  lastpos;
4470  double    delta_min, delta_max;
4471  Vector    x;
4472
4473  const double  eps = 1e-4;
4474
4475  delta_min = 0.0; 
4476  delta_max = positive ? fabs(delta_RF) : -fabs(delta_RF);
4477 
4478  while (fabs(delta_max-delta_min) > eps) 
4479  {
4480    delta = (delta_max+delta_min)/2.0; /* binary serach*/
4481   
4482    // propagate initial conditions
4483    CopyVec(6, globval.CODvect, x); Cell_Pass(0, k, x, lastpos); 
4484    // generate Touschek event
4485    x[delta_] += delta;
4486    // complete one turn
4487    Cell_Pass(k, globval.Cell_nLoc, x, lastpos); 
4488   
4489    if (lastpos < globval.Cell_nLoc)
4490      // particle lost
4491      delta_max = delta;
4492    else { 
4493      // track
4494      for(j = 0; j < n_turn; j++) {
4495        Cell_Pass(0, globval.Cell_nLoc, x, lastpos); 
4496
4497        if ((delta_max > delta_RF) || (lastpos < globval.Cell_nLoc)) { 
4498          // particle lost
4499          delta_max = delta; 
4500          break; 
4501        }
4502      }
4503     
4504      if ((delta_max <= delta_RF) && (lastpos == globval.Cell_nLoc))
4505        // particle not lost
4506        delta_min = delta; 
4507    }
4508  } 
4509}
4510
4511/****************************************************************************/
4512/* double Touschek(const double Qb, const double delta_RF,const bool consistent,
4513                const double eps_x, const double eps_y,
4514                const double sigma_delta, double sigma_s,
4515                const int n_turn, const bool aper_on,
4516                double sum_delta[][2], double sum2_delta[][2])
4517
4518   Purpose:
4519        Calculate the Touschek lifetime
4520        using the normal theoretical formula
4521        the momentum acceptance is determined by the RF acceptance delta_RF
4522        and the momentum aperture at each element location
4523        which is tracked over n turns,
4524        the vacuum chamber limitation is read from a outside file 
4525     
4526   Parameters:   
4527             Qb  electric charge per bunch
4528       delta_RF  RF momentum acceptance
4529     consistent
4530          eps_x  emittance x
4531          eps_y  emittance y
4532    sigma_delta  longitudinal beam size
4533        sigma_s
4534         n_turn  number of turns for tracking 
4535        aper_on   must be TRUE
4536                 to indicate vacuum chamber setting is read from a outside file         
4537 sum_delta[][2]
4538sum2_delta[][2]         
4539                 
4540   
4541   Input:
4542       none
4543               
4544   Output:
4545       
4546
4547   Return:
4548       Touschek lifetime [second]
4549
4550   Global variables:
4551       
4552
4553   Specific functions:
4554      mom_aper
4555
4556   Comments:
4557       none
4558
4559****************************************************************************/
4560double Touschek(const double Qb, const double delta_RF,const bool consistent,
4561                const double eps_x, const double eps_y,
4562                const double sigma_delta, double sigma_s,
4563                const int n_turn, const bool aper_on,
4564                double sum_delta[][2], double sum2_delta[][2])
4565{
4566  bool      cav, aper;
4567  long int  k;
4568  double    rate, delta_p, delta_m, curly_H0, curly_H1, L;
4569  double    sigma_x, sigma_y, sigma_xp;
4570
4571  const bool  prt = false;
4572
4573  //  const char  file_name[] = "Touschek.out";
4574  const double  eps = 1e-5, gamma = 1e9*globval.Energy/m_e, N_e = Qb/q_e;
4575
4576  cav = globval.Cavity_on; aper = globval.Aperture_on;
4577
4578  globval.Cavity_on = true;
4579
4580  Ring_GetTwiss(true, 0.0);
4581
4582  globval.Aperture_on = aper_on;
4583
4584  printf("\n");
4585  printf("Qb = %4.2f nC, delta_RF = %4.2f%%"
4586         ", eps_x = %9.3e m.rad, eps_y = %9.3e m.rad\n",
4587         1e9*Qb, 1e2*delta_RF, eps_x, eps_y);
4588  printf("sigma_delta = %8.2e, sigma_s = %4.2f mm\n",
4589         sigma_delta, 1e3*sigma_s);
4590
4591  printf("\n");
4592  printf("Momentum aperture:" );
4593  printf("\n");
4594
4595  delta_p = delta_RF; mom_aper(delta_p, delta_RF, 0, n_turn, true);
4596  delta_m = -delta_RF; mom_aper(delta_m, delta_RF, 0, n_turn, false);
4597  delta_p = min(delta_RF, delta_p); delta_m = max(-delta_RF, delta_m);
4598  sum_delta[0][0] += delta_p; sum_delta[0][1] += delta_m;
4599  sum2_delta[0][0] += sqr(delta_p); sum2_delta[0][1] += sqr(delta_m);
4600 
4601  rate = 0.0; curly_H0 = -1e30;
4602  for (k = 1; k <= globval.Cell_nLoc; k++) {
4603    L = Cell[k].S - Cell[k-1].S;
4604
4605    curly_H1 = get_curly_H(Cell[k].Alpha[X_], Cell[k].Beta[X_],
4606                           Cell[k].Eta[X_], Cell[k].Etap[X_]);
4607
4608    if (fabs(curly_H0-curly_H1) > eps) {
4609      mom_aper(delta_p, delta_RF, k, n_turn, true);
4610      delta_m = -delta_p; mom_aper(delta_m, delta_RF, k, n_turn, false);
4611      delta_p = min(delta_RF, delta_p); delta_m = max(-delta_RF, delta_m);
4612      printf("%4ld %6.2f %3.2lf%% %3.2lf%%\n",
4613             k, Cell[k].S, 1e2*delta_p, 1e2*delta_m);
4614      curly_H0 = curly_H1;
4615    }
4616
4617    sum_delta[k][0] += delta_p; 
4618    sum_delta[k][1] += delta_m;
4619    sum2_delta[k][0] += sqr(delta_p); 
4620    sum2_delta[k][1] += sqr(delta_m);
4621   
4622    if (prt)
4623      printf("%4ld %6.2f %3.2lf %3.2lf\n",
4624             k, Cell[k].S, 1e2*delta_p, 1e2*delta_m);
4625
4626    if (!consistent) {
4627      // compute estimated beam sizes for given
4628      // hor.,  ver. emittance, sigma_s, and sigma_delta (x ~ 0)
4629      sigma_x = sqrt(Cell[k].Beta[X_]*eps_x+sqr(sigma_delta*Cell[k].Eta[X_]));
4630      sigma_y = sqrt(Cell[k].Beta[Y_]*eps_y);
4631      sigma_xp = (eps_x/sigma_x)*sqrt(1.0+curly_H1*sqr(sigma_delta)/eps_x);
4632    } else {
4633      // use self-consistent beam sizes
4634      sigma_x = sqrt(Cell[k].sigma[x_][x_]);
4635      sigma_y = sqrt(Cell[k].sigma[y_][y_]);
4636      sigma_xp = sqrt(Cell[k].sigma[px_][px_]);
4637      sigma_s = sqrt(Cell[k].sigma[ct_][ct_]);
4638    }
4639
4640    u_Touschek = sqr(delta_p/(gamma*sigma_xp));
4641
4642    rate += qromb(f_int_Touschek, 0.0, 1.0)
4643            /(sigma_x*sigma_xp*sigma_y*sqr(delta_p))*L;
4644
4645    fflush(stdout);
4646  }
4647
4648  rate *= N_e*sqr(r_e)*c0/(8.0*M_PI*cube(gamma)*sigma_s)
4649         /Cell[globval.Cell_nLoc].S;
4650
4651  printf("\n");
4652  printf("Touschek lifetime [hrs]: %4.2f\n", 1.0/(3600.0*rate));
4653
4654  globval.Cavity_on = cav; globval.Aperture_on = aper;
4655
4656  return 1/rate;
4657}
4658
4659
4660double f_IBS(const double chi_m)
4661{
4662  // Interpolated
4663
4664  double  f, ln_chi_m;
4665
4666  const double A = 1.0, B = 0.579, C = 0.5;
4667
4668  ln_chi_m = log(chi_m); f = A + B*ln_chi_m + C*sqr(ln_chi_m);
4669
4670  return f;
4671}
4672
4673
4674float f_int_IBS(const float chi)
4675{
4676  // Integrated
4677
4678  double  f;
4679
4680  f = exp(-chi)*log(chi/chi_m)/chi;
4681
4682  return f;
4683}
4684
4685
4686void IBS(const double Qb,
4687         const double eps_SR[], double eps[],
4688         const double alpha_z, const double beta_z)
4689{
4690  /* The equilibrium emittance is given by:
4691
4692       sigma_delta^2 = tau_delta*(D_delta_IBS + D_delta_SR)
4693
4694       D_x = <D_delta*curly_H> =>
4695
4696       eps_x = eps_x_SR + eps_x_IBS,    D_x_IBS ~ 1/eps_x
4697
4698     i.e., the 2nd term is essentially constant.  From
4699
4700       eps_x^2 = eps_x*eps_x_SR + eps_x*eps_x_IBS
4701
4702     one obtains the IBS limited emittance (eps_x_SR = 0, eps_x^2=eps_x_IBS*tau_x*D_x_IBS->cst):
4703
4704       eps_x,IBS = sqrt(tau_x*(eps_x*D_x_IBS)
4705
4706     and by solving for eps_x
4707
4708       eps_x = eps_x_SR/2 + sqrt(eps_x_SR^2/4+eps_x_IBS^2)
4709
4710  */
4711
4712  long int  k;
4713  double    D_x, D_delta, b_max, L, gamma_z;
4714  double    sigma_x, sigma_xp, sigma_y, sigma_p, sigma_s, sigma_delta;
4715  double    incr, C, curly_H, eps_IBS[3], eps_x_tot, sigma_E;
4716  double    sigma_s_SR, sigma_delta_SR, sigma_delta_IBS;
4717
4718  const bool    integrate = false;
4719  const double  gamma = 1e9*globval.Energy/m_e, N_e = Qb/q_e;
4720
4721  // bunch size
4722  gamma_z = (1.0+sqr(alpha_z))/beta_z;
4723  sigma_s_SR = sqrt(beta_z*eps_SR[Z_]);
4724  sigma_delta_SR = sqrt(gamma_z*eps_SR[Z_]);
4725
4726  sigma_s = sqrt(beta_z*eps[Z_]); sigma_delta = sqrt(gamma_z*eps[Z_]);
4727
4728  printf("\n");
4729  printf("Qb             = %4.2f nC\n", 1e9*Qb);
4730  printf("eps_x_SR       = %9.3e m.rad, eps_x       = %9.3e m.rad\n",
4731         eps_SR[X_], eps[X_]);
4732  printf("eps_y_SR       = %9.3e m.rad, eps_y       = %9.3e m.rad\n",
4733         eps_SR[Y_], eps[Y_]);
4734  printf("eps_z_SR       = %9.3e,       eps_z       = %9.3e\n",
4735         eps_SR[Z_], eps[Z_]);
4736  printf("alpha_z        = %9.3e,       beta_z      = %9.3e\n",
4737         alpha_z, beta_z);
4738  printf("sigma_s_SR     = %9.3e mm,    sigma_s     = %9.3e mm\n",
4739         1e3*sigma_s_SR, 1e3*sigma_s);
4740  printf("sigma_delta_SR = %9.3e,       sigma_delta = %9.3e\n",
4741         sigma_delta_SR, sigma_delta);
4742
4743  D_delta = 0.0; D_x = 0.0;
4744  for(k = 1; k <= globval.Cell_nLoc; k++) {
4745    L = Cell[k].Elem.PL;
4746
4747    curly_H = get_curly_H(Cell[k].Alpha[X_], Cell[k].Beta[X_],
4748                          Cell[k].Eta[X_], Cell[k].Etap[X_]);
4749
4750    // compute estimated beam sizes for given
4751    // hor.,  ver. emittance, sigma_s, and sigma_delta (x ~ 0)
4752    sigma_x = sqrt(Cell[k].Beta[X_]*eps[X_]+sqr(sigma_delta*Cell[k].Eta[X_])); 
4753    sigma_y = sqrt(Cell[k].Beta[Y_]*eps[Y_]);
4754    sigma_xp = (eps[X_]/sigma_x)*sqrt(1.0+curly_H*sqr(sigma_delta)/eps[X_]);
4755
4756    sigma_E = gamma*m_e*q_e*sigma_delta;
4757    sigma_p = gamma*m_e*q_e*sigma_xp;
4758
4759//    b_max = sqrt(2.0*M_PI)/pow(N_e/(sigma_x*sigma_y*sigma_s), 1.0/3.0);
4760    b_max = 2.0*sqrt(M_PI)/pow(N_e/(sigma_x*sigma_y*sigma_s), 1.0/3.0);
4761
4762    chi_m = r_e*sqr(m_e*q_e/sigma_E)/b_max;
4763//    chi_m = r_e*sqr(m_e*q_e/sigma_p)/b_max;
4764
4765    if (!integrate)
4766      incr = f_IBS(chi_m)/sigma_y*L;
4767    else
4768      incr = qromb(f_int_IBS, chi_m, 1e5*chi_m)/sigma_y*L;
4769
4770    D_delta += incr; D_x += incr*curly_H;
4771  }
4772
4773  C = Cell[globval.Cell_nLoc].S;
4774
4775  if (true) {
4776    eps_x_tot = 1.0;
4777    D_delta *= N_e*sqr(r_e)*c0/(32.0*M_PI*cube(gamma)*eps_x_tot*sigma_s*C);
4778    D_x *= N_e*sqr(r_e)*c0/(32.0*M_PI*cube(gamma)*eps_x_tot*sigma_s*C);
4779    eps_IBS[X_] = sqrt(globval.tau[X_]*D_x);
4780    eps_x_tot = eps_SR[X_]*(1.0+sqrt(1.0+4.0*sqr(eps_IBS[X_]/eps_SR[X_])))/2.0;
4781    D_delta /= eps_x_tot; D_x /= eps_x_tot;
4782    sigma_delta_IBS = sqrt(globval.tau[Z_]*D_delta);
4783    eps[X_] = eps_x_tot;
4784
4785    eps[Y_] = eps_SR[Y_]/eps_SR[X_]*eps[X_];
4786    sigma_delta = sqrt(sqr(sigma_delta_SR)+sqr(sigma_delta_IBS));
4787    eps_IBS[Z_] = sigma_s*sigma_delta;
4788    eps[Z_] = eps_SR[Z_] + eps_IBS[Z_];
4789  } else {
4790    D_delta *= N_e*sqr(r_e)*c0/(32.0*M_PI*cube(gamma)*eps[X_]*sigma_s*C);
4791    D_x *= N_e*sqr(r_e)*c0/(32.0*M_PI*cube(gamma)*eps[X_]*sigma_s*C);
4792    eps_IBS[X_] = globval.tau[X_]*D_x;
4793    eps_IBS[Z_] = globval.tau[Z_]*D_delta;
4794    eps[X_] = eps_SR[X_] + eps_IBS[X_];
4795    eps[Y_] = eps_SR[Y_]/eps_SR[X_]*eps[X_];
4796    sigma_delta_IBS = sqrt(globval.tau[Z_]*D_delta);
4797    sigma_delta = sqrt(sqr(sigma_delta_SR)+sqr(sigma_delta_IBS));
4798    // approx.
4799    eps_IBS[Z_] = sigma_s*sigma_delta;
4800    eps[Z_] = eps_SR[Z_] + eps_IBS[Z_];
4801  }
4802
4803  // bunch size
4804  sigma_s = sqrt(beta_z*eps[Z_]);
4805  sigma_delta = sqrt(gamma_z*eps[Z_]);
4806
4807  printf("\n");
4808  printf("D_x,IBS         = %9.3e\n", D_x);
4809  printf("eps_x,IBS       = %5.3f nm.rad\n", 1e9*eps_IBS[X_]);
4810  printf("eps_x,tot       = %5.3f nm.rad\n", 1e9*eps[X_]);
4811  printf("eps_y,tot       = %5.3f nm.rad\n", 1e9*eps[Y_]);
4812  printf("eps_z,tot       = %9.3e\n", eps[Z_]);
4813  printf("D_delta,IBS     = %9.3e\n", D_delta);
4814  printf("sigma_delta,IBS = %9.3e\n", sigma_delta_IBS);
4815  printf("sigma_delta,tot = %9.3e\n", sigma_delta);
4816  printf("sigma_s,tot     = %9.3e\n", sigma_s);
4817}
4818
4819
4820void rm_space(char *name)
4821{
4822  int  i, k;
4823
4824  i = 0;
4825  while (name[i] == ' ')
4826    i++;
4827
4828  for (k = i; k <= (int)strlen(name)+1; k++)
4829    name[k-i] = name[k];
4830}
4831
4832
4833void get_bn(char file_name[], int n, const bool prt)
4834{
4835  char      line[max_str], str[max_str], str1[max_str], *token, *name;
4836  int       n_prm, Fnum, Knum, order;
4837  double    bnL, bn, C, L;
4838  FILE      *inf, *fp_lat;
4839
4840  inf = file_read(file_name); fp_lat = file_write("get_bn.lat");
4841
4842  // if n = 0: go to last data set
4843  if (n == 0) {
4844    while (fgets(line, max_str, inf) != NULL )
4845      if (strstr(line, "n = ") != NULL) sscanf(line, "n = %d", &n);
4846
4847    fclose(inf); inf = file_read(file_name);
4848  }
4849
4850  if (prt) {
4851    printf("\n");
4852    printf("reading values (n=%d): %s\n", n, file_name);
4853    printf("\n");
4854  }
4855
4856  sprintf(str, "n = %d", n);
4857  do
4858    fgets(line, max_str, inf);
4859  while (strstr(line, str) == NULL);
4860
4861  fprintf(fp_lat, "\n");
4862  n_prm = 0;
4863  while (fgets(line, max_str, inf) != NULL) {
4864    if (strcmp(line, "\n") == 0) break;
4865    n_prm++;
4866    name = strtok(line, "(");
4867    rm_space(name);
4868    strcpy(str, name); Fnum = ElemIndex(str);
4869    strcpy(str1, name); upr_case(str1);
4870    token = strtok(NULL, ")"); sscanf(token, "%d", &Knum);
4871    strtok(NULL, "="); token = strtok(NULL, "\n");
4872    sscanf(token, "%lf %d", &bnL, &order);
4873    if (prt) printf("%6s(%2d) = %10.6f %d\n", name, Knum, bnL, order);
4874   if (Fnum != 0) {
4875      if (order == 0)
4876        SetL(Fnum, bnL);
4877      else
4878        SetbnL(Fnum, order, bnL);
4879   
4880      L = GetL(Fnum, 1);
4881      if (Knum == 1)
4882        if (order == 0)
4883          fprintf(fp_lat, "%s: Drift, L = %8.6f;\n", str1, bnL);
4884        else {
4885          bn = (L != 0.0)? bnL/L : bnL;
4886          if (order == Quad)
4887            fprintf(fp_lat, "%s: Quadrupole, L = %8.6f, K = %10.6f, N = Nquad"
4888                    ", Method = Meth;\n", str1, L, bn);
4889          else if (order == Sext)
4890            fprintf(fp_lat, "%s: Sextupole, L = %8.6f, K = %10.6f"
4891                    ", N = 1, Method = Meth;\n", str1, L, bn);
4892          else {
4893            fprintf(fp_lat, "%s: Multipole, L = %8.6f"
4894                    ", N = 1, Method = Meth,\n", str1, L);
4895            fprintf(fp_lat, "     HOM = (%d, %13.6f, %3.1f);\n",
4896                    order, bn, 0.0);
4897          }
4898        }
4899    } else {
4900      printf("element %s not found\n", name);
4901      exit_(1);
4902    }
4903  }
4904  if (prt) printf("\n");
4905
4906  C = Cell[globval.Cell_nLoc].S; recalc_S();
4907  if (prt)
4908    printf("New Cell Length: %5.3f (%5.3f)\n", Cell[globval.Cell_nLoc].S, C);
4909
4910  fclose(inf); fclose(fp_lat);
4911}
4912
4913
4914double get_dynap(const double delta)
4915{
4916  char      str[max_str];
4917  int       i;
4918  double    x_aper[n_aper], y_aper[n_aper], DA;
4919  FILE      *fp;
4920
4921  const int  prt = true;
4922
4923  fp = file_write("dynap.out");
4924  dynap(fp, 5e-3, 0.0, 0.1e-3, n_aper, n_track, x_aper, y_aper, false, prt);
4925  fclose(fp);
4926  DA = get_aper(n_aper, x_aper, y_aper);
4927
4928  if (true) {
4929    sprintf(str, "dynap_dp%3.1f.out", 1e2*delta);
4930    fp = file_write(str);
4931    dynap(fp, 5e-3, delta, 0.1e-3, n_aper, n_track,
4932          x_aper, y_aper, false, prt);
4933    fclose(fp);
4934    DA += get_aper(n_aper, x_aper, y_aper);
4935
4936    for (i = 0; i < nv_; i++)
4937      globval.CODvect[i] = 0.0;
4938    sprintf(str, "dynap_dp%3.1f.out", -1e2*delta);
4939    fp = file_write(str);
4940    dynap(fp, 5e-3, -delta, 0.1e-3, n_aper,
4941          n_track, x_aper, y_aper, false, prt);
4942    fclose(fp);
4943    DA += get_aper(n_aper, x_aper, y_aper);
4944  }
4945
4946  return DA/3.0;
4947}
4948
4949
4950double get_chi2(long int n, double x[], double y[], long int m, Vector b)
4951{
4952  /* Compute chi2 for polynomial fit */
4953
4954  int     i, j;
4955  double  sum, z;
4956 
4957  sum = 0.0;
4958  for (i = 0; i < n; i++) {
4959    z = 0.0;
4960    for (j = 0; j <= m; j++)
4961      z += b[j]*pow(x[i], j);
4962    sum += pow(y[i]-z, 2);
4963  }
4964  return(sum/n);
4965}
4966
4967
4968void pol_fit(int n, double x[], double y[], int order, Vector &b,
4969             double &sigma)
4970{
4971  /* Polynomial fit by linear chi-square */
4972
4973  int     i, j, k;
4974  Matrix  T1;
4975  bool   prt = false;
4976 
4977  const double sigma_k = 1.0, chi2 = 4.0;
4978
4979  /* initialization */
4980  for (i = 0; i <= order; i++) {
4981    b[i] = 0.0;
4982    for (j = 0; j <= order; j++)
4983      T1[i][j] = 0.0;
4984  }
4985
4986  /* polynomal fiting */
4987  for (i = 0; i < n; i++)
4988    for (j = 0; j <= order; j++) {
4989      b[j] += y[i]*pow(x[i], j)/pow(sigma_k, 2);
4990      for (k = 0; k <= order; k++)
4991        T1[j][k] += pow(x[i], j+k)/pow(sigma_k, 2);
4992    }
4993
4994  if (InvMat(order+1, T1)) {
4995    LinTrans(order+1, T1, b); sigma = get_chi2(n, x, y, order, b);
4996    if(prt){
4997      printf("\n");
4998      printf("  n    Coeff.\n");
4999      for (i = 0; i <= order; i++)
5000        printf("%3d %10.3e+/-%8.2e\n", i, b[i], sigma*sqrt(chi2*T1[i][i]));
5001    }   
5002  } else {
5003    printf("pol_fit: Matrix is singular\n");
5004    exit_(1);
5005  }
5006}
5007
5008
5009void get_ksi2(const double d_delta)
5010{
5011  const int     n_points = 20, order = 5;
5012
5013  int       i, n;
5014  double    delta[2*n_points+1], nu[2][2*n_points+1], sigma;
5015  Vector    b;
5016  FILE      *fp;
5017 
5018  fp = file_write("chrom2.out");
5019  n = 0;
5020  for (i = -n_points; i <= n_points; i++) {
5021    n++; delta[n-1] = i*(double)d_delta/(double)n_points;
5022    Ring_GetTwiss(false, delta[n-1]);
5023    nu[0][n-1] = globval.TotalTune[X_]; nu[1][n-1] = globval.TotalTune[Y_];
5024    fprintf(fp, "%5.2f %8.5f %8.5f\n", 1e2*delta[n-1], nu[0][n-1], nu[1][n-1]);
5025  }
5026  printf("\n");
5027  printf("Horizontal chromaticity:\n");
5028  pol_fit(n, delta, nu[0], order, b, sigma);
5029  printf("\n");
5030  printf("Vertical chromaticity:\n");
5031  pol_fit(n, delta, nu[1], order, b, sigma);
5032  printf("\n");
5033  fclose(fp);
5034}
5035
5036
5037bool find_nu(const int n, const double nus[], const double eps, double &nu)
5038{
5039  bool  lost;
5040  int   k;
5041
5042  const bool  prt = false;
5043
5044  if (prt) printf("\n");
5045  k = 0;
5046  while ((k < n) && (fabs(nus[k]-nu) > eps)) {
5047    if (prt) printf("nu = %7.5f(%7.5f) eps = %7.1e\n", nus[k], nu, eps);
5048    k++;
5049  }
5050
5051  if (k < n) {
5052    if (prt) printf("nu = %7.5f(%7.5f)\n", nus[k], nu);
5053    nu = nus[k]; lost = false;
5054  } else {
5055    if (prt) printf("lost\n");
5056    lost = true;
5057  }
5058
5059  return !lost;
5060}
5061
5062
5063bool get_nu(const double Ax, const double Ay, const double delta,
5064            const double eps, double &nu_x, double &nu_y)
5065{
5066  const int  n_turn = 512, n_peaks = 5;
5067
5068  bool      lost, ok_x, ok_y;
5069  char      str[max_str];
5070  int       i;
5071  long int  lastpos, lastn, n;
5072  double    x[n_turn], px[n_turn], y[n_turn], py[n_turn];
5073  double    nu[2][n_peaks], A[2][n_peaks];
5074  Vector    x0;
5075  FILE     *fp;
5076
5077  const bool   prt = false;
5078  const char   file_name[] = "track.out";
5079
5080  // complex FFT in Floquet space
5081  x0[x_] = Ax; x0[px_] = 0.0; x0[y_] = Ay; x0[py_] = 0.0;
5082  LinTrans(4, globval.Ascrinv, x0);
5083  track(file_name, x0[x_], x0[px_], x0[y_], x0[py_], delta,
5084        n_turn, lastn, lastpos, 1, 0.0);
5085  if (lastn == n_turn) {
5086    GetTrack(file_name, &n, x, px, y, py);
5087    sin_FFT((int)n, x, px); sin_FFT((int)n, y, py);
5088
5089    if (prt) {
5090      strcpy(str, file_name); strcat(str, ".fft");
5091      fp = file_write(str);
5092      for (i = 0; i < n; i++)
5093        fprintf(fp, "%5.3f %12.5e %8.5f %12.5e %8.5f\n",
5094                (double)i/(double)n, x[i], px[i], y[i], py[i]);
5095      fclose(fp);
5096    }
5097   
5098    GetPeaks1(n, x, n_peaks, nu[X_], A[X_]);
5099    GetPeaks1(n, y, n_peaks, nu[Y_], A[Y_]);
5100
5101    ok_x = find_nu(n_peaks, nu[X_], eps, nu_x);
5102    ok_y = find_nu(n_peaks, nu[Y_], eps, nu_y);
5103
5104    lost = !ok_x || !ok_y;
5105  } else
5106    lost = true;
5107 
5108  return !lost;
5109}
5110
5111
5112// tune shift with amplitude
5113void dnu_dA(const double Ax_max, const double Ay_max, const double delta)
5114{
5115  bool      ok;
5116  int       i;
5117  double    nu_x, nu_y, Ax, Ay, Jx, Jy;
5118  Vector    ps;
5119  FILE      *fp;
5120
5121  const int     n_ampl = 25;
5122  const double  A_min  = 1e-3;
5123//  const double  eps0   = 0.04, eps   = 0.02;
5124//  const double  eps0   = 0.025, eps   = 0.02;
5125  const double  eps0   = 0.04, eps   = 0.015;
5126
5127  Ring_GetTwiss(false, 0.0);
5128
5129  ok = false;
5130  fp = file_write("dnu_dAx.out");
5131  fprintf(fp, "#   x[mm]        y[mm]        J_x        J_y      nu_x    nu_y\n");
5132  fprintf(fp, "#\n");
5133  for (i = -n_ampl; i <= n_ampl; i++) {
5134    Ax = i*Ax_max/n_ampl;
5135    if (Ax == 0.0) Ax = A_min;
5136    Ay = A_min;
5137    ps[x_] = Ax; ps[px_] = 0.0; ps[y_] = Ay; ps[py_] = 0.0; getfloqs(ps);
5138    Jx = (sqr(ps[x_])+sqr(ps[px_]))/2.0; Jy = (sqr(ps[y_])+sqr(ps[py_]))/2.0;
5139    if (!ok) {
5140      nu_x = fract(globval.TotalTune[X_]); nu_y = fract(globval.TotalTune[Y_]);
5141      ok = get_nu(Ax, Ay, delta, eps0, nu_x, nu_y);
5142    } else
5143      ok = get_nu(Ax, Ay, delta, eps, nu_x, nu_y);
5144    if (ok)
5145      fprintf(fp, "%10.3e %10.3e %10.3e %10.3e %7.5f %7.5f\n",
5146              1e3*Ax, 1e3*Ay, 1e6*Jx, 1e6*Jy, fract(nu_x), fract(nu_y));
5147    else
5148      fprintf(fp, "# %10.3e %10.3e particle lost\n", 1e3*Ax, 1e3*Ay);
5149  }
5150  fclose(fp);
5151
5152  ok = false;
5153  fp = file_write("dnu_dAy.out");
5154  fprintf(fp, "#   x[mm]        y[mm]      J_x        J_y      nu_x    nu_y\n");
5155  fprintf(fp, "#\n");
5156  for (i = -n_ampl; i <= n_ampl; i++) {
5157    Ax = A_min; Ay = i*Ay_max/n_ampl;
5158    Jx = pow(Ax, 2)/(2.0*Cell[globval.Cell_nLoc].Beta[X_]);
5159    if (Ay == 0.0) Ay = A_min;
5160    Jy = pow(Ay, 2)/(2.0*Cell[globval.Cell_nLoc].Beta[Y_]);
5161//    if (i == 1) exit_(0);
5162    if (!ok) {
5163      nu_x = fract(globval.TotalTune[X_]); nu_y = fract(globval.TotalTune[Y_]);
5164      ok = get_nu(Ax, Ay, delta, eps0, nu_x, nu_y);
5165    } else
5166      ok = get_nu(Ax, Ay, delta, eps, nu_x, nu_y);
5167    if (ok)
5168      fprintf(fp, "%10.3e %10.3e %10.3e %10.3e %7.5f %7.5f\n",
5169              1e3*Ax, 1e3*Ay, 1e6*Jx, 1e6*Jy, fract(nu_x), fract(nu_y));
5170    else
5171      fprintf(fp, "# %10.3e %10.3e particle lost\n", 1e3*Ax, 1e3*Ay);
5172  }
5173  fclose(fp);
5174}
5175
5176
5177bool orb_corr(const int n_orbit)
5178{
5179  bool      cod = false;
5180  int       i;
5181  long      lastpos;
5182  Vector2   xmean, xsigma, xmax;
5183
5184  printf("\n");
5185
5186//  FitTune(qf, qd, nu_x, nu_y);
5187//  printf("\n");
5188//  printf("  qf = %8.5f qd = %8.5f\n",
5189//         GetKpar(qf, 1, Quad), GetKpar(qd, 1, Quad));
5190
5191  globval.CODvect.zero();
5192  for (i = 1; i <= n_orbit; i++) {
5193    cod = getcod(0.0, lastpos);
5194    if (cod) {
5195      codstat(xmean, xsigma, xmax, globval.Cell_nLoc, false);
5196      printf("\n");
5197      printf("RMS orbit [mm]: (%8.1e+/-%7.1e, %8.1e+/-%7.1e)\n", 
5198             1e3*xmean[X_], 1e3*xsigma[X_], 1e3*xmean[Y_], 1e3*xsigma[Y_]);
5199      lsoc(1, globval.bpm, globval.hcorr, 1);
5200      lsoc(1, globval.bpm, globval.vcorr, 2);
5201      cod = getcod(0.0, lastpos);
5202      if (cod) {
5203        codstat(xmean, xsigma, xmax, globval.Cell_nLoc, false);
5204        printf("RMS orbit [mm]: (%8.1e+/-%7.1e, %8.1e+/-%7.1e)\n", 
5205               1e3*xmean[X_], 1e3*xsigma[X_], 1e3*xmean[Y_], 1e3*xsigma[Y_]);
5206      } else
5207        printf("orb_corr: failed\n");
5208    } else
5209      printf("orb_corr: failed\n");
5210  }
5211 
5212  prt_cod("orb_corr.out", globval.bpm, true);
5213
5214  return cod;
5215}
5216
5217
5218double get_code(CellType &Cell)
5219{
5220  double  code = 0.0;
5221
5222  switch (Cell.Elem.Pkind) {
5223  case drift:
5224    code = 0.0;
5225    break;
5226  case Mpole:
5227    if (Cell.Elem.M->Pirho != 0.0) {
5228      code = 0.5;
5229    } else if (Cell.Elem.M->PBpar[Quad+HOMmax] != 0)
5230      code = sgn(Cell.Elem.M->PBpar[Quad+HOMmax]);
5231    else if (Cell.Elem.M->PBpar[Sext+HOMmax] != 0)
5232      code = 1.5*sgn(Cell.Elem.M->PBpar[Sext+HOMmax]);
5233    else if (Cell.Fnum == globval.bpm)
5234      code = 2.0;
5235    else
5236      code = 0.0;
5237    break;
5238  default:
5239    code = 0.0;
5240    break;
5241  }
5242
5243  return code;
5244}
5245
5246
5247void get_alphac(void)
5248{
5249  CellType  Cell;
5250
5251  getelem(globval.Cell_nLoc, &Cell);
5252  globval.Alphac = globval.OneTurnMat[ct_][delta_]/Cell.S;
5253}
5254
5255/******************************************************************
5256void get_alphac2(void)
5257
5258  Purpose:
5259           Calculate momentum compact factor up to third order.
5260           Method:
5261           track the orbit offset c*tau at the first element,
5262           at 11 different energy offset(-10-3 to 10-3), then
5263           use polynomal to fit the momentum compact faction factor
5264           up to 3rd order.
5265           The initial coorinates are (x_co,px_co,y_co,py_co,delta,0).
5266           
5267 Input:
5268       none
5269
5270   Output:
5271       none
5272
5273   Return:
5274      alphac
5275
5276   specific functions:
5277      none
5278
5279   Comments:
5280       change the initial coorinates from (0 0 0 0 delta 0) to
5281              (x_co,px_co,y_co,py_co,delta,0).  i.e., track
5282              around the off-momentum close orbit but not zero orbit.     
5283
5284*******************************************************************/
5285
5286void get_alphac2(void)
5287{
5288  /* Note, do not extract from M[5][4], i.e. around delta
5289     dependent fixed point.                                */
5290
5291  const int     n_points = 5;
5292  const double  d_delta  = 1e-3;
5293 //const double  d_delta  = 1e-4;
5294
5295  int       i, n;
5296  long int  lastpos;  // last tracking position when the particle is stable
5297  double    delta[2*n_points+1], alphac[2*n_points+1], sigma;
5298  Vector    x, b;
5299  CellType  Cell;
5300  bool      cod;
5301 
5302  globval.pathlength = false;
5303  getelem(globval.Cell_nLoc, &Cell); 
5304  n = 0;
5305 
5306  for (i = -n_points; i <= n_points; i++) {
5307    n++; 
5308    delta[n-1] = i*(double)d_delta/(double)n_points;
5309   
5310   // for (j = 0; j < nv_; j++)
5311    //  x[j] = 0.0;
5312    //   x[delta_] = delta[n-1];
5313   
5314    cod = getcod(delta[n-1], lastpos);   
5315    x =  globval.CODvect;
5316     
5317   
5318      Cell_Pass(0, globval.Cell_nLoc, x, lastpos);
5319      alphac[n-1] = x[ct_]/Cell.S; // this is not mcf but DT/T
5320  }
5321  pol_fit(n, delta, alphac, 3, b, sigma);
5322  printf("Momentum compaction factor:\n");
5323  printf("   alpha1 = %10.4e  alpha2 = %10.4e  alpha3 = %10.4e\n",
5324         b[1], b[2], b[3]);
5325}
5326
5327
5328float f_bend(float b0L[])
5329{
5330  long int  lastpos;
5331  Vector    ps;
5332
5333  const int   n_prt = 10;
5334
5335  n_iter_Cart++;
5336
5337  SetbnL_sys(Fnum_Cart, Dip, b0L[1]);
5338
5339  ps.zero();
5340  Cell_Pass(Elem_GetPos(Fnum_Cart, 1)-1, Elem_GetPos(Fnum_Cart, 1),
5341            ps, lastpos);
5342
5343  if (n_iter_Cart % n_prt == 0)
5344    cout << scientific << setprecision(3)
5345         << setw(4) << n_iter_Cart
5346         << setw(11) << ps[x_] << setw(11) << ps[px_] << setw(11) << ps[ct_]
5347         << setw(11) << b0L[1] << endl;
5348
5349  return sqr(ps[x_]) + sqr(ps[px_]);
5350}
5351
5352
5353void bend_cal_Fam(const int Fnum)
5354{
5355  /* Adjusts b1L_sys to zero the orbit for a given gradient. */
5356  const int  n_prm = 1;
5357
5358  int       iter;
5359  float     *b0L, **xi, fret;
5360  Vector    ps;
5361
5362  const float  ftol = 1e-15;
5363
5364  b0L = vector(1, n_prm); xi = matrix(1, n_prm, 1, n_prm);
5365
5366  cout << endl;
5367  cout << "bend_cal: " << ElemFam[Fnum-1].ElemF.PName << ":" << endl;
5368
5369  Fnum_Cart = Fnum;  b0L[1] = 0.0; xi[1][1] = 1e-3;
5370
5371  cout << endl;
5372  n_iter_Cart = 0;
5373  powell(b0L, xi, n_prm, ftol, &iter, &fret, f_bend);
5374
5375  free_vector(b0L, 1, n_prm); free_matrix(xi, 1, n_prm, 1, n_prm);
5376}
5377
5378
5379void bend_cal(void)
5380{
5381  long int  k;
5382
5383  for (k = 1; k <= globval.Elem_nFam; k++)
5384    if ((ElemFam[k-1].ElemF.Pkind == Mpole) &&
5385        (ElemFam[k-1].ElemF.M->Pirho != 0.0) &&
5386        (ElemFam[k-1].ElemF.M->PBpar[Quad+HOMmax] != 0.0))
5387      if (ElemFam[k-1].nKid > 0) bend_cal_Fam(k);
5388}
5389
5390
5391double h_ijklm(const tps &h, const int i, const int j, const int k,
5392               const int l, const int m)
5393{
5394  int      i1;
5395  iVector  jj;
5396
5397  for (i1 = 0; i1 < nv_tps; i1++)
5398    jj[i1] = 0;
5399  jj[x_] = i; jj[px_] = j; jj[y_] = k; jj[py_] = l; jj[delta_] = m;
5400  return h[jj];
5401}
5402
5403
5404ss_vect<tps> get_A(const double alpha[], const double beta[],
5405                   const double eta[], const double etap[])
5406{
5407  ss_vect<tps>  A, Id;
5408
5409  Id.identity();
5410
5411  A.identity();
5412  A[x_]  = sqrt(beta[X_])*Id[x_];
5413  A[px_] = -alpha[X_]/sqrt(beta[X_])*Id[x_] + 1.0/sqrt(beta[X_])*Id[px_];
5414  A[y_]  = sqrt(beta[Y_])*Id[y_];
5415  A[py_] = -alpha[Y_]/sqrt(beta[Y_])*Id[y_] + 1.0/sqrt(beta[Y_])*Id[py_];
5416
5417  A[x_] += eta[X_]*Id[delta_]; A[px_] += etap[X_]*Id[delta_];
5418
5419  return A;
5420}
5421
5422
5423void get_ab(const ss_vect<tps> &A,
5424            double alpha[], double beta[], double eta[], double etap[])
5425{
5426  int           k;
5427  ss_vect<tps>  A_Atp;
5428
5429  A_Atp = A*tp_S(2, A);
5430
5431  for (k = 0; k <= 1; k++) {
5432      eta[k] = A[2*k][delta_]; 
5433     etap[k] = A[2*k+1][delta_];
5434    alpha[k] = -A_Atp[2*k][2*k+1]; 
5435     beta[k] = A_Atp[2*k][2*k];
5436  }
5437}
5438
5439
5440void set_tune(const char file_name1[], const char file_name2[], const int n)
5441{
5442  const int  n_b2 = 8;
5443
5444  char      line[max_str], names[n_b2][max_str];
5445  int       j, k, Fnum;
5446  double    b2s[n_b2], nu[2];
5447  ifstream  inf1, inf2;
5448  FILE      *fp_lat;
5449
5450  file_rd(inf1, file_name1); file_rd(inf2, file_name2);
5451
5452  // skip 1st and 2nd line
5453  inf1.getline(line, max_str); inf2.getline(line, max_str);
5454  inf1.getline(line, max_str); inf2.getline(line, max_str);
5455 
5456  inf1.getline(line, max_str);
5457  sscanf(line, "%*s %*s %*s %s %s %s %s",
5458         names[0], names[1], names[2], names[3]);
5459  inf2.getline(line, max_str);
5460  sscanf(line, "%*s %*s %*s %s %s %s %s",
5461         names[4], names[5], names[6], names[7]);
5462
5463  printf("\n");
5464  printf("quads:");
5465  for (k = 0; k <  n_b2; k++) {
5466    lwr_case(names[k]); rm_space(names[k]);
5467    printf(" %s", names[k]);
5468  }
5469  printf("\n");
5470
5471  // skip 4th line
5472  inf1.getline(line, max_str); inf2.getline(line, max_str);
5473
5474  fp_lat = file_write("set_tune.lat");
5475  while (inf1.getline(line, max_str) != NULL) {
5476    sscanf(line, "%d%lf %lf %lf %lf %lf %lf",
5477           &j, &nu[X_], &nu[Y_], &b2s[0], &b2s[1], &b2s[2], &b2s[3]);
5478    inf2.getline(line, max_str);
5479    sscanf(line, "%d%lf %lf %lf %lf %lf %lf",
5480           &j, &nu[X_], &nu[Y_], &b2s[4], &b2s[5], &b2s[6], &b2s[7]);
5481
5482    if (j == n) {
5483      printf("\n");
5484      printf("%3d %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
5485             n,
5486             b2s[0], b2s[1], b2s[2], b2s[3], b2s[4], b2s[5], b2s[6], b2s[7]);
5487
5488      for (k = 0; k <  n_b2; k++) {
5489        Fnum = ElemIndex(names[k]);
5490        set_bn_design_fam(Fnum, Quad, b2s[k], 0.0);
5491
5492        fprintf(fp_lat, "%s: Quadrupole, L = %8.6f, K = %10.6f, N = Nquad"
5493                ", Method = Meth;\n",
5494                names[k], ElemFam[Fnum-1].ElemF.PL, b2s[k]);
5495      }
5496      break;
5497    }
5498  }
5499
5500  fclose(fp_lat);
5501}
5502
5503
5504
5505
5506
Note: See TracBrowser for help on using the repository browser.