FLAME  devel
 All Classes Functions Variables Typedefs Enumerations Pages
moment.cpp
1 
2 #include <fstream>
3 
4 #include <limits>
5 
6 #include <boost/lexical_cast.hpp>
7 
8 #include <boost/numeric/ublas/vector_proxy.hpp>
9 #include <boost/numeric/ublas/matrix_proxy.hpp>
10 
11 #include "flame/constants.h"
12 
13 #include "flame/moment.h"
14 #include "flame/moment_sup.h"
15 #include "flame/rf_cavity.h"
16 #include "flame/chg_stripper.h"
17 
18 #define sqr(x) ((x)*(x))
19 #define cube(x) ((x)*(x)*(x))
20 
21 namespace ub = boost::numeric::ublas;
22 
23 namespace {
24 
25 // ARR should be an array-like object (std::vector or or ublas vector or matrix storage)
26 // fill 'to' using config.get<>(name)
27 // 'T' selects throw error (true) or return boolean (false)
28 template<typename ARR>
29 bool load_storage(ARR& to, const Config& conf, const std::string& name, bool T=true)
30 {
31  try{
32  const std::vector<double>& val(conf.get<std::vector<double> >(name));
33  if(to.size()!=val.size()) {
34  throw std::invalid_argument(SB()<<"Array "<<name<<" must have "<<to.size()<<" elements, not "<<val.size());
35  }
36  std::copy(val.begin(), val.end(), to.begin());
37 
38  return true;
39  }catch(std::exception&){
40  if(T)
41  throw std::invalid_argument(name+" not defined or has wrong type (must be vector)");
42  else
43  return false;
44  // default to identity
45  }
46 }
47 
48 } // namespace
49 
50 std::ostream& operator<<(std::ostream& strm, const Particle& P)
51 {
52  strm
53  <<"IonZ="<<std::scientific << std::setprecision(10)<<P.IonZ
54  <<" IonQ="<<P.IonQ
55  <<" IonEs="<<P.IonEs
56  <<" IonEk="<<P.IonEk
57  <<" SampleIonK="<<P.SampleIonK
58  <<" phis="<<P.phis
59  <<" IonW="<<P.IonW
60  <<" gamma="<<P.gamma
61  <<" beta="<<P.beta
62  <<" bg="<<P.bg
63  ;
64  return strm;
65 }
66 
67 MomentState::MomentState(const Config& c)
68  :StateBase(c)
69  ,ref()
70  ,real()
71  ,moment0_env(maxsize, 0e0)
72  ,moment0_rms(maxsize, 0e0)
73  ,moment1_env(boost::numeric::ublas::identity_matrix<double>(maxsize))
74 {
75  // hack. getArray() promises that returned pointers will remain valid for our lifetime.
76  // This may not be true if std::vectors are resized.
77  // Reserve for up to 10 states and hope for the best...
78  // either need to provide real limit to max. states, or change getArray() iface
79  real.reserve(10);
80 
81  double icstate_f = 0.0;
82  bool have_cstate = c.tryGet<double>("cstate", icstate_f);
83  size_t icstate = (size_t)icstate_f;
84 
85  std::string vectorname(c.get<std::string>("vector_variable", "moment0"));
86  std::string matrixname(c.get<std::string>("matrix_variable", "initial"));
87 
88  std::vector<double> ics, nchg;
89  bool have_ics = c.tryGet<std::vector<double> >("IonChargeStates", ics);
90 
91  ref.IonEs = c.get<double>("IonEs", 0e0),
92  ref.IonEk = c.get<double>("IonEk", 0e0),
93  ref.SampleFreq = c.get<double>("SampleFreq", SampleFreqDefault);
94  ref.recalc();
95 
96  if(!have_ics) {
97  ref.IonZ = c.get<double>("IonZ", 0e0);
98  ref.IonQ = c.get<double>("IonQ", 1e0);
99 
100  ics.push_back(ref.IonZ);
101  nchg.push_back(ref.IonQ);
102 
103  } else {
104  if(ics.empty())
105  throw std::invalid_argument("IonChargeStates w/ length 0");
106  if(icstate>=ics.size())
107  throw std::invalid_argument("IonChargeStates[cstate] is out of bounds");
108 
109  nchg = c.get<std::vector<double> >("NCharge");
110  if(nchg.size()!=ics.size())
111  throw std::invalid_argument("NCharge[] and IonChargeStates[] must have equal length");
112 
113  bool have_ionz = c.tryGet<double>("IonZ", ref.IonZ);
114  if(!have_ionz) {
115  ref.IonZ = ics[0];
116  }
117 
118  bool have_ionq = c.tryGet<double>("IonQ", ref.IonQ);
119  if(!have_ionq) {
120  ref.IonQ = nchg[0];
121  }
122  }
123 
124  /* Possible configurations
125  * 1. Neither 'cstate' nor 'IonChargeStates' defined (empty Config).
126  * No charge states, must go through source element to be useful
127  * 2. 'IonChargeStates' defined, but not 'cstate'.
128  * Load all charge states
129  * 3. 'cstate' and 'IonChargeStates' defined.
130  * Load a single charge state
131  */
132  if(!have_cstate && !have_ics) {
133  // no-op
134  } else if(!have_cstate && have_ics) {
135  // many charge states
136  icstate = 0;
137 
138  } else if(have_cstate && have_ics) {
139  // single charge state
140 
141  // drop other than selected state
142  ics[0] = ics[icstate];
143  nchg[0] = nchg[icstate];
144  ics.resize(1);
145  nchg.resize(1);
146 
147  } else {
148  throw std::invalid_argument("MomentState: must define IonChargeStates and NCharge when cstate is set");
149  }
150 
151  if(have_ics) {
152  real.resize(ics.size());
153  moment0.resize(ics.size());
154  moment1.resize(ics.size());
155  transmat.resize(ics.size());
156 
157  for(size_t i=0; i<ics.size(); i++) {
158  std::string num(boost::lexical_cast<std::string>(icstate+i));
159 
160  moment0[i].resize(maxsize);
161  moment1[i].resize(maxsize, maxsize);
162  moment1[i] = boost::numeric::ublas::identity_matrix<double>(maxsize);
163  transmat[i].resize(maxsize, maxsize);
164  transmat[i] = boost::numeric::ublas::identity_matrix<double>(maxsize);
165 
166  load_storage(moment0[i].data(), c, vectorname+num);
167  load_storage(moment1[i].data(), c, matrixname+num);
168 
169  real[i] = ref;
170 
171  real[i].IonZ = ics[i];
172  real[i].IonQ = nchg[i];
173 
174  real[i].phis = moment0[i][PS_S];
175  real[i].IonEk += moment0[i][PS_PS]*MeVtoeV;
176 
177  real[i].recalc();
178  }
179  } else {
180  real.resize(1); // hack, ensure at least one element so getArray() can return some pointer
181  real[0] = ref;
182 
183  moment0.resize(1);
184  moment1.resize(1);
185  transmat.resize(1);
186  moment0[0].resize(maxsize);
187  moment1[0].resize(maxsize, maxsize);
188  transmat[0].resize(maxsize, maxsize);
189 
190  load_storage(moment0[0].data(), c, vectorname, false);
191  load_storage(moment1[0].data(), c, matrixname, false);
192  transmat[0] = boost::numeric::ublas::identity_matrix<double>(maxsize);
193  }
194 
195  last_caviphi0 = 0e0;
196  calc_rms();
197 }
198 
199 MomentState::~MomentState() {}
200 
201 void MomentState::calc_rms()
202 {
203  assert(real.size()>0);
204  assert(moment0_env.size()==maxsize);
205  assert(moment0_rms.size()==maxsize);
206  assert(moment1_env.size1()==maxsize);
207  assert(moment1_env.size2()==maxsize);
208 
209  double totQ = 0.0;
210  for(size_t n=0; n<real.size(); n++) {
211  totQ += real[n].IonQ;
212  if(n==0)
213  noalias(moment0_env) = moment0[n]*real[n].IonQ;
214  else
215  noalias(moment0_env) += moment0[n]*real[n].IonQ;
216  }
217  moment0_env /= totQ;
218 
219  // Zero orbit terms.
220  moment1_env = ub::zero_matrix<double>(state_t::maxsize);
221  boost::numeric::ublas::slice S(0, 1, 6);
222  for(size_t n=0; n<real.size(); n++) {
223  const double Q = real[n].IonQ;
224  state_t::vector_t m0diff(moment0[n]-moment0_env);
225 
226  ub::project(moment1_env, S, S) += ub::project(Q*(moment1[n]+ub::outer_prod(m0diff, m0diff)), S, S);
227  }
228  moment1_env /= totQ;
229 
230  for(size_t j=0; j<maxsize; j++) {
231  moment0_rms[j] = sqrt(moment1_env(j,j));
232  }
233 }
234 
235 MomentState::MomentState(const MomentState& o, clone_tag t)
236  :StateBase(o, t)
237  ,ref(o.ref)
238  ,real(o.real)
239  ,moment0(o.moment0)
240  ,moment1(o.moment1)
241  ,transmat(o.transmat)
242  ,moment0_env(o.moment0_env)
243  ,moment0_rms(o.moment0_rms)
244  ,moment1_env(o.moment1_env)
245  ,last_caviphi0(o.last_caviphi0)
246 {}
247 
248 void MomentState::assign(const StateBase& other)
249 {
250  const MomentState *O = dynamic_cast<const MomentState*>(&other);
251  if(!O)
252  throw std::invalid_argument("Can't assign State: incompatible types");
253  ref = O->ref;
254  real = O->real;
255  moment0 = O->moment0;
256  moment1 = O->moment1;
257  transmat = O->transmat;
258  moment0_env = O->moment0_env;
259  moment0_rms = O->moment0_rms;
260  moment1_env = O->moment1_env;
261  last_caviphi0 = O->last_caviphi0;
262  StateBase::assign(other);
263 }
264 
265 void MomentState::show(std::ostream& strm, int level) const
266 {
267  if(real.empty()) {
268  strm<<"State: empty";
269  return;
270  }
271 
272  if(level<=0) {
273  strm<<"State: moment0 mean="<<moment0_env;
274  }
275  if(level>=1) {
276  strm << std::scientific << std::setprecision(8)
277  << "\nState:\n energy [eV] =\n" << std::setw(20) << real[0].IonEk << "\n moment0 mean =\n ";
278  for (size_t k = 0; k < MomentState::maxsize; k++)
279  strm << std::scientific << std::setprecision(10) << std::setw(18) << moment0_env(k) << ",";
280  strm << std::scientific << std::setprecision(10)
281  << "\n moment0 rms =\n ";
282  for (size_t k = 0; k < MomentState::maxsize; k++)
283  strm << std::scientific << std::setprecision(10) << std::setw(18) << moment0_rms(k) << ",";
284  strm << "\n moment1 mean =\n";
285  for (size_t j = 0; j < MomentState::maxsize; j++) {
286  strm << " ";
287  for (size_t k = 0; k < MomentState::maxsize; k++) {
288  strm << std::scientific << std::setprecision(10) << std::setw(18) << moment1_env(j, k) << ",";
289  }
290  if (j < MomentState::maxsize-1) strm << "\n";
291  }
292  }
293  if(level>=2) {
294  strm<< "\n Reference state:\n "<<ref<<"\n";
295  for(size_t k=0; k<size(); k++) {
296  strm<<" Charge state "<<k<<"\n"
297  " "<<real[k]<<"\n"
298  " moment0 "<<moment0[k]<<"\n";
299  if(level>=3)
300  strm<<
301  " moment1 "<<moment1[k]<<"\n";
302  }
303  }
304 }
305 
306 bool MomentState::getArray(unsigned idx, ArrayInfo& Info) {
307  unsigned I=0;
308  if(idx==I++) {
309  Info.name = "moment1_env";
310  Info.ptr = &moment1_env(0,0);
311  Info.type = ArrayInfo::Double;
312  Info.ndim = 2;
313  Info.dim[0] = moment1_env.size1();
314  Info.dim[1] = moment1_env.size2();
315  Info.stride[0] = sizeof(double)*moment1_env.size2();
316  Info.stride[1] = sizeof(double);
317  return true;
318  } else if(idx==I++) {
319  /* Slight evilness here
320  * moment0 is vector of ublas::matrix
321  * We assume ublax::matrix uses storage bounded_array<>, and that this storage
322  * is really a C array, which means that everything is part of one big allocation.
323  * Further we assume that all entries in the vector have the same shape.
324  * If this isn't the case, then SIGSEGV here we come...
325  */
326  static_assert(sizeof(moment1[0])>=sizeof(double)*maxsize*maxsize,
327  "storage assumption violated");
328  Info.name = "moment1";
329  Info.ptr = &moment1[0](0,0);
330  Info.type = ArrayInfo::Double;
331  Info.ndim = 3;
332  Info.dim[0] = moment1[0].size1();
333  Info.dim[1] = moment1[0].size2();
334  Info.dim[2] = moment1.size();
335  Info.stride[0] = sizeof(double)*moment1_env.size2();
336  Info.stride[1] = sizeof(double);
337  Info.stride[2] = sizeof(moment1[0]);
338  return true;
339  } else if(idx==I++) {
340  static_assert(sizeof(transmat[0])>=sizeof(double)*maxsize*maxsize,
341  "storage assumption violated");
342  Info.name = "transmat";
343  Info.ptr = &transmat[0](0,0);
344  Info.type = ArrayInfo::Double;
345  Info.ndim = 3;
346  Info.dim[0] = transmat[0].size1();
347  Info.dim[1] = transmat[0].size2();
348  Info.dim[2] = transmat.size();
349  Info.stride[0] = sizeof(double)*moment1_env.size2();
350  Info.stride[1] = sizeof(double);
351  Info.stride[2] = sizeof(transmat[0]);
352  return true;
353  } else if(idx==I++) {
354  Info.name = "moment0_env";
355  Info.ptr = &moment0_env(0);
356  Info.type = ArrayInfo::Double;
357  Info.ndim = 1;
358  Info.dim[0] = moment0_env.size();
359  Info.stride[0] = sizeof(double);
360  return true;
361  } else if(idx==I++) {
362  Info.name = "moment0_rms";
363  Info.ptr = &moment0_rms(0);
364  Info.type = ArrayInfo::Double;
365  Info.ndim = 1;
366  Info.dim[0] = moment0_rms.size();
367  Info.stride[0] = sizeof(double);
368  return true;
369  } else if(idx==I++) {
370  // more evilness here, see above
371  static_assert(sizeof(moment0[0])>=sizeof(double)*maxsize,
372  "storage assumption violated");
373  Info.name = "moment0";
374  Info.ptr = &moment0[0][0];
375  Info.type = ArrayInfo::Double;
376  Info.ndim = 2;
377  Info.dim[0] = moment0[0].size();
378  Info.dim[1] = moment0.size();
379  Info.stride[0] = sizeof(double);
380  Info.stride[1] = sizeof(moment0[0]);
381  return true;
382  } else if(idx==I++) {
383  Info.name = "ref_IonZ";
384  Info.ptr = &ref.IonZ;
385  Info.type = ArrayInfo::Double;
386  Info.ndim = 0;
387  return true;
388  } else if(idx==I++) {
389  Info.name = "ref_IonQ";
390  Info.ptr = &ref.IonQ;
391  Info.type = ArrayInfo::Double;
392  Info.ndim = 0;
393  return true;
394  } else if(idx==I++) {
395  Info.name = "ref_IonEs";
396  Info.ptr = &ref.IonEs;
397  Info.type = ArrayInfo::Double;
398  Info.ndim = 0;
399  return true;
400  } else if(idx==I++) {
401  Info.name = "ref_IonW";
402  Info.ptr = &ref.IonW;
403  Info.type = ArrayInfo::Double;
404  Info.ndim = 0;
405  return true;
406  } else if(idx==I++) {
407  Info.name = "ref_gamma";
408  Info.ptr = &ref.gamma;
409  Info.type = ArrayInfo::Double;
410  Info.ndim = 0;
411  return true;
412  } else if(idx==I++) {
413  Info.name = "ref_beta";
414  Info.ptr = &ref.beta;
415  Info.type = ArrayInfo::Double;
416  Info.ndim = 0;
417  return true;
418  } else if(idx==I++) {
419  Info.name = "ref_bg";
420  Info.ptr = &ref.bg;
421  Info.type = ArrayInfo::Double;
422  Info.ndim = 0;
423  return true;
424  } else if(idx==I++) {
425  Info.name = "ref_SampleFreq";
426  Info.ptr = &ref.SampleFreq;
427  Info.type = ArrayInfo::Double;
428  Info.ndim = 0;
429  return true;
430  } else if(idx==I++) {
431  Info.name = "ref_SampleIonK";
432  Info.ptr = &ref.SampleIonK;
433  Info.type = ArrayInfo::Double;
434  Info.ndim = 0;
435  return true;
436  } else if(idx==I++) {
437  Info.name = "ref_phis";
438  Info.ptr = &ref.phis;
439  Info.type = ArrayInfo::Double;
440  Info.ndim = 0;
441  return true;
442  } else if(idx==I++) {
443  Info.name = "ref_IonEk";
444  Info.ptr = &ref.IonEk;
445  Info.type = ArrayInfo::Double;
446  Info.ndim = 0;
447  return true;
448  } else if(idx==I++) {
449  Info.name = "IonZ";
450  Info.ptr = &real[0].IonZ;
451  Info.type = ArrayInfo::Double;
452  Info.ndim = 1;
453  Info.dim [0] = real.size();
454  Info.stride[0] = sizeof(real[0]);
455  // Note: this array is discontigious as we reference a single member from a Particle[]
456  return true;
457  } else if(idx==I++) {
458  Info.name = "IonEs";
459  Info.ptr = &real[0].IonEs;
460  Info.type = ArrayInfo::Double;
461  Info.ndim = 1;
462  Info.dim [0] = real.size();
463  Info.stride[0] = sizeof(real[0]);
464  return true;
465  } else if(idx==I++) {
466  Info.name = "IonW";
467  Info.ptr = &real[0].IonW;
468  Info.type = ArrayInfo::Double;
469  Info.ndim = 1;
470  Info.dim [0] = real.size();
471  Info.stride[0] = sizeof(real[0]);
472  return true;
473  } else if(idx==I++) {
474  Info.name = "gamma";
475  Info.ptr = &real[0].gamma;
476  Info.type = ArrayInfo::Double;
477  Info.ndim = 1;
478  Info.dim [0] = real.size();
479  Info.stride[0] = sizeof(real[0]);
480  return true;
481  } else if(idx==I++) {
482  Info.name = "beta";
483  Info.ptr = &real[0].beta;
484  Info.type = ArrayInfo::Double;
485  Info.ndim = 1;
486  Info.dim [0] = real.size();
487  Info.stride[0] = sizeof(real[0]);
488  return true;
489  } else if(idx==I++) {
490  Info.name = "bg";
491  Info.ptr = &real[0].bg;
492  Info.type = ArrayInfo::Double;
493  Info.ndim = 1;
494  Info.dim [0] = real.size();
495  Info.stride[0] = sizeof(real[0]);
496  return true;
497  } else if(idx==I++) {
498  Info.name = "SampleFreq";
499  Info.ptr = &real[0].SampleFreq;
500  Info.type = ArrayInfo::Double;
501  Info.ndim = 1;
502  Info.dim [0] = real.size();
503  Info.stride[0] = sizeof(real[0]);
504  return true;
505  } else if(idx==I++) {
506  Info.name = "SampleIonK";
507  Info.ptr = &real[0].SampleIonK;
508  Info.type = ArrayInfo::Double;
509  Info.ndim = 1;
510  Info.dim [0] = real.size();
511  Info.stride[0] = sizeof(real[0]);
512  return true;
513  } else if(idx==I++) {
514  Info.name = "phis";
515  Info.ptr = &real[0].phis;
516  Info.type = ArrayInfo::Double;
517  Info.ndim = 1;
518  Info.dim [0] = real.size();
519  Info.stride[0] = sizeof(real[0]);
520  return true;
521  } else if(idx==I++) {
522  Info.name = "IonEk";
523  Info.ptr = &real[0].IonEk;
524  Info.type = ArrayInfo::Double;
525  Info.ndim = 1;
526  Info.dim [0] = real.size();
527  Info.stride[0] = sizeof(real[0]);
528  return true;
529  } else if(idx==I++) {
530  Info.name = "IonQ";
531  Info.ptr = &real[0].IonQ;
532  Info.type = ArrayInfo::Double;
533  Info.ndim = 1;
534  Info.dim [0] = real.size();
535  Info.stride[0] = sizeof(real[0]);
536  // Note: this array is discontigious as we reference a single member from a Particle[]
537  return true;
538  } else if(idx==I++) {
539  Info.name = "last_caviphi0";
540  Info.ptr = &last_caviphi0;
541  Info.type = ArrayInfo::Double;
542  Info.ndim = 0;
543  // driven phase [degree]
544  return true;
545  }
546  return StateBase::getArray(idx-I, Info);
547 }
548 
549 MomentElementBase::MomentElementBase(const Config& c)
550  :ElementVoid(c)
551  ,dx (c.get<double>("dx", 0e0)*MtoMM)
552  ,dy (c.get<double>("dy", 0e0)*MtoMM)
553  ,pitch(c.get<double>("pitch", 0e0))
554  ,yaw (c.get<double>("yaw", 0e0))
555  ,roll (c.get<double>("roll", 0e0))
556  ,skipcache(c.get<double>("skipcache", 0.0)!=0.0)
557  ,scratch(state_t::maxsize, state_t::maxsize)
558 {
559 }
560 
561 MomentElementBase::~MomentElementBase() {}
562 
563 void MomentElementBase::assign(const ElementVoid *other)
564 {
565  const MomentElementBase *O = static_cast<const MomentElementBase*>(other);
566  last_real_in = O->last_real_in;
567  last_real_out= O->last_real_out;
568  last_ref_in = O->last_ref_in;
569  last_ref_out = O->last_ref_out;
570  transfer = O->transfer;
571  misalign = O->misalign;
572  misalign_inv = O->misalign_inv;
573  dx = O->dx;
574  dy = O->dy;
575  pitch = O->pitch;
576  yaw = O->yaw;
577  roll = O->roll;
578  skipcache = O->skipcache;
579  ElementVoid::assign(other);
580 }
581 
582 void MomentElementBase::show(std::ostream& strm, int level) const
583 {
584  using namespace boost::numeric::ublas;
585  ElementVoid::show(strm, level);
586  /*
587  strm<<"Length "<<length<<"\n"
588  "Transfer: "<<transfer<<"\n"
589  "Mis-align: "<<misalign<<"\n";
590  */
591 }
592 
593 void MomentElementBase::get_misalign(const state_t &ST, const Particle &real, value_t &M, value_t &IM) const
594 {
595  state_t::matrix_t R,
596  scl = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize),
597  scl_inv = scl,
598  T = scl,
599  T_inv = scl,
600  R_inv = scl;
601 
602  scl(state_t::PS_S, state_t::PS_S) /= -real.SampleIonK;
603  scl(state_t::PS_PS, state_t::PS_PS) /= sqr(real.beta)*real.gamma*ST.ref.IonEs/MeVtoeV;
604 
605  inverse(scl_inv, scl);
606 
607  // Translate to center of element.
608  T(state_t::PS_S, 6) = -length/2e0*MtoMM;
609  T(state_t::PS_PS, 6) = 1e0;
610  inverse(T_inv, T);
611 
612  RotMat(dx, dy, pitch, yaw, roll, R);
613 
614  M = prod(T, scl);
615  M = prod(R, M);
616  M = prod(T_inv, M);
617  M = prod(scl_inv, M);
618 
619  inverse(R_inv, R);
620 
621  // Translate to center of element.
622  T(state_t::PS_S, 6) = length/2e0*MtoMM;
623  T(state_t::PS_PS, 6) = 1e0;
624  inverse(T_inv, T);
625 
626  IM = prod(T, scl);
627  IM = prod(R_inv, IM);
628  IM = prod(T_inv, IM);
629  IM = prod(scl_inv, IM);
630 }
631 
632 unsigned MomentElementBase::get_flag(const Config& c, const std::string& name, const unsigned& def_value)
633 {
634  unsigned read_value;
635  double check_value;
636 
637  try {
638  check_value = boost::lexical_cast<double>(c.get<std::string>(name));
639  }catch (std::exception&){
640  try {
641  check_value = c.get<double>(name);
642  }catch (std::exception&){
643  check_value = boost::lexical_cast<double>(def_value);
644  }
645  }
646 
647  //Check the value is an unsigned integer
648  try {
649  read_value = boost::lexical_cast<unsigned>(check_value);
650  if (boost::lexical_cast<double>(read_value) != check_value)
651  throw std::runtime_error(SB()<< name << " must be an unsigned integer");
652  }catch (std::exception&){
653  throw std::runtime_error(SB()<< name << " must be an unsigned integer");
654  }
655  return read_value;
656 }
657 
658 void MomentElementBase::advance(StateBase& s)
659 {
660  state_t& ST = static_cast<state_t&>(s);
661  using namespace boost::numeric::ublas;
662 
663  // IonEk is Es + E_state; the latter is set by user.
664  ST.recalc();
665 
666  if(!check_cache(ST)){
667  // need to re-calculate energy dependent terms
668  last_ref_in = ST.ref;
669  last_real_in = ST.real;
670  resize_cache(ST);
671 
672  recompute_matrix(ST); // updates transfer and last_Kenergy_out
673 
674  ST.recalc();
675 
676  if(!ST.retreat){
677  ST.ref.phis += ST.ref.SampleIonK*length*MtoMM;
678  for(size_t k=0; k<last_real_in.size(); k++)
679  ST.real[k].phis += ST.real[k].SampleIonK*length*MtoMM;
680  } else {
681  ST.ref.phis -= ST.ref.SampleIonK*length*MtoMM;
682  for(size_t k=0; k<last_real_in.size(); k++)
683  ST.real[k].phis -= ST.real[k].SampleIonK*length*MtoMM;
684  }
685 
686  last_ref_out = ST.ref;
687  last_real_out = ST.real;
688  } else {
689  ST.ref = last_ref_out;
690  assert(last_real_out.size()==ST.real.size()); // should be true if check_cache() -> true
691  std::copy(last_real_out.begin(),
692  last_real_out.end(),
693  ST.real.begin());
694  }
695 
696  if(!ST.retreat){
697  // Forward propagation
698  ST.pos += length;
699 
700  for(size_t k=0; k<last_real_in.size(); k++) {
701  ST.moment0[k] = prod(transfer[k], ST.moment0[k]);
702 
703  scratch = prod(transfer[k], ST.moment1[k]);
704  ST.moment1[k] = prod(scratch, trans(transfer[k]));
705 
706  ST.transmat[k] = transfer[k];
707  }
708  } else {
709  // Backward propagation
710  ST.pos -= length;
711 
712  value_t invmat = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
713  for(size_t k=0; k<last_real_in.size(); k++) {
714  inverse(invmat, transfer[k]);
715 
716  ST.moment0[k] = prod(invmat, ST.moment0[k]);
717 
718  scratch = prod(invmat, ST.moment1[k]);
719  ST.moment1[k] = prod(scratch, trans(invmat));
720 
721  ST.transmat[k] = invmat;
722  }
723  }
724 
725  ST.calc_rms();
726 }
727 
728 bool MomentElementBase::check_cache(const state_t& ST) const
729 {
730  return !skipcache
731  && last_real_in.size()==ST.size()
732  && last_ref_in==ST.ref
733  && std::equal(last_real_in.begin(),
734  last_real_in.end(),
735  ST.real.begin());
736 }
737 
738 bool MomentElementBase::check_backward(const state_t& ST) const
739 {
740  bool reals = true;
741  if (last_real_out.size()==ST.size()) {
742  reals = last_ref_out<=ST.ref;
743  for(size_t k=0; k<last_real_out.size(); k++) {
744  reals &= last_real_out[k]<=ST.real[k];
745  }
746  } else {
747  reals = false;
748  }
749 
750  return reals;
751 }
752 
753 void MomentElementBase::resize_cache(const state_t& ST)
754 {
755  transfer.resize(ST.real.size(), boost::numeric::ublas::identity_matrix<double>(state_t::maxsize));
756  misalign.resize(ST.real.size(), boost::numeric::ublas::identity_matrix<double>(state_t::maxsize));
757  misalign_inv.resize(ST.real.size(), boost::numeric::ublas::identity_matrix<double>(state_t::maxsize));
758 }
759 
760 void MomentElementBase::recompute_matrix(state_t& ST)
761 {
762  // Default, initialize as no-op
763 
764  for(size_t k=0; k<last_real_in.size(); k++) {
765  transfer[k] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
766  }
767 }
768 
769 namespace {
770 
771 struct ElementSource : public MomentElementBase
772 {
773  typedef ElementSource self_t;
774  typedef MomentElementBase base_t;
775  typedef typename base_t::state_t state_t;
776 
777  ElementSource(const Config& c): base_t(c), istate(c) {}
778 
779  virtual void advance(StateBase& s)
780  {
781  state_t& ST = static_cast<state_t&>(s);
782  if (!ST.retreat)
783  // Replace state with our initial values
784  ST.assign(istate);
785  }
786 
787  virtual void show(std::ostream& strm, int level) const
788  {
789  ElementVoid::show(strm, level);
790  strm<<"Initial: "<<istate.moment0_env<<"\n";
791  }
792 
793  state_t istate;
794  // note that 'transfer' is not used by this element type
795 
796  virtual ~ElementSource() {}
797 
798  virtual const char* type_name() const {return "source";}
799 
800  virtual void assign(const ElementVoid *other) {
801  base_t::assign(other);
802  const self_t* O=static_cast<const self_t*>(other);
803  istate.assign(O->istate);
804  }
805 };
806 
807 struct ElementMark : public MomentElementBase
808 {
809  // Transport (identity) matrix for Marker.
810  typedef ElementMark self_t;
811  typedef MomentElementBase base_t;
812  typedef typename base_t::state_t state_t;
813 
814  ElementMark(const Config& c): base_t(c) {length = 0e0;}
815  virtual ~ElementMark() {}
816  virtual const char* type_name() const {return "marker";}
817 
818  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
819 };
820 
821 struct ElementBPM : public MomentElementBase
822 {
823  // Transport (identity) matrix for BPM.
824  typedef ElementBPM self_t;
825  typedef MomentElementBase base_t;
826  typedef typename base_t::state_t state_t;
827 
828  ElementBPM(const Config& c): base_t(c) {length = 0e0;}
829  virtual ~ElementBPM() {}
830  virtual const char* type_name() const {return "bpm";}
831 
832  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
833 };
834 
835 struct ElementDrift : public MomentElementBase
836 {
837  // Transport matrix for Drift.
838  typedef ElementDrift self_t;
839  typedef MomentElementBase base_t;
840  typedef typename base_t::state_t state_t;
841 
842  ElementDrift(const Config& c) : base_t(c) {}
843  virtual ~ElementDrift() {}
844  virtual const char* type_name() const {return "drift";}
845 
846  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
847 
848  virtual void recompute_matrix(state_t& ST)
849  {
850  // Re-initialize transport matrix.
851 
852  const double L = length*MtoMM; // Convert from [m] to [mm].
853 
854  for(size_t i=0; i<last_real_in.size(); i++) {
855  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
856  transfer[i](state_t::PS_X, state_t::PS_PX) = L;
857  transfer[i](state_t::PS_Y, state_t::PS_PY) = L;
858  transfer[i](state_t::PS_S, state_t::PS_PS) =
859  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*L;
860  }
861  }
862 };
863 
864 struct ElementOrbTrim : public MomentElementBase
865 {
866  // Transport matrix for Orbit Trim.
867  typedef ElementOrbTrim self_t;
868  typedef MomentElementBase base_t;
869  typedef typename base_t::state_t state_t;
870 
871  ElementOrbTrim(const Config& c) : base_t(c) {length = 0e0;}
872  virtual ~ElementOrbTrim() {}
873  virtual const char* type_name() const {return "orbtrim";}
874 
875  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
876 
877  virtual void recompute_matrix(state_t& ST)
878  {
879  // Re-initialize transport matrix.
880  double theta_x = conf().get<double>("theta_x", 0e0),
881  theta_y = conf().get<double>("theta_y", 0e0);
882  const double tm_xkick = conf().get<double>("tm_xkick", 0e0),
883  tm_ykick = conf().get<double>("tm_ykick", 0e0),
884  xyrotate = conf().get<double>("xyrotate", 0e0)*M_PI/180e0;
885  const bool realpara = conf().get<double>("realpara", 0e0) == 1e0;
886 
887  if (realpara) {
888  double ecpi = ST.ref.IonZ*C0/sqrt(sqr(ST.ref.IonW) - sqr(ST.ref.IonEs));
889  theta_x = tm_xkick*ecpi;
890  theta_y = tm_ykick*ecpi;
891  }
892 
893  for(size_t i=0; i<last_real_in.size(); i++) {
894  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
895  transfer[i](state_t::PS_PX, 6) = theta_x*ST.real[i].IonZ/ST.ref.IonZ;
896  transfer[i](state_t::PS_PY, 6) = theta_y*ST.real[i].IonZ/ST.ref.IonZ;
897 
898  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
899 
900  noalias(scratch) = prod(transfer[i], misalign[i]);
901  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
902 
903  if (xyrotate != 0e0) {
904  state_t::matrix_t R;
905  RotMat(0e0, 0e0, 0e0, 0e0, xyrotate, R);
906  noalias(scratch) = transfer[i];
907  noalias(transfer[i]) = prod(scratch, R);
908  }
909 
910  }
911  }
912 };
913 
914 struct ElementSBend : public MomentElementBase
915 {
916  // Transport matrix for Gradient Sector Bend; with edge focusing (cylindrical coordinates).
917  // Note, TLM only includes energy offset for the orbit; not the transport matrix.
918 
919  typedef ElementSBend self_t;
920  typedef MomentElementBase base_t;
921  typedef typename base_t::state_t state_t;
922 
923  unsigned HdipoleFitMode;
924 
925  ElementSBend(const Config& c) : base_t(c), HdipoleFitMode(0) {
926 
927  HdipoleFitMode = get_flag(c, "HdipoleFitMode", 1);
928  if (HdipoleFitMode != 0 && HdipoleFitMode != 1)
929  throw std::runtime_error(SB()<< "Undefined HdipoleFitMode: " << HdipoleFitMode);
930  }
931  virtual ~ElementSBend() {}
932  virtual const char* type_name() const {return "sbend";}
933 
934  virtual void assign(const ElementVoid *other) {
935  base_t::assign(other);
936  const self_t* O=static_cast<const self_t*>(other);
937  HdipoleFitMode = O->HdipoleFitMode;
938  }
939 
940  virtual void advance(StateBase& s)
941  {
942  state_t& ST = static_cast<state_t&>(s);
943  using namespace boost::numeric::ublas;
944 
945  // IonEk is Es + E_state; the latter is set by user.
946  ST.recalc();
947 
948  if(!check_cache(ST)) {
949  // need to re-calculate energy dependent terms
950  last_ref_in = ST.ref;
951  last_real_in = ST.real;
952  resize_cache(ST);
953 
954  recompute_matrix(ST); // updates transfer and last_Kenergy_out
955 
956  ST.recalc();
957  last_ref_out = ST.ref;
958  last_real_out = ST.real;
959  } else {
960  ST.ref = last_ref_out;
961  assert(last_real_out.size()==ST.real.size()); // should be true if check_cache() -> true
962  std::copy(last_real_out.begin(),
963  last_real_out.end(),
964  ST.real.begin());
965  }
966 
967  if(!ST.retreat){
968  // Forward propagation
969  ST.pos += length;
970  ST.ref.phis += ST.ref.SampleIonK*length*MtoMM;
971 
972  for(size_t i=0; i<last_real_in.size(); i++) {
973  double phis_temp = ST.moment0[i][state_t::PS_S];
974 
975  ST.moment0[i] = prod(transfer[i], ST.moment0[i]);
976 
977  noalias(scratch) = prod(transfer[i], ST.moment1[i]);
978  noalias(ST.moment1[i]) = prod(scratch, trans(transfer[i]));
979 
980  double dphis_temp = ST.moment0[i][state_t::PS_S] - phis_temp;
981 
982  ST.real[i].phis += ST.real[i].SampleIonK*length*MtoMM + dphis_temp;
983  ST.transmat[i] = transfer[i];
984  }
985  } else {
986  // Backward propagation
987  ST.pos -= length;
988  ST.ref.phis -= ST.ref.SampleIonK*length*MtoMM;
989 
990  value_t invmat = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
991  for(size_t i=0; i<last_real_in.size(); i++) {
992  double phis_temp = ST.moment0[i][state_t::PS_S];
993 
994  inverse(invmat, transfer[i]);
995  ST.moment0[i] = prod(invmat, ST.moment0[i]);
996 
997  noalias(scratch) = prod(invmat, ST.moment1[i]);
998  noalias(ST.moment1[i]) = prod(scratch, trans(invmat));
999 
1000  double dphis_temp = ST.moment0[i][state_t::PS_S] - phis_temp;
1001 
1002  ST.real[i].phis -= ST.real[i].SampleIonK*length*MtoMM - dphis_temp;
1003  ST.transmat[i] = invmat;
1004  }
1005  }
1006 
1007  ST.calc_rms();
1008  }
1009 
1010  virtual void recompute_matrix(state_t& ST)
1011  {
1012  // Re-initialize transport matrix.
1013  bool ver = conf().get<double>("ver", 0e0) == 1.0;
1014  double L = conf().get<double>("L")*MtoMM,
1015  phi = conf().get<double>("phi")*M_PI/180e0,
1016  phi1 = conf().get<double>("phi1")*M_PI/180e0,
1017  phi2 = conf().get<double>("phi2")*M_PI/180e0,
1018  dphi1 = conf().get<double>("dphi1", 0e0)*M_PI/180e0,
1019  dphi2 = conf().get<double>("dphi2", 0e0)*M_PI/180e0,
1020  K = conf().get<double>("K", 0e0)/sqr(MtoMM);
1021 
1022  unsigned EFcorrection = get_flag(conf(), "EFcorrection", 0);
1023 
1024  if (EFcorrection != 0 && EFcorrection != 1)
1025  throw std::runtime_error(SB()<< "Undefined EFcorrection: " << EFcorrection);
1026 
1027  for(size_t i=0; i<last_real_in.size(); i++) {
1028 
1029  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1030 
1031  if (L != 0.0) {
1032  if (!HdipoleFitMode) {
1033  double dip_bg = conf().get<double>("bg", ST.ref.bg),
1034  dip_IonZ = conf().get<double>("ref_IonZ", ST.ref.IonZ),
1035  qmrel = (ST.real[i].IonZ-dip_IonZ)/dip_IonZ,
1036  // Dipole reference energy.
1037  dip_Ek = (sqrt(sqr(dip_bg)+1e0)-1e0)*ST.ref.IonEs,
1038  dip_gamma = (dip_Ek+ST.ref.IonEs)/ST.ref.IonEs,
1039  dip_beta = sqrt(1e0-1e0/sqr(dip_gamma)),
1040  d = (ST.ref.gamma-dip_gamma)/(sqr(dip_beta)*dip_gamma) - qmrel,
1041  dip_IonK = 2e0*M_PI/(dip_beta*ST.ref.SampleLambda);
1042 
1043  GetSBendMatrix(L, phi, phi1, phi2, K, ST.ref.IonEs, ST.ref.gamma, qmrel,
1044  dphi1, dphi2, EFcorrection, dip_beta, dip_gamma, d, dip_IonK, transfer[i]);
1045  } else {
1046  double qmrel = (ST.real[i].IonZ-ST.ref.IonZ)/ST.ref.IonZ;
1047  GetSBendMatrix(L, phi, phi1, phi2, K, ST.ref.IonEs, ST.ref.gamma, qmrel,
1048  dphi1, dphi2, EFcorrection, ST.ref.beta, ST.ref.gamma, - qmrel,
1049  ST.ref.SampleIonK, transfer[i]);
1050  }
1051 
1052  if (ver) {
1053  // Rotate transport matrix by 90 degrees.
1054  value_t
1055  R = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1056  R(state_t::PS_X, state_t::PS_X) = 0e0;
1057  R(state_t::PS_PX, state_t::PS_PX) = 0e0;
1058  R(state_t::PS_Y, state_t::PS_Y) = 0e0;
1059  R(state_t::PS_PY, state_t::PS_PY) = 0e0;
1060  R(state_t::PS_X, state_t::PS_Y) = -1e0;
1061  R(state_t::PS_PX, state_t::PS_PY) = -1e0;
1062  R(state_t::PS_Y, state_t::PS_X) = 1e0;
1063  R(state_t::PS_PY, state_t::PS_PX) = 1e0;
1064 
1065  noalias(scratch) = prod(R, transfer[i]);
1066  noalias(transfer[i]) = prod(scratch, trans(R));
1067  //TODO: no-op code? results are unconditionally overwritten
1068  }
1069 
1070  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1071 
1072  noalias(scratch) = prod(transfer[i], misalign[i]);
1073  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1074  }
1075  }
1076  }
1077 };
1078 
1079 struct ElementQuad : public MomentElementBase
1080 {
1081  // Transport matrix for Quadrupole; K = B2/Brho.
1082  typedef ElementQuad self_t;
1083  typedef MomentElementBase base_t;
1084  typedef typename base_t::state_t state_t;
1085 
1086  ElementQuad(const Config& c) : base_t(c) {}
1087  virtual ~ElementQuad() {}
1088  virtual const char* type_name() const {return "quadrupole";}
1089 
1090  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
1091 
1092  virtual void recompute_matrix(state_t& ST)
1093  {
1094  const double L = conf().get<double>("L")*MtoMM;
1095  const unsigned ncurve = get_flag(conf(), "ncurve", 0);
1096 
1097  if (ncurve != 0) {
1098  std::vector<std::vector<double> > Curves;
1099  std::vector<double> Scales;
1100  GetCurveData(conf(), ncurve, Scales, Curves);
1101 
1102  for(size_t i=0; i<last_real_in.size(); i++) {
1103  double K;
1104  double dL = L/double(Curves[0].size()),
1105  Brho = ST.real[i].Brho();
1106  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1107  for (size_t j=0; j<Curves[0].size(); j++){
1108  K = 0.0;
1109  for (size_t n=0; n<Curves.size(); n++) K += Scales[n]*Curves[n][j]/Brho/sqr(MtoMM);
1110  value_t tmstep = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1111  GetQuadMatrix(dL, K, (unsigned)state_t::PS_X, tmstep);
1112  GetQuadMatrix(dL, -K, (unsigned)state_t::PS_Y, tmstep);
1113 
1114  tmstep(state_t::PS_S, state_t::PS_PS) =
1115  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*dL;
1116 
1117  transfer[i] = prod(tmstep, transfer[i]);
1118  }
1119  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1120  noalias(scratch) = prod(transfer[i], misalign[i]);
1121  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1122  }
1123 
1124  } else {
1125  const double B2= conf().get<double>("B2");
1126  for(size_t i=0; i<last_real_in.size(); i++) {
1127  // Re-initialize transport matrix.
1128  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1129 
1130  double Brho = ST.real[i].Brho(),
1131  K = B2/Brho/sqr(MtoMM);
1132 
1133  // Horizontal plane.
1134  GetQuadMatrix(L, K, (unsigned)state_t::PS_X, transfer[i]);
1135  // Vertical plane.
1136  GetQuadMatrix(L, -K, (unsigned)state_t::PS_Y, transfer[i]);
1137  // Longitudinal plane.
1138 
1139  transfer[i](state_t::PS_S, state_t::PS_PS) =
1140  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*L;
1141 
1142  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1143  noalias(scratch) = prod(transfer[i], misalign[i]);
1144  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1145  }
1146  }
1147  }
1148 };
1149 
1150 struct ElementSext : public MomentElementBase
1151 {
1152  // Transport matrix for Sextupole; K = B3/Brho.
1153  typedef ElementSext self_t;
1154  typedef MomentElementBase base_t;
1155  typedef typename base_t::state_t state_t;
1156 
1157  ElementSext(const Config& c) : base_t(c) {}
1158 
1159  virtual ~ElementSext() {}
1160  virtual const char* type_name() const {return "sextupole";}
1161 
1162  virtual void assign(const ElementVoid *other) {base_t::assign(other); }
1163 
1164  virtual void advance(StateBase& s)
1165  {
1166  const double B3= conf().get<double>("B3"),
1167  L = conf().get<double>("L")*MtoMM;
1168  const int step = conf().get<double>("step", 1.0);
1169  const bool thinlens = conf().get<double>("thinlens", 0.0) == 1.0,
1170  dstkick = conf().get<double>("dstkick", 1.0) == 1.0;
1171 
1172  state_t& ST = static_cast<state_t&>(s);
1173  using namespace boost::numeric::ublas;
1174  value_t invmat = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1175 
1176  ST.recalc();
1177 
1178  last_ref_in = ST.ref;
1179  last_real_in = ST.real;
1180  resize_cache(ST);
1181 
1182  //if(ST.retreat) throw std::runtime_error(SB()<<
1183  // "Backward propagation error: Backward propagation does not support sextupole.");
1184 
1185  const double dL = L/step;
1186 
1187  for(size_t k=0; k<last_real_in.size(); k++) {
1188 
1189  transfer[k] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1190  ST.transmat[k] = transfer[k];
1191 
1192  double Brho = ST.real[k].Brho(),
1193  K = B3/Brho/cube(MtoMM);
1194 
1195  get_misalign(ST, ST.real[k], misalign[k], misalign_inv[k]);
1196 
1197  if(!ST.retreat){
1198  ST.moment0[k] = prod(misalign[k], ST.moment0[k]);
1199  noalias(scratch) = prod(misalign[k], ST.moment1[k]);
1200  ST.moment1[k] = prod(scratch, trans(misalign[k]));
1201  } else {
1202  inverse(invmat, misalign_inv[k]);
1203  ST.moment0[k] = prod(invmat, ST.moment0[k]);
1204  noalias(scratch) = prod(invmat, ST.moment1[k]);
1205  ST.moment1[k] = prod(scratch, trans(invmat));
1206  }
1207 
1208  for(int i=0; i<step; i++){
1209  double Dx = ST.moment0[k][state_t::PS_X],
1210  Dy = ST.moment0[k][state_t::PS_Y],
1211  D2x = ST.moment1[k](state_t::PS_X, state_t::PS_X),
1212  D2y = ST.moment1[k](state_t::PS_Y, state_t::PS_Y),
1213  D2xy = ST.moment1[k](state_t::PS_X, state_t::PS_Y);
1214 
1215  GetSextMatrix(dL, K, Dx, Dy, D2x, D2y, D2xy, thinlens, dstkick, transfer[k]);
1216 
1217  transfer[k](state_t::PS_S, state_t::PS_PS) =
1218  -2e0*M_PI/(ST.real[k].SampleLambda*ST.real[k].IonEs/MeVtoeV*cube(ST.real[k].bg))*dL;
1219 
1220  if(!ST.retreat){
1221  ST.moment0[k] = prod(transfer[k], ST.moment0[k]);
1222  noalias(scratch) = prod(transfer[k], ST.moment1[k]);
1223  ST.moment1[k] = prod(scratch, trans(transfer[k]));
1224  ST.transmat[k] = prod(transfer[k], ST.transmat[k]);
1225  } else {
1226  inverse(invmat, transfer[k]);
1227  ST.moment0[k] = prod(invmat, ST.moment0[k]);
1228  noalias(scratch) = prod(invmat, ST.moment1[k]);
1229  ST.moment1[k] = prod(scratch, trans(invmat));
1230  ST.transmat[k] = prod(invmat, ST.transmat[k]);
1231  }
1232  }
1233 
1234  if(!ST.retreat){
1235  ST.moment0[k] = prod(misalign_inv[k], ST.moment0[k]);
1236  noalias(scratch) = prod(misalign_inv[k], ST.moment1[k]);
1237  ST.moment1[k] = prod(scratch, trans(misalign_inv[k]));
1238  noalias(scratch) = prod(ST.transmat[k], misalign[k]);
1239  ST.transmat[k] = prod(misalign_inv[k], scratch);
1240  } else {
1241  inverse(invmat, misalign_inv[k]);
1242  noalias(scratch) = prod(ST.transmat[k], invmat);
1243  inverse(invmat, misalign[k]);
1244  ST.transmat[k] = prod(invmat, scratch);
1245  ST.moment0[k] = prod(invmat, ST.moment0[k]);
1246  noalias(scratch) = prod(invmat, ST.moment1[k]);
1247  ST.moment1[k] = prod(scratch, trans(invmat));
1248  }
1249  }
1250 
1251  ST.recalc();
1252 
1253  if(!ST.retreat){
1254  for(size_t k=0; k<last_real_in.size(); k++)
1255  ST.real[k].phis += ST.real[k].SampleIonK*length*MtoMM;
1256  ST.ref.phis += ST.ref.SampleIonK*length*MtoMM;
1257  ST.pos += length;
1258  } else {
1259  for(size_t k=0; k<last_real_in.size(); k++)
1260  ST.real[k].phis -= ST.real[k].SampleIonK*length*MtoMM;
1261  ST.ref.phis -= ST.ref.SampleIonK*length*MtoMM;
1262  ST.pos -= length;
1263  }
1264 
1265  last_ref_out = ST.ref;
1266  last_real_out = ST.real;
1267 
1268  ST.calc_rms();
1269  }
1270 };
1271 
1272 struct ElementSolenoid : public MomentElementBase
1273 {
1274  // Transport (identity) matrix for a Solenoid; K = B/(2 Brho).
1275  typedef ElementSolenoid self_t;
1276  typedef MomentElementBase base_t;
1277  typedef typename base_t::state_t state_t;
1278 
1279  ElementSolenoid(const Config& c) : base_t(c) {}
1280  virtual ~ElementSolenoid() {}
1281  virtual const char* type_name() const {return "solenoid";}
1282 
1283  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
1284 
1285  virtual void recompute_matrix(state_t& ST)
1286  {
1287  const double L = conf().get<double>("L")*MtoMM; // Convert from [m] to [mm].
1288  const unsigned ncurve = get_flag(conf(), "ncurve", 0);
1289 
1290  if (ncurve != 0) {
1291  std::vector<std::vector<double> > Curves;
1292  std::vector<double> Scales;
1293  GetCurveData(conf(), ncurve, Scales, Curves);
1294 
1295  for(size_t i=0; i<last_real_in.size(); i++) {
1296  double K;
1297  double dL = L/double(Curves[0].size()),
1298  Brho = ST.real[i].Brho();
1299  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1300  for (size_t j=0; j<Curves[0].size(); j++){
1301  K = 0.0;
1302  for (size_t n=0; n<Curves.size(); n++) K += Scales[n]*Curves[n][j]/(2e0*Brho)/MtoMM;
1303  value_t tmstep = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1304 
1305  GetSolMatrix(dL, K, tmstep);
1306 
1307  tmstep(state_t::PS_S, state_t::PS_PS) =
1308  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*dL;
1309 
1310  transfer[i] = prod(tmstep, transfer[i]);
1311  }
1312  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1313  noalias(scratch) = prod(transfer[i], misalign[i]);
1314  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1315  }
1316  } else {
1317  const double B = conf().get<double>("B");
1318  for(size_t i=0; i<last_real_in.size(); i++) {
1319  // Re-initialize transport matrix.
1320  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1321 
1322  double Brho = ST.real[i].Brho(),
1323  K = B/(2e0*Brho)/MtoMM;
1324 
1325  GetSolMatrix(L, K, transfer[i]);
1326 
1327  transfer[i](state_t::PS_S, state_t::PS_PS) =
1328  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*L;
1329 
1330  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1331 
1332  noalias(scratch) = prod(transfer[i], misalign[i]);
1333  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1334  }
1335  }
1336  }
1337 };
1338 
1339 struct ElementEDipole : public MomentElementBase
1340 {
1341  // Transport matrix for Electrostatic Dipole with edge focusing.
1342  typedef ElementEDipole self_t;
1343  typedef MomentElementBase base_t;
1344  typedef typename base_t::state_t state_t;
1345 
1346  ElementEDipole(const Config& c) : base_t(c) {}
1347  virtual ~ElementEDipole() {}
1348  virtual const char* type_name() const {return "edipole";}
1349 
1350  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
1351 
1352  virtual void recompute_matrix(state_t& ST)
1353  {
1354  // Re-initialize transport matrix.
1355 
1356  //value_mat R;
1357 
1358  bool ver = conf().get<double>("ver", 0e0) == 1.0;
1359  double L = conf().get<double>("L")*MtoMM,
1360  phi = conf().get<double>("phi")*M_PI/180e0,
1361  // fit to TLM unit.
1362  fringe_x = conf().get<double>("fringe_x", 0e0)/MtoMM,
1363  fringe_y = conf().get<double>("fringe_y", 0e0)/MtoMM,
1364  kappa = conf().get<double>("asym_fac", 0e0),
1365  // spher: cylindrical - 0, spherical - 1.
1366  spher = conf().get<double>("spher"),
1367  // magnetic - 0, electrostatic - 1.
1368  h = 1e0,
1369  dip_beta = conf().get<double>("beta", ST.ref.beta);
1370 
1371  unsigned HdipoleFitMode = get_flag(conf(), "HdipoleFitMode", 1);
1372 
1373  if (HdipoleFitMode != 0 && HdipoleFitMode != 1)
1374  throw std::runtime_error(SB()<< "Undefined HdipoleFitMode: " << HdipoleFitMode);
1375 
1376  if (HdipoleFitMode) dip_beta = ST.ref.beta;
1377 
1378  for(size_t i=0; i<last_real_in.size(); i++) {
1379  double eta0 = (1e0/sqrt(1e0 - sqr(dip_beta)) - 1e0)/2e0,
1380  Erho = sqr(ST.real[i].beta)/ST.real[i].IonZ,
1381  Erho0 = sqr(dip_beta)/ST.ref.IonZ,
1382  eL = Erho/Erho0*L,
1383  rho = eL/phi,
1384  Kx = (1e0-spher+sqr(1e0+2e0*eta0))/sqr(rho),
1385  Ky = spher/sqr(rho),
1386  delta_K = (Erho/Erho0 - 1e0) - (ST.real[i].IonEk - ST.ref.IonEk)/ST.real[i].IonEk,
1387  delta_KZ = ST.ref.IonZ/ST.real[i].IonZ - 1e0,
1388  SampleIonK = 2e0*M_PI/(ST.real[i].beta*ST.real[i].SampleLambda);
1389 
1390  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1391 
1392  if (L != 0e0) {
1393  GetEBendMatrix(eL, phi, fringe_x, fringe_y, kappa, Kx, Ky, ST.ref.IonEs, ST.real[i].gamma,
1394  eta0, h, delta_K, delta_KZ, SampleIonK, transfer[i]);
1395 
1396  if (ver) {
1397  // Rotate transport matrix by 90 degrees.
1398  value_t
1399  R = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1400  R(state_t::PS_X, state_t::PS_X) = 0e0;
1401  R(state_t::PS_PX, state_t::PS_PX) = 0e0;
1402  R(state_t::PS_Y, state_t::PS_Y) = 0e0;
1403  R(state_t::PS_PY, state_t::PS_PY) = 0e0;
1404  R(state_t::PS_X, state_t::PS_Y) = -1e0;
1405  R(state_t::PS_PX, state_t::PS_PY) = -1e0;
1406  R(state_t::PS_Y, state_t::PS_X) = 1e0;
1407  R(state_t::PS_PY, state_t::PS_PX) = 1e0;
1408 
1409  noalias(scratch) = prod(R, transfer[i]);
1410  noalias(transfer[i]) = prod(scratch, trans(R));
1411  //TODO: no-op code? results are unconditionally overwritten
1412  }
1413 
1414  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1415 
1416  noalias(scratch) = prod(transfer[i], misalign[i]);
1417  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1418  }
1419  }
1420  }
1421 };
1422 
1423 struct ElementEQuad : public MomentElementBase
1424 {
1425  // Transport matrix for Electrostatic Quadrupole.
1426  typedef ElementEQuad self_t;
1427  typedef MomentElementBase base_t;
1428  typedef typename base_t::state_t state_t;
1429 
1430  ElementEQuad(const Config& c) : base_t(c) {}
1431  virtual ~ElementEQuad() {}
1432  virtual const char* type_name() const {return "equad";}
1433 
1434  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
1435 
1436  virtual void recompute_matrix(state_t& ST)
1437  {
1438  const double L = conf().get<double>("L")*MtoMM;
1439  const unsigned ncurve = get_flag(conf(), "ncurve", 0);
1440 
1441  if (ncurve != 0) {
1442  std::vector<std::vector<double> > Curves;
1443  std::vector<double> Scales;
1444  GetCurveData(conf(), ncurve, Scales, Curves);
1445 
1446  for(size_t i=0; i<last_real_in.size(); i++) {
1447  double K;
1448  double dL = L/double(Curves[0].size()),
1449  Brho = ST.real[i].Brho();
1450  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1451  for (size_t j=0; j<Curves[0].size(); j++){
1452  K = 0.0;
1453  for (size_t n=0; n<Curves.size(); n++) K += 2e0*Scales[n]*Curves[n][j]/(C0*ST.real[i].beta)/Brho/sqr(MtoMM);
1454  value_t tmstep = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1455  GetQuadMatrix(dL, K, (unsigned)state_t::PS_X, tmstep);
1456  GetQuadMatrix(dL, -K, (unsigned)state_t::PS_Y, tmstep);
1457 
1458  tmstep(state_t::PS_S, state_t::PS_PS) =
1459  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*dL;
1460 
1461  transfer[i] = prod(tmstep, transfer[i]);
1462  }
1463  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1464  noalias(scratch) = prod(transfer[i], misalign[i]);
1465  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1466  }
1467 
1468  } else {
1469  const double V0 = conf().get<double>("V"),
1470  R = conf().get<double>("radius");
1471 
1472  for(size_t i=0; i<last_real_in.size(); i++) {
1473  // Re-initialize transport matrix.
1474  // V0 [V] electrode voltage and R [m] electrode half-distance.
1475  transfer[i] = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
1476 
1477  double Brho = ST.real[i].Brho(),
1478  K = 2e0*V0/(C0*ST.real[i].beta*sqr(R))/Brho/sqr(MtoMM);
1479 
1480  // Horizontal plane.
1481  GetQuadMatrix(L, K, (unsigned)state_t::PS_X, transfer[i]);
1482  // Vertical plane.
1483  GetQuadMatrix(L, -K, (unsigned)state_t::PS_Y, transfer[i]);
1484  // Longitudinal plane.
1485  // transfer(state_t::PS_S, state_t::PS_S) = L;
1486 
1487  transfer[i](state_t::PS_S, state_t::PS_PS) =
1488  -2e0*M_PI/(ST.real[i].SampleLambda*ST.real[i].IonEs/MeVtoeV*cube(ST.real[i].bg))*L;
1489 
1490  get_misalign(ST, ST.real[i], misalign[i], misalign_inv[i]);
1491 
1492  noalias(scratch) = prod(transfer[i], misalign[i]);
1493  noalias(transfer[i]) = prod(misalign_inv[i], scratch);
1494  }
1495  }
1496  }
1497 };
1498 
1499 struct ElementTMatrix : public MomentElementBase
1500 {
1501  // Transport matrix by user input
1502  typedef ElementTMatrix self_t;
1503  typedef MomentElementBase base_t;
1504  typedef typename base_t::state_t state_t;
1505 
1506  ElementTMatrix(const Config& c) : base_t(c) {}
1507  virtual ~ElementTMatrix() {}
1508  virtual const char* type_name() const {return "tmatrix";}
1509 
1510  virtual void assign(const ElementVoid *other) { base_t::assign(other); }
1511 
1512  virtual void recompute_matrix(state_t& ST)
1513  {
1514  for(size_t i=0; i<last_real_in.size(); i++) {
1515  load_storage(transfer[i].data(), conf(), "matrix");
1516  }
1517  }
1518 };
1519 
1520 } // namespace
1521 
1522 void registerMoment()
1523 {
1524  Machine::registerState<MomentState>("MomentMatrix");
1525 
1526  Machine::registerElement<ElementSource >("MomentMatrix", "source");
1527 
1528  Machine::registerElement<ElementMark >("MomentMatrix", "marker");
1529 
1530  Machine::registerElement<ElementBPM >("MomentMatrix", "bpm");
1531 
1532  Machine::registerElement<ElementDrift >("MomentMatrix", "drift");
1533 
1534  Machine::registerElement<ElementOrbTrim >("MomentMatrix", "orbtrim");
1535 
1536  Machine::registerElement<ElementSBend >("MomentMatrix", "sbend");
1537 
1538  Machine::registerElement<ElementQuad >("MomentMatrix", "quadrupole");
1539 
1540  Machine::registerElement<ElementSext >("MomentMatrix", "sextupole");
1541 
1542  Machine::registerElement<ElementSolenoid >("MomentMatrix", "solenoid");
1543 
1544  Machine::registerElement<ElementRFCavity >("MomentMatrix", "rfcavity");
1545 
1546  Machine::registerElement<ElementStripper >("MomentMatrix", "stripper");
1547 
1548  Machine::registerElement<ElementEDipole >("MomentMatrix", "edipole");
1549 
1550  Machine::registerElement<ElementEQuad >("MomentMatrix", "equad");
1551 
1552  Machine::registerElement<ElementTMatrix >("MomentMatrix", "tmatrix");
1553 }
Particle::gamma
double gamma
Gamma for ion. (dependent)
Definition: moment.h:22
Particle::IonQ
double IonQ
Ion Charge.
Definition: moment.h:22
MomentElementBase::check_cache
virtual bool check_cache(const state_t &S) const
Definition: moment.cpp:728
StateBase::retreat
bool retreat
retreat (backward) simulation flag
Definition: base.h:39
MomentState::getArray
virtual bool getArray(unsigned idx, ArrayInfo &Info)
Introspect named parameter of the derived class.
Definition: moment.cpp:306
ElementVoid
Base class for all simulated elements.
Definition: base.h:165
ElementVoid::length
double length
Longitudual length of this element (added to StateBase::pos)
Definition: base.h:195
Particle::IonW
double IonW
Total energy. (dependent)
Definition: moment.h:22
Particle::IonEk
double IonEk
Kinetic energy.
Definition: moment.h:22
MomentElementBase::show
virtual void show(std::ostream &strm, int level) const
Definition: moment.cpp:582
ElementVoid::type_name
virtual const char * type_name() const =0
Particle::IonZ
double IonZ
Charge state.
Definition: moment.h:22
StateBase
The abstract base class for all simulation state objects.
Definition: base.h:29
MomentElementBase::advance
virtual void advance(StateBase &s)
Propogate the given State through this Element.
Definition: moment.cpp:658
Particle::phis
double phis
Absolute synchrotron phase [rad].
Definition: moment.h:22
ElementVoid::assign
virtual void assign(const ElementVoid *other)=0
Definition: base.cpp:74
MomentElementBase::check_backward
virtual bool check_backward(const state_t &S) const
Check input state for backward propagation.
Definition: moment.cpp:738
Particle::IonEs
double IonEs
Rest energy.
Definition: moment.h:22
Particle::SampleFreq
double SampleFreq
Sampling frequency [Hz].
Definition: moment.h:22
Config
Associative configuration container.
Definition: config.h:66
MomentElementBase::recompute_matrix
virtual void recompute_matrix(state_t &ST)
recalculate &#39;transfer&#39; taking into consideration the provided input state
Definition: moment.cpp:760
StateBase::ArrayInfo
Used with StateBase::getArray() to describe a single parameter.
Definition: base.h:51
StateBase::assign
virtual void assign(const StateBase &other)=0
Definition: base.cpp:32
StateBase::pos
double pos
absolute longitudinal position at end of Element
Definition: base.h:37
MomentElementBase::assign
virtual void assign(const ElementVoid *other)=0
Definition: moment.cpp:563
MomentState::assign
virtual void assign(const StateBase &other)
Definition: moment.cpp:248
ElementVoid::name
const std::string name
Name of this element (unique in its Machine)
Definition: base.h:192
StateBase::ArrayInfo::dim
size_t dim[maxdims]
Array dimensions in elements.
Definition: base.h:69
ElementVoid::conf
const Config & conf() const
The Config used to construct this element.
Definition: base.h:190
MomentElementBase::dx
double dx
constituents of misalign
Definition: moment.h:185
MomentElementBase::skipcache
bool skipcache
If set, check_cache() will always return false.
Definition: moment.h:188
MomentState::size
size_t size() const
of charge states
Definition: moment.h:138
StateBase::ArrayInfo::stride
size_t stride[maxdims]
Array strides in bytes.
Definition: base.h:71
MomentElementBase::transfer
std::vector< value_t > transfer
final transfer matricies
Definition: moment.h:181
Particle::beta
double beta
Beta for ion. (dependent)
Definition: moment.h:22
Particle::bg
double bg
Beta*gamma. (dependent)
Definition: moment.h:22
MomentState::show
virtual void show(std::ostream &strm, int level) const
Definition: moment.cpp:265
Config::tryGet
bool tryGet(const std::string &name, T &val) const
Definition: config.h:156
StateBase::getArray
virtual bool getArray(unsigned index, ArrayInfo &Info)
Introspect named parameter of the derived class.
Definition: base.cpp:37
Config::get
detail::RT< T >::type get(const std::string &name) const
Definition: config.h:123
Particle::SampleIonK
double SampleIonK
Sample rate; different RF Cavity due to RF frequenies. (dependent)
Definition: moment.h:22
ElementVoid::show
virtual void show(std::ostream &, int level) const
Definition: base.cpp:69
MomentElementBase::resize_cache
void resize_cache(const state_t &ST)
Definition: moment.cpp:753
StateBase::ArrayInfo::name
const char * name
The parameter name.
Definition: base.h:55
MomentElementBase
An Element which propagates the statistical moments of a bunch.
Definition: moment.h:146
StateBase::ArrayInfo::ndim
unsigned ndim
Definition: base.h:67
Particle::SampleLambda
double SampleLambda
Sampling distance [m].
Definition: moment.h:22
Generated on Thu Apr 17 2025 01:29:13 for FLAME by   doxygen 1.8.5