FLAME  devel
 All Classes Functions Variables Typedefs Enumerations Pages
moment_sup.cpp
1 
2 #include <boost/numeric/ublas/lu.hpp>
3 
4 #include "flame/constants.h"
5 #include "flame/moment.h"
6 
7 #define sqr(x) ((x)*(x))
8 #define cube(x) ((x)*(x)*(x))
9 
10 // http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?LU_Matrix_Inversion
11 // by LU-decomposition.
12 void inverse(MomentElementBase::value_t& out, const MomentElementBase::value_t& in)
13 {
14  using boost::numeric::ublas::permutation_matrix;
15  using boost::numeric::ublas::lu_factorize;
16  using boost::numeric::ublas::lu_substitute;
17  using boost::numeric::ublas::identity_matrix;
18 
19  MomentElementBase::value_t scratch(in); // copy
20  permutation_matrix<size_t> pm(scratch.size1());
21  if(lu_factorize(scratch, pm)!=0)
22  throw std::runtime_error("Failed to invert matrix");
23  out.assign(identity_matrix<double>(scratch.size1()));
24  //out = identity_matrix<double>(scratch.size1());
25  lu_substitute(scratch, pm, out);
26 }
27 
28 void RotMat(const double dx, const double dy,
29  const double theta_x, const double theta_y, const double theta_z,
30  typename MomentElementBase::value_t &R)
31 {
32  typedef typename MomentElementBase::state_t state_t;
33 
34  MomentState::matrix_t T = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
35 
36  R = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
37 
38  // Left-handed coordinate system => theta_y -> -theta_y.
39 
40  double m11 = cos(-theta_y)*cos(theta_z),
41  m12 = sin(theta_x)*sin(-theta_y)*cos(theta_z) + cos(theta_x)*sin(theta_z),
42  m13 = -cos(theta_x)*sin(-theta_y)*cos(theta_z) + sin(theta_x)*sin(theta_z),
43 
44  m21 = -cos(-theta_y)*sin(theta_z),
45  m22 = -sin(theta_x)*sin(-theta_y)*sin(theta_z) + cos(theta_x)*cos(theta_z),
46  m23 = cos(theta_x)*sin(-theta_y)*sin(theta_z) + sin(theta_x)*cos(theta_z),
47 
48  m31 = sin(-theta_y),
49  m32 = -sin(theta_x)*cos(-theta_y),
50  m33 = cos(theta_x)*cos(-theta_y);
51 
52  R(0, 0) = m11, R(0, 2) = m12, R(0, 4) = m13;
53  R(2, 0) = m21, R(2, 2) = m22, R(2, 4) = m23;
54  R(4, 0) = m31, R(4, 2) = m32, R(4, 4) = m33;
55 
56  R(1, 1) = m11, R(1, 3) = m12, R(1, 5) = m13;
57  R(3, 1) = m21, R(3, 3) = m22, R(3, 5) = m23;
58  R(5, 1) = m31, R(5, 3) = m32, R(5, 5) = m33;
59 
60  T(0, 6) = -dx, T(2, 6) = -dy;
61 
62  R = prod(R, T);
63 }
64 
65 
66 void GetQuadMatrix(const double L, const double K, const unsigned ind, typename MomentElementBase::value_t &M)
67 {
68  // 2D quadrupole transport matrix.
69  double sqrtK,
70  psi,
71  cs,
72  sn;
73 
74  if (K > 0e0) {
75  // Focusing.
76  sqrtK = sqrt(K);
77  psi = sqrtK*L;
78  cs = ::cos(psi);
79  sn = ::sin(psi);
80 
81  M(ind, ind) = M(ind+1, ind+1) = cs;
82  if (sqrtK != 0e0)
83  M(ind, ind+1) = sn/sqrtK;
84  else
85  M(ind, ind+1) = L;
86  if (sqrtK != 0e0)
87  M(ind+1, ind) = -sqrtK*sn;
88  else
89  M(ind+1, ind) = 0e0;
90  } else {
91  // Defocusing.
92  sqrtK = sqrt(-K);
93  psi = sqrtK*L;
94  cs = ::cosh(psi);
95  sn = ::sinh(psi);
96 
97  M(ind, ind) = M(ind+1, ind+1) = cs;
98  if (sqrtK != 0e0)
99  M(ind, ind+1) = sn/sqrtK;
100  else
101  M(ind, ind+1) = L;
102  if (sqrtK != 0e0)
103  M(ind+1, ind) = sqrtK*sn;
104  else
105  M(ind+1, ind) = 0e0;
106  }
107 }
108 
109 
110 void GetEdgeMatrix(const double rho, const double phi, typename MomentElementBase::value_t &M)
111 {
112  typedef typename MomentElementBase::state_t state_t;
113 
114  M = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
115 
116  M(state_t::PS_PX, state_t::PS_X) = tan(phi)/rho;
117  M(state_t::PS_PY, state_t::PS_Y) = -tan(phi)/rho;
118 }
119 
120 
121 void GetEEdgeMatrix(const double fringe_x, const double fringe_y, const double kappa, typename MomentElementBase::value_t &M)
122 {
123  // Edge focusing for electrostatic dipole.
124  typedef typename MomentElementBase::state_t state_t;
125 
126  M = boost::numeric::ublas::identity_matrix<double>(state_t::maxsize);
127 
128  M(state_t::PS_PX, state_t::PS_X) = fringe_x;
129  M(state_t::PS_PX, state_t::PS_PX) = sqrt(1e0+kappa);
130  M(state_t::PS_PY, state_t::PS_Y) = fringe_y;
131  M(state_t::PS_PS, state_t::PS_PS) = 1e0+kappa;
132 }
133 
134 
135 void GetSBendMatrix(const double L, const double phi, const double phi1, const double phi2, const double K,
136  const double IonEs, const double ref_gamma, const double qmrel,
137  const double dip_beta, const double dip_gamma, const double d, const double dip_IonK, typename MomentElementBase::value_t &M)
138 {
139  typedef typename MomentElementBase::state_t state_t;
140 
141  MomentState::matrix_t edge1, edge2, R;
142 
143  double rho = L/phi,
144  Kx = K + 1e0/sqr(rho),
145  Ky = -K,
146  dx = 0e0,
147  sx = 0e0;
148 
149  // Horizontal plane.
150  GetQuadMatrix(L, Kx, (unsigned)state_t::PS_X, M);
151  // Vertical plane.
152  GetQuadMatrix(L, Ky, (unsigned)state_t::PS_Y, M);
153 
154  // Include dispersion.
155  if (Kx == 0e0) {
156  dx = sqr(L)/(2e0*rho);
157  sx = L/rho;
158  } else if (Kx > 0e0) {
159  dx = (1e0-cos(sqrt(Kx)*L))/(rho*Kx);
160  sx = sin(sqrt(Kx)*L)/(rho*sqrt(Kx));
161  } else {
162  dx = (1e0-cosh(sqrt(-Kx)*L))/(rho*Kx);
163  sx = sin(sqrt(Kx)*L)/(rho*sqrt(Kx));
164  }
165 
166  M(state_t::PS_X, state_t::PS_PS) = dx/(sqr(dip_beta)*dip_gamma*IonEs/MeVtoeV);
167  M(state_t::PS_PX, state_t::PS_PS) = sx/(sqr(dip_beta)*dip_gamma*IonEs/MeVtoeV);
168 
169  M(state_t::PS_S, state_t::PS_X) = sx*dip_IonK;
170  M(state_t::PS_S, state_t::PS_PX) = dx*dip_IonK;
171  // Low beta approximation.
172  M(state_t::PS_S, state_t::PS_PS) =
173  ((L/rho-sx)/(Kx*rho)-L/sqr(ref_gamma))*dip_IonK
174  /(sqr(dip_beta)*dip_gamma*IonEs/MeVtoeV);
175 
176  // Add dipole terms.
177  M(state_t::PS_S, 6) = ((L/rho-sx)/(Kx*rho)*d-L/sqr(ref_gamma)*(d+qmrel))*dip_IonK;
178  M(state_t::PS_X, 6) = dx*d;
179  M(state_t::PS_PX, 6) = sx*d;
180 
181  // Edge focusing.
182  GetEdgeMatrix(rho, phi1, edge1);
183  GetEdgeMatrix(rho, phi2, edge2);
184 
185  M = prod(M, edge1);
186  M = prod(edge2, M);
187 
188  // Longitudinal plane.
189  // For total path length.
190  // M(state_t::PS_S, state_t::PS_S) = L;
191 }
192 
193 
194 void GetSolMatrix(const double L, const double K, typename MomentElementBase::value_t &M)
195 {
196  typedef typename MomentElementBase::state_t state_t;
197 
198  double C = ::cos(K*L),
199  S = ::sin(K*L);
200 
201  M(state_t::PS_X, state_t::PS_X)
202  = M(state_t::PS_PX, state_t::PS_PX)
203  = M(state_t::PS_Y, state_t::PS_Y)
204  = M(state_t::PS_PY, state_t::PS_PY)
205  = sqr(C);
206 
207  if (K != 0e0)
208  M(state_t::PS_X, state_t::PS_PX) = S*C/K;
209  else
210  M(state_t::PS_X, state_t::PS_PX) = L;
211  M(state_t::PS_X, state_t::PS_Y) = S*C;
212  if (K != 0e0)
213  M(state_t::PS_X, state_t::PS_PY) = sqr(S)/K;
214  else
215  M(state_t::PS_X, state_t::PS_PY) = 0e0;
216 
217  M(state_t::PS_PX, state_t::PS_X) = -K*S*C;
218  M(state_t::PS_PX, state_t::PS_Y) = -K*sqr(S);
219  M(state_t::PS_PX, state_t::PS_PY) = S*C;
220 
221  M(state_t::PS_Y, state_t::PS_X) = -S*C;
222  if (K != 0e0)
223  M(state_t::PS_Y, state_t::PS_PX) = -sqr(S)/K;
224  else
225  M(state_t::PS_Y, state_t::PS_PX) = 0e0;
226  if (K != 0e0)
227  M(state_t::PS_Y, state_t::PS_PY) = S*C/K;
228  else
229  M(state_t::PS_Y, state_t::PS_PY) = L;
230 
231  M(state_t::PS_PY, state_t::PS_X) = K*sqr(S);
232  M(state_t::PS_PY, state_t::PS_PX) = -S*C;
233  M(state_t::PS_PY, state_t::PS_Y) = -K*S*C;
234 
235  // Longitudinal plane.
236  // For total path length.
237 // M(state_t::PS_S, state_t::PS_S) = L;
238 }
239 
240 
241 void GetEBendMatrix(const double L, const double phi, const double fringe_x, const double fringe_y, const double kappa, const double Kx, const double Ky,
242  const double IonEs, const double ref_beta, const double real_gamma, const double eta0, const double h,
243  const double dip_beta, const double dip_gamma, const double delta_KZ, const double SampleIonK, typename MomentElementBase::value_t &M)
244 {
245  typedef typename MomentElementBase::state_t state_t;
246 
247  MomentState::matrix_t edge;
248 
249  double rho = L/phi,
250  scl = (real_gamma - 1e0)*IonEs/MeVtoeV,
251  dx = 0e0,
252  sx = 0e0;
253 
254  // Horizontal plane.
255  GetQuadMatrix(L, Kx, (unsigned)state_t::PS_X, M);
256  // Vertical plane.
257  GetQuadMatrix(L, Ky, (unsigned)state_t::PS_Y, M);
258 
259  // Include dispersion.
260  if (Kx == 0e0) {
261  dx = 0e0;
262  sx = 0e0;
263  } else if (Kx > 0e0) {
264  dx = (1e0-cos(sqrt(Kx)*L))/(rho*Kx);
265  sx = sin(sqrt(Kx)*L)/sqrt(Kx);
266  } else {
267  dx = (1e0-cosh(sqrt(-Kx)*L))/(rho*Kx);
268  sx = sin(sqrt(Kx)*L)/sqrt(Kx);
269  }
270 
271  double Nk = (sqr(1e0+2e0*eta0)+h)/(2e0*(1e0+eta0)*(1e0+2e0*eta0)),
272  Nt = 1e0 + h/sqr(1e0+2e0*eta0),
273  CorT = -real_gamma/(1e0+real_gamma),
274  tx = sx/rho*Nt*CorT,
275  txp = dx*Nt*CorT,
276  tzp = (-L/(2e0*(1e0+eta0)*(1e0+2e0*eta0))+((L-sx)/Kx/sqr(rho))*Nk*Nt)*CorT;
277 
278  M(state_t::PS_X, state_t::PS_PS) = dx*Nk/scl;
279  M(state_t::PS_PX, state_t::PS_PS) = sx/rho*Nk/scl;
280 
281  M(state_t::PS_S, state_t::PS_X) = -tx*SampleIonK;
282  M(state_t::PS_S, state_t::PS_PX) = -txp*SampleIonK;
283  // Low beta approximation.
284  M(state_t::PS_S, state_t::PS_PS) = -tzp*SampleIonK/scl;
285 
286  // Add dipole terms.
287  double delta_K = sqr(ref_beta/dip_beta) - 1e0;
288 
289 // M(state_t::PS_X, 6) = dx*Nk*(delta_K+delta_KZ);
290 // M(state_t::PS_PX, 6) = sx/rho*Nk*(delta_K+delta_KZ);
291 
292  double Nkz = (1e0+2e0*eta0-h)/(2e0*(1e0+eta0)*(1e0+2e0*eta0));
293  M(state_t::PS_X, 6) = dx*(Nk*delta_K+Nkz*delta_KZ);
294  M(state_t::PS_PX, 6) = sx/rho*(Nk*delta_K+Nkz*delta_KZ);
295 
296  // Edge focusing.
297  GetEEdgeMatrix(fringe_x, fringe_y ,kappa, edge);
298 
299  M = prod(M, edge);
300  M = prod(edge, M);
301 
302  // Longitudinal plane.
303  // For total path length.
304  // M(state_t::PS_S, state_t::PS_S) = L;
305 }
Generated on Mon Jul 18 2016 13:55:44 for FLAME by   doxygen 1.8.8