source: TRACY3/trunk/tracy/tracy/inc/field.h @ 32

Last change on this file since 32 was 32, checked in by zhangj, 10 years ago

active the transport of the twiss functions and orbits of the transfer line.
  • Property svn:executable set to *
File size: 12.1 KB
Line 
1 /* Author:      Johan Bengtsson
2
3   Definitions:  Polymorphic number class.              */
4
5
6const int  max_str = 132;
7
8const int  n_m2    = 21;  // no of 2nd moments
9
10#if ORDER == 1
11const int  ss_dim  = 6;   // state space dimension
12#else
13const int  ss_dim  = 6+1; /* state space dimension:
14                             phase space and parameter dependence */
15#endif
16
17// spatial components
18enum spatial_index { X_ = 0, Y_ = 1, Z_ = 2 };
19
20// phase space components
21// (Note, e.g. spin components should be added here)
22enum ps_index { x_ = 0, px_ = 1, y_ = 2, py_ = 3, delta_ = 4, ct_ = 5 };
23
24#if ORDER == 1
25typedef double  tps_buf[ss_dim+1]; // const. and linear terms
26#endif
27
28template<typename T> class ss_vect;
29
30typedef int  iVector[ss_dim];
31
32// Polymorphic class for floating point and TPSA
33
34struct MNF_struct;
35
36class tps {
37 public:
38  tps(void);
39  tps(const double);
40  tps(const double, const int);
41  tps(const tps &);
42  ~tps(void);
43
44  // initialize TPSA library
45   friend void TPSAEps(const double);
46  // trace level for TPSALib and LieLib
47  friend void idprset(const int);
48
49  const double cst(void) const;
50  double operator[](const int) const;
51  double operator[](const int []) const;
52  void pook(const int [], const double);
53
54  tps& operator=(const double);
55  tps& operator+=(const double);
56  tps& operator-=(const double);
57  tps& operator*=(const double);
58  tps& operator/=(const double);
59
60  tps& operator=(const tps &);
61  tps& operator+=(const tps &);
62  tps& operator-=(const tps &);
63  tps& operator*=(const tps &);
64  tps& operator/=(const tps &);
65
66  friend double getmat(const ss_vect<tps> &map, const int i, const int j);
67  friend void putmat(ss_vect<tps> &map, const int i, const int j,
68                     const double r);
69#if ORDER == 1
70  friend void dacct_(const ss_vect<tps> &x, const int i,
71                     const ss_vect<tps> &y, const int j,
72                     ss_vect<tps> &z, const int k);
73#endif
74
75  friend istream& operator>>(istream &, tps &);
76  friend ostream& operator<<(ostream &, const tps &);
77
78  friend double abs(const tps &);
79  friend double abs2(const tps &);
80  friend tps sqrt(const tps &);
81  friend tps sqr(const tps &);
82  friend tps pow(const tps &, const int);
83  friend tps exp(const tps &);
84  friend tps log(const tps &);
85  friend tps sin(const tps &);
86  friend tps cos(const tps &);
87  friend tps tan(const tps &);
88  friend tps ctan(const tps &);
89  friend tps asin(const tps &);
90  friend tps acos(const tps &);
91  friend tps atan(const tps &);
92  friend tps actan(const tps &);
93  friend tps sinh(const tps &);
94  friend tps cosh(const tps &);
95
96  friend tps Der(const tps &, const int);
97  friend tps LieExp(const tps &, const tps &);
98  friend tps PB(const tps &, const tps &);
99  friend tps Take(const tps &, const int);
100
101  // R(nd2, nv) = P(nd2, nd2)*Q(nd2, nv)
102  friend ss_vect<tps> operator*(const ss_vect<tps> &, const ss_vect<tps> &);
103  // R(nv, nv) = P(nv, nv)*Q(nv, nv)
104  friend void CCT(const tps x[], const int n_x, const tps y[], const int n_y,
105                  tps z[], const int n_z);
106  friend tps operator*(const tps &, const ss_vect<tps> &);
107
108  friend ss_vect<tps> FExpo(const tps &, const ss_vect<tps> &,
109                            const int, const int, const int);
110  friend ss_vect<tps> LieExp(const tps &, const ss_vect<tps> &);
111  // Q(nv, nv) = P(nd2, nd2)^-1
112  friend ss_vect<tps> Inv(const ss_vect<tps> &);
113  // Q(nv, nv) = P(nv, nv)^-1
114  friend ss_vect<tps> Inv_Ext(const ss_vect<tps> &);
115  friend ss_vect<tps> PInv(const ss_vect<tps> &, const iVector &);
116  friend void GoFix(const ss_vect<tps> &, ss_vect<tps> &,
117                    ss_vect<tps> &, const int);
118  friend MNF_struct MapNorm(const ss_vect<tps> &, const int);
119  friend ss_vect<tps> MapNormF(const ss_vect<tps> &, ss_vect<tps> &,
120                               ss_vect<tps> &, ss_vect<tps> &,
121                               ss_vect<tps> &, const int, const int);
122  friend ss_vect<tps> dHdJ(const tps &);
123  friend void CtoR(const tps &, tps &, tps &);
124  friend tps RtoC(const tps &, const tps &);
125  friend tps LieFact_DF(const ss_vect<tps> &, ss_vect<tps> &);
126  friend ss_vect<tps> FlowFact(const ss_vect<tps> &);
127  friend tps Intd(const ss_vect<tps> &, const double);
128  friend ss_vect<tps> Difd(const tps &, const double);
129  friend ss_vect<tps> Taked(const ss_vect<tps> &, const int);
130 private:
131#if ORDER == 1
132  tps_buf  ltps;  // linear TPSA
133#else
134  int     intptr; // index used by Fortran implementation
135#endif
136  double  r;      // floating-point calc. if intptr = 0
137};
138
139
140// Class for single particle phase space dynamics
141
142// pre-declare template friend functions: f<>()
143template<typename T>
144ss_vect<T> operator+(const ss_vect<T> &);
145template<typename T>
146ss_vect<T> operator-(const ss_vect<T> &);
147
148template<typename T> class ss_vect {
149 public:
150  typedef T value_type;
151
152  ss_vect(void);
153// Let's the compiler synthetize the copy constructor
154//  ss_vect(const T &a) { }
155//  ss_vect(const ss_vect<T> &a) { }
156  template<typename U>
157    ss_vect(const U &);
158  template<typename U>
159    ss_vect(const ss_vect<U> &);
160
161  ss_vect<double> cst(void) const;
162  T& operator[](const int i) { return ss[i]; }
163  const T& operator[](const int i) const { return ss[i]; }
164
165  ss_vect<T>& operator*=(const double);
166  ss_vect<T>& operator*=(const tps &);
167
168  ss_vect<T>& operator=(const ss_vect<T> &);
169  ss_vect<T>& operator+=(const ss_vect<T> &);
170  ss_vect<T>& operator-=(const ss_vect<T> &);
171
172  friend ss_vect<T> operator+<>(const ss_vect<T> &);
173  friend ss_vect<T> operator-<>(const ss_vect<T> &);
174
175  friend ss_vect<double> operator+(const ss_vect<double> &,
176                                   const ss_vect<double> &);
177  friend ss_vect<tps> operator+(const ss_vect<tps> &, const ss_vect<double> &);
178  friend ss_vect<tps> operator+(const ss_vect<double> &, const ss_vect<tps> &);
179  friend ss_vect<tps> operator+(const ss_vect<tps> &, const ss_vect<tps> &);
180
181  friend ss_vect<double> operator-(const ss_vect<double> &,
182                                   const ss_vect<double> &);
183  friend ss_vect<tps> operator-(const ss_vect<tps> &, const ss_vect<double> &);
184  friend ss_vect<tps> operator-(const ss_vect<double> &, const ss_vect<tps> &);
185  friend ss_vect<tps> operator-(const ss_vect<tps> &, const ss_vect<tps> &);
186
187//  friend ss_vect<double> operator*(const ss_vect<tps> &,
188//                                 const ss_vect<double> &);
189  // R(nd2, nv) = P(nd2, nd2)*Q(nd2, nv)
190  friend ss_vect<tps> operator*(const ss_vect<tps> &, const ss_vect<tps> &);
191  // R(nv, nv) = P(nv, nv)*Q(nv, nv)
192  friend ss_vect<tps> CCT(const ss_vect<tps> &, const ss_vect<tps> &);
193  friend tps operator*(const tps &, const ss_vect<tps> &);
194
195  template<typename CharT, class Traits>
196    friend basic_istream<CharT, Traits>&
197    operator>>(basic_istream<CharT, Traits> &, ss_vect<tps> &);
198
199  template<typename CharT, class Traits>
200    friend basic_ostream<CharT, Traits>&
201    operator<<(basic_ostream<CharT, Traits> &, const ss_vect<T> &);
202
203  ss_vect<T> zero(void);
204  ss_vect<tps> identity(void);
205
206  friend ss_vect<tps> FExpo(const tps &, const ss_vect<tps> &,
207                            const int, const int, const int);
208  friend ss_vect<tps> LieExp(const tps &, const ss_vect<tps> &);
209  // Q(nv, nv) = P(nd2, nd2)^-1
210  friend ss_vect<tps> Inv(const ss_vect<tps> &);
211  // Q(nv, nv) = P(nv, nv)^-1
212  friend ss_vect<tps> Inv_Ext(const ss_vect<tps> &);
213  friend ss_vect<tps> PInv(const ss_vect<tps> &, const iVector &);
214  friend void GoFix(const ss_vect<tps> &, ss_vect<tps> &,
215                    ss_vect<tps> &, const int);
216  friend MNF_struct MapNorm(const ss_vect<tps> &, const int);
217  friend ss_vect<tps> MapNormF(const ss_vect<tps> &, ss_vect<tps> &,
218                               ss_vect<tps> &, ss_vect<tps> &,
219                               ss_vect<tps> &, const int, const int);
220  friend void dHdJ(const tps &, ss_vect<tps> &);
221  friend void CtoR(const tps &, tps &, tps &);
222  friend tps RtoC(const tps &, const tps &);
223  friend tps LieFact_DF(const ss_vect<tps> &, ss_vect<tps> &);
224  friend tps LieFact(const ss_vect<tps> &);
225  friend ss_vect<tps> FlowFact(const ss_vect<tps> &);
226  friend tps Intd(const ss_vect<tps> &, const double);
227  friend ss_vect<tps> Difd(const tps &, const double);
228  friend ss_vect<tps> Taked(const ss_vect<tps> &, const int);
229 private:
230  // (Note, e.g. spin components should be added here)
231  T  ss[ss_dim];
232};
233
234
235typedef ss_vect<double>  Vector;                 /*6*1 vector*/
236typedef Vector           Matrix[ss_dim];         /*6*6 matrix*/
237
238
239typedef struct MNF_struct
240{
241  tps           K;       // new effective Hamiltonian
242  ss_vect<tps>  A0;      // transformation to fixed point
243  ss_vect<tps>  A1;      // (linear) transformation to Floquet space
244  tps           g;       /* generator for nonlinear transformation to
245                            Floquet space */
246  ss_vect<tps>  map_res; // residual map
247} MNF_struct;
248
249
250// partial template-class specialization
251// primary version: is_double<>
252template<typename T>
253class is_double { };
254
255// partial specialization
256template<>
257class is_double<double> {
258 public:
259  static inline double cst(const double x) { return x; }
260};
261
262// partial specialization
263template<>
264class is_double<tps> {
265 public:
266  static inline double cst(const tps &x) { return x.cst(); }
267};
268
269
270template<typename T>
271inline T sqr(const T &a)
272{ return a*a; }
273
274//overloading operator functions, specific for
275//the user-defined data type tps
276inline tps operator+(const tps &a, const tps &b)
277{ return tps(a) += b; }
278
279inline tps operator+(const tps &a, const double b)
280{ return tps(a) += b; }
281
282inline tps operator+(const double a, const tps &b)
283{ return tps(a) += b; }
284
285inline tps operator-(const tps &a, const tps &b)
286{ return tps(a) -= b; }
287
288inline tps operator-(const tps &a, const double b)
289{ return tps(a) -= b; }
290
291inline tps operator-(const double a, const tps &b)
292{ return tps(a) -= b; }
293
294inline tps operator*(const tps &a, const tps &b)
295{ return tps(a) *= b; }
296
297inline tps operator*(const tps &a, const double b)
298{ return tps(a) *= b; }
299
300inline tps operator*(const double a, const tps &b)
301{ return tps(a) *= b; }
302
303inline tps operator/(const tps &a, const tps &b)
304{ return tps(a) /= b; }
305
306inline tps operator/(const tps &a, const double b)
307{ return tps(a) /= b; }
308
309inline tps operator/(const double a, const tps &b)
310{ return tps(a) /= b; }
311
312
313inline tps operator+(const tps &x)
314{ return tps(x); }
315
316inline tps operator-(const tps &x)
317{ return tps(x) *= -1.0; }
318
319
320inline bool operator>(const tps &a, const tps &b)
321{ return a.cst() > b.cst(); }
322
323inline bool operator>(const tps &a, const double b)
324{ return a.cst() > b; }
325
326inline bool operator>(const double a, const tps &b)
327{ return a > b.cst(); }
328
329
330inline bool operator<(const tps &a, const tps &b)
331{ return a.cst() < b.cst(); }
332
333inline bool operator<(const tps &a, const double b)
334{ return a.cst() < b; }
335
336inline bool operator<(const double a, const tps &b)
337{ return a < b.cst(); }
338
339
340inline bool operator>=(const tps &a, const tps &b)
341{ return a.cst() >= b.cst(); }
342
343inline bool operator>=(const tps &a, const double b)
344{ return a.cst() >= b; }
345
346inline bool operator>=(const double a, const tps &b)
347{ return a >= b.cst(); }
348
349
350inline bool operator<=(const tps &a, const tps &b)
351{ return a.cst() <= b.cst(); }
352
353inline bool operator<=(const tps &a, const double b)
354{ return a.cst() <= b; }
355
356inline bool operator<=(const double a, const tps &b)
357{ return a <= b.cst(); }
358
359
360inline bool operator==(const tps &a, const tps &b)
361{ return a.cst() == b.cst(); }
362
363inline bool operator==(const tps &a, const double b)
364{ return a.cst() == b; }
365
366inline bool operator==(const double a, const tps &b)
367{ return a == b.cst(); }
368
369
370inline bool operator!=(const tps &a, const tps &b)
371{ return a.cst() != b.cst(); }
372
373inline bool operator!=(const tps &a, const double b)
374{ return a.cst() != b; }
375
376inline bool operator!=(const double a, const tps &b)
377{ return a != b.cst(); }
378
379/******************************************************
380template<typename T>
381inline ss_vect<T>::ss_vect(void)
382
383  Purpose:
384     Initialize 6-D vector(T:double) or 6*6 matrix(T:tps) ss, based on
385     the type of T, if ORDER = 1
386      Initialize 7-D vector or 7*7 matrix ss, based on
387     the type of T, if ORDER != 1
388     
389   
390     
391 
392 
393  Parameters:
394      ss_dim:  6      if ORDER = 1
395      ss_dim:  6+1    if ORDER != 1 
396     
397******************************************************/
398
399template<typename T>
400inline ss_vect<T>::ss_vect(void)
401{
402  int  i;
403
404  for (i = 0; i < ss_dim; i++)
405    ss[i] = T();
406}
407
408template<typename T>
409template<typename U>
410inline ss_vect<T>::ss_vect(const U &a) : ss(a) { }
411
412template<typename T>
413template<typename U>
414inline ss_vect<T>::ss_vect(const ss_vect<U> &a)
415{
416  int              i;
417
418  for (i = 0; i < ss_dim; i++)
419    ss[i] = a[i];
420 }
Note: See TracBrowser for help on using the repository browser.