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