FLAME  devel
 All Classes Functions Variables Typedefs Enumerations Pages
rf_cavity.cpp
1 
2 #include <fstream>
3 
4 #include <boost/lexical_cast.hpp>
5 #include <boost/filesystem.hpp>
6 
7 #include "flame/constants.h"
8 #include "flame/moment.h"
9 #include "flame/moment_sup.h"
10 #include "flame/rf_cavity.h"
11 #include "flame/config.h"
12 
13 #define sqr(x) ((x)*(x))
14 #define cube(x) ((x)*(x)*(x))
15 
16 #ifdef DEFPATH
17  #define defpath DEFPATH
18 #else
19  #define defpath "."
20 #endif
21 
22 std::map<std::string,boost::shared_ptr<Config> > ElementRFCavity::CavConfMap;
23 
24 // RF Cavity beam dynamics functions.
25 
26 static double ipow(double base, int exp)
27 {
28  double result = 1.0;
29  while (exp)
30  {
31  if (exp & 1)
32  result *= base;
33  exp >>= 1;
34  base *= base;
35  }
36 
37  return result;
38 }
39 
40 
41 static
42 // Evaluate the beam energy and phase in the acceleration gap.
43 void EvalGapModel(const double dis, const double IonW0, const Particle &real, const double IonFy0,
44  const double k, const double Lambda, const double Ecen,
45  const double T, const double S, const double Tp, const double Sp, const double V0,
46  double &IonW_f, double &IonFy_f);
47 
48 
49 static
50 // Calculate driven phase from synchronous phase which defined by sinusoidal fitting.
51 double GetCavPhase(const int cavi, const Particle& ref, const double IonFys, const double multip, const std::vector<double>& P);
52 
53 
54 static
55 // Calculate driven phase from synchronous phase which defined by complex fitting (e.g. peak-base model).
56 double GetCavPhaseComplex(const Particle& ref, const double IonFys, const double scale,
57  const double multip, const std::vector<double>& P);
58 
59 
60 void CavDataType::show(std::ostream& strm, const int k) const
61 {
62  strm << std::scientific << std::setprecision(5)
63  << std::setw(13) << s[k] << std::setw(13) << Elong[k] << "\n";
64 }
65 
66 
67 void CavDataType::show(std::ostream& strm) const
68 {
69  for (unsigned int k = 0; k < s.size(); k++)
70  show(strm, k);
71 }
72 
73 
74 void CavTLMLineType::clear(void)
75 {
76  s.clear(); Elem.clear(); E0.clear();
77  T.clear(); S.clear(); Accel.clear();
78 }
79 
80 
81 void CavTLMLineType::set(const double s, const std::string &Elem, const double E0,
82  const double T, const double S, const double Accel)
83 {
84  this->s.push_back(s); this->Elem.push_back(Elem); this->E0.push_back(E0);
85  this->T.push_back(T); this->S.push_back(S); this->Accel.push_back(Accel);
86 }
87 
88 
89 void CavTLMLineType::show(const int k) const
90 {
91  FLAME_LOG(FINE) << std::fixed << std::setprecision(5)
92  << std::setw(9) << s[k] << std::setw(10) << Elem[k]
93  << std::scientific << std::setprecision(10)
94  << std::setw(18) << T[k] << std::setw(18) << S[k]
95  << std::setw(18) << Accel[k] << "\n";
96 }
97 
98 
99 void CavTLMLineType::show() const
100 {
101  for (unsigned int k = 0; k < s.size(); k++)
102  show(k);
103 }
104 
105 
106 int get_column(const std::string &str)
107 {
108  int column;
109 
110  if (str == "CaviMlp_EFocus1")
111  column = 3;
112  else if (str == "CaviMlp_EFocus2")
113  column = 4;
114  else if (str == "CaviMlp_EDipole")
115  column = 2;
116  else if (str == "CaviMlp_EQuad")
117  column = 5;
118  else if (str == "CaviMlp_HMono")
119  column = 7;
120  else if (str == "CaviMlp_HDipole")
121  column = 6;
122  else if (str == "CaviMlp_HQuad")
123  column = 8;
124  else {
125  throw std::runtime_error(SB()<<"get_column: undef. column: " << str);
126  }
127 
128  return column;
129 }
130 
131 
132 void calTransfac(const numeric_table& fldmap, int column_no, const int gaplabel, const double IonK, const bool half,
133  double &Ecenter, double &T, double &Tp, double &S, double &Sp, double &V0)
134 {
135  // Compute electric center, amplitude, and transit time factors [T, Tp, S, Sp] for RF cavity mode.
136  int n, k;
137  double len, dz, eml, em_mom;
138  std::vector<double> z, EM;
139 
140  column_no--; // F*** fortran
141  assert(column_no>0);
142 
143  n = fldmap.table.size1();
144 
145  if(n<=0 || (size_t)column_no>=fldmap.table.size2())
146  throw std::runtime_error("field map size invalid");
147 
148  z.resize(n);
149  EM.resize(n);
150 
151  std::copy(fldmap.table.find1(2, 0, 0),
152  fldmap.table.find1(2, n, 0),
153  z.begin());
154  std::copy(fldmap.table.find1(2, 0, column_no),
155  fldmap.table.find1(2, n, column_no),
156  EM.begin());
157 
158  // Used at end of function.
159  len = z[n-1];
160 
161  if (half) n = (int)round((n-1)/2e0);
162 
163  dz = (z[n-1]-z[0])/(n-1);
164 
165 // prt_data(z, EM);
166 
167  // Start at zero.
168  for (k = n-1; k >= 0; k--)
169  z[k] -= z[0];
170 
171  eml = 0e0, em_mom = 0e0;
172  for (k = 0; k < n-1; k++) {
173  eml += (fabs(EM[k])+fabs(EM[k+1]))/2e0*dz;
174  em_mom += (z[k]+z[k+1])/2e0*(fabs(EM[k])+fabs(EM[k+1]))/2e0*dz;
175  }
176  Ecenter = em_mom/eml;
177 
178  for (k = 0; k < n; k++)
179  z[k] -= Ecenter;
180 
181  T = 0;
182  for (k = 0; k < n-1; k++)
183  T += (EM[k]+EM[k+1])/2e0*cos(IonK*(z[k]+z[k+1])/2e0)*dz;
184  T /= eml;
185 
186  Tp = 0;
187  for (k = 0; k < n-1; k++)
188  Tp -= ((z[k]+z[k+1])/2e0)*(EM[k]+EM[k+1])/2e0*sin(IonK*(z[k]+z[k+1])/2e0)*dz;
189  Tp /= eml;
190 
191  S = 0e0;
192  for (k = 0; k < n-1; k++)
193  S += (EM[k]+EM[k+1])/2e0*sin(IonK*(z[k]+z[k+1])/2e0)*(z[k+1]-z[k]);
194  S /= eml;
195 
196  Sp = 0e0;
197  for (k = 0; k < n-1; k++) {
198  Sp += (z[k]+z[k+1])/2e0*(EM[k]+EM[k+1])/2e0*cos(IonK*(z[k]+z[k+1])/2e0)*dz;
199  }
200  Sp /= eml;
201 
202  V0 = eml/MeVtoeV/MtoMM;
203 
204  if (gaplabel == 2) {
205  // Second gap.
206  Ecenter = len - Ecenter;
207  T = -T;
208  Tp = -Tp;
209  }
210 }
211 
212 
213 double PwrSeries(const double beta,
214  const double a0, const double a1, const double a2, const double a3,
215  const double a4, const double a5, const double a6, const double a7,
216  const double a8, const double a9)
217 {
218  int k;
219  double f;
220 
221  const int n = 10;
222  const double a[] = {a0, a1, a2, a3, a4, a5, a6, a7, a8, a9};
223 
224  f = a[0];
225  for (k = 1; k < n; k++)
226  f += a[k]*pow(beta, k);
227 
228  return f;
229 }
230 
231 
232 void ElementRFCavity::TransFacts(const int cavilabel, double beta, const double CaviIonK, const int gaplabel, const double EfieldScl,
233  double &Ecen, double &T, double &Tp, double &S, double &Sp, double &V0) const
234 {
235  // Evaluate Electric field center, transit factors [T, T', S, S'] and cavity field.
236  std::ostringstream strm;
237 
238  // For debugging of TTF function.
239  if (forcettfcalc) {
240  calTransfac(CavData, 2, gaplabel, CaviIonK, true, Ecen, T, Tp, S, Sp, V0);
241  V0 *= EfieldScl;
242  return;
243  }
244 
245  switch (cavilabel) {
246  case 41:
247  if (beta < 0.025 || beta > 0.08) {
248  FLAME_LOG(DEBUG) << "*** TransFacts: CaviIonK out of Range " << cavilabel << "\n";
249  calTransfac(CavData, 2, gaplabel, CaviIonK, true, Ecen, T, Tp, S, Sp, V0);
250  V0 *= EfieldScl;
251  return;
252  }
253  switch (gaplabel) {
254  case 0:
255  // One gap evaluation.
256  Ecen = 120.0; // [mm].
257  T = 0.0;
258  Tp = 0.0;
259  S = PwrSeries(beta, -4.109, 399.9, -1.269e4, 1.991e5, -1.569e6, 4.957e6, 0.0, 0.0, 0.0, 0.0);
260  Sp = PwrSeries(beta, 61.98, -1.073e4, 4.841e5, 9.284e6, 8.379e7, -2.926e8, 0.0, 0.0, 0.0, 0.0);
261  V0 = 0.98477*EfieldScl;
262  break;
263  case 1:
264  // Two gap calculation, first gap.
265  Ecen = 0.0006384*pow(beta, -1.884) + 86.69;
266  T = PwrSeries(beta, 0.9232, -123.2, 3570, -5.476e4, 4.316e5, -1.377e6, 0.0, 0.0, 0.0, 0.0);
267  Tp = PwrSeries(beta, 1.699, 924.7, -4.062e4, 7.528e5, -6.631e6, 2.277e7, 0.0, 0.0, 0.0, 0.0);
268  S = 0.0;
269  Sp = PwrSeries(beta, -1.571, 25.59, 806.6, -2.98e4, 3.385e5, -1.335e6, 0.0, 0.0, 0.0, 0.0);
270  V0 = 0.492385*EfieldScl;
271  break;
272  case 2:
273  // Two gap calculation, second gap.
274  Ecen = -0.0006384*pow(beta, -1.884) + 33.31;
275  T = PwrSeries(beta, -0.9232, 123.2, -3570, 5.476e4, -4.316e5, 1.377e6, 0.0, 0.0, 0.0, 0.0);
276  Tp = PwrSeries(beta, -1.699, -924.7, 4.062e4, -7.528e5, 6.631e6, -2.277e7, 0.0, 0.0, 0.0, 0.0);
277  S = 0.0;
278  Sp = PwrSeries(beta, -1.571, 25.59, 806.6, -2.98e4, 3.385e5, -1.335e6, 0.0, 0.0, 0.0, 0.0);
279  V0 = 0.492385*EfieldScl;
280  break;
281  default:
282  strm << "*** GetTransitFac: undef. number of gaps " << gaplabel << "\n";
283  throw std::runtime_error(strm.str());
284  }
285  break;
286  case 85:
287  if (beta < 0.05 || beta > 0.25) {
288  FLAME_LOG(DEBUG) << "*** TransFacts: CaviIonK out of Range " << cavilabel << "\n";
289  calTransfac(CavData, 2, gaplabel, CaviIonK, true, Ecen, T, Tp, S, Sp, V0);
290  V0 *= EfieldScl;
291  return;
292  }
293  switch (gaplabel) {
294  case 0:
295  Ecen = 150.0; // [mm].
296  T = 0.0;
297  Tp = 0.0;
298  S = PwrSeries(beta, -6.811, 343.9, -6385, 6.477e4, -3.914e5, 1.407e6, -2.781e6, 2.326e6, 0.0, 0.0);
299  Sp = PwrSeries(beta, 162.7, -1.631e4, 4.315e5, -5.344e6, 3.691e7, -1.462e8, 3.109e8, -2.755e8, 0.0, 0.0);
300  V0 = 1.967715*EfieldScl;
301  break;
302  case 1:
303  Ecen = 0.0002838*pow(beta, -2.13) + 76.5;
304  T = 0.0009467*pow(beta, -1.855) - 1.002;
305  Tp = PwrSeries(beta, 24.44, -334, 2468, -1.017e4, 2.195e4, -1.928e4, 0.0, 0.0, 0.0, 0.0);
306  S = 0.0;
307  Sp = -0.0009751*pow(beta, -1.898) + 0.001568;
308  V0 = 0.9838574*EfieldScl;
309  break;
310  case 2:
311  Ecen = -0.0002838*pow(beta, -2.13) + 73.5;
312  T = -0.0009467*pow(beta, -1.855) + 1.002;
313  Tp = PwrSeries(beta, -24.44, 334, -2468, 1.017e4, -2.195e4, 1.928e4, 0.0, 0.0, 0.0, 0.0);
314  S = 0.0;
315  Sp = -0.0009751*pow(beta, -1.898) + 0.001568;
316  V0 = 0.9838574*EfieldScl;
317  break;
318  default:
319  strm << "*** GetTransitFac: undef. number of gaps " << gaplabel << "\n";
320  throw std::runtime_error(strm.str());
321  }
322  break;
323  case 29:
324  if (beta < 0.15 || beta > 0.4) {
325  FLAME_LOG(DEBUG) << "*** TransFacts: CaviIonK out of Range " << cavilabel << "\n";
326  calTransfac(CavData, 2, gaplabel, CaviIonK, true, Ecen, T, Tp, S, Sp, V0);
327  V0 *= EfieldScl;
328  return;
329  }
330  switch (gaplabel) {
331  case 0:
332  Ecen = 150.0; // [mm].
333  T = 0.0;
334  Tp = 0.0;
335  S = PwrSeries(beta, -4.285, 58.08, -248, 486, -405.6, 76.54, 0.0, 0.0, 0.0, 0.0);
336  Sp = PwrSeries(beta, 888, -2.043e4, 1.854e5, -9.127e5, 2.695e6, -4.791e6, 4.751e6, -2.025e6, 0.0, 0.0);
337  V0 = 2.485036*EfieldScl;
338  break;
339  case 1:
340  Ecen = 0.01163*pow(beta, -2.001) + 91.77;
341  T = 0.02166*pow(beta, -1.618) - 1.022;
342  Tp = PwrSeries(beta, -11.25, 534.7, -3917, 1.313e4, -2.147e4, 1.389e4, 0.0, 0.0, 0.0, 0.0);
343  S = 0.0;
344  Sp = PwrSeries(beta, -0.8283, -4.409, 78.77, -343.9, 645.1, -454.4, 0.0, 0.0, 0.0, 0.0);
345  V0 = 1.242518*EfieldScl;
346  break;
347  case 2:
348  Ecen =-0.01163*pow(beta, -2.001) + 58.23;
349  T =-0.02166*pow(beta, -1.618) + 1.022;
350  Tp = PwrSeries(beta, 11.25, -534.7, 3917, -1.313e4, 2.147e4, -1.389e4, 0.0, 0.0, 0.0, 0.0);
351  S = 0.0;
352  Sp = PwrSeries(beta, -0.8283, -4.409, 78.77, -343.9, 645.1, -454.4, 0.0, 0.0, 0.0, 0.0);
353  V0 = 1.242518*EfieldScl;
354  break;
355  default:
356  strm << "*** GetTransitFac: undef. number of gaps " << gaplabel << "\n";
357  throw std::runtime_error(strm.str());
358  }
359  break;
360  case 53:
361  if (beta < 0.3 || beta > 0.6) {
362  FLAME_LOG(DEBUG) << "*** TransFacts: CaviIonK out of Range " << cavilabel << "\n";
363  calTransfac(CavData, 2, gaplabel, CaviIonK, true, Ecen, T, Tp, S, Sp, V0);
364  V0 *= EfieldScl;
365  return;
366  }
367  switch (gaplabel) {
368  case 0:
369  Ecen = 250.0; // [mm].
370  T = 0.0;
371  Tp = 0.0;
372  S = PwrSeries(beta, -4.222, 26.64, -38.49, -17.73, 84.12, -52.93, 0.0, 0.0, 0.0, 0.0);
373  Sp = PwrSeries(beta, -1261, -1.413e4, 5.702e4, -1.111e5, 1.075e5, -4.167e4, 0.0, 0.0 , 0.0, 0.0);
374  V0 = 4.25756986*EfieldScl;
375  break;
376  case 1:
377  Ecen = 0.01219*pow(beta, -2.348) + 137.8;
378  T = 0.04856*pow(beta, -1.68) - 1.018;
379  Tp = PwrSeries(beta, -3.612, 422.8, -1973, 4081, -4109, 1641, 0.0, 0.0, 0.0, 0.0);
380  S = 0.0;
381  Sp = -0.03969*pow(beta, -1.775) +0.009034;
382  V0 = 2.12878493*EfieldScl;
383  break;
384  case 2:
385  Ecen = -0.01219*pow(beta, -2.348) + 112.2;
386  T = -0.04856*pow(beta, -1.68) + 1.018;
387  Tp = PwrSeries(beta, 3.612, -422.8, 1973, -4081, 4109, -1641, 0.0, 0.0, 0.0, 0.0);
388  S = 0.0;
389  Sp = -0.03969*pow(beta, -1.775) +0.009034;
390  V0 = 2.12878493*EfieldScl;
391  break;
392  default:
393  strm << "*** GetTransitFac: undef. number of gaps " << gaplabel << "\n";
394  throw std::runtime_error(strm.str());
395  }
396  break;
397  default:
398  strm << "*** GetTransitFac: undef. cavity type" << "\n";
399  throw std::runtime_error(strm.str());
400  }
401 
402  // Convert from [mm] to [m].
403 // Ecen /= MtoMM;
404 }
405 
406 
407 void ElementRFCavity::TransitFacMultipole(const int cavi, const std::string &flabel, const double CaviIonK,
408  double &T, double &S) const
409 {
410  double Ecen, Tp, Sp, V0;
411 
412  // For debugging of TTF function.
413  if (forcettfcalc) {
414  calTransfac(mlptable, get_column(flabel), 0, CaviIonK, false, Ecen, T, Tp, S, Sp, V0);
415  return;
416  }
417 
418  if (((cavi == 1) && (CaviIonK < 0.025 || CaviIonK > 0.055)) ||
419  ((cavi == 2) && (CaviIonK < 0.006 || CaviIonK > 0.035)) ||
420  ((cavi == 3) && (CaviIonK < 0.01687155 || CaviIonK > 0.0449908)) ||
421  ((cavi == 4) && (CaviIonK < 0.0112477 || CaviIonK > 0.0224954))) {
422  FLAME_LOG(DEBUG) << "*** TransitFacMultipole: CaviIonK out of Range" << "\n";
423  calTransfac(mlptable, get_column(flabel), 0, CaviIonK, false, Ecen, T, Tp, S, Sp, V0);
424  return;
425  }
426 
427  if (flabel == "CaviMlp_EFocus1") {
428  switch (cavi) {
429  case 1:
430  T = PwrSeries(CaviIonK, 1.256386e+02, -3.108322e+04, 3.354464e+06, -2.089452e+08, 8.280687e+09, -2.165867e+11,
431  3.739846e+12, -4.112154e+13, 2.613462e14, -7.316972e14);
432  S = PwrSeries(CaviIonK, 1.394183e+02, -3.299673e+04, 3.438044e+06, -2.070369e+08, 7.942886e+09, -2.013750e+11,
433  3.374738e+12, -3.605780e+13, 2.229446e+14, -6.079177e+14);
434  break;
435  case 2:
436  T = PwrSeries(CaviIonK, -9.450041e-01, -3.641390e+01, 9.926186e+03, -1.449193e+06, 1.281752e+08, -7.150297e+09,
437  2.534164e+11, -5.535252e+12, 6.794778e+13, -3.586197e+14);
438  S = PwrSeries(CaviIonK, 9.928055e-02, -5.545119e+01, 1.280168e+04, -1.636888e+06, 1.279801e+08, -6.379800e+09,
439  2.036575e+11, -4.029152e+12, 4.496323e+13, -2.161712e+14);
440  break;
441  case 3:
442  T = PwrSeries(CaviIonK, -1.000000e+00, 2.778823e-07, 6.820327e+01, 4.235106e-03, -1.926935e+03, 1.083516e+01,
443  2.996807e+04, 6.108642e+03, -3.864554e+05, 6.094390e+05);
444  S = PwrSeries(CaviIonK, -4.530303e-10, 1.625011e-07, -2.583224e-05, 2.478684e+01, -1.431967e-01, -1.545412e+03,
445  -1.569820e+02, 3.856713e+04, -3.159828e+04, -2.700076e+05);
446 
447  break;
448  case 4:
449  T = PwrSeries(CaviIonK, -1.000000e+00, -2.406447e-07, 9.480040e+01, -7.659927e-03, -4.926996e+03, -3.504383e+01,
450  1.712590e+05, -1.964643e+04, -4.142976e+06, 6.085390e+06);
451  S = PwrSeries(CaviIonK, 3.958048e-11, -2.496811e-08, 7.027794e-06, -8.662787e+01, 1.246098e-01, 9.462491e+03,
452  4.481784e+02, -4.552412e+05, 3.026543e+05, 8.798256e+06);
453  break;
454  }
455  } else if (flabel == "CaviMlp_EFocus2") {
456  switch (cavi) {
457  case 1:
458  T = PwrSeries(CaviIonK, 1.038803e+00, -9.121320e+00, 8.943931e+02, -5.619149e+04, 2.132552e+06, -5.330725e+07,
459  8.799404e+08, -9.246033e+09, 5.612073e+10, -1.499544e+11);
460  S = PwrSeries(CaviIonK, 1.305154e-02, -2.585211e+00, 2.696971e+02, -1.488249e+04, 5.095765e+05, -1.154148e+07,
461  1.714580e+08, -1.604935e+09, 8.570757e+09, -1.983302e+10);
462  break;
463  case 2:
464  T = PwrSeries(CaviIonK, 9.989307e-01, 7.299233e-01, -2.932580e+02, 3.052166e+04, -2.753614e+06, 1.570331e+08,
465  -5.677804e+09, 1.265012e+11, -1.584238e+12, 8.533351e+12);
466  S = PwrSeries(CaviIonK, -3.040839e-03, 2.016667e+00, -4.313590e+02, 5.855139e+04, -4.873584e+06, 2.605444e+08,
467  -8.968899e+09, 1.923697e+11, -2.339920e+12, 1.233014e+13);
468  break;
469  case 3:
470  T = PwrSeries(CaviIonK, 1.000000e+00, -4.410575e-06, -8.884752e+01, -7.927594e-02, 4.663277e+03, -2.515405e+02,
471  -1.797134e+05, -1.904305e+05, 8.999378e+06, -2.951362e+07);
472  S = PwrSeries(CaviIonK, 6.387813e-08, -2.300899e-05, 3.676251e-03, -1.703282e+02, 2.066461e+01, 1.704569e+04,
473  2.316653e+04, -1.328926e+06, 4.853676e+06, 1.132796e+06);
474 
475  break;
476  case 4:
477  T = PwrSeries(CaviIonK, 1.000000e+00, -5.025186e-06, -1.468976e+02, -2.520376e-01,2.048799e+04, -2.224267e+03,
478  -2.532091e+06, -4.613480e+06, 3.611911e+08, -1.891951e+09);
479  S = PwrSeries(CaviIonK, -1.801149e-08, 1.123280e-05, -3.126902e-03, 4.655245e+02, -5.431878e+01, -1.477730e+05,
480  -1.922110e+05, 2.795761e+07, -1.290046e+08, -4.656951e+08);
481  break;
482  }
483  } else if (flabel == "CaviMlp_EDipole") {
484  switch (cavi) {
485  case 1:
486  T = PwrSeries(CaviIonK, -1.005885e+00, 1.526489e+00, -1.047651e+02, 1.125013e+04, -4.669147e+05, 1.255841e+07,
487  -2.237287e+08, 2.535541e+09, -1.656906e+10, 4.758398e+10);
488  S = PwrSeries(CaviIonK, -2.586200e-02, 5.884367e+00, -6.407538e+02, 3.888964e+04, -1.488484e+06, 3.782592e+07,
489  -6.361033e+08, 6.817810e+09, -4.227114e+10, 1.155597e+11);
490  break;
491  case 2:
492  T = PwrSeries(CaviIonK, -9.999028e-01, -6.783669e-02, 1.415756e+02, -2.950990e+03, 2.640980e+05, -1.570742e+07,
493  5.770450e+08, -1.303686e+10, 1.654958e+11, -9.030017e+11);
494  S = PwrSeries(CaviIonK, 2.108581e-04, -3.700608e-01, 2.851611e+01, -3.502994e+03, 2.983061e+05, -1.522679e+07,
495  4.958029e+08, -1.002040e+10, 1.142835e+11, -5.617061e+11);
496  break;
497  default:
498  std::ostringstream strm;
499  strm << "*** 0.29 HWR and 0.53HWR havr no dipole term\n";
500  throw std::runtime_error(strm.str());
501  break;
502  }
503  } else if (flabel == "CaviMlp_EQuad") {
504  switch (cavi) {
505  case 1:
506  T = PwrSeries(CaviIonK, 1.038941e+00, -9.238897e+00, 9.127945e+02, -5.779110e+04, 2.206120e+06, -5.544764e+07,
507  9.192347e+08, -9.691159e+09, 5.896915e+10, -1.578312e+11);
508  S = PwrSeries(CaviIonK, 1.248096e-01, -2.923507e+01, 3.069331e+03, -1.848380e+05, 7.094882e+06, -1.801113e+08,
509  3.024208e+09, -3.239241e+10, 2.008767e+11, -5.496217e+11);
510  break;
511  case 2:
512  T = PwrSeries(CaviIonK, 1.000003e+00, -1.015639e-03, -1.215634e+02, 1.720764e+01, 3.921401e+03, 2.674841e+05,
513  -1.236263e+07, 3.128128e+08, -4.385795e+09, 2.594631e+10);
514  S = PwrSeries(CaviIonK, -1.756250e-05, 2.603597e-01, -2.551122e+00, -4.840638e+01, -2.870201e+04, 1.552398e+06,
515  -5.135200e+07, 1.075958e+09, -1.277425e+10, 6.540748e+10);
516  break;
517  case 3:
518  T = PwrSeries(CaviIonK, 1.000000e+00, 6.239107e-06, -1.697479e+02, 3.444883e-02, 1.225241e+04, -1.663533e+02,
519  -5.526645e+05, -3.593353e+05, 2.749580e+07, -9.689870e+07);
520  S = PwrSeries(CaviIonK, 2.128708e-07, -7.985618e-05, 1.240259e-02, -3.211339e+02, 7.098731e+01, 3.474652e+04,
521  8.187145e+04, -3.731688e+06, 1.802053e+07, -1.819958e+07);
522  break;
523  case 4:
524  T = PwrSeries(CaviIonK, 9.998746e-01, -2.431292e-05, -5.019138e+02, -1.176338e+00, 1.006054e+05, -9.908805e+03,
525  -1.148028e+07, -1.922707e+07, 1.432258e+09, -7.054482e+09);
526  S = PwrSeries(CaviIonK, 6.003340e-08, -1.482633e-02, 1.037590e-02, -2.235440e+03, 1.790006e+02, 6.456882e+05,
527  6.261020e+05, -1.055477e+08, 4.110502e+08, 2.241301e+09);
528  break;
529  }
530  } else if (flabel == "CaviMlp_HMono") {
531  switch (cavi) {
532  case 1:
533  T = PwrSeries(CaviIonK, 1.703336e+00, -1.671357e+02, 1.697657e+04, -9.843253e+05, 3.518178e+07, -8.043084e+08,
534  1.165760e+10, -1.014721e+11, 4.632851e+11, -7.604796e+11);
535  S = PwrSeries(CaviIonK, 1.452657e+01, -3.409550e+03, 3.524921e+05, -2.106663e+07, 8.022856e+08,
536  -2.019481e+10, 3.360597e+11, -3.565836e+12, 2.189668e+13, -5.930241e+13);
537  break;
538  case 2:
539  T = PwrSeries(CaviIonK, 1.003228e+00, -1.783406e+00, 1.765330e+02, -5.326467e+04, 4.242623e+06, -2.139672e+08,
540  6.970488e+09, -1.411958e+11, 1.617248e+12, -8.000662e+12);
541  S = PwrSeries(CaviIonK, -1.581533e-03, 1.277444e+00, -2.742508e+02, 3.966879e+04, -3.513478e+06, 1.962939e+08,
542  -6.991916e+09, 1.539708e+11, -1.910236e+12, 1.021016e+13);
543  break;
544  case 3:
545  T = PwrSeries(CaviIonK, 9.999990e-01, 3.477993e-04, -2.717994e+02, 4.554376e+00, 3.083481e+04, 8.441315e+03,
546  -2.439843e+06, 1.322379e+06, 1.501750e+08, -6.822135e+08);
547  S = PwrSeries(CaviIonK, 1.709084e-06, -6.240506e-04, 1.013278e-01, -2.649338e+02, 5.944163e+02,
548  4.588900e+04, 7.110518e+05, -2.226574e+07, 1.658524e+08, -3.976459e+08);
549  break;
550  case 4:
551  T = PwrSeries(CaviIonK, 1.000000e+00, -4.358956e-05, -7.923870e+02, -2.472669e+00, 2.241378e+05, -2.539286e+04,
552  -3.385480e+07, -6.375134e+07, 5.652166e+09, -3.355877e+10);
553  S = PwrSeries(CaviIonK, 1.163146e-07, -7.302018e-05, 2.048587e-02, -3.689694e+02, 3.632907e+02, 1.757838e+05,
554  1.327057e+06, -9.520645e+07, 9.406709e+08, -2.139562e+09);
555  break;
556  }
557  } else if (flabel == "CaviMlp_HDipole") {
558  switch (cavi) {
559  case 1:
560  T = PwrSeries(CaviIonK, 6.853803e-01, 7.075414e+01, -7.117391e+03, 3.985674e+05, -1.442888e+07, 3.446369e+08,
561  -5.420826e+09, 5.414689e+10, -3.116216e+11, 7.869717e+11);
562  S = PwrSeries(CaviIonK, 1.021102e+00, -2.441117e+02, 2.575274e+04, -1.569273e+06, 6.090118e+07, -1.562284e+09,
563  2.649289e+10, -2.864139e+11, 1.791634e+12, -4.941947e+12);
564  break;
565  case 2:
566  T = PwrSeries(CaviIonK, 1.014129e+00, -8.016304e+00, 1.631339e+03, -2.561826e+05, 2.115355e+07, -1.118723e+09,
567  3.821029e+10, -8.140248e+11, 9.839613e+12, -5.154137e+13);
568  S = PwrSeries(CaviIonK, -4.688714e-03, 3.299051e+00, -8.101936e+02, 1.163814e+05, -1.017331e+07, 5.607330e+08,
569  -1.967300e+10, 4.261388e+11, -5.194592e+12, 2.725370e+13);
570  break;
571  default:
572  std::ostringstream strm;
573  strm << "*** 0.29 HWR and 0.53HWR have no dipole term\n";
574  throw std::runtime_error(strm.str());
575  break;
576  }
577  } else if (flabel == "CaviMlp_HQuad") {
578  switch (cavi) {
579  case 1:
580  T = PwrSeries(CaviIonK, -1.997432e+00, 2.439177e+02, -2.613724e+04, 1.627837e+06, -6.429625e+07, 1.676173e+09,
581  -2.885455e+10, 3.163675e+11, -2.005326e+12, 5.600545e+12);
582  S = PwrSeries(CaviIonK, -2.470704e+00, 5.862902e+02, -6.135071e+04, 3.711527e+06, -1.431267e+08, 3.649414e+09,
583  -6.153570e+10, 6.617859e+11, -4.119861e+12, 1.131390e+13);
584  break;
585  case 2:
586  T = PwrSeries(CaviIonK, -1.000925e+00, 5.170302e-01, 9.311761e+01, 1.591517e+04, -1.302247e+06, 6.647808e+07,
587  -2.215417e+09, 4.603390e+10, -5.420873e+11, 2.764042e+12);
588  S = PwrSeries(CaviIonK, 3.119419e-04, -4.540868e-01, 5.433028e+01, -7.571946e+03, 6.792565e+05, -3.728390e+07,
589  1.299263e+09, -2.793705e+10, 3.377097e+11, -1.755126e+12);
590  break;
591  case 3:
592  T = PwrSeries(CaviIonK, -9.999997e-01, -1.049624e-04, 2.445420e+02, -1.288731e+00, -2.401575e+04, -1.972894e+03,
593  1.494708e+06, 2.898145e+05, -8.782506e+07, 3.566907e+08);
594  S = PwrSeries(CaviIonK, -7.925695e-07, 2.884963e-04, -4.667266e-02, 2.950936e+02, -2.712131e+02, -4.260259e+04,
595  -3.199682e+05, 1.103376e+07, -7.304474e+07, 1.479036e+08);
596  break;
597  case 4:
598  T = PwrSeries(CaviIonK, -1.000000e+00, 4.357777e-05, 7.605879e+02, 2.285787e+00, -2.009415e+05, 2.149581e+04,
599  2.773856e+07, 4.886782e+07, -4.127019e+09, 2.299278e+10);
600  S = PwrSeries(CaviIonK, -1.483304e-07, 9.278457e-05, -2.592071e-02, 1.690272e+03, -4.545599e+02, -6.192487e+05,
601  -1.632321e+06, 1.664856e+08, -1.124066e+09, -3.121299e+08);
602  break;
603  }
604  } else {
605  std::ostringstream strm;
606  strm << "*** TransitFacMultipole: undef. multipole type " << flabel << "\n";
607  throw std::runtime_error(strm.str());
608  }
609 }
610 
611 
612 static
613 double GetCavPhase(const int cavi, const Particle& ref, const double IonFys, const double multip, const std::vector<double>& P)
614 {
615  /* If the cavity is not at full power, the method gives synchrotron
616  * phase slightly different from the nominal value. */
617 
618  double IonEk, Fyc;
619 
620  IonEk = (ref.IonW-ref.IonEs)/MeVtoeV;
621 
622  switch (cavi) {
623  case 1:
624  Fyc = 4.394*pow(IonEk, -0.4965) - 4.731;
625  break;
626  case 2:
627  Fyc = 5.428*pow(IonEk, -0.5008) + 1.6;
628  break;
629  case 3:
630  Fyc = 22.35*pow(IonEk, -0.5348) + 2.026;
631  break;
632  case 4:
633  Fyc = 41.43*pow(IonEk, -0.5775) + 2.59839;
634  break;
635  case 5:
636  Fyc = 5.428*pow(IonEk, -0.5008) + 1.6;
637  break;
638  case 0:
639  Fyc = P[0]*pow(IonEk,P[1])- P[2];
640  break;
641  default:
642  std::ostringstream strm;
643  strm << "*** GetCavPhase: undef. cavity type" << "\n";
644  throw std::runtime_error(strm.str());
645  }
646 
647  return IonFys - Fyc - ref.phis*multip;
648 }
649 
650 
651 static
652 double GetCavPhaseComplex(const Particle& ref, const double IonFys, const double scale,
653  const double multip, const std::vector<double>& P)
654 {
655  double Ek, Fyc;
656  size_t npoly = P.size()/5;
657 
658  if (P.size()%5 != 0) {
659  throw std::runtime_error(SB()<<"*** Size of fitting parameter for the synchronous phase must be 5*n.");
660  }
661 
662  Ek = (ref.IonW-ref.IonEs)/MeVtoeV;
663 
664  Fyc = 0.0;
665  for (size_t i=0; i<npoly; i++){
666  Fyc += (P[5*i]*pow(Ek,P[5*i+1]) + P[5*i+2]*log(Ek) + P[5*i+3]*exp(-Ek) + P[5*i+4])
667  *ipow(scale, static_cast<int>(i));
668  }
669 
670  return IonFys - Fyc - ref.phis*multip;
671 }
672 
673 
674 static
675 void EvalGapModel(const double dis, const double IonW0, const Particle &real, const double IonFy0,
676  const double k, const double Lambda, const double Ecen,
677  const double T, const double S, const double Tp, const double Sp, const double V0,
678  double &IonW_f, double &IonFy_f)
679 {
680  double Iongamma_f, IonBeta_f, k_f;
681 
682  IonW_f = IonW0 + real.IonZ*V0*T*cos(IonFy0+k*Ecen)*MeVtoeV - real.IonZ*V0*S*sin(IonFy0+k*Ecen)*MeVtoeV;
683  Iongamma_f = IonW_f/real.IonEs;
684  IonBeta_f = sqrt(1e0-1e0/sqr(Iongamma_f));
685  k_f = 2e0*M_PI/(IonBeta_f*Lambda);
686 
687  IonFy_f = IonFy0 + k*Ecen + k_f*(dis-Ecen)
688  + real.IonZ*V0*k*(Tp*sin(IonFy0+k*Ecen)+Sp*cos(IonFy0+k*Ecen))/(2e0*(IonW0-real.IonEs)/MeVtoeV);
689 }
690 
691 
692 ElementRFCavity::ElementRFCavity(const Config& c)
693  :base_t(c)
694 {
695  ElementRFCavity::LoadCavityFile(c);
696 }
697 
698 
699 void ElementRFCavity::LoadCavityFile(const Config& c)
700 {
701  fRF = c.get<double>("f");
702  IonFys = c.get<double>("phi")*M_PI/180e0;
703  phi_ref = std::numeric_limits<double>::quiet_NaN();
704  cRm = c.get<double>("Rm", 0.0);
705  forcettfcalc = c.get<double>("forcettfcalc", 0.0)!=0.0;
706  MpoleLevel = get_flag(c, "MpoleLevel", 2);
707  EmitGrowth = get_flag(c, "EmitGrowth", 0);
708 
709  if (MpoleLevel != 0 && MpoleLevel != 1 && MpoleLevel != 2)
710  throw std::runtime_error(SB()<< "Undefined MpoleLevel: " << MpoleLevel);
711 
712  if (EmitGrowth != 0 && EmitGrowth != 1)
713  throw std::runtime_error(SB()<< "Undefined EmitGrowth: " << EmitGrowth);
714 
715  CavType = c.get<std::string>("cavtype");
716  std::string cavfile(c.get<std::string>("Eng_Data_Dir", defpath)),
717  fldmap(cavfile),
718  mlpfile(cavfile);
719 
720  if (CavType == "0.041QWR") {
721  cavi = 1;
722  fldmap += "/axisData_41.txt";
723  cavfile += "/Multipole41/thinlenlon_41.txt";
724  mlpfile += "/Multipole41/CaviMlp_41.txt";
725  } else if (CavType == "0.085QWR") {
726  cavi = 2;
727  fldmap += "/axisData_85.txt";
728  cavfile += "/Multipole85/thinlenlon_85.txt";
729  mlpfile += "/Multipole85/CaviMlp_85.txt";
730  } else if (CavType == "0.29HWR") {
731  cavi = 3;
732  fldmap += "/axisData_29.txt";
733  cavfile += "/Multipole29/thinlenlon_29.txt";
734  mlpfile += "/Multipole29/CaviMlp_29.txt";
735  } else if (CavType == "0.53HWR") {
736  cavi = 4;
737  fldmap += "/axisData_53.txt";
738  cavfile += "/Multipole53/thinlenlon_53.txt";
739  mlpfile += "/Multipole53/CaviMlp_53.txt";
740  } else if (CavType == "Generic") {
741  DataFile = cavfile+"/"+c.get<std::string>("datafile");
742  cavi = 0;
743  } else {
744  throw std::runtime_error(SB()<<"*** InitRFCav: undef. cavity type: " << CavType);
745  }
746 
747  if (cavi != 0)
748  {
749  numeric_table_cache *cache = numeric_table_cache::get();
750 
751  try{
752  numeric_table_cache::table_pointer ent = cache->fetch(fldmap);
753  CavData = *ent;
754  if(CavData.table.size1()==0 || CavData.table.size2()<2)
755  throw std::runtime_error("field map needs 2+ columns");
756  }catch(std::exception& e){
757  throw std::runtime_error(SB()<<"Error parsing '"<<fldmap<<"' : "<<e.what());
758  }
759 
760  try{
761  numeric_table_cache::table_pointer ent = cache->fetch(mlpfile);
762  mlptable = *ent;
763  if(mlptable.table.size1()==0 || mlptable.table.size2()<7)
764  throw std::runtime_error("CaviMlp needs 7+ columns");
765  }catch(std::exception& e){
766  throw std::runtime_error(SB()<<"Error parsing '"<<mlpfile<<"' : "<<e.what());
767  }
768 
769  {
770  std::ifstream fstrm(cavfile.c_str());
771 
772  std::string rawline;
773  unsigned line=0;
774  while(std::getline(fstrm, rawline)) {
775  line++;
776 
777  size_t cpos = rawline.find_first_not_of(" \t");
778  if(cpos==rawline.npos || rawline[cpos]=='%')
779  continue; // skip blank and comment lines
780 
781  cpos = rawline.find_last_not_of("\r\n");
782  if(cpos!=rawline.npos)
783  rawline = rawline.substr(0, cpos+1);
784 
785  std::istringstream lstrm(rawline);
786  RawParams params;
787  lstrm >> params.type >> params.name >> params.length >> params.aperature;
788  bool needE0 = params.type!="drift" && params.type!="AccGap";
789  if(needE0)
790  lstrm >> params.E0;
791  else
792  params.E0 = 0.0;
793 
794  if(lstrm.fail() && !lstrm.eof()) {
795  throw std::runtime_error(SB()<<"Error parsing line '"<<line<<"' in '"<<cavfile<<"'");
796  }
797  lattice.push_back(params);
798  }
799 
800  if(fstrm.fail() && !fstrm.eof()) {
801  throw std::runtime_error(SB()<<"Error, extra chars at end of file (line "<<line<<") in '"<<cavfile<<"'");
802  }
803  }
804  }
805  else
806  {
807  boost::shared_ptr<Config> conf;
808  std::string key(SB()<<DataFile<<"|"<<boost::filesystem::last_write_time(DataFile));
809  if ( CavConfMap.find(key) == CavConfMap.end() ) {
810  // not found in CavConfMap
811  try {
812  try {
813  GLPSParser P;
814  conf.reset(P.parse_file(DataFile.c_str()));
815  }catch(std::exception& e){
816  std::cerr<<"Parse error: "<<e.what()<<"\n";
817  }
818 
819  }catch(std::exception& e){
820  std::cerr<<"Error: "<<e.what()<<"\n";
821  }
822  CavConfMap.insert(std::make_pair(key, conf));
823  } else {
824  // found in CavConfMap
825  conf=CavConfMap[key];
826  }
827 
828  typedef Config::vector_t elements_t;
829  elements_t Es(conf->get<elements_t>("elements"));
830  for(elements_t::iterator it=Es.begin(), end=Es.end(); it!=end; ++it)
831  {
832  const Config& EC = *it;
833  const std::string& etype(EC.get<std::string>("type"));
834  const double elength(EC.get<double>("L"));
835  // fill in the lattice
836  RawParams params;
837  params.type = etype;
838  params.length = elength;
839  std::vector<double> attrs;
840  bool notdrift = etype!="drift";
841  if(notdrift)
842  {
843  const double eV0(EC.get<double>("V0"));
844  params.E0 = eV0;
845  EC.tryGet<std::vector<double> >("attr", attrs);
846  for(int i=0; i<10; i++)
847  {
848  params.Tfit.push_back(attrs[i]);
849  }
850  for(int i=0; i<10; i++)
851  {
852  params.Sfit.push_back(attrs[i+10]);
853  }
854  }
855  bool needSynAccTab = params.type=="AccGap";
856  // SynAccTab should only update once
857  if(needSynAccTab && SynAccTab.size()==0)
858  {
859  for(int i=0; i<3; i++)
860  {
861  SynAccTab.push_back(attrs[i+20]);
862  }
863  }
864  lattice.push_back(params);
865  }
866 
867  std::vector<double> Ez;
868  bool checker = conf->tryGet<std::vector<double> >("Ez", Ez);
869  if (!checker) throw std::runtime_error(SB()<<"'Ez' is missing in RF cavity file.\n");
870  CavData.readvec(Ez,2);
871  conf->tryGet<double>("Rm",cRm);
872 
873  // Get extra parameters for complex synchronous phase definition
874  have_RefNrm = conf->tryGet<double>("RefNorm", RefNrm);
875  have_SynComplex = conf->tryGet<std::vector<double> >("SyncFit", SynComplex);
876  have_EkLim = conf->tryGet<std::vector<double> >("EnergyLimit", EkLim);
877  have_NrLim = conf->tryGet<std::vector<double> >("NormLimit", NrLim);
878  }
879 }
880 
881 void ElementRFCavity::GetCavMatParams(const int cavi, const double beta_tab[], const double gamma_tab[], const double CaviIonK[],
882  CavTLMLineType& lineref) const
883 {
884  if(lattice.empty())
885  throw std::runtime_error("empty RF cavity lattice");
886 
887  lineref.clear();
888 
889  size_t i;
890  double s=CavData.table(0,0);
891  for(i=0; i<lattice.size(); i++) {
892  const RawParams& P = lattice[i];
893  {
894  double E0=0.0, T=0.0, S=0.0, Accel=0.0;
895 
896  if ((P.type != "drift") && (P.type != "AccGap"))
897  E0 = P.E0;
898 
899  s+=lattice[i].length;
900 
901  if (P.type == "drift") {
902  } else if (P.type == "EFocus1") {
903  if (s < 0e0) {
904  // First gap. By reflection 1st Gap EFocus1 is 2nd gap EFocus2.
905  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EFocus2", CaviIonK[0], T, S);
906  // First gap *1, transverse E field the same.
907  S = -S;
908  } else {
909  // Second gap.
910  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EFocus1", CaviIonK[1], T, S);
911  }
912  } else if (P.type == "EFocus2") {
913  if (s < 0e0) {
914  // First gap.
915  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EFocus1", CaviIonK[0], T, S);
916  S = -S;
917  } else {
918  // Second gap.
919  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EFocus2", CaviIonK[1], T, S);
920  }
921  } else if (P.type == "EDipole") {
922  if (MpoleLevel >= 1) {
923  if (s < 0e0) {
924  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EDipole", CaviIonK[0], T, S);
925  // First gap *1, transverse E field the same.
926  S = -S;
927  } else {
928  // Second gap.
929  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EDipole", CaviIonK[1], T, S);
930  }
931  }
932  } else if (P.type == "EQuad") {
933  if (MpoleLevel >= 2) {
934  if (s < 0e0) {
935  // First gap.
936  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EQuad", CaviIonK[0], T, S);
937  S = -S;
938  } else {
939  // Second Gap
940  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_EQuad", CaviIonK[1], T, S);
941  }
942  }
943  } else if (P.type == "HMono") {
944  if (MpoleLevel >= 2) {
945  if (s < 0e0) {
946  // First gap.
947  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_HMono", CaviIonK[0], T, S);
948  T = -T;
949  } else {
950  // Second Gap
951  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_HMono", CaviIonK[1], T, S);
952  }
953  }
954  } else if (P.type == "HDipole") {
955  if (MpoleLevel >= 1) {
956  if (s < 0e0) {
957  // First gap.
958  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_HDipole", CaviIonK[0], T, S);
959  T = -T;
960  } else {
961  // Second gap.
962  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_HDipole", CaviIonK[1], T, S);
963  }
964  }
965  } else if (P.type == "HQuad") {
966  if (MpoleLevel >= 2) {
967  if (s < 0e0) {
968  // First gap.
969  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_HQuad", CaviIonK[0], T, S);
970  T = -T;
971  } else {
972  // Second gap.
973  ElementRFCavity::TransitFacMultipole(cavi, "CaviMlp_HQuad", CaviIonK[1], T, S);
974  }
975  }
976  } else if (P.type == "AccGap") {
977  if (s < 0e0) {
978  // First gap.
979  Accel = (beta_tab[0]*gamma_tab[0])/((beta_tab[1]*gamma_tab[1]));
980  } else {
981  // Second gap.
982  Accel = (beta_tab[1]*gamma_tab[1])/((beta_tab[2]*gamma_tab[2]));
983  }
984  } else {
985  std::ostringstream strm;
986  strm << "*** GetCavMatParams: undef. multipole element " << P.type << "\n";
987  throw std::runtime_error(strm.str());
988  }
989 
990  lineref.set(s, P.type, E0, T, S, Accel);
991  }
992  }
993 
994  if (FLAME_LOG_CHECK(DEBUG)) {
995  std::cout << "\n";
996  lineref.show();
997  }
998 }
999 
1000 
1001 void ElementRFCavity::GenCavMat2(const int cavi, const double dis, const double EfieldScl, const double TTF_tab[],
1002  const double beta_tab[], const double gamma_tab[], const double Lambda,
1003  Particle &real, const double IonFys[], const double Rm, state_t::matrix_t &M,
1004  const CavTLMLineType& linetab) const
1005 {
1006  /* RF cavity model, transverse only defocusing.
1007  * 2-gap matrix model. */
1008 
1009  int seg;
1010  double k_s[3];
1011  double Ecens[2], Ts[2], Ss[2], V0s[2], ks[2], L1, L2, L3;
1012  double beta, gamma, kfac, V0, T, S, kfdx, kfdy, dpy, Accel, IonFy;
1013  state_t::matrix_t Idmat, Mlon_L1, Mlon_K1, Mlon_L2;
1014  state_t::matrix_t Mlon_K2, Mlon_L3, Mlon, Mtrans, Mprob;
1015 
1016  // fetch the log level once to speed our loop
1017  const bool logme = FLAME_LOG_CHECK(DEBUG);
1018 
1019  const double IonA = 1e0;
1020 
1021  using boost::numeric::ublas::prod;
1022 
1023  Idmat = boost::numeric::ublas::identity_matrix<double>(PS_Dim);
1024 
1025  k_s[0] = 2e0*M_PI/(beta_tab[0]*Lambda);
1026  k_s[1] = 2e0*M_PI/(beta_tab[1]*Lambda);
1027  k_s[2] = 2e0*M_PI/(beta_tab[2]*Lambda);
1028 
1029  // Longitudinal model: Drift-Kick-Drift, dis: total lenghth centered at 0,
1030  // Ecens[0] & Ecens[1]: Electric Center position where accel kick applies, Ecens[0] < 0
1031  // TTFtab: 2*6 vector, Ecens, T Tp S Sp, V0;
1032 
1033  Ecens[0] = TTF_tab[0];
1034  Ts[0] = TTF_tab[1];
1035  Ss[0] = TTF_tab[3];
1036  V0s[0] = TTF_tab[5];
1037  ks[0] = 0.5*(k_s[0]+k_s[1]);
1038  L1 = dis + Ecens[0]; //try change dis/2 to dis 14/12/12
1039 
1040  Mlon_L1 = Idmat;
1041  Mlon_K1 = Idmat;
1042  // Pay attention, original is -
1043  Mlon_L1(4, 5) = -2e0*M_PI/Lambda*(1e0/cube(beta_tab[0]*gamma_tab[0])*MeVtoeV/real.IonEs*L1);
1044  // Pay attention, original is -k1-k2
1045  Mlon_K1(5, 4) = -real.IonZ*V0s[0]*Ts[0]*sin(IonFys[0]+ks[0]*L1)-real.IonZ*V0s[0]*Ss[0]*cos(IonFys[0]+ks[0]*L1);
1046 
1047  Ecens[1] = TTF_tab[6];
1048  Ts[1] = TTF_tab[7];
1049  Ss[1] = TTF_tab[9];
1050  V0s[1] = TTF_tab[11];
1051  ks[1] = 0.5*(k_s[1]+k_s[2]);
1052  L2 = Ecens[1] - Ecens[0];
1053 
1054  Mlon_L2 = Idmat;
1055  Mlon_K2 = Idmat;
1056 
1057  Mlon_L2(4, 5) = -2e0*M_PI/Lambda*(1e0/cube(beta_tab[1]*gamma_tab[1])*MeVtoeV/real.IonEs*L2); //Problem is Here!!
1058  Mlon_K2(5, 4) = -real.IonZ*V0s[1]*Ts[1]*sin(IonFys[1]+ks[1]*Ecens[1])-real.IonZ*V0s[1]*Ss[1]*cos(IonFys[1]+ks[1]*Ecens[1]);
1059 
1060  L3 = dis - Ecens[1]; //try change dis/2 to dis 14/12/12
1061 
1062  Mlon_L3 = Idmat;
1063  Mlon_L3(4, 5) = -2e0*M_PI/Lambda*(1e0/cube(beta_tab[2]*gamma_tab[2])*MeVtoeV/real.IonEs*L3);
1064 
1065  Mlon = prod(Mlon_K1, Mlon_L1);
1066  Mlon = prod(Mlon_L2, Mlon);
1067  Mlon = prod(Mlon_K2, Mlon);
1068  Mlon = prod(Mlon_L3, Mlon);
1069 // std::cout<<__FUNCTION__<<" Mlon "<<Mlon<<"\n";
1070 
1071  // Transverse model
1072  // Drift-FD-Drift-LongiKick-Drift-FD-Drift-0-Drift-FD-Drift-LongiKick-Drift-FD-Drift
1073 
1074  seg = 0;
1075 
1076  Mtrans = Idmat;
1077  Mprob = Idmat;
1078 
1079  beta = beta_tab[0];
1080  gamma = gamma_tab[0];
1081  IonFy = IonFys[0];
1082  kfac = k_s[0];
1083 
1084  V0 = 0e0, T = 0e0, S = 0e0, kfdx = 0e0, kfdy = 0e0, dpy = 0e0;
1085  size_t n;
1086  double s=CavData.table(0,0);
1087  for(n=0; n<lattice.size(); n++) {
1088  const RawParams& P = lattice[n];
1089 
1090  s+=lattice[n].length;
1091 
1092  if (false)
1093  printf("%9.5f %8s %8s %9.5f %9.5f %9.5f\n",
1094  s, P.type.c_str(), P.name.c_str(), P.length, P.aperature, P.E0);
1095 
1096  Mprob = Idmat;
1097  if (P.type == "drift") {
1098  IonFy = IonFy + kfac*P.length;
1099 
1100  Mprob(0, 1) = P.length;
1101  Mprob(2, 3) = P.length;
1102  Mtrans = prod(Mprob, Mtrans);
1103  } else if (P.type == "EFocus1") {
1104  V0 = linetab.E0[n]*EfieldScl;
1105  T = linetab.T[n];
1106  S = linetab.S[n];
1107  kfdx = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy))/Rm;
1108  kfdy = kfdx;
1109  if(logme) {
1110  FLAME_LOG(FINE)<<" X EFocus1 kfdx="<<kfdx<<"\n"
1111  <<" Y "<<linetab.E0[n]<<" "<<EfieldScl<<" "<<beta
1112  <<" "<<gamma<<" "<<IonFy<<" "<<Rm<<"\n Z "<<T<<" "<<S<<"\n";
1113  }
1114 
1115  Mprob(1, 0) = kfdx;
1116  Mprob(3, 2) = kfdy;
1117  Mtrans = prod(Mprob, Mtrans);
1118  } else if (P.type == "EFocus2") {
1119  V0 = linetab.E0[n]*EfieldScl;
1120  T = linetab.T[n];
1121  S = linetab.S[n];
1122  kfdx = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy))/Rm;
1123  kfdy = kfdx;
1124  if(logme) FLAME_LOG(FINE)<<" X EFocus2 kfdx="<<kfdx<<"\n";
1125 
1126  Mprob(1, 0) = kfdx;
1127  Mprob(3, 2) = kfdy;
1128  Mtrans = prod(Mprob, Mtrans);
1129  } else if (P.type == "EDipole") {
1130  if (MpoleLevel >= 1) {
1131  V0 = linetab.E0[n]*EfieldScl;
1132  T = linetab.T[n];
1133  S = linetab.S[n];
1134  dpy = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy));
1135  if(logme) FLAME_LOG(FINE)<<" X EDipole dpy="<<dpy<<"\n";
1136 
1137  Mprob(3, 6) = dpy;
1138  Mtrans = prod(Mprob, Mtrans);
1139  }
1140  } else if (P.type == "EQuad") {
1141  if (MpoleLevel >= 2) {
1142  V0 = linetab.E0[n]*EfieldScl;
1143  T = linetab.T[n];
1144  S = linetab.S[n];
1145  kfdx = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy))/Rm;
1146  kfdy = -kfdx;
1147  if(logme) FLAME_LOG(FINE)<<" X EQuad kfdx="<<kfdx<<"\n";
1148 
1149  Mprob(1, 0) = kfdx;
1150  Mprob(3, 2) = kfdy;
1151  Mtrans = prod(Mprob, Mtrans);
1152  }
1153  } else if (P.type == "HMono") {
1154  if (MpoleLevel >= 2) {
1155  V0 = linetab.E0[n]*EfieldScl;
1156  T = linetab.T[n];
1157  S = linetab.S[n];
1158  kfdx = -MU0*C0*real.IonZ*V0/beta/gamma/IonA/AU*(T*cos(IonFy+M_PI/2e0)-S*sin(IonFy+M_PI/2e0))/Rm;
1159  kfdy = kfdx;
1160  if(logme) FLAME_LOG(FINE)<<" X HMono kfdx="<<kfdx<<"\n";
1161 
1162  Mprob(1, 0) = kfdx;
1163  Mprob(3, 2) = kfdy;
1164  Mtrans = prod(Mprob, Mtrans);
1165  }
1166  } else if (P.type == "HDipole") {
1167  if (MpoleLevel >= 1) {
1168  V0 = linetab.E0[n]*EfieldScl;
1169  T = linetab.T[n];
1170  S = linetab.S[n];
1171  dpy = -MU0*C0*real.IonZ*V0/beta/gamma/IonA/AU*(T*cos(IonFy+M_PI/2e0)-S*sin(IonFy+M_PI/2e0));
1172  if(logme) FLAME_LOG(FINE)<<" X HDipole dpy="<<dpy<<"\n";
1173 
1174  Mprob(3, 6) = dpy;
1175  Mtrans = prod(Mprob, Mtrans);
1176  }
1177  } else if (P.type == "HQuad") {
1178  if (MpoleLevel >= 2) {
1179  if (s < 0e0) {
1180  // First gap.
1181  beta = (beta_tab[0]+beta_tab[1])/2e0;
1182  gamma = (gamma_tab[0]+gamma_tab[1])/2e0;
1183  } else {
1184  beta = (beta_tab[1]+beta_tab[2])/2e0;
1185  gamma = (gamma_tab[1]+gamma_tab[2])/2e0;
1186  }
1187  V0 = linetab.E0[n]*EfieldScl;
1188  T = linetab.T[n];
1189  S = linetab.S[n];
1190  kfdx = -MU0*C0*real.IonZ*V0/beta/gamma/IonA/AU*(T*cos(IonFy+M_PI/2e0)-S*sin(IonFy+M_PI/2e0))/Rm;
1191  kfdy = -kfdx;
1192  if(logme) FLAME_LOG(FINE)<<" X HQuad kfdx="<<kfdx<<"\n";
1193 
1194  Mprob(1, 0) = kfdx;
1195  Mprob(3, 2) = kfdy;
1196  Mtrans = prod(Mprob, Mtrans);
1197  }
1198  } else if (P.type == "AccGap") {
1199  //IonFy = IonFy + real.IonZ*V0s[0]*kfac*(TTF_tab[2]*sin(IonFy)
1200  // + TTF_tab[4]*cos(IonFy))/2/((gamma-1)*real.IonEs/MeVtoeV); //TTF_tab[2]~Tp
1201  seg = seg + 1;
1202  beta = beta_tab[seg];
1203  gamma = gamma_tab[seg];
1204  kfac = 2e0*M_PI/(beta*Lambda);
1205  Accel = linetab.Accel[n];
1206  if(logme) FLAME_LOG(FINE)<<" X AccGap Accel="<<Accel<<"\n";
1207 
1208  Mprob(1, 1) = Accel;
1209  Mprob(3, 3) = Accel;
1210  Mtrans = prod(Mprob, Mtrans);
1211  } else {
1212  std::ostringstream strm;
1213  strm << "*** GetCavMat: undef. multipole type " << P.type << "\n";
1214  throw std::runtime_error(strm.str());
1215  }
1216  // FLAME_LOG(FINE) << Elem << "\n";
1217  // PrtMat(Mprob);
1218 
1219  if(logme) FLAME_LOG(FINE)<<"Elem "<<P.name<<":"<<P.type<<"\n Mtrans "<<Mtrans<<"\nMprob "<<Mprob<<"\n";
1220  }
1221 
1222  M = Mtrans;
1223 
1224  M(4, 4) = Mlon(4, 4);
1225  M(4, 5) = Mlon(4, 5);
1226  M(5, 4) = Mlon(5, 4);
1227  M(5, 5) = Mlon(5, 5);
1228 }
1229 
1230 
1231 void ElementRFCavity::GetCavMat(const int cavi, const int cavilabel, const double Rm, Particle &real,
1232  const double EfieldScl, const double IonFyi_s,
1233  const double IonEk_s, state_t::matrix_t &M,
1234  CavTLMLineType &linetab) const
1235 {
1236  double CaviLambda, Ecen[2], T[2], Tp[2], S[2], Sp[2], V0[2];
1237  double dis, IonW_s[3], IonFy_s[3], gamma_s[3], beta_s[3], CaviIonK_s[3];
1238  double CaviIonK[2];
1239 
1240  CaviLambda = C0/fRF*MtoMM;
1241 
1242  IonW_s[0] = IonEk_s + real.IonEs;
1243  IonFy_s[0] = IonFyi_s;
1244  gamma_s[0] = IonW_s[0]/real.IonEs;
1245  beta_s[0] = sqrt(1e0-1e0/sqr(gamma_s[0]));
1246  CaviIonK_s[0] = 2e0*M_PI/(beta_s[0]*CaviLambda);
1247 
1248  size_t n = CavData.table.size1();
1249  assert(n>0);
1250  dis = (CavData.table(n-1,0)-CavData.table(0,0))/2e0;
1251 
1252  ElementRFCavity::TransFacts(cavilabel, beta_s[0], CaviIonK_s[0], 1, EfieldScl,
1253  Ecen[0], T[0], Tp[0], S[0], Sp[0], V0[0]);
1254  EvalGapModel(dis, IonW_s[0], real, IonFy_s[0], CaviIonK_s[0], CaviLambda,
1255  Ecen[0], T[0], S[0], Tp[0], Sp[0], V0[0], IonW_s[1], IonFy_s[1]);
1256  gamma_s[1] = IonW_s[1]/real.IonEs;
1257  beta_s[1] = sqrt(1e0-1e0/sqr(gamma_s[1]));
1258  CaviIonK_s[1] = 2e0*M_PI/(beta_s[1]*CaviLambda);
1259 
1260  ElementRFCavity::TransFacts(cavilabel, beta_s[1], CaviIonK_s[1], 2, EfieldScl,
1261  Ecen[1], T[1], Tp[1], S[1], Sp[1], V0[1]);
1262  EvalGapModel(dis, IonW_s[1], real, IonFy_s[1], CaviIonK_s[1], CaviLambda,
1263  Ecen[1], T[1], S[1], Tp[1], Sp[1], V0[1], IonW_s[2], IonFy_s[2]);
1264  gamma_s[2] = IonW_s[2]/real.IonEs;
1265  beta_s[2] = sqrt(1e0-1e0/sqr(gamma_s[2]));
1266  CaviIonK_s[2] = 2e0*M_PI/(beta_s[2]*CaviLambda);
1267 
1268  Ecen[0] = Ecen[0] - dis;
1269 
1270  double TTF_tab[] = {Ecen[0], T[0], Tp[0], S[0], Sp[0], V0[0], Ecen[1], T[1], Tp[1], S[1], Sp[1], V0[1]};
1271  CaviIonK[0] = (CaviIonK_s[0]+CaviIonK_s[1])/2e0;
1272  CaviIonK[1] = (CaviIonK_s[1]+CaviIonK_s[2])/2e0;
1273 
1274  if (false) {
1275  printf("\n GetCavMat:\n");
1276  printf("CaviIonK: %15.10f %15.10f %15.10f\n", CaviIonK_s[0], CaviIonK_s[1], CaviIonK_s[2]);
1277  printf("CaviIonK: %15.10f %15.10f\n", CaviIonK[0], CaviIonK[1]);
1278  printf("beta: %15.10f %15.10f %15.10f\n", beta_s[0], beta_s[1], beta_s[2]);
1279  printf("gamma: %15.10f %15.10f %15.10f\n", gamma_s[0], gamma_s[1], gamma_s[2]);
1280  printf("Ecen: %15.10f %15.10f\n", Ecen[0], Ecen[1]);
1281  printf("T: %15.10f %15.10f\n", T[0], T[1]);
1282  printf("Tp: %15.10f %15.10f\n", Tp[0], Tp[1]);
1283  printf("S: %15.10f %15.10f\n", S[0], S[1]);
1284  printf("Sp: %15.10f %15.10f\n", Sp[0], Sp[1]);
1285  printf("V0: %15.10f %15.10f\n", V0[0], V0[1]);
1286  }
1287 
1288  ElementRFCavity::GetCavMatParams(cavi, beta_s, gamma_s, CaviIonK, linetab);
1289  ElementRFCavity::GenCavMat2(cavi, dis, EfieldScl, TTF_tab, beta_s, gamma_s, CaviLambda, real, IonFy_s, Rm, M, linetab);
1290 }
1291 
1292 
1293 void ElementRFCavity::GetCavMatGeneric(Particle &real, const double EfieldScl, const double IonFyi_s,
1294  const double IonEk_s, state_t::matrix_t &M, CavTLMLineType &linetab) const
1295 {
1296 
1297  const double IonA = 1e0;
1298  const bool logme = FLAME_LOG_CHECK(DEBUG);
1299  state_t::matrix_t Idmat, Mlon, Mtrans,Mprob,MprobLon;
1300  Idmat = boost::numeric::ublas::identity_matrix<double>(PS_Dim);
1301  Mtrans = Idmat;
1302  Mlon = Idmat;
1303 
1304  double IonW0 = IonEk_s + real.IonEs;
1305  double gamma0 = IonW0/real.IonEs;
1306  double beta0 = sqrt(1e0-1e0/sqr(gamma0));
1307  double IonFy0 = IonFyi_s;
1308  double CaviLambda = C0/fRF*MtoMM;
1309  double kfac0 = 2e0*M_PI/(beta0*CaviLambda);
1310  double dis = 0.0, V0 = 0.0, T = 0.0, S = 0.0, kfdx = 0.0, kfdy = 0.0, Accel = 0.0;
1311  double dpy = 0.0;
1312 
1313  double IonW=IonW0,gamma=gamma0,beta=beta0,IonFy=IonFy0,kfac=kfac0;
1314 
1315  assert(cRm>0);
1316 
1317  for(unsigned n=0; n<lattice.size(); n++) {
1318  const RawParams& P = lattice[n];
1319 
1320  dis+=lattice[n].length;
1321 
1322  if (false)
1323  printf("%9.5f %8s %8s %9.5f %9.5f %9.5f\n",
1324  dis, P.type.c_str(), P.name.c_str(), P.length, P.aperature, P.E0);
1325 
1326  Mprob = Idmat;
1327  MprobLon = Idmat;
1328  if (P.type == "drift") {
1329  IonFy = IonFy + kfac*P.length;
1330 
1331  Mprob(0, 1) = P.length;
1332  Mprob(2, 3) = P.length;
1333  Mtrans = prod(Mprob, Mtrans);
1334 
1335  // Pay attention, original is -
1336  MprobLon(4, 5) = -2e0*M_PI/CaviLambda*(1e0/cube(beta*gamma)*MeVtoeV/real.IonEs*P.length);
1337  Mlon = prod(MprobLon, Mlon);
1338 
1339  } else if (P.type == "EFocus") {
1340  V0 = P.E0*EfieldScl;
1341  T = calFitPow(kfac,P.Tfit);
1342  S = calFitPow(kfac,P.Sfit);
1343  kfdx = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy))/cRm;
1344  kfdy = kfdx;
1345  if(logme) {
1346  FLAME_LOG(FINE)<<" X EFocus1 kfdx="<<kfdx<<"\n"
1347  <<" Y "<<linetab.E0[n]<<" "<<EfieldScl<<" "<<beta
1348  <<" "<<gamma<<" "<<IonFy<<" "<<cRm<<"\n Z "<<T<<" "<<S<<"\n";
1349  }
1350 
1351  Mprob(1, 0) = kfdx;
1352  Mprob(3, 2) = kfdy;
1353  Mtrans = prod(Mprob, Mtrans);
1354 
1355  } else if (P.type == "EDipole") {
1356  if (MpoleLevel >= 1) {
1357  V0 = P.E0*EfieldScl;
1358  T = calFitPow(kfac,P.Tfit);
1359  S = calFitPow(kfac,P.Sfit);
1360  dpy = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy));
1361  if(logme) FLAME_LOG(FINE)<<" X EDipole dpy="<<dpy<<"\n";
1362 
1363  Mprob(3, 6) = dpy;
1364  Mtrans = prod(Mprob, Mtrans);
1365  }
1366  } else if (P.type == "EQuad") {
1367  if (MpoleLevel >= 2) {
1368  V0 = P.E0*EfieldScl;
1369  T = calFitPow(kfac,P.Tfit);
1370  S = calFitPow(kfac,P.Sfit);
1371  kfdx = real.IonZ*V0/sqr(beta)/gamma/IonA/AU*(T*cos(IonFy)-S*sin(IonFy))/cRm;
1372  kfdy = -kfdx;
1373  if(logme) FLAME_LOG(FINE)<<" X EQuad kfdx="<<kfdx<<"\n";
1374 
1375  Mprob(1, 0) = kfdx;
1376  Mprob(3, 2) = kfdy;
1377  Mtrans = prod(Mprob, Mtrans);
1378  }
1379  } else if (P.type == "HMono") {
1380  if (MpoleLevel >= 2) {
1381  V0 = P.E0*EfieldScl;
1382  T = calFitPow(kfac,P.Tfit);
1383  S = calFitPow(kfac,P.Sfit);
1384  kfdx = -MU0*C0*real.IonZ*V0/beta/gamma/IonA/AU*(T*cos(IonFy+M_PI/2e0)-S*sin(IonFy+M_PI/2e0))/cRm;
1385  kfdy = kfdx;
1386  if(logme) FLAME_LOG(FINE)<<" X HMono kfdx="<<kfdx<<"\n";
1387 
1388  Mprob(1, 0) = kfdx;
1389  Mprob(3, 2) = kfdy;
1390  Mtrans = prod(Mprob, Mtrans);
1391  }
1392  } else if (P.type == "HDipole") {
1393  if (MpoleLevel >= 1) {
1394  V0 = P.E0*EfieldScl;
1395  T = calFitPow(kfac,P.Tfit);
1396  S = calFitPow(kfac,P.Sfit);
1397  dpy = -MU0*C0*real.IonZ*V0/beta/gamma/IonA/AU*(T*cos(IonFy+M_PI/2e0)-S*sin(IonFy+M_PI/2e0));
1398  if(logme) FLAME_LOG(FINE)<<" X HDipole dpy="<<dpy<<"\n";
1399 
1400  Mprob(3, 6) = dpy;
1401  Mtrans = prod(Mprob, Mtrans);
1402  }
1403  } else if (P.type == "HQuad") {
1404  if (MpoleLevel >= 2) {
1405  V0 = P.E0*EfieldScl;
1406  T = calFitPow(kfac,P.Tfit);
1407  S = calFitPow(kfac,P.Sfit);
1408  kfdx = -MU0*C0*real.IonZ*V0/beta/gamma/IonA/AU*(T*cos(IonFy+M_PI/2e0)-S*sin(IonFy+M_PI/2e0))/cRm;
1409  kfdy = -kfdx;
1410  if(logme) FLAME_LOG(FINE)<<" X HQuad kfdx="<<kfdx<<"\n";
1411 
1412  Mprob(1, 0) = kfdx;
1413  Mprob(3, 2) = kfdy;
1414  Mtrans = prod(Mprob, Mtrans);
1415  }
1416  } else if (P.type == "AccGap") {
1417  V0 = P.E0*EfieldScl;
1418  T = calFitPow(kfac,P.Tfit);
1419  S = calFitPow(kfac,P.Sfit);
1420  IonW=IonW+real.IonZ*V0*MeVtoeV*T*cos(IonFy)-real.IonZ*V0*MeVtoeV*S*sin(IonFy);
1421  double gamma_f=IonW/real.IonEs;
1422  double beta_f=sqrt(1.0-1.0/(gamma_f*gamma_f));
1423  kfac = 2e0*M_PI/(beta_f*CaviLambda);
1424  Accel = (beta*gamma)/((beta_f*gamma_f));
1425  if(logme) FLAME_LOG(FINE)<<" X AccGap Accel="<<Accel<<"\n";
1426 
1427  Mprob(1, 1) = Accel;
1428  Mprob(3, 3) = Accel;
1429  Mtrans = prod(Mprob, Mtrans);
1430 
1431  MprobLon(5, 4) = -real.IonZ*V0*T*sin(IonFy)-real.IonZ*V0*S*cos(IonFy);
1432  Mlon = prod(MprobLon, Mlon);
1433 
1434  beta=beta_f;
1435  gamma=gamma_f;
1436 
1437  } else {
1438  std::ostringstream strm;
1439  strm << "*** GetCavMat: undef. multipole type " << P.type << "\n";
1440  throw std::runtime_error(strm.str());
1441  }
1442  // FLAME_LOG(FINE) << Elem << "\n";
1443  // PrtMat(Mprob);
1444 
1445  if(logme) FLAME_LOG(FINE)<<"Elem "<<P.name<<":"<<P.type<<"\n Mtrans "<<Mtrans<<"\nMprob "<<Mprob<<"\n";
1446  }
1447 
1448  M = Mtrans;
1449 
1450  M(4, 4) = Mlon(4, 4);
1451  M(4, 5) = Mlon(4, 5);
1452  M(5, 4) = Mlon(5, 4);
1453  M(5, 5) = Mlon(5, 5);
1454 }
1455 
1456 double ElementRFCavity::calFitPow(double kfac,const std::vector<double>& Tfit) const
1457 {
1458  int order=Tfit.size();
1459  double res=0.0;
1460  for (int ii=0; ii<order; ii++)
1461  {
1462  res=res+Tfit[ii]*ipow(kfac,order-ii-1);
1463  }
1464  return res;
1465 }
1466 
1467 void ElementRFCavity::GetCavBoost(const numeric_table &CavData, Particle &state, const double IonFy0,
1468  const double EfieldScl, double &IonFy) const
1469 {
1470  size_t n = CavData.table.size1();
1471 
1472  assert(n>1);
1473 
1474  const bool logme = FLAME_LOG_CHECK(DEBUG);
1475 
1476  const double dis = CavData.table(n-1, 0) - CavData.table(0, 0),
1477  dz = dis/(n-1),
1478  CaviLambda = C0/fRF*MtoMM;
1479  // assumes dz is constant even though CavData contains actual z positions are available
1480 
1481  FLAME_LOG(DEBUG)<<__FUNCTION__
1482  <<" IonFy0="<<IonFy0
1483  <<" fRF="<<fRF
1484  <<" EfieldScl="<<EfieldScl
1485  <<" state="<<state
1486  <<"\n";
1487  IonFy = IonFy0;
1488  // Sample rate is different for RF Cavity; due to different RF frequencies.
1489  // IonK = state.SampleIonK;
1490  double CaviIonK = 2e0*M_PI*fRF/(state.beta*C0*MtoMM);
1491  for (size_t k = 0; k < n-1; k++) {
1492  double IonFylast = IonFy;
1493  IonFy += CaviIonK*dz;
1494  state.IonW += state.IonZ*EfieldScl*(CavData.table(k,1)+CavData.table(k+1,1))/2e0
1495  *cos((IonFylast+IonFy)/2e0)*dz/MtoMM;
1496  double IonGamma = state.IonW/state.IonEs;
1497  double IonBeta = sqrt(1e0-1e0/sqr(IonGamma));
1498  if ((state.IonW-state.IonEs) < 0e0) {
1499  state.IonW = state.IonEs;
1500  IonBeta = 0e0;
1501  //TODO: better handling of this error?
1502  // will be divide by zero (NaN)
1503  }
1504  CaviIonK = 2e0*M_PI/(IonBeta*CaviLambda);
1505  if(logme) FLAME_LOG(DEBUG)<<" "<<k<<" CaviIonK="<<CaviIonK<<" IonW="<<state.IonW<<"\n";
1506  }
1507 }
1508 
1509 
1510 void ElementRFCavity::PropagateLongRFCav(Particle &ref, double& phi_ref) const
1511 {
1512  double multip, EfieldScl, caviFy, IonFy_i, IonFy_o;
1513  double fsync = conf().get<double>("syncflag", 1.0);
1514 
1515  multip = fRF/ref.SampleFreq;
1516  EfieldScl = conf().get<double>("scl_fac"); // Electric field scale factor.
1517 
1518  if (cavi == 0 && have_EkLim) {
1519  if (ref.IonEk/MeVtoeV < EkLim[0] || ref.IonEk/MeVtoeV > EkLim[1])
1520  FLAME_LOG(WARN)<< "Warning: RF cavity incident energy (" << ref.IonEk/MeVtoeV
1521  << " [MeV]) is out of range (" << EkLim[0] << " ~ " << EkLim[1] << ").\n";
1522  }
1523 
1524  if (fsync >= 1.0) {
1525  if (cavi == 0 && have_RefNrm && have_SynComplex && fsync == 1.0) {
1526  // Get driven phase from synchronous phase based on peak position
1527  double NormScl = EfieldScl*ref.IonZ/RefNrm;
1528  if (have_NrLim) {
1529  if (NormScl < NrLim[0] || NormScl > NrLim[1])
1530  FLAME_LOG(WARN)<< "Warning: RF cavity normalized scale (" << NormScl
1531  << ") is out of range (" << NrLim[0] << " ~ " << NrLim[1] << ").\n";
1532  }
1533  caviFy = GetCavPhaseComplex(ref, IonFys, NormScl, multip, SynComplex);
1534  } else {
1535  // Get driven phase from synchronous phase based on sin fit model
1536  caviFy = GetCavPhase(cavi, ref, IonFys, multip, SynAccTab);
1537  }
1538  } else {
1539  caviFy = conf().get<double>("phi")*M_PI/180e0;
1540  }
1541 
1542  IonFy_i = multip*ref.phis + caviFy;
1543  phi_ref = caviFy;
1544  FLAME_LOG(DEBUG)<<"RF long phase"
1545  " caviFy="<<caviFy
1546  <<" multip="<<multip
1547  <<" phis="<<ref.phis
1548  <<" IonFy_i="<<IonFy_i
1549  <<" EfieldScl="<<EfieldScl
1550  <<"\n";
1551 
1552  // For the reference particle, evaluate the change of:
1553  // kinetic energy, absolute phase, beta, and gamma.
1554  GetCavBoost(CavData, ref, IonFy_i, EfieldScl, IonFy_o);
1555 
1556  ref.IonEk = ref.IonW - ref.IonEs;
1557  ref.recalc();
1558  ref.phis += (IonFy_o-IonFy_i)/multip;
1559 }
1560 
1561 
1562 void ElementRFCavity::calRFcaviEmitGrowth(const state_t::matrix_t &matIn, Particle &state, const int n, const double betaf, const double gamaf,
1563  const double aveX2i, const double cenX, const double aveY2i, const double cenY,
1564  state_t::matrix_t &matOut)
1565 {
1566  // Evaluate emittance growth.
1567  int k;
1568  double ionLamda, E0TL, DeltaPhi, kpX, fDeltaPhi, f2DeltaPhi, gPhisDeltaPhi, deltaAveXp2f, XpIncreaseFactor;
1569  double kpY, deltaAveYp2f, YpIncreaseFactor, kpZ, ionK, aveZ2i, deltaAveZp2, longiTransFactor, ZpIncreaseFactor;
1570 
1571  matOut = matIn;
1572 
1573  ionLamda = C0/fRF*MtoMM;
1574 
1575  // safe to look at last_real_out[] here as we are called (from advance() ) after it is updated
1576  const double accIonW = last_real_out[n].IonW -last_real_in[n].IonW,
1577  ave_beta = (last_real_out[n].beta +last_real_in[n].beta)/2.0,
1578  ave_gamma = (last_real_out[n].gamma+last_real_in[n].gamma)/2.0;
1579 
1580  E0TL = accIonW/cos(IonFys)/state.IonZ;
1581 
1582  // for rebuncher, because no acceleration, E0TL would be wrong when cos(ionFys) is devided.
1583  if (cos(IonFys) > -0.0001 && cos(IonFys) < 0.0001) E0TL = 0e0;
1584  DeltaPhi = sqrt(matIn(4, 4));
1585  // ionLamda in m, kpX in 1/mm
1586  kpX = -M_PI*fabs(state.IonZ)*E0TL/state.IonEs/sqr(ave_beta*ave_gamma)/betaf/gamaf/ionLamda;
1587  fDeltaPhi = 15e0/sqr(DeltaPhi)*(3e0/sqr(DeltaPhi)*(sin(DeltaPhi)/DeltaPhi-cos(DeltaPhi))-(sin(DeltaPhi)/DeltaPhi));
1588  f2DeltaPhi = 15e0/sqr(2e0*DeltaPhi)*(3e0/sqr(2e0*DeltaPhi)
1589  *(sin(2e0*DeltaPhi)/(2e0*DeltaPhi)-cos(2e0*DeltaPhi))-(sin(2e0*DeltaPhi)/(2e0*DeltaPhi)));
1590  gPhisDeltaPhi = 0.5e0*(1+(sqr(sin(IonFys))-sqr(cos(IonFys)))*f2DeltaPhi);
1591  deltaAveXp2f = kpX*kpX*(gPhisDeltaPhi-sqr(sin(IonFys)*fDeltaPhi))*(aveX2i+cenX*cenX);
1592  XpIncreaseFactor = 1e0;
1593 
1594  if (deltaAveXp2f+matIn(1, 1) > 0e0) XpIncreaseFactor = sqrt((deltaAveXp2f+matIn(1, 1))/matIn(1, 1));
1595 
1596  // ionLamda in m
1597  kpY = -M_PI*fabs(state.IonZ)*E0TL/state.IonEs/sqr(ave_beta*ave_gamma)/betaf/gamaf/ionLamda;
1598  deltaAveYp2f = sqr(kpY)*(gPhisDeltaPhi-sqr(sin(IonFys)*fDeltaPhi))*(aveY2i+sqr(cenY));
1599  YpIncreaseFactor = 1.0;
1600  if (deltaAveYp2f+matIn(3, 3)>0) {
1601  YpIncreaseFactor = sqrt((deltaAveYp2f+matIn(3, 3))/matIn(3, 3));
1602  }
1603 
1604  kpZ = -2e0*kpX*sqr(ave_gamma);
1605  //unit: 1/mm
1606  ionK = 2e0*M_PI/(ave_beta*ionLamda);
1607  aveZ2i = sqr(DeltaPhi)/sqr(ionK);
1608  deltaAveZp2 = sqr(kpZ*DeltaPhi)*aveZ2i*(sqr(cos(IonFys))/8e0+DeltaPhi*sin(IonFys)/576e0);
1609  longiTransFactor = 1e0/(ave_gamma-1e0)/state.IonEs*MeVtoeV;
1610  ZpIncreaseFactor = 1e0;
1611  if (deltaAveZp2+matIn(5, 5)*sqr(longiTransFactor) > 0e0)
1612  ZpIncreaseFactor = sqrt((deltaAveZp2+matIn(5, 5)*sqr(longiTransFactor))/(matIn(5, 5)*sqr(longiTransFactor)));
1613 
1614  for (k = 0; k < PS_Dim; k++) {
1615  matOut(1, k) *= XpIncreaseFactor;
1616  matOut(k, 1) *= XpIncreaseFactor;
1617  matOut(3, k) *= YpIncreaseFactor;
1618  matOut(k, 3) *= YpIncreaseFactor;
1619  matOut(5, k) *= ZpIncreaseFactor;
1620  matOut(k, 5) *= ZpIncreaseFactor;
1621  }
1622 }
1623 
1624 
1625 void ElementRFCavity::InitRFCav(Particle &real, state_t::matrix_t &M, CavTLMLineType &linetab)
1626 {
1627  int cavilabel;
1628  double Rm, multip, IonFy_i, Ek_i, EfieldScl, IonFy_o;
1629 
1630  FLAME_LOG(DEBUG)<<"RF recompute start "<<real<<"\n";
1631 
1632  if (cavi == 1) {
1633  cavilabel = 41;
1634  Rm = 17e0;
1635  } else if (cavi== 2) {
1636  cavilabel = 85;
1637  Rm = 17e0;
1638  } else if (cavi== 3) {
1639  cavilabel = 29;
1640  Rm = 20e0;
1641  } else if (cavi== 4) {
1642  cavilabel = 53;
1643  Rm = 20e0;
1644  } else if (cavi== 5) {
1645  // 5 Cell elliptical.
1646  cavilabel = 53;
1647  Rm = 20e0;
1648  } else if (cavi == 0) {
1649  // Generic2
1650  cavilabel = 0;
1651  Rm = 0e0; //this is dummy.
1652  } else {
1653  throw std::logic_error(SB()<<"*** InitRFCav: undef. cavity type: after ctor");
1654  }
1655 
1656  multip = fRF/real.SampleFreq;
1657 
1658  IonFy_i = multip*real.phis + phi_ref;
1659  Ek_i = real.IonEk;
1660  real.IonW = real.IonEk + real.IonEs;
1661 
1662  EfieldScl = conf().get<double>("scl_fac"); // Electric field scale factor.
1663  ElementRFCavity::GetCavBoost(CavData, real, IonFy_i, EfieldScl, IonFy_o); // updates IonW
1664 
1665  real.IonEk = real.IonW - real.IonEs;
1666  real.recalc();
1667  real.phis += (IonFy_o-IonFy_i)/multip;
1668 
1669  FLAME_LOG(DEBUG)<<"RF recompute before "<<real
1670  <<" cavi="<<cavi
1671  <<" cavilabel="<<cavilabel
1672  <<" Rm="<<Rm
1673  <<" EfieldScl="<<EfieldScl
1674  <<" IonFy_i="<<IonFy_i
1675  <<" Ek_i="<<Ek_i
1676  <<" fRF="<<fRF
1677  <<"\n";
1678 
1679  if (cavi> 0)
1680  {
1681  GetCavMat(cavi, cavilabel, Rm, real, EfieldScl, IonFy_i, Ek_i, M, linetab);
1682  }
1683  else
1684  {
1685  GetCavMatGeneric(real, EfieldScl, IonFy_i, Ek_i, M, linetab);
1686  }
1687 
1688 
1689  //Wrapper for fequency jump in rf cavity
1690  if (multip != 1) {
1691  for (int i=0; i < state_t::maxsize; i++){
1692  M(i,4) *= multip;
1693  M(4,i) /= multip;
1694  }
1695  }
1696 
1697  FLAME_LOG(DEBUG)<<"RF recompute after "<<real<<"\n"
1698  <<" YY "<<M<<"\n"
1699  ;
1700 }
GLPSParser::parse_file
Config * parse_file(const char *fname, const bool have_lattice=true)
Open and parse a file.
Definition: config.cpp:364
Particle::IonW
double IonW
Total energy. (dependent)
Definition: moment.h:22
GLPSParser
Interface to lattice file parser.
Definition: config.h:240
Particle::IonEk
double IonEk
Kinetic energy.
Definition: moment.h:22
Particle::IonZ
double IonZ
Charge state.
Definition: moment.h:22
Particle::phis
double phis
Absolute synchrotron phase [rad].
Definition: moment.h:22
Particle::IonEs
double IonEs
Rest energy.
Definition: moment.h:22
Particle::recalc
void recalc()
Definition: moment.h:47
Particle::SampleFreq
double SampleFreq
Sampling frequency [Hz].
Definition: moment.h:22
Config
Associative configuration container.
Definition: config.h:66
Config::tryGet
bool tryGet(const std::string &name, T &val) const
Definition: config.h:156
Config::get
detail::RT< T >::type get(const std::string &name) const
Definition: config.h:123
numeric_table
Definition: util.h:72
Generated on Thu Apr 17 2025 01:29:13 for FLAME by   doxygen 1.8.5