OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
Twiss.cpp
Go to the documentation of this file.
1 // ------------------------------------------------------------------------
2 // $RCSfile: Twiss.cpp,v $
3 // ------------------------------------------------------------------------
4 // $Revision: 1.4.2.3 $
5 // ------------------------------------------------------------------------
6 // Copyright: see Copyright.readme
7 // ------------------------------------------------------------------------
8 //
9 // Class: Twiss
10 // The abstract class Twiss implements the interface for a table buffer
11 // holding lattice function.
12 //
13 // ------------------------------------------------------------------------
14 //
15 // $Date: 2004/11/12 20:10:11 $
16 // $Author: adelmann $
17 //
18 // ------------------------------------------------------------------------
19 
20 #include "Tables/Twiss.h"
21 #include "AbstractObjects/BeamSequence.h"
22 #include "AbstractObjects/Element.h"
23 #include "AbstractObjects/Expressions.h"
24 #include "AbstractObjects/OpalData.h"
25 #include "AbstractObjects/PlaceRep.h"
26 #include "AbstractObjects/RangeRep.h"
27 #include "Algorithms/ThickMapper.h"
28 #include "Algorithms/LinearMapper.h"
29 #include "Algorithms/ThinMapper.h"
30 #include "Attributes/Attributes.h"
31 #include "FixedAlgebra/FMatrix.h"
32 #include "FixedAlgebra/FVector.h"
33 #include "Physics/Physics.h"
34 #include "Structure/Beam.h"
35 #include "Tables/Flatten.h"
36 #include "Utilities/OpalException.h"
37 #include "Utilities/Options.h"
38 #include "Utilities/Round.h"
39 
40 #include <cmath>
41 #include <iomanip>
42 #include <iostream>
43 #include <sstream>
44 
45 // Local structures.
46 // ------------------------------------------------------------------------
47 
48 namespace {
49 
50  struct ColDesc {
51 
52  // Column name.
53  const char *colName;
54 
55  // Pointer to the row method which returns the column value.
56  double(Twiss::*get)(const Twiss::Row &, int, int) const;
57 
58  // Default field width and precision.
59  int printWidth, printPrecision;
60 
61  // Indices to be given to get().
62  int ind_1, ind_2;
63  };
64 
65 
66  // The complete column entries table.
67  const ColDesc allColumns[] = {
68 
69  // Arc length.
70  { "S", &Twiss::getS, 14, 6, 0, 0 },
71 
72  // "Naive" Twiss.
73  { "ALFX", &Twiss::getALFi, 10, 6, 0, 0 },
74  { "ALFY", &Twiss::getALFi, 10, 6, 1, 0 },
75  { "ALFT", &Twiss::getALFi, 10, 6, 2, 0 },
76 
77  { "BETX", &Twiss::getBETi, 10, 3, 0, 0 },
78  { "BETY", &Twiss::getBETi, 10, 3, 1, 0 },
79  { "BETT", &Twiss::getBETi, 10, 3, 2, 0 },
80 
81  { "MUX", &Twiss::getMUi, 10, 6, 0, 0 },
82  { "MUY", &Twiss::getMUi, 10, 6, 1, 0 },
83  { "MUT", &Twiss::getMUi, 10, 6, 2, 0 },
84 
85  // Mais-Ripken functions.
86  { "ALFX1", &Twiss::getALFik, 8, 3, 0, 0 },
87  { "ALFX2", &Twiss::getALFik, 8, 3, 0, 1 },
88  { "ALFX3", &Twiss::getALFik, 8, 3, 0, 2 },
89  { "ALFY1", &Twiss::getALFik, 8, 3, 1, 0 },
90  { "ALFY2", &Twiss::getALFik, 8, 3, 1, 1 },
91  { "ALFY3", &Twiss::getALFik, 8, 3, 1, 2 },
92  { "ALFT1", &Twiss::getALFik, 8, 3, 2, 0 },
93  { "ALFT2", &Twiss::getALFik, 8, 3, 2, 1 },
94  { "ALFT3", &Twiss::getALFik, 8, 3, 2, 2 },
95 
96  { "BETX1", &Twiss::getBETik, 10, 3, 0, 0 },
97  { "BETX2", &Twiss::getBETik, 10, 3, 0, 1 },
98  { "BETX3", &Twiss::getBETik, 10, 3, 0, 2 },
99  { "BETY1", &Twiss::getBETik, 10, 3, 1, 0 },
100  { "BETY2", &Twiss::getBETik, 10, 3, 1, 1 },
101  { "BETY3", &Twiss::getBETik, 10, 3, 1, 2 },
102  { "BETT1", &Twiss::getBETik, 10, 3, 2, 0 },
103  { "BETT2", &Twiss::getBETik, 10, 3, 2, 1 },
104  { "BETT3", &Twiss::getBETik, 10, 3, 2, 2 },
105 
106  { "GAMX1", &Twiss::getGAMik, 10, 3, 0, 0 },
107  { "GAMX2", &Twiss::getGAMik, 10, 3, 0, 1 },
108  { "GAMX3", &Twiss::getGAMik, 10, 3, 0, 2 },
109  { "GAMY1", &Twiss::getGAMik, 10, 3, 1, 0 },
110  { "GAMY2", &Twiss::getGAMik, 10, 3, 1, 1 },
111  { "GAMY3", &Twiss::getGAMik, 10, 3, 1, 2 },
112  { "GAMT1", &Twiss::getGAMik, 10, 3, 2, 0 },
113  { "GAMT2", &Twiss::getGAMik, 10, 3, 2, 1 },
114  { "GAMT3", &Twiss::getGAMik, 10, 3, 2, 2 },
115 
116  // Closed orbit.
117  { "XC", &Twiss::getCO, 10, 6, 0, 0 },
118  { "PXC", &Twiss::getCO, 10, 6, 1, 0 },
119  { "YC", &Twiss::getCO, 10, 6, 2, 0 },
120  { "PYC", &Twiss::getCO, 10, 6, 3, 0 },
121  { "TC", &Twiss::getCO, 10, 6, 4, 0 },
122  { "PTC", &Twiss::getCO, 10, 6, 5, 0 },
123 
124  // Dispersion.
125  { "DX", &Twiss::getDisp, 10, 6, 0, 0 },
126  { "DPX", &Twiss::getDisp, 10, 6, 1, 0 },
127  { "DY", &Twiss::getDisp, 10, 6, 2, 0 },
128  { "DPY", &Twiss::getDisp, 10, 6, 3, 0 },
129  { "DT", &Twiss::getDisp, 10, 6, 4, 0 },
130  { "DPT", &Twiss::getDisp, 10, 6, 5, 0 },
131 
132  // Eigenvectors.
133  { "E11", &Twiss::getEigen, 10, 6, 0, 0 },
134  { "E21", &Twiss::getEigen, 10, 6, 1, 0 },
135  { "E31", &Twiss::getEigen, 10, 6, 2, 0 },
136  { "E41", &Twiss::getEigen, 10, 6, 3, 0 },
137  { "E51", &Twiss::getEigen, 10, 6, 4, 0 },
138  { "E61", &Twiss::getEigen, 10, 6, 5, 0 },
139  { "E12", &Twiss::getEigen, 10, 6, 0, 1 },
140  { "E22", &Twiss::getEigen, 10, 6, 1, 1 },
141  { "E32", &Twiss::getEigen, 10, 6, 2, 1 },
142  { "E42", &Twiss::getEigen, 10, 6, 3, 1 },
143  { "E52", &Twiss::getEigen, 10, 6, 4, 1 },
144  { "E62", &Twiss::getEigen, 10, 6, 5, 1 },
145  { "E13", &Twiss::getEigen, 10, 6, 0, 2 },
146  { "E23", &Twiss::getEigen, 10, 6, 1, 2 },
147  { "E33", &Twiss::getEigen, 10, 6, 2, 2 },
148  { "E43", &Twiss::getEigen, 10, 6, 3, 2 },
149  { "E53", &Twiss::getEigen, 10, 6, 4, 2 },
150  { "E63", &Twiss::getEigen, 10, 6, 5, 2 },
151  { "E14", &Twiss::getEigen, 10, 6, 0, 3 },
152  { "E24", &Twiss::getEigen, 10, 6, 1, 3 },
153  { "E34", &Twiss::getEigen, 10, 6, 2, 3 },
154  { "E44", &Twiss::getEigen, 10, 6, 3, 3 },
155  { "E54", &Twiss::getEigen, 10, 6, 4, 3 },
156  { "E64", &Twiss::getEigen, 10, 6, 5, 3 },
157  { "E15", &Twiss::getEigen, 10, 6, 0, 4 },
158  { "E25", &Twiss::getEigen, 10, 6, 1, 4 },
159  { "E35", &Twiss::getEigen, 10, 6, 2, 4 },
160  { "E45", &Twiss::getEigen, 10, 6, 3, 4 },
161  { "E55", &Twiss::getEigen, 10, 6, 4, 4 },
162  { "E65", &Twiss::getEigen, 10, 6, 5, 4 },
163  { "E16", &Twiss::getEigen, 10, 6, 0, 5 },
164  { "E26", &Twiss::getEigen, 10, 6, 1, 5 },
165  { "E36", &Twiss::getEigen, 10, 6, 2, 5 },
166  { "E46", &Twiss::getEigen, 10, 6, 3, 5 },
167  { "E56", &Twiss::getEigen, 10, 6, 4, 5 },
168  { "E66", &Twiss::getEigen, 10, 6, 5, 5 },
169 
170  // Sigma matrix.
171  { "S11", &Twiss::getSigma, 10, 6, 0, 0 },
172  { "S21", &Twiss::getSigma, 10, 6, 1, 0 },
173  { "S31", &Twiss::getSigma, 10, 6, 2, 0 },
174  { "S41", &Twiss::getSigma, 10, 6, 3, 0 },
175  { "S51", &Twiss::getSigma, 10, 6, 4, 0 },
176  { "S61", &Twiss::getSigma, 10, 6, 5, 0 },
177  { "S12", &Twiss::getSigma, 10, 6, 0, 1 },
178  { "S22", &Twiss::getSigma, 10, 6, 1, 1 },
179  { "S32", &Twiss::getSigma, 10, 6, 2, 1 },
180  { "S42", &Twiss::getSigma, 10, 6, 3, 1 },
181  { "S52", &Twiss::getSigma, 10, 6, 4, 1 },
182  { "S62", &Twiss::getSigma, 10, 6, 5, 1 },
183  { "S13", &Twiss::getSigma, 10, 6, 0, 2 },
184  { "S23", &Twiss::getSigma, 10, 6, 1, 2 },
185  { "S33", &Twiss::getSigma, 10, 6, 2, 2 },
186  { "S43", &Twiss::getSigma, 10, 6, 3, 2 },
187  { "S53", &Twiss::getSigma, 10, 6, 4, 2 },
188  { "S63", &Twiss::getSigma, 10, 6, 5, 2 },
189  { "S14", &Twiss::getSigma, 10, 6, 0, 3 },
190  { "S24", &Twiss::getSigma, 10, 6, 1, 3 },
191  { "S34", &Twiss::getSigma, 10, 6, 2, 3 },
192  { "S44", &Twiss::getSigma, 10, 6, 3, 3 },
193  { "S54", &Twiss::getSigma, 10, 6, 4, 3 },
194  { "S64", &Twiss::getSigma, 10, 6, 5, 3 },
195  { "S15", &Twiss::getSigma, 10, 6, 0, 4 },
196  { "S25", &Twiss::getSigma, 10, 6, 1, 4 },
197  { "S35", &Twiss::getSigma, 10, 6, 2, 4 },
198  { "S45", &Twiss::getSigma, 10, 6, 3, 4 },
199  { "S55", &Twiss::getSigma, 10, 6, 4, 4 },
200  { "S65", &Twiss::getSigma, 10, 6, 5, 4 },
201  { "S16", &Twiss::getSigma, 10, 6, 0, 5 },
202  { "S26", &Twiss::getSigma, 10, 6, 1, 5 },
203  { "S36", &Twiss::getSigma, 10, 6, 2, 5 },
204  { "S46", &Twiss::getSigma, 10, 6, 3, 5 },
205  { "S56", &Twiss::getSigma, 10, 6, 4, 5 },
206  { "S66", &Twiss::getSigma, 10, 6, 5, 5 },
207 
208  // Transfer matrix.
209  { "R11", &Twiss::getMatrix, 10, 6, 0, 0 },
210  { "R21", &Twiss::getMatrix, 10, 6, 1, 0 },
211  { "R31", &Twiss::getMatrix, 10, 6, 2, 0 },
212  { "R41", &Twiss::getMatrix, 10, 6, 3, 0 },
213  { "R51", &Twiss::getMatrix, 10, 6, 4, 0 },
214  { "R61", &Twiss::getMatrix, 10, 6, 5, 0 },
215  { "R12", &Twiss::getMatrix, 10, 6, 0, 1 },
216  { "R22", &Twiss::getMatrix, 10, 6, 1, 1 },
217  { "R32", &Twiss::getMatrix, 10, 6, 2, 1 },
218  { "R42", &Twiss::getMatrix, 10, 6, 3, 1 },
219  { "R52", &Twiss::getMatrix, 10, 6, 4, 1 },
220  { "R62", &Twiss::getMatrix, 10, 6, 5, 1 },
221  { "R13", &Twiss::getMatrix, 10, 6, 0, 2 },
222  { "R23", &Twiss::getMatrix, 10, 6, 1, 2 },
223  { "R33", &Twiss::getMatrix, 10, 6, 2, 2 },
224  { "R43", &Twiss::getMatrix, 10, 6, 3, 2 },
225  { "R53", &Twiss::getMatrix, 10, 6, 4, 2 },
226  { "R63", &Twiss::getMatrix, 10, 6, 5, 2 },
227  { "R14", &Twiss::getMatrix, 10, 6, 0, 3 },
228  { "R24", &Twiss::getMatrix, 10, 6, 1, 3 },
229  { "R34", &Twiss::getMatrix, 10, 6, 2, 3 },
230  { "R44", &Twiss::getMatrix, 10, 6, 3, 3 },
231  { "R54", &Twiss::getMatrix, 10, 6, 4, 3 },
232  { "R64", &Twiss::getMatrix, 10, 6, 5, 3 },
233  { "R15", &Twiss::getMatrix, 10, 6, 0, 4 },
234  { "R25", &Twiss::getMatrix, 10, 6, 1, 4 },
235  { "R35", &Twiss::getMatrix, 10, 6, 2, 4 },
236  { "R45", &Twiss::getMatrix, 10, 6, 3, 4 },
237  { "R55", &Twiss::getMatrix, 10, 6, 4, 4 },
238  { "R65", &Twiss::getMatrix, 10, 6, 5, 4 },
239  { "R16", &Twiss::getMatrix, 10, 6, 0, 5 },
240  { "R26", &Twiss::getMatrix, 10, 6, 1, 5 },
241  { "R36", &Twiss::getMatrix, 10, 6, 2, 5 },
242  { "R46", &Twiss::getMatrix, 10, 6, 3, 5 },
243  { "R56", &Twiss::getMatrix, 10, 6, 4, 5 },
244  { "R66", &Twiss::getMatrix, 10, 6, 5, 5 },
245 
246  { 0, 0, 0, 0, 0, 0 }
247  };
248 
249 
250  // The default column entries table
251  const ColDesc defaultColumns[] = {
252  { "S", &Twiss::getS, 14, 6, 0, 0 },
253 
254  { "MUX", &Twiss::getMUi, 10, 6, 0, 0 },
255  { "BETX", &Twiss::getBETi, 10, 3, 0, 0 },
256  { "ALFX", &Twiss::getALFi, 10, 6, 0, 0 },
257 
258  { "MUY", &Twiss::getMUi, 10, 6, 1, 0 },
259  { "BETY", &Twiss::getBETi, 10, 3, 1, 0 },
260  { "ALFY", &Twiss::getALFi, 10, 6, 1, 0 },
261 
262  { "XC", &Twiss::getCO, 10, 6, 0, 0 },
263  { "YC", &Twiss::getCO, 10, 6, 2, 0 },
264  { "DX", &Twiss::getDisp, 10, 6, 0, 0 },
265  { "DPX", &Twiss::getDisp, 10, 6, 1, 0 },
266 
267  { 0, 0, 0, 0, 0, 0 }
268  };
269 
270 
271  const ColDesc *findCol(const Twiss &table, const std::string &colName) {
272  for(const ColDesc *col = allColumns; col->colName; ++col) {
273  if(colName == col->colName) {
274  return col;
275  }
276  }
277 
278  throw OpalException("Twiss::findCol()",
279  "Twiss table \"" + table.getOpalName() +
280  "\" has no column named \"" + colName + "\".");
281  }
282 
283 
284  // Local class Column.
285  // Returns the value for a given column in the current row of a given
286  // table.
287  class Column: public Expressions::Scalar<double> {
288 
289  public:
290 
291  //: Constructor.
292  // Identify the table by its name [b]tab[/b], and the column by its
293  // name [b]col[/b] and the function [b]col[/b].
294  // The row is specified as the ``current'' row of the table.
295  Column(const Twiss &tab, const std::string &colName, const ColDesc &desc);
296 
297  Column(const Column &);
298  virtual ~Column();
299 
300  //: Make clone.
301  virtual Expressions::Scalar<double> *clone() const;
302 
303  //: Evaluate.
304  virtual double evaluate() const;
305 
306  //: Print expression.
307  virtual void print(std::ostream &os, int precedence = 99) const;
308 
309  private:
310 
311  // Not implemented.
312  Column();
313  const Column &operator=(const Column &);
314 
315  // The Table referred.
316  const Twiss &itsTable;
317 
318  // Column name.
319  std::string colName;
320 
321  // The function returning the column value.
322  double(Twiss::*get)(const Twiss::Row &, int, int) const;
323 
324  // The indices to be transmitted to get().
325  int ind_1, ind_2;
326  };
327 
328 
329  // Implementation.
330  // ------------------------------------------------------------------------
331 
332  Column::Column(const Twiss &tab, const std::string &colName, const ColDesc &desc):
333  itsTable(tab), colName(colName),
334  get(desc.get), ind_1(desc.ind_1), ind_2(desc.ind_2)
335  {}
336 
337 
338  Column::Column(const Column &rhs):
339  Scalar<double>(rhs),
340  itsTable(rhs.itsTable), colName(rhs.colName),
341  get(rhs.get), ind_1(rhs.ind_1), ind_2(rhs.ind_2)
342  {}
343 
344 
345  Column::~Column()
346  {}
347 
348 
349  Expressions::Scalar<double> *Column::clone() const {
350  return new Column(*this);
351  }
352 
353 
354  double Column::evaluate() const {
355  return (itsTable.*get)(itsTable.getCurrent(), ind_1, ind_2);
356  }
357 
358 
359  void Column::print(std::ostream &os, int) const {
360  os << colName;
361  }
362 
363 };
364 
365 
366 // Class Twiss::Row
367 // ------------------------------------------------------------------------
368 
369 Twiss::Row::Row(ElementBase *elem, int occur):
370  FlaggedElmPtr(elem) {
371  setCounter(occur);
372 }
373 
374 
375 Twiss::Row::Row(const FlaggedElmPtr &rhs):
376  FlaggedElmPtr(rhs)
377 {}
378 
379 
380 Twiss::Row::~Row()
381 {}
382 
383 
384 const FVector<double, 6> &Twiss::Row::getCO() const {
385  return orbit;
386 }
387 
388 
389 const FMatrix<double, 6, 6> &Twiss::Row::getMatrix() const {
390  return matrix;
391 }
392 
393 
394 double Twiss::Row::getS() const {
395  return arc;
396 }
397 
398 
399 double Twiss::Row::getMUi(int i) const {
400  return mu[i];
401 }
402 
403 
404 // Class Twiss
405 // ------------------------------------------------------------------------
406 
407 Twiss::Twiss(int size, const char *name, const char *help):
408  Table(size, name, help), itsTable(0), itsMapper(0) {
409  itsAttr[LINE] = Attributes::makeString
410  ("LINE", "The beam line use for filling");
411  itsAttr[BEAM] = Attributes::makeString
412  ("BEAM", "The beam to be used", "UNNAMED_BEAM");
413  itsAttr[RANGE] = Attributes::makeRange
414  ("RANGE", "The range in the lattice");
415  itsAttr[STATIC] = Attributes::makeBool
416  ("STATIC", "If true, the table is not recalculated at each iteration");
417  itsAttr[ORDER] = Attributes::makeReal
418  ("ORDER", "The order of the calculation", 2);
419  itsAttr[METHOD] = Attributes::makeString
420  ("METHOD", "the algorithm used for filling:\n"
421  "\t\t\t\"LINEAR\", \"THIN\", \"THICK\", or \"TRANSPORT\"\n"
422  "\t\t\tDefault value is \"LINEAR\"", "LINEAR");
423  itsAttr[REVBEAM] = Attributes::makeBool
424  ("REVBEAM", "Set true to run beam backwards through lattice");
425  itsAttr[REVTRACK] = Attributes::makeBool
426  ("REVTRACK", "Set true to track against the beam");
427 
428  // READ ONLY.
429  itsAttr[BETXMAX] = Attributes::makeReal
430  ("BETXMAX", "Maximum horizontal beta in m");
431  itsAttr[BETXMAX].setReadOnly(true);
432 
433  itsAttr[BETYMAX] = Attributes::makeReal
434  ("BETYMAX", "Maximum vertical beta in m");
435  itsAttr[BETYMAX].setReadOnly(true);
436 
437  itsAttr[XCMAX] = Attributes::makeReal
438  ("XCMAX", "Maximum horizontal closed orbit in m");
439  itsAttr[XCMAX].setReadOnly(true);
440 
441  itsAttr[YCMAX] = Attributes::makeReal
442  ("YCMAX", "Maximum vertical closed orbit in m");
443  itsAttr[YCMAX].setReadOnly(true);
444 
445  itsAttr[XCRMS] = Attributes::makeReal
446  ("XCRMS", "R.m.s. horizontal closed orbit in m");
447  itsAttr[XCRMS].setReadOnly(true);
448 
449  itsAttr[YCRMS] = Attributes::makeReal
450  ("YCRMS", "R.m.s. vertical closed orbit in m");
451  itsAttr[YCRMS].setReadOnly(true);
452 
453  itsAttr[DXMAX] = Attributes::makeReal
454  ("DXMAX", "Maximum horizontal dispersion in m");
455  itsAttr[DXMAX].setReadOnly(true);
456 
457  itsAttr[DYMAX] = Attributes::makeReal
458  ("DYMAX", "Maximum vertical dispersion in m");
459  itsAttr[DYMAX].setReadOnly(true);
460 
461  itsAttr[DXRMS] = Attributes::makeReal
462  ("DXRMS", "R.m.s. horizontal dispersion in m");
463  itsAttr[DXRMS].setReadOnly(true);
464 
465  itsAttr[DYRMS] = Attributes::makeReal
466  ("DYRMS", "R.m.s. vertical dispersion in m");
467  itsAttr[DYRMS].setReadOnly(true);
468 }
469 
470 
471 Twiss::Twiss(const std::string &name, Twiss *parent):
472  Table(name, parent), itsTable(new TLine(name)), itsMapper(0)
473 {}
474 
475 
476 Twiss::~Twiss() {
477  delete itsTable;
478  delete itsMapper;
479 }
480 
481 /*
482 void Twiss::doomPut(DoomWriter &writer) const {
483  // Write the definition for this table.
484  Object::doomPut(writer);
485 
486  // Make sure the table is up to date.
487  // Cast away const, to allow logically constant table to update.
488  const_cast<Twiss *>(this)->fill();
489 
490  // Make the table header.
491  const std::string &tableName = getOpalName();
492  {
493  static const int headKeyList[] = { 1, TABLE_HEAD };
494  DoomWriter headWriter(tableName, headKeyList);
495  headWriter.setParentName(itsLine);
496  headWriter.setTypeName("TABLE_HEADER");
497 
498  static const char *columnName[numColumns] = {
499  "ALFX", "ALFY", "BETX", "BETY", "DPX", "DPY", "DX", "DY", "MUX", "MUY", "S"
500  };
501  static const int columnNumber[numColumns] = {
502  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
503  };
504  for(int i = 0; i < numColumns; ++i) {
505  headWriter.putInt(i, columnNumber[i]);
506  headWriter.putString(i, columnName[i]);
507  }
508 
509  int numRows = itsTable->size();
510  headWriter.putInt(numColumns + 0, 1); // 1 delta value.
511  headWriter.putInt(numColumns + 1, numRows); // number of rows.
512  headWriter.putInt(numColumns + 3, 3); // position at end.
513  // Header object will be written when headWriter goes out of scope.
514  }
515 
516  // Write the table row by row.
517  {
518  static const int rowKeyList[] = { 2, TABLE_BODY, 0 };
519  DoomWriter rowWriter(tableName, rowKeyList);
520  rowWriter.setParentName(itsLine);
521  rowWriter.setTypeName("SUB_TABLE");
522  int lineCount = 0;
523  int realCount = 0;
524  for(TLine::const_iterator row = begin(); row != end(); ++row) {
525  rowWriter.putInt(lineCount, row->getCounter());
526  rowWriter.putString(lineCount, row->getElement()->getName());
527  rowWriter.putReal(realCount + 0, getALFi(*row, 0, 0));
528  rowWriter.putReal(realCount + 1, getALFi(*row, 1, 0));
529  rowWriter.putReal(realCount + 2, getBETi(*row, 0, 0));
530  rowWriter.putReal(realCount + 3, getBETi(*row, 1, 0));
531  rowWriter.putReal(realCount + 4, getDisp(*row, 1, 0));
532  rowWriter.putReal(realCount + 5, getDisp(*row, 3, 0));
533  rowWriter.putReal(realCount + 6, getDisp(*row, 0, 0));
534  rowWriter.putReal(realCount + 7, getDisp(*row, 2, 0));
535  rowWriter.putReal(realCount + 8, getMUi(*row, 0, 0));
536  rowWriter.putReal(realCount + 9, getMUi(*row, 1, 0));
537  rowWriter.putReal(realCount + 10, getS(*row, 0, 0));
538  lineCount++;
539  realCount += numColumns;
540  }
541  // Body object will be written when rowWriter goes out of scope.
542  }
543 }
544 */
545 
546 Twiss::TLine::iterator Twiss::begin() {
547  return itsTable->begin();
548 }
549 
550 
551 Twiss::TLine::const_iterator Twiss::begin() const {
552  return itsTable->begin();
553 }
554 
555 
556 Twiss::TLine::iterator Twiss::end() {
557  return itsTable->end();
558 }
559 
560 
561 Twiss::TLine::const_iterator Twiss::end() const {
562  return itsTable->end();
563 }
564 
565 
566 void Twiss::execute() {
567  //std::cerr << "In Twiss::execute()" << std::endl;
568  // Find Table definition.
569  itsLine = Attributes::getString(itsAttr[LINE]);
570  BeamSequence *use = BeamSequence::find(itsLine);
571 
572  // Find Beam data.
573  const std::string &beamName = Attributes::getString(itsAttr[BEAM]);
574  beam = Beam::find(beamName);
575  reference = &beam->getReference();
576 
577  // Get the algorithm name.
578  std::string method = Attributes::getString(itsAttr[METHOD]);
579 
580  // Make sure all is up-to-date.
581  OpalData::getInstance()->update();
582 
583  // Create flat list with space for data storage.
584  RangeRep range = Attributes::getRange(itsAttr[RANGE]);
585  Flatten<Row> flattener(*use->fetchLine(), *itsTable, range);
586  flattener.execute();
587 
588  if(itsTable->empty() && Options::warn) {
589  std::cerr << "\n### Warning ### Lattice function table \""
590  << getOpalName() << "\" contains no elements.\n" << std::endl;
591  } else {
592  itsTable->front().setSelectionFlag(true);
593  itsTable->back().setSelectionFlag(true);
594  }
595 
596  // Assign the correct mapper.
597  revBeam = Attributes::getBool(itsAttr[REVBEAM]);
598  revTrack = Attributes::getBool(itsAttr[REVTRACK]);
599  revPath = ( revBeam && !revTrack ) || ( !revBeam && revTrack );
600  order = int(Round(Attributes::getReal(itsAttr[ORDER])));
601 
602  if(method == "THICK") {
603  //std::cerr << " method == \"THICK\"" << std::endl;
604  itsMapper = new ThickMapper(*itsTable, *reference, revBeam, revTrack);
605  } else if(method == "THIN") {
606  //std::cerr << " method == \"THIN\"" << std::endl;
607  itsMapper = new ThinMapper(*itsTable, *reference, revBeam, revTrack);
608  } else if(method == "LINEAR") {
609  //std::cerr << " method == \"LINEAR\"" << std::endl;
610  itsMapper = new LinearMapper(*itsTable, *reference, revBeam, revTrack);
611  } else {
612  throw OpalException("Twiss::execute()",
613  "Method name \"" + method + "\" is unknown.");
614  }
615 
616  // Fill the table.
617  fill();
618 
619  // Set static flag.
620  if(Attributes::getBool(itsAttr[STATIC])) dynamic = false;
621  printTable(std::cout, getDefault());
622  //std::cerr << "Leaving Twiss::execute()" << std::endl;
623 }
624 
625 
626 double Twiss::getCell(const PlaceRep &place, const std::string &colName) {
627  Row &row = findRow(place);
628  const ColDesc *col = findCol(*this, colName);
629  return (this->*(col->get))(row, col->ind_1, col->ind_2);
630 }
631 
632 
633 Table::CellArray Twiss::getDefault() const {
634  CellArray columns;
635  for(const ColDesc *col = defaultColumns; col->colName; ++col) {
636  Expressions::PtrToScalar<double> expr =
637  new Column(*this, col->colName, *col);
638  columns.push_back(Cell(expr, col->printWidth, col->printPrecision));
639  }
640  return columns;
641 }
642 
643 
644 std::vector<double>
645 Twiss::getColumn(const RangeRep &rng, const std::string &colName) {
646  // Find proper column function.
647  const ColDesc *col = findCol(*this, colName);
648  RangeRep range(rng);
649  range.initialize();
650  std::vector<double> column;
651 
652  for(TLine::const_iterator row = begin(); row != end(); ++row) {
653  range.enter(*row);
654  if(range.isActive()) {
655  column.push_back((this->*(col->get))(*row, col->ind_1, col->ind_2));
656  }
657  range.leave(*row);
658  }
659 
660  return column;
661 }
662 
663 
664 const Twiss::Row &Twiss::getCurrent() const {
665  return *current;
666 }
667 
668 
669 FMatrix<double, 6, 6> Twiss::getCurlyA() const {
670  return curly_A;
671 }
672 
673 
674 FMatrix<double, 6, 6> Twiss::getCurlyA(const Row &row) const {
675  return row.matrix * curly_A;
676 }
677 
678 
679 double Twiss::getEX() const {
680  return beam->getEX();
681 }
682 
683 
684 double Twiss::getEY() const {
685  return beam->getEY();
686 }
687 
688 
689 double Twiss::getET() const {
690  return beam->getET();
691 }
692 
693 
694 double Twiss::getLength() {
695  return itsTable->getElementLength();
696 }
697 
698 
699 const Beamline *Twiss::getLine() const {
700  return itsTable;
701 }
702 
703 
704 FMatrix<double, 6, 6> Twiss::getMatrix(const Row &row) const {
705  return row.matrix;
706 }
707 
708 
709 FVector<double, 6> Twiss::getOrbit(const Row &row) const {
710  return row.orbit;
711 }
712 
713 
714 FVector<double, 6> Twiss::getOrbit() const {
715  return orbit;
716 }
717 
718 
719 FMatrix<double, 6, 6> Twiss::getSigma(const Row &row) const {
720  double E1 = getEX();
721  double E2 = getEY();
722  double E3 = getET();
723  const FMatrix<double, 6, 6> &eigen = getCurlyA(row);
724  FMatrix<double, 6, 6> sigma;
725  for(int i = 0; i < 6; ++i) {
726  for(int j = 0; j <= i; ++j) {
727  sigma[i][j] = sigma[j][i] =
728  E1 * (eigen[i][0] * eigen[j][0] + eigen[i][1] * eigen[j][1]) +
729  E2 * (eigen[i][2] * eigen[j][2] + eigen[i][3] * eigen[j][3]) +
730  E3 * (eigen[i][4] * eigen[j][4] + eigen[i][5] * eigen[j][5]);
731  }
732  }
733  return sigma;
734 }
735 
736 
737 std::vector<double>
738 Twiss::getRow(const PlaceRep &pos, const std::vector<std::string> &cols) {
739  Row &row = findRow(pos);
740  std::vector<double> result;
741 
742  if(cols.empty()) {
743  // Standard column selection.
744  for(const ColDesc *col = defaultColumns; col->colName != 0; ++col) {
745  result.push_back((this->*col->get)(row, col->ind_1, col->ind_2));
746  }
747  } else {
748  // User column selection.
749  for(std::vector<std::string>::const_iterator iter = cols.begin();
750  iter != cols.end(); ++iter) {
751  const ColDesc *col = findCol(*this, *iter);
752  result.push_back((this->*(col->get))(row, col->ind_1, col->ind_2));
753  }
754  }
755 
756  return result;
757 }
758 
759 
760 bool Twiss::isDependent(const std::string &name) const {
761  // Test if name refers to USE attribute.
762  if(itsLine == name) return true;
763 
764  // Test if name occurs in table.
765  for(TLine::const_iterator row = begin(); row != end(); ++row) {
766  if(row->getElement()->getName() == name) return true;
767  }
768 
769  // Otherwise replacement is not required.
770  return false;
771 }
772 
773 
774 Expressions::PtrToScalar<double>
775 Twiss::makeColumnExpression(const std::string &colName) const {
776  const ColDesc *col = findCol(*this, colName);
777  return new Column(*this, colName, *col);
778 }
779 
780 
781 bool Twiss::matches(Table *rhs) const {
782  return dynamic_cast<Twiss *>(rhs) != 0;
783 }
784 
785 
786 void Twiss::printTableBody(std::ostream &os, const CellArray &cells) const {
787  // Find line length.
788  int lineLength = 16;
789  for(CellArray::const_iterator cell = cells.begin();
790  cell < cells.end(); ++cell) {
791  lineLength += cell->printWidth;
792  }
793 
794  // Write table header.
795  os << std::string(lineLength, '-') << '\n';
796  os << "Element ";
797  for(CellArray::const_iterator cell = cells.begin();
798  cell < cells.end(); ++cell) {
799  std::ostringstream ss;
800  cell->itsExpr->print(ss, 0);
801  ss << std::ends;
802  std::string image = ss.str();
803 
804  if(int(image.length()) < cell->printWidth) {
805  // Right adjust the column header.
806  os << std::string(cell->printWidth - image.length(), ' ') << image;
807  } else {
808  // Truncate the column header.
809  os << ' ' << std::string(image, 0, cell->printWidth - 3) << "..";
810  }
811  }
812  os << '\n';
813  os << std::string(lineLength, '-') << '\n';
814 
815  // Write table body.
816  for(current = begin(); current != end(); ++current) {
817  if(current->getSelectionFlag()) {
818  std::string name = current->getElement()->getName();
819  if(int occur = current->getCounter()) {
820  std::ostringstream tos;
821  tos << name << '[' << occur << ']' << std::ends;
822  name = tos.str();
823  }
824 
825  if(name.length() > 16) {
826  // Truncate the element name.
827  os << std::string(name, 0, 13) << ".. ";
828  } else {
829  // Left adjust the element name.
830  os << name << std::string(16 - name.length(), ' ');
831  }
832 
833  for(CellArray::const_iterator cell = cells.begin();
834  cell != cells.end(); ++cell) {
835  os << std::setw(cell->printWidth)
836  << std::setprecision(cell->printPrecision)
837  << cell->itsExpr->evaluate();
838  }
839  os << '\n';
840  }
841  }
842 
843  os << std::string(lineLength, '-') << '\n';
844 }
845 
846 
847 void Twiss::printTableTitle(std::ostream &os, const char *title) const {
848  OpalData::getInstance()->printTitle(os);
849  os << '\n' << title
850  << ", LINE: " << itsAttr[LINE]
851  << ", BEAM: " << itsAttr[BEAM]
852  << ", RANGE: " << itsAttr[RANGE]
853  << ", METHOD: " << itsAttr[METHOD]
854  << ", ORDER: " << order << ".\n";
855 }
856 
857 
858 // Protected methods.
859 // ------------------------------------------------------------------------
860 
861 Twiss::Row &Twiss::findRow(const PlaceRep &place) {
862  PlaceRep row(place);
863  row.initialize();
864 
865  for(TLine::iterator i = begin(); i != end(); ++i) {
866  row.enter(*i);
867  if(row.isActive()) return *i;
868  row.leave(*i);
869  }
870 
871  std::ostringstream os;
872  os << row << std::ends;
873  throw OpalException("Twiss::findRow()", "Row \"" + os.str() +
874  "\" not found in twiss table \"" + getOpalName() + "\".");
875 }
876 
877 
878 void Twiss::put() {
879  //std::cerr << "==> In Twiss::put()..." << std::endl;
880  // Initial (final) phase angles.
881  static const double epsilon = 1.0e-8;
882  static const LinearFun<double, 6> nrgy = LinearFun<double, 6>::makeVariable(5);
883  double arc = 0.0;
884  LinearMap<double, 6> map;
885  itsMapper->getMap(map);
886  bool staticQ = (map[5] == map[5][0] + nrgy);
887 
888  //std::cerr << " [Twiss::put] map =\n" << map << std::endl;
889  //std::cerr << " [Twiss::put] curly_A =\n" << curly_A << std::endl;
890  int imux = 0;
891  int imuy = 0;
892  int imut = 0;
893  double mux = atan2(curly_A[0][1], curly_A[0][0]) / Physics::two_pi;
894  double muy = atan2(curly_A[2][3], curly_A[2][2]) / Physics::two_pi;
895  double mut = atan2(curly_A[4][5], curly_A[4][4]) / Physics::two_pi;
896  //std::cerr << " tunes = (" << mux << ", " << muy << ", " << mut << ")" << std::endl;
897 
898  if(revPath) {
899  //std::cerr << "Twiss::put(): doing reverse path ..." << std::endl;
900  for(TLine::reverse_iterator row = itsTable->rbegin();
901  row != itsTable->rend(); ++row) {
902  // Store values at end of element to the table.
903  row->orbit = map.constantTerm();
904  row->matrix = map.linearTerms();
905  row->arc = arc;
906  row->mu[0] = mux + double(imux);
907  row->mu[1] = muy + double(imuy);
908  row->mu[2] = mut + double(imut);
909 
910  // Traverse element.
911  row->accept(*itsMapper);
912 
913  // Update values for beginning of element.
914  ElementBase &elem = *row->getElement();
915  if(! dynamic_cast<Beamline *>(&elem)) {
916  arc -= elem.getElementLength();
917  }
918 
919  itsMapper->getMap(map);
920  double A11 = 0.0;
921  double A12 = 0.0;
922  double A33 = 0.0;
923  double A34 = 0.0;
924  double A55 = 0.0;
925  double A56 = 0.0;
926  FMatrix<double, 6, 6> matrix = map.linearTerms();
927  const double *row1 = matrix[0];
928  const double *row3 = matrix[2];
929  const double *row5 = matrix[4];
930 
931  for(int i = 0; i < 6; ++i) {
932  A11 += row1[i] * curly_A[i][0];
933  A12 += row1[i] * curly_A[i][1];
934  A33 += row3[i] * curly_A[i][2];
935  A34 += row3[i] * curly_A[i][3];
936  A55 += row5[i] * curly_A[i][4];
937  A56 += row5[i] * curly_A[i][5];
938  }
939 
940  double mu = atan2(A12, A11) / Physics::two_pi;
941  if(mu - epsilon > mux) --imux;
942  mux = mu;
943  mu = atan2(A34, A33) / Physics::two_pi;
944  if(mu - epsilon > muy) --imuy;
945  muy = mu;
946  mu = atan2(A56, A55) / Physics::two_pi;
947  if(mu - epsilon > mut) --imut;
948  mut = mu;
949 
950  // FMatrix<double,6,6> A_scr = map.linearTerms() * curly_A;
951  // double mu;
952  // mu = atan2(A_scr[0][1], A_scr[0][0]) / Physics::two_pi;
953  // if (mu - epsilon > mux) --imux;
954  // mux = mu;
955  // mu = atan2(A_scr[2][3], A_scr[2][2]) / Physics::two_pi;
956  // if (mu - epsilon > muy) --imuy;
957  // muy = mu;
958  // mu = atan2(A_scr[4][5], A_scr[4][4]) / Physics::two_pi;
959  // if (mu - epsilon > mut) --imut;
960  // mut = mu;
961  }
962 
963  double arc1 = arc;
964  arc = - arc;
965  double mux1 = mux + double(imux);
966  double muy1 = muy + double(imuy);
967  double mut1 = mut + double(imut);
968 
969  // Change arc length and phases to their values from the origin.
970  for(TLine::iterator row = begin(); row != end(); ++row) {
971  row->arc -= arc1;
972  row->mu[0] -= mux1;
973  row->mu[1] -= muy1;
974  row->mu[2] -= mut1;
975  }
976  } else {
977  //std::cerr << "Twiss::put(): doing forward path ..." << std::endl;
978  double arc1 = 0.0;
979  double mux1 = mux;
980  double muy1 = muy;
981  double mut1 = mut;
982 
983  if(staticQ) { // reset 'T' component of map
984  //std::cerr << "Twiss::put(): resetting 'T' component of map ..." << std::endl;
985  itsMapper->getMap(map);
986  map[4][0] = 0.0;
987  itsMapper->setMap(map);
988  }
989 
990  for(TLine::iterator row = begin(); row != end(); ++row) {
991  // Traverse element.
992  row->accept(*itsMapper);
993 
994  // Update values for end of element.
995  ElementBase &elem = *row->getElement();
996  if(! dynamic_cast<Beamline *>(&elem)) {
997  arc += elem.getElementLength();
998  }
999 
1000  // Store values at end of element to the table.
1001  itsMapper->getMap(map);
1002  row->orbit = map.constantTerm();
1003  //std::cerr << "Twiss::put(): row->orbit =\n" << row->orbit << std::endl;
1004  row->matrix = map.linearTerms();
1005  double A11 = 0.0;
1006  double A12 = 0.0;
1007  double A33 = 0.0;
1008  double A34 = 0.0;
1009  double A55 = 0.0;
1010  double A56 = 0.0;
1011  FMatrix<double, 6, 6> &matrix = row->matrix;
1012  const double *row1 = matrix[0];
1013  const double *row3 = matrix[2];
1014  const double *row5 = matrix[4];
1015 
1016  //std::cerr << " matrix =\n" << matrix << std::endl;
1017  //std::cerr << " curly_A =\n" << curly_A << std::endl;
1018 
1019  for(int i = 0; i < 6; ++i) {
1020  A11 += row1[i] * curly_A[i][0];
1021  A12 += row1[i] * curly_A[i][1];
1022  A33 += row3[i] * curly_A[i][2];
1023  A34 += row3[i] * curly_A[i][3];
1024  A55 += row5[i] * curly_A[i][4];
1025  A56 += row5[i] * curly_A[i][5];
1026  }
1027 
1028  double mu = atan2(A12, A11) / Physics::two_pi;
1029  if(mu + epsilon < mux) ++imux;
1030  mux = mu;
1031  mu = atan2(A34, A33) / Physics::two_pi;
1032  if(mu + epsilon < muy) ++imuy;
1033  muy = mu;
1034  mu = atan2(A56, A55) / Physics::two_pi;
1035  if(mu + epsilon < mut) ++imut;
1036  mut = mu;
1037 
1038  // FMatrix<double,6,6> A_scr = row->matrix * curly_A;
1039  // double mu;
1040  // mu = atan2(A_scr[0][1], A_scr[0][0]) / Physics::two_pi;
1041  // if (mu + epsilon < mux) ++imux;
1042  // mux = mu;
1043  // mu = atan2(A_scr[2][3], A_scr[2][2]) / Physics::two_pi;
1044  // if (mu + epsilon < muy) ++imuy;
1045  // muy = mu;
1046  // mu = atan2(A_scr[4][5], A_scr[4][4]) / Physics::two_pi;
1047  // if (mu + epsilon < mut) ++imut;
1048  // mut = mu;
1049 
1050  row->arc = arc - arc1;
1051  row->mu[0] = mux + double(imux) - mux1;
1052  row->mu[1] = muy + double(imuy) - muy1;
1053  row->mu[2] = mut + double(imut) - mut1;
1054  }
1055  }
1056 
1057  // Compute statistic quantities.
1058  double betxmax = 0.0;
1059  double betymax = 0.0;
1060  double xmax = 0.0;
1061  double ymax = 0.0;
1062  double xrms = 0.0;
1063  double yrms = 0.0;
1064  double dxmax = 0.0;
1065  double dymax = 0.0;
1066  double dxrms = 0.0;
1067  double dyrms = 0.0;
1068 
1069  for(TLine::iterator row = begin(); row != end(); ++row) {
1070  betxmax = std::max(std::abs(getBETi(*row, 0, 0)), betxmax);
1071  betymax = std::max(std::abs(getBETi(*row, 1, 0)), betymax);
1072  double x = getCO(*row, 0, 0);
1073  double y = getCO(*row, 2, 0);
1074  double dx = getDisp(*row, 0, 0);
1075  double dy = getDisp(*row, 2, 0);
1076  xmax = std::max(std::abs(x), xmax);
1077  ymax = std::max(std::abs(y), ymax);
1078  xrms += x * x;
1079  yrms += y * y;
1080  dxmax = std::max(std::abs(dx), dxmax);
1081  dymax = std::max(std::abs(dy), dymax);
1082  dxrms += dx * dx;
1083  dyrms += dy * dy;
1084  }
1085 
1086  double size = itsTable->size();
1087  Attributes::setReal(itsAttr[BETXMAX], betxmax);
1088  Attributes::setReal(itsAttr[BETYMAX], betymax);
1089  Attributes::setReal(itsAttr[XCMAX], xmax);
1090  Attributes::setReal(itsAttr[YCMAX], ymax);
1091  Attributes::setReal(itsAttr[DXMAX], dxmax);
1092  Attributes::setReal(itsAttr[DYMAX], dymax);
1093  Attributes::setReal(itsAttr[XCRMS], sqrt(xrms / size));
1094  Attributes::setReal(itsAttr[YCRMS], sqrt(yrms / size));
1095  Attributes::setReal(itsAttr[DXRMS], sqrt(dxrms / size));
1096  Attributes::setReal(itsAttr[DYRMS], sqrt(dyrms / size));
1097  // std::cerr << "==> Leaving Twiss::put()" << std::endl;
1098 }
1099 
1100 
1101 double Twiss::getS(const Twiss::Row &row, int, int) const {
1102  return row.getS();
1103 }
1104 
1105 
1106 double Twiss::getMUi(const Twiss::Row &row, int i1, int) const {
1107  return row.getMUi(i1);
1108 }
1109 
1110 
1111 double Twiss::getBETi(const Twiss::Row &row, int i1, int) const {
1112  const double *row1 = row.matrix[2*i1];
1113  const double *row2 = row.matrix[2*i1+1];
1114  double r11 = 0.0;
1115  double r12 = 0.0;
1116  double r21 = 0.0;
1117  double r22 = 0.0;
1118 
1119  for(int i = 0; i < 6; ++i) {
1120  double t1 = row1[i];
1121  r11 += t1 * curly_A[i][2*i1];
1122  r12 += t1 * curly_A[i][2*i1+1];
1123  double t2 = row2[i];
1124  r21 += t2 * curly_A[i][2*i1];
1125  r22 += t2 * curly_A[i][2*i1+1];
1126  }
1127 
1128  // const FMatrix<double,6,6> eigen = row.matrix * curly_A;
1129  // const double r11 = eigen[2*i1][2*i1];
1130  // const double r12 = eigen[2*i1][2*i1+1];
1131  // const double r21 = eigen[2*i1+1][2*i1];
1132  // const double r22 = eigen[2*i1+1][2*i1+1];
1133  return (r11 * r11 + r12 * r12) / (r11 * r22 - r12 * r21);
1134 }
1135 
1136 
1137 double Twiss::getALFi(const Twiss::Row &row, int i1, int) const {
1138  const double *row1 = row.matrix[2*i1];
1139  const double *row2 = row.matrix[2*i1+1];
1140  double r11 = 0.0;
1141  double r12 = 0.0;
1142  double r21 = 0.0;
1143  double r22 = 0.0;
1144 
1145  for(int i = 0; i < 6; ++i) {
1146  double t1 = row1[i];
1147  r11 += t1 * curly_A[i][2*i1];
1148  r12 += t1 * curly_A[i][2*i1+1];
1149  double t2 = row2[i];
1150  r21 += t2 * curly_A[i][2*i1];
1151  r22 += t2 * curly_A[i][2*i1+1];
1152  }
1153 
1154  // ada Mon Mar 27 23:45:26 CEST 2000
1155  // const FMatrix<double,6,6> eigen = row.matrix * curly_A;
1156  // const double r11 = eigen[2*i1][2*i1];
1157  // const double r12 = eigen[2*i1][2*i1+1];
1158  // const double r21 = eigen[2*i1+1][2*i1];
1159  // const double r22 = eigen[2*i1+1][2*i1+1];
1160  return - (r11 * r21 + r12 * r22) / (r11 * r22 - r12 * r21);
1161 }
1162 
1163 
1164 double Twiss::getBETik(const Twiss::Row &row, int i1, int i2) const {
1165  const double *row1 = row.matrix[2*i1];
1166  double r11 = 0.0;
1167  double r12 = 0.0;
1168 
1169  for(int i = 0; i < 6; ++i) {
1170  double t1 = row1[i];
1171  r11 += t1 * curly_A[i][2*i2];
1172  r12 += t1 * curly_A[i][2*i2+1];
1173  }
1174 
1175  // ada Mon Mar 27 23:45:26 CEST 2000
1176  // const FMatrix<double,6,6> eigen = row.matrix * curly_A;
1177  // const double r11 = eigen[2*i1][2*i2];
1178  // const double r12 = eigen[2*i1][2*i2+1];
1179  return (r11 * r11 + r12 * r12);
1180 }
1181 
1182 
1183 double Twiss::getALFik(const Twiss::Row &row, int i1, int i2) const {
1184  const double *row1 = row.matrix[2*i1];
1185  const double *row2 = row.matrix[2*i1+1];
1186  double r11 = 0.0;
1187  double r12 = 0.0;
1188  double r21 = 0.0;
1189  double r22 = 0.0;
1190 
1191  for(int i = 0; i < 6; ++i) {
1192  double t1 = row1[i];
1193  r11 += t1 * curly_A[i][2*i2];
1194  r12 += t1 * curly_A[i][2*i2+1];
1195  double t2 = row2[i];
1196  r21 += t2 * curly_A[i][2*i2];
1197  r22 += t2 * curly_A[i][2*i2+1];
1198  }
1199 
1200  // ada Mon Mar 27 23:45:26 CEST 2000
1201  // const FMatrix<double,6,6> eigen = row.matrix * curly_A;
1202  // const double r11 = eigen[2*i1][2*i2];
1203  // const double r12 = eigen[2*i1][2*i2+1];
1204  // const double r21 = eigen[2*i1+1][2*i2];
1205  // const double r22 = eigen[2*i1+1][2*i2+1];
1206  return (r11 * r21 + r12 * r22);
1207 }
1208 
1209 
1210 double Twiss::getGAMik(const Twiss::Row &row, int i1, int i2) const {
1211  const double *row1 = row.matrix[2*i1+1];
1212  double r21 = 0.0;
1213  double r22 = 0.0;
1214 
1215  for(int i = 0; i < 6; ++i) {
1216  double t1 = row1[i];
1217  r21 += t1 * curly_A[i][2*i2];
1218  r22 += t1 * curly_A[i][2*i2+1];
1219  }
1220 
1221  // ada Mon Mar 27 23:45:26 CEST 2000
1222  // const FMatrix<double,6,6> eigen = row.matrix * curly_A;
1223  // const double r21 = eigen[2*i1+1][2*i2];
1224  // const double r22 = eigen[2*i1+1][2*i2+1];
1225  return (r21 * r21 + r22 * r22);
1226 }
1227 
1228 
1229 double Twiss::getCO(const Twiss::Row &row, int i1, int) const {
1230  return row.getCO()[i1];
1231 }
1232 
1233 
1234 double Twiss::getDisp(const Twiss::Row &row, int i1, int) const {
1235  // FMatrix<double,6,6> matrix = row.matrix;
1236  const double *row1 = row.matrix[i1];
1237  double result = 0.0;
1238 
1239  for(int i = 0; i < 6; ++i) {
1240  // result += matrix[i1][i] * curly_A[i][5];
1241  result += row1[i] * curly_A[i][5];
1242  }
1243 
1244  return result;
1245 }
1246 
1247 
1248 double Twiss::getEigen(const Twiss::Row &row, int i1, int i2) const {
1249  // FMatrix<double,6,6> matrix = row.matrix;
1250  const double *row1 = row.matrix[i1];
1251  double result = 0.0;
1252 
1253  for(int i = 0; i < 6; ++i) {
1254  // result += matrix[i1][i] * curly_A[i][i2];
1255  result += row1[i] * curly_A[i][i2];
1256  }
1257 
1258  return result;
1259 }
1260 
1261 /*
1262 1334 double Twiss::getEigen(const Twiss::Row &row, int i1, int i2) const
1263 1335 {
1264 1336 FMatrix<double,6,6> matrix = row.matrix;
1265 1337 const double *row1 = row.matrix[i1];
1266 1338 double result = 0.0;
1267 1339
1268 1340 for (int i = 0; i < 6; ++i) {
1269 1341 result += matrix[i1][i] * curly_A[i][i2];
1270 1342 //result += row1[i] * curly_A[i][i2];
1271 1343 }
1272 1344
1273 1345 return result;
1274 1346 }
1275 */
1276 
1277 
1278 
1279 
1280 double Twiss::getSigma(const Twiss::Row &row, int i1, int i2) const {
1281  const FMatrix<double, 6, 6> eigen = row.matrix * curly_A;
1282  const double E1 = getEX();
1283  const double E2 = getEY();
1284  const double E3 = getET();
1285  return
1286  (E1 * (eigen[i1][0] * eigen[i2][0] + eigen[i1][1] * eigen[i2][1]) +
1287  E2 * (eigen[i1][2] * eigen[i2][2] + eigen[i1][3] * eigen[i2][3]) +
1288  E3 * (eigen[i1][4] * eigen[i2][4] + eigen[i1][5] * eigen[i2][5]));
1289 }
1290 
1291 
1292 double Twiss::getMatrix(const Twiss::Row &row, int i1, int i2) const {
1293  return row.matrix[i1][i2];
1294 }
Twiss::~Twiss
virtual ~Twiss()
Definition: Twiss.cpp:476
Attributes::setReal
void setReal(Attribute &attr, double val)
Set real value.
Definition: Attributes.cpp:236
Round
double Round(double value)
Round the double argument.
Definition: Round.cpp:23
Twiss::revBeam
bool revBeam
Definition: Twiss.h:305
Table::CellArray
std::vector< Cell > CellArray
An array of cell descriptors.
Definition: Table.h:63
Twiss::matches
virtual bool matches(Table *rhs) const
Check compatibility.
Definition: Twiss.cpp:781
Twiss::isDependent
virtual bool isDependent(const std::string &name) const
Check dependency.
Definition: Twiss.cpp:760
Expressions::Scalar
A scalar expression.
Definition: Expressions.h:79
Table::fill
virtual void fill()=0
Refill the buffer.
abs
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
Definition: IpplExpressions.h:62
Twiss::makeColumnExpression
virtual Expressions::PtrToScalar< double > makeColumnExpression(const std::string &colName) const
Return column expression.
Definition: Twiss.cpp:775
Expressions::Scalar::print
virtual void print(std::ostream &, int precedence=99) const =0
Print expression.
Twiss::current
TLine::const_iterator current
Definition: Twiss.h:287
ElementBase
Interface for basic beam line object.
Definition: ElementBase.h:128
Expressions::Scalar::operator=
void operator=(const Scalar &)
Twiss::getET
double getET() const
Return emittance for mode 3.
Definition: Twiss.cpp:689
Twiss::end
TLine::const_iterator end() const
Access to last row.
Definition: Twiss.cpp:561
Attributes::getRange
RangeRep getRange(const Attribute &attr)
Get range value.
Definition: Attributes.cpp:177
Twiss::getEY
double getEY() const
Return emittance for mode 2.
Definition: Twiss.cpp:684
Twiss::getSigma
FMatrix< double, 6, 6 > getSigma() const
Initial envelope (Sigma) matrix.
Twiss::getGAMik
double getGAMik(const Row &, int i1, int i2) const
Mais-Ripken gamma functions.
Definition: Twiss.cpp:1210
Twiss::revPath
bool revPath
Definition: Twiss.h:307
Twiss::Row::orbit
FVector< double, 6 > orbit
The closed orbit after the element.
Definition: Twiss.h:76
Twiss::begin
TLine::const_iterator begin() const
Access to first row.
Definition: Twiss.cpp:551
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:123
RangeRep::leave
void leave(const FlaggedElmPtr &) const
Leave an element or line.
Definition: RangeRep.cpp:81
OpalException
The base class for all OPAL exceptions.
Definition: OpalException.h:28
Twiss::getRow
virtual std::vector< double > getRow(const PlaceRep &, const std::vector< std::string > &)
Return a table row, possible user-defined.
Definition: Twiss.cpp:738
Twiss::Row::Row
Row(ElementBase *, int)
Definition: Twiss.cpp:369
Physics::two_pi
constexpr double two_pi
The value of .
Definition: Physics.h:34
Twiss::getBETik
double getBETik(const Row &, int i1, int i2) const
Mais-Ripken beta functions.
Definition: Twiss.cpp:1164
Twiss::itsLine
std::string itsLine
Definition: Twiss.h:310
TBeamline::getElementLength
virtual double getElementLength() const
Get design length.
Definition: TBeamline.h:311
RangeRep::enter
void enter(const FlaggedElmPtr &) const
Enter an element or line.
Definition: RangeRep.cpp:71
Twiss::Row::getCO
const FVector< double, 6 > & getCO() const
Closed orbit.
Definition: Twiss.cpp:384
Twiss::findRow
Row & findRow(const PlaceRep &row)
Definition: Twiss.cpp:861
Attributes::makeRange
Attribute makeRange(const std::string &name, const std::string &help)
Create a range attribute.
Definition: Attributes.cpp:171
Beam::getET
double getET() const
Return emittance for mode 3.
Definition: Beam.cpp:178
Flatten
Flatten a beamline.
Definition: Flatten.h:42
OpalData::update
void update()
Update all objects.
Definition: OpalData.cpp:723
Object::itsAttr
std::vector< Attribute > itsAttr
The object attributes (see Attribute.hh).
Definition: Object.h:214
LinearMapper
Build a map using a linear map for each element.
Definition: LinearMapper.h:68
PlaceRep::isActive
bool isActive() const
Return status.
Definition: PlaceRep.cpp:98
Table::dynamic
bool dynamic
Flag dynamic table.
Definition: Table.h:153
LinearFun::makeVariable
static LinearFun makeVariable(int var)
Make variable.
Definition: LinearFun.hpp:91
Options::warn
bool warn
Warn flag.
Definition: Options.cpp:10
Beam::getEY
double getEY() const
Return emittance for mode 2.
Definition: Beam.cpp:173
Twiss::getEigen
double getEigen(const Row &, int i1, int i2) const
Eigenvectors.
Definition: Twiss.cpp:1248
Twiss::beam
const Beam * beam
Definition: Twiss.h:296
ThinMapper
Construct thin lens map.
Definition: ThinMapper.h:43
Twiss::itsTable
TLine * itsTable
Definition: Twiss.h:290
Twiss::getCurrent
const Row & getCurrent() const
Return current table row in iteration.
Definition: Twiss.cpp:664
Attributes::getBool
bool getBool(const Attribute &attr)
Return logical value.
Definition: Attributes.cpp:66
Twiss::printTableBody
void printTableBody(std::ostream &, const CellArray &) const
Print the body to this TWISS table.
Definition: Twiss.cpp:786
RangeRep::isActive
bool isActive() const
Test for active range.
Definition: RangeRep.cpp:66
OpalData::printTitle
void printTitle(std::ostream &)
Print the page title.
Definition: OpalData.cpp:705
PlaceRep
Representation of a place within a beam line or sequence.
Definition: PlaceRep.h:41
Twiss::curly_A
FMatrix< double, 6, 6 > curly_A
The initial curly A matrix.
Definition: Twiss.h:272
OpalData::getInstance
static OpalData * getInstance()
Definition: OpalData.cpp:209
Twiss::Twiss
Twiss(int size, const char *name, const char *help)
Exemplar constructor.
Definition: Twiss.cpp:407
Object::getOpalName
const std::string & getOpalName() const
Return object name.
Definition: Object.cpp:284
RangeRep::initialize
void initialize()
Initialise data for search.
Definition: RangeRep.cpp:55
BeamSequence::find
static BeamSequence * find(const std::string &name)
Find a BeamSequence by name.
Definition: BeamSequence.cpp:36
Twiss::getMUi
double getMUi(const Row &, int i1, int=0) const
Three modes, &quot;naive&quot; Twiss functions.
Definition: Twiss.cpp:1106
Twiss::Row::matrix
FMatrix< double, 6, 6 > matrix
The transfer matrix up to and including the element.
Definition: Twiss.h:79
LinearMap::constantTerm
FVector< T, N > constantTerm(const FVector< T, N > &P) const
Evaluate map at point [b]P[/b].
Definition: LinearMap.hpp:220
ElementBase::getElementLength
virtual double getElementLength() const
Get design length.
Definition: ElementBase.h:511
PlaceRep::enter
void enter(const FlaggedElmPtr &) const
Enter an element or line.
Definition: PlaceRep.cpp:70
FlaggedElmPtr::setCounter
void setCounter(int) const
Set clone counter.
Definition: FlaggedElmPtr.h:111
Beam::find
static Beam * find(const std::string &name)
Find named BEAM.
Definition: Beam.cpp:150
Twiss::getDefault
virtual CellArray getDefault() const
Return the default print columns.
Definition: Twiss.cpp:633
BeamSequence::fetchLine
virtual Beamline * fetchLine() const =0
Return the embedded CLASSIC beam line.
Twiss::getColumn
virtual std::vector< double > getColumn(const RangeRep &range, const std::string &col)
Return column [b]col[/b] of this table, limited by [b]range[/b].
Definition: Twiss.cpp:645
Twiss::getLength
virtual double getLength()
Return the length of the table.
Definition: Twiss.cpp:694
Expressions::Scalar::clone
virtual Scalar< T > * clone() const =0
Copy scalar expression.
Twiss::getCell
virtual double getCell(const PlaceRep &row, const std::string &col)
Return a selected value in a selected row.
Definition: Twiss.cpp:626
Twiss::execute
virtual void execute()
Check validity of the table definition.
Definition: Twiss.cpp:566
atan2
PETE_TBTree< FnArcTan2, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > atan2(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
Definition: IpplExpressions.h:244
Twiss::getBETi
double getBETi(const Row &, int i1, int=0) const
Definition: Twiss.cpp:1111
Twiss::orbit
FVector< double, 6 > orbit
The initial closed orbit.
Definition: Twiss.h:269
AbstractMapper::setMap
virtual void setMap(const LinearMap< double, 6 > &)=0
Reset the linear part of the accumulated map for restart.
Flatten::execute
virtual void execute()
Apply the algorithm to the top-level beamline.
Definition: Flatten.h:92
PlaceRep::initialize
void initialize()
Initialise data for search.
Definition: PlaceRep.cpp:64
RangeRep
Representation of a range within a beam line or sequence.
Definition: RangeRep.h:34
PlaceRep::leave
void leave(const FlaggedElmPtr &) const
Leave an element or line.
Definition: PlaceRep.cpp:84
Twiss
Class Twiss.
Definition: Twiss.h:41
Beamline
An abstract sequence of beam line components.
Definition: Beamline.h:37
Twiss::getDisp
double getDisp(const Row &, int i1, int=0) const
Dispersion.
Definition: Twiss.cpp:1234
Table::printTable
virtual void printTable(std::ostream &, const CellArray &) const =0
Print list for the table.
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
Twiss::Row::getMatrix
const FMatrix< double, 6, 6 > & getMatrix() const
Transfer matrix.
Definition: Twiss.cpp:389
BeamSequence
The base class for all OPAL beam lines and sequences.
Definition: BeamSequence.h:32
Beam::getEX
double getEX() const
Return emittance for mode 1.
Definition: Beam.cpp:168
Twiss::Row
Structure for a row of the Twiss table.
Definition: Twiss.h:49
Twiss::getEX
double getEX() const
Return emittance for mode 1.
Definition: Twiss.cpp:679
Beam::getReference
const PartData & getReference() const
Return the embedded CLASSIC PartData.
Definition: Beam.cpp:183
Twiss::reference
const PartData * reference
Definition: Twiss.h:299
Twiss::itsMapper
AbstractMapper * itsMapper
Definition: Twiss.h:293
Twiss::getCO
double getCO(const Row &, int i1, int=0) const
Closed orbit.
Definition: Twiss.cpp:1229
Attributes::makeBool
Attribute makeBool(const std::string &name, const std::string &help)
Make logical attribute.
Definition: Attributes.cpp:56
name
const std::string name
Definition: MaxNormRadialPeak.cpp:3
Table::Cell
Descriptor for printing a table cell.
Definition: Table.h:47
LinearMap::linearTerms
FMatrix< T, N, N > linearTerms() const
Extract linear terms at origin.
Definition: LinearMap.hpp:243
Twiss::Row::getS
double getS() const
Arc length.
Definition: Twiss.cpp:394
Twiss::order
int order
Definition: Twiss.h:302
Twiss::printTableTitle
void printTableTitle(std::ostream &, const char *title) const
Print standard information about the TWISS table.
Definition: Twiss.cpp:847
Twiss::getMatrix
FMatrix< double, 6, 6 > getMatrix(const Row &) const
Accumulated transfer map.
Definition: Twiss.cpp:704
Expressions::Scalar::evaluate
virtual T evaluate() const =0
Evaluate.
Twiss::getCurlyA
FMatrix< double, 6, 6 > getCurlyA() const
Return initial curly A matrix.
Definition: Twiss.cpp:669
Twiss::getLine
virtual const Beamline * getLine() const
Return embedded CLASSIC beamline.
Definition: Twiss.cpp:699
mslang::iterator
std::string::iterator iterator
Definition: MSLang.h:16
LinearFun
Linear function in N variables of type T.
Definition: LinearFun.h:39
Twiss::getALFik
double getALFik(const Row &, int i1, int i2) const
Mais-Ripken alpha functions.
Definition: Twiss.cpp:1183
Twiss::revTrack
bool revTrack
Definition: Twiss.h:306
Attributes::getReal
double getReal(const Attribute &attr)
Return real value.
Definition: Attributes.cpp:217
Twiss::getOrbit
FVector< double, 6 > getOrbit() const
Return initial closed orbit.
Definition: Twiss.cpp:714
Table
The base class for all OPAL tables.
Definition: Table.h:42
Attributes::makeString
Attribute makeString(const std::string &name, const std::string &help)
Make string attribute.
Definition: Attributes.cpp:296
Twiss::getALFi
double getALFi(const Row &, int i1, int=0) const
Definition: Twiss.cpp:1137
Attributes::makeReal
Attribute makeReal(const std::string &name, const std::string &help)
Make real attribute.
Definition: Attributes.cpp:205
ThickMapper
Build a map using a finite-length lens for each element.
Definition: ThickMapper.h:84
Twiss::getS
double getS(const Row &, int=0, int=0) const
Arc length for given row.
Definition: Twiss.cpp:1101
AbstractMapper::getMap
virtual void getMap(LinearMap< double, 6 > &) const =0
Return the linear part of the accumulated map.
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
Twiss::put
void put()
Definition: Twiss.cpp:878
FlaggedElmPtr
A section of a beam line.
Definition: FlaggedElmPtr.h:36
Attributes::getString
std::string getString(const Attribute &attr)
Get string value.
Definition: Attributes.cpp:307
Twiss::Row::getMUi
double getMUi(int i) const
Phase for mode i.
Definition: Twiss.cpp:399
Generated on Sat Jan 4 2020 14:48:07 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5