[3] | 1 | /*Tracy-3 |
---|
| 2 | |
---|
| 3 | J. Bengtsson, CBP, LBL 1990 - 1994 Pascal version |
---|
| 4 | SLS, PSI 1995 - 1997 |
---|
| 5 | M. Boege SLS, PSI 1998 C translation |
---|
| 6 | L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps |
---|
| 7 | J. Bengtsson NSLS-II, BNL 2004 - |
---|
| 8 | |
---|
| 9 | Element propagators. */ |
---|
| 10 | |
---|
| 11 | |
---|
| 12 | double c_1, d_1, c_2, d_2, cl_rad, q_fluct; |
---|
| 13 | double I2, I4, I5, dcurly_H, dI4; |
---|
| 14 | ElemFamType ElemFam[Elem_nFamMax]; |
---|
| 15 | CellType Cell[Cell_nLocMax+1]; |
---|
| 16 | |
---|
| 17 | // for IBS |
---|
| 18 | int i_, j_; |
---|
| 19 | double **C_; |
---|
| 20 | |
---|
| 21 | // Dynamic model |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | /****************************************************************************/ |
---|
| 26 | /* void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, |
---|
| 27 | const double c0, const double c1, const double s1) |
---|
| 28 | |
---|
| 29 | Purpose: |
---|
| 30 | Global to local coordinates |
---|
| 31 | |
---|
| 32 | ****************************************************************************/ |
---|
| 33 | template<typename T> |
---|
| 34 | void GtoL(ss_vect<T> &X, Vector2 &S, Vector2 &R, |
---|
| 35 | const double c0, const double c1, const double s1) |
---|
| 36 | { |
---|
| 37 | ss_vect<T> x1; |
---|
| 38 | |
---|
| 39 | /* Simplified rotated p_rot */ |
---|
| 40 | X[px_] += c1; X[3] += s1; |
---|
| 41 | /* Translate */ |
---|
| 42 | X[x_] -= S[X_]; X[y_] -= S[Y_]; |
---|
| 43 | /* Rotate */ |
---|
| 44 | x1 = X; |
---|
| 45 | X[x_] = R[X_]*x1[x_] + R[Y_]*x1[y_]; |
---|
| 46 | X[px_] = R[X_]*x1[px_] + R[Y_]*x1[py_]; |
---|
| 47 | X[y_] = -R[Y_]*x1[x_] + R[X_]*x1[y_]; |
---|
| 48 | X[py_] = -R[Y_]*x1[px_] + R[X_]*x1[py_] ; |
---|
| 49 | /* Simplified p_rot */ |
---|
| 50 | X[px_] -= c0; |
---|
| 51 | } |
---|
| 52 | |
---|
| 53 | /****************************************************************************/ |
---|
| 54 | /* void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R, |
---|
| 55 | double c0, double c1, double s1) |
---|
| 56 | |
---|
| 57 | Purpose: |
---|
| 58 | Local to global coordinates |
---|
| 59 | |
---|
| 60 | ****************************************************************************/ |
---|
| 61 | template<typename T> |
---|
| 62 | void LtoG(ss_vect<T> &X, Vector2 &S, Vector2 &R, |
---|
| 63 | double c0, double c1, double s1) |
---|
| 64 | { |
---|
| 65 | ss_vect<T> x1; |
---|
| 66 | |
---|
| 67 | /* Simplified p_rot */ |
---|
| 68 | X[px_] -= c0; |
---|
| 69 | /* Rotate */ |
---|
| 70 | x1 = X; |
---|
| 71 | X[x_] = R[X_]*x1[x_] - R[Y_]*x1[y_]; |
---|
| 72 | X[px_] = R[X_]*x1[px_] - R[Y_]*x1[py_]; |
---|
| 73 | X[y_] = R[Y_]*x1[x_] + R[X_]*x1[y_]; |
---|
| 74 | X[py_] = R[Y_]*x1[px_] + R[X_]*x1[py_]; |
---|
| 75 | /* Translate */ |
---|
| 76 | X[x_] += S[X_]; |
---|
| 77 | X[y_] += S[Y_]; |
---|
| 78 | /* p_rot rotated */ |
---|
| 79 | X[px_] += c1; |
---|
| 80 | X[py_] += s1; |
---|
| 81 | } |
---|
| 82 | /**********************************************************/ |
---|
| 83 | /* |
---|
| 84 | |
---|
| 85 | Purpose: |
---|
| 86 | Get the longitudinal momentum ps; for the exact/approximation Hamiltonian |
---|
| 87 | |
---|
| 88 | for approximation Hamitonian: |
---|
| 89 | ps = 1+delta |
---|
| 90 | |
---|
[10] | 91 | |
---|
[3] | 92 | For exact hamitonian: |
---|
| 93 | |
---|
[10] | 94 | (1) ps = sqrt((1+delta)^2 - px_^2 - py_^2) |
---|
| 95 | |
---|
| 96 | (2) using the check of TEAPOT to check ps < 0 |
---|
[3] | 97 | **********************************************************/ |
---|
| 98 | template<typename T> |
---|
| 99 | inline T get_p_s(ss_vect<T> &x) |
---|
| 100 | { |
---|
| 101 | T p_s, p_s2; |
---|
| 102 | |
---|
| 103 | if (!globval.H_exact) |
---|
| 104 | p_s = 1.0+x[delta_]; |
---|
| 105 | else { |
---|
| 106 | p_s2 = sqr(1.0+x[delta_]) - sqr(x[px_]) - sqr(x[py_]); |
---|
| 107 | if (p_s2 >= 0.0) |
---|
| 108 | p_s = sqrt(p_s2); |
---|
| 109 | else { |
---|
| 110 | // avoid compile warning |
---|
| 111 | //p_s = 0.0; |
---|
| 112 | p_s = 1.0e-20; |
---|
| 113 | cout << " " << 1+x[delta_] << " " << x[px_] << " " << x[py_] << endl; |
---|
| 114 | printf("get_p_s: *** Speed of light exceeded!\n"); |
---|
[11] | 115 | // exit_(1); |
---|
[3] | 116 | |
---|
| 117 | } |
---|
| 118 | } |
---|
| 119 | return(p_s); |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | /****************************************************************************/ |
---|
| 123 | /* Drift(double L, ss_vect<T> &x) |
---|
| 124 | |
---|
| 125 | Purpose: |
---|
[11] | 126 | Exact Drift pass in the curvilinear coordinate for the dipole |
---|
[3] | 127 | |
---|
[11] | 128 | See Eqn. (10.26) on page 307 in E. Forest's beam dynamics, a new attitude |
---|
| 129 | |
---|
| 130 | H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta |
---|
| 131 | |
---|
| 132 | Input: |
---|
| 133 | L: drift length |
---|
| 134 | &x: pointer to initial vector: x |
---|
| 135 | |
---|
| 136 | Output: |
---|
| 137 | none |
---|
| 138 | |
---|
| 139 | Return: |
---|
| 140 | none |
---|
| 141 | |
---|
| 142 | Global variables: |
---|
| 143 | |
---|
| 144 | |
---|
| 145 | Specific functions: |
---|
| 146 | |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | Comments: |
---|
| 150 | |
---|
| 151 | Test version....rewritten by Jianfeng zhang @ LAL, 31/07/2013 |
---|
| 152 | |
---|
| 153 | ****************************************************************************/ |
---|
| 154 | |
---|
| 155 | |
---|
| 156 | template<typename T> |
---|
| 157 | void Drift(double L,double h_bend, ss_vect<T> &x) { |
---|
| 158 | T u; |
---|
| 159 | ss_vect<T> x0; |
---|
| 160 | |
---|
| 161 | x0=x; |
---|
| 162 | |
---|
| 163 | if (globval.H_exact) { |
---|
| 164 | if(get_p_s(x) == 0) |
---|
| 165 | return; |
---|
| 166 | |
---|
| 167 | x[x_] = x0[x_]/cos(L*h_bend)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend)); |
---|
| 168 | x[px_]= x0[px_]*cos(L*h_bend) + get_p_s(x)*sin(L*h_bend); |
---|
| 169 | x[y_] = x0[y_]+x0[px_]*x0[x_]*tan(L*h_bend)/get_p_s(x)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend)); |
---|
| 170 | x[py_]= x0[py_]; |
---|
| 171 | x[ct_]= x0[ct_]+((1+x0[delta_])*x0[x_]*tan(L*h_bend))/get_p_s(x)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend)); |
---|
| 172 | |
---|
| 173 | |
---|
| 174 | } |
---|
| 175 | } |
---|
| 176 | |
---|
| 177 | /****************************************************************************/ |
---|
| 178 | /* Drift(double L, ss_vect<T> &x) |
---|
| 179 | |
---|
| 180 | Purpose: |
---|
| 181 | Drift pass in the Cartesian coordinate |
---|
| 182 | |
---|
[3] | 183 | If H_exact = false, using approximation Hamiltonian (only valid for big ring): |
---|
| 184 | |
---|
| 185 | px^2+py^2 |
---|
| 186 | H(x,y,cT,px,py,delta) = ------------- |
---|
| 187 | 2(1+delta) |
---|
| 188 | |
---|
| 189 | |
---|
| 190 | Otherwise, using the exact Hamiltonian (valid for small ring) |
---|
| 191 | (J. Bengtsson: Tracy-2 user's mannual): |
---|
| 192 | |
---|
| 193 | H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta |
---|
| 194 | |
---|
| 195 | Input: |
---|
| 196 | L: drift length |
---|
| 197 | &x: pointer to initial vector: x |
---|
| 198 | |
---|
| 199 | Output: |
---|
| 200 | none |
---|
| 201 | |
---|
| 202 | Return: |
---|
| 203 | none |
---|
| 204 | |
---|
| 205 | Global variables: |
---|
| 206 | |
---|
| 207 | |
---|
| 208 | Specific functions: |
---|
| 209 | |
---|
| 210 | ****************************************************************************/ |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | template<typename T> |
---|
| 214 | void Drift(double L, ss_vect<T> &x) { |
---|
| 215 | T u; |
---|
| 216 | |
---|
| 217 | if (globval.H_exact) { |
---|
| 218 | if(get_p_s(x) == 0) |
---|
| 219 | return; |
---|
| 220 | |
---|
| 221 | u = L/get_p_s(x); |
---|
| 222 | x[ct_] += u*(1.0+x[delta_]) - L; |
---|
| 223 | } else { |
---|
| 224 | u = L/(1.0+x[delta_]); |
---|
| 225 | x[ct_] += u*(sqr(x[px_])+sqr(x[py_]))/(2.0*(1.0+x[delta_])); |
---|
| 226 | } |
---|
| 227 | |
---|
| 228 | x[x_] += x[px_] * u; |
---|
| 229 | x[y_] += x[py_] * u; |
---|
| 230 | if (globval.pathlength) |
---|
| 231 | x[ct_] += L; |
---|
| 232 | } |
---|
| 233 | |
---|
| 234 | template<typename T> |
---|
| 235 | void Drift_Pass(CellType &Cell, ss_vect<T> &x) { |
---|
| 236 | Drift(Cell.Elem.PL, x); |
---|
| 237 | } |
---|
| 238 | |
---|
| 239 | /****************************************************************************/ |
---|
| 240 | /* zero_mat(const int n, double** A) |
---|
| 241 | |
---|
| 242 | Purpose: |
---|
| 243 | Initionize n x n matrix A with 0 |
---|
| 244 | |
---|
| 245 | |
---|
| 246 | ****************************************************************************/ |
---|
| 247 | void zero_mat(const int n, double** A) { |
---|
| 248 | int i, j; |
---|
| 249 | |
---|
| 250 | for (i = 1; i <= n; i++) |
---|
| 251 | for (j = 1; j <= n; j++) |
---|
| 252 | A[i][j] = 0.0; |
---|
| 253 | } |
---|
| 254 | |
---|
| 255 | /****************************************************************************/ |
---|
| 256 | /* void identity_mat(const int n, double** A) |
---|
| 257 | |
---|
| 258 | Purpose: |
---|
| 259 | generate n x n identity matrix A |
---|
| 260 | |
---|
| 261 | |
---|
| 262 | ****************************************************************************/ |
---|
| 263 | void identity_mat(const int n, double** A) { |
---|
| 264 | int i, j; |
---|
| 265 | |
---|
| 266 | for (i = 1; i <= n; i++) |
---|
| 267 | for (j = 1; j <= n; j++) |
---|
| 268 | A[i][j] = (i == j) ? 1.0 : 0.0; |
---|
| 269 | } |
---|
| 270 | |
---|
| 271 | /****************************************************************************/ |
---|
| 272 | /* void det_mat(const int n, double **A) |
---|
| 273 | |
---|
| 274 | Purpose: |
---|
| 275 | get the determinant of n x n matrix A |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | ****************************************************************************/ |
---|
| 279 | double det_mat(const int n, double **A) { |
---|
| 280 | int i, *indx; |
---|
| 281 | double **U, d; |
---|
| 282 | |
---|
| 283 | indx = ivector(1, n); |
---|
| 284 | U = dmatrix(1, n, 1, n); |
---|
| 285 | |
---|
| 286 | dmcopy(A, n, n, U); |
---|
| 287 | dludcmp(U, n, indx, &d); |
---|
| 288 | |
---|
| 289 | for (i = 1; i <= n; i++) |
---|
| 290 | d *= U[i][i]; |
---|
| 291 | |
---|
| 292 | free_dmatrix(U, 1, n, 1, n); |
---|
| 293 | free_ivector(indx, 1, n); |
---|
| 294 | |
---|
| 295 | return d; |
---|
| 296 | } |
---|
| 297 | |
---|
| 298 | /****************************************************************************/ |
---|
| 299 | /* void trace_mat(const int n, double **A) |
---|
| 300 | |
---|
| 301 | Purpose: |
---|
| 302 | get the trace of n x n matrix A |
---|
| 303 | |
---|
| 304 | |
---|
| 305 | ****************************************************************************/ |
---|
| 306 | double trace_mat(const int n, double **A) { |
---|
| 307 | int i; |
---|
| 308 | double d; |
---|
| 309 | |
---|
| 310 | d = 0.0; |
---|
| 311 | for (i = 1; i <= n; i++) |
---|
| 312 | d += A[i][i]; |
---|
| 313 | |
---|
| 314 | return d; |
---|
| 315 | } |
---|
| 316 | |
---|
| 317 | float K_fun(float lambda) { |
---|
| 318 | double **Id, **Lambda, **Lambda_inv, **U, **V, **K_int; |
---|
| 319 | |
---|
| 320 | Id = dmatrix(1, DOF, 1, DOF); |
---|
| 321 | Lambda = dmatrix(1, DOF, 1, DOF); |
---|
| 322 | Lambda_inv = dmatrix(1, DOF, 1, DOF); |
---|
| 323 | U = dmatrix(1, DOF, 1, DOF); |
---|
| 324 | V = dmatrix(1, DOF, 1, DOF); |
---|
| 325 | K_int = dmatrix(1, DOF, 1, DOF); |
---|
| 326 | |
---|
| 327 | identity_mat(DOF, Id); |
---|
| 328 | |
---|
| 329 | dmsmy(Id, DOF, DOF, lambda, U); |
---|
| 330 | // dmsub(C_, DOF, DOF, U, Lambda); |
---|
| 331 | dmadd(C_, DOF, DOF, U, Lambda); |
---|
| 332 | dinverse(Lambda, DOF, Lambda_inv); |
---|
| 333 | |
---|
| 334 | dmsmy(Id, DOF, DOF, trace_mat(DOF, Lambda_inv), U); |
---|
| 335 | dmsmy(Lambda_inv, DOF, DOF, 3.0, V); |
---|
| 336 | dmsub(U, DOF, DOF, V, K_int); |
---|
| 337 | dmsmy(K_int, DOF, DOF, 2.0 * sqr(M_PI) |
---|
| 338 | * sqrt(lambda / det_mat(DOF, Lambda)), K_int); |
---|
| 339 | |
---|
| 340 | free_dmatrix(Id, 1, DOF, 1, DOF); |
---|
| 341 | free_dmatrix(Lambda, 1, DOF, 1, DOF); |
---|
| 342 | free_dmatrix(Lambda_inv, 1, DOF, 1, DOF); |
---|
| 343 | free_dmatrix(U, 1, DOF, 1, DOF); |
---|
| 344 | free_dmatrix(V, 1, DOF, 1, DOF); |
---|
| 345 | free_dmatrix(K_int, 1, DOF, 1, DOF); |
---|
| 346 | |
---|
| 347 | return K_int[i_][j_]; |
---|
| 348 | } |
---|
| 349 | |
---|
| 350 | // partial template-class specialization |
---|
| 351 | // primary version |
---|
| 352 | template<typename T> |
---|
| 353 | class is_tps { |
---|
| 354 | }; |
---|
| 355 | |
---|
| 356 | // partial specialization |
---|
| 357 | template<> |
---|
| 358 | class is_tps<double> { |
---|
| 359 | public: |
---|
| 360 | static inline void get_ps(const ss_vect<double> &x, CellType &Cell) { |
---|
| 361 | Cell.BeamPos = x; |
---|
| 362 | } |
---|
| 363 | |
---|
| 364 | static inline double set_prm(const int k) { |
---|
| 365 | return 1.0; |
---|
| 366 | } |
---|
| 367 | |
---|
| 368 | static inline double get_curly_H(const ss_vect<tps> &x) { |
---|
| 369 | cout << "get_curly_H: operation not defined for double" << endl; |
---|
| 370 | exit_(1); |
---|
| 371 | return 0.0; |
---|
| 372 | } |
---|
| 373 | |
---|
| 374 | static inline double get_dI4(const double h, const double b2, |
---|
| 375 | const double L, const ss_vect<tps> &x) { |
---|
| 376 | cout << "get_dI4: operation not defined for double" << endl; |
---|
| 377 | exit_(1); |
---|
| 378 | return 0.0; |
---|
| 379 | } |
---|
| 380 | |
---|
| 381 | static inline void emittance(const double B2, const double u, |
---|
| 382 | const double ps0, const ss_vect<double> &xp) { |
---|
| 383 | } |
---|
| 384 | |
---|
| 385 | static inline void do_IBS(const double L, const ss_vect<double> &A_tps) { |
---|
| 386 | } |
---|
| 387 | |
---|
| 388 | static inline void diff_mat(const double B2, const double u, |
---|
| 389 | const double ps0, const ss_vect<double> &xp) { |
---|
| 390 | } |
---|
| 391 | |
---|
| 392 | }; |
---|
| 393 | |
---|
| 394 | // partial specialization |
---|
| 395 | template<> |
---|
| 396 | class is_tps<tps> { |
---|
| 397 | public: |
---|
| 398 | static inline void get_ps(const ss_vect<tps> &x, CellType &Cell) { |
---|
| 399 | getlinmat(6, x, Cell.A); |
---|
| 400 | } |
---|
| 401 | |
---|
| 402 | static inline tps set_prm(const int k) { |
---|
| 403 | return tps(0.0, k); |
---|
| 404 | } |
---|
| 405 | |
---|
| 406 | static inline double get_curly_H(const ss_vect<tps> &A) { |
---|
| 407 | int j; |
---|
| 408 | double curly_H[2]; |
---|
| 409 | ss_vect<double> eta; |
---|
| 410 | |
---|
| 411 | eta.zero(); |
---|
| 412 | for (j = 0; j < 4; j++) |
---|
| 413 | eta[j] = A[j][delta_]; |
---|
| 414 | |
---|
| 415 | get_twoJ(2, eta, A, curly_H); |
---|
| 416 | |
---|
| 417 | return curly_H[X_]; |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | static inline double get_dI4(const ss_vect<tps> &A) { |
---|
| 421 | return A[x_][delta_]; |
---|
| 422 | } |
---|
| 423 | |
---|
| 424 | static inline void emittance(const tps &B2, const tps &ds, const tps &ps0, |
---|
| 425 | const ss_vect<tps> &A) { |
---|
| 426 | int j; |
---|
| 427 | double B_66; |
---|
| 428 | ss_vect<tps> A_inv; |
---|
| 429 | |
---|
| 430 | if (B2 > 0.0) { |
---|
| 431 | B_66 = (q_fluct * pow(B2.cst(), 1.5) * pow(ps0, 4) * ds).cst(); |
---|
| 432 | A_inv = Inv(A); |
---|
| 433 | // D_11 = D_22 = curly_H_x,y * B_66 / 2, |
---|
| 434 | // curly_H_x,y = eta_Fl^2 + etap_Fl^2 |
---|
| 435 | for (j = 0; j < 3; j++) |
---|
| 436 | globval.D_rad[j] += (sqr(A_inv[j * 2][delta_]) + sqr(A_inv[j |
---|
| 437 | * 2 + 1][delta_])) * B_66 / 2.0; |
---|
| 438 | } |
---|
| 439 | } |
---|
| 440 | |
---|
| 441 | static void do_IBS(const double L, const ss_vect<tps> &A_tps) { |
---|
| 442 | /* A is passed, compute the invariants and emittances, |
---|
| 443 | The invariants for the uncoupled case are: |
---|
| 444 | |
---|
| 445 | [gamma alpha] |
---|
| 446 | Sigma ^-1 = [ ] |
---|
| 447 | x,y,z [alpha beta ] |
---|
| 448 | |
---|
| 449 | Note, ps = [x, y, ct, p_x, p_y, delta] */ |
---|
| 450 | |
---|
| 451 | int i, j, k; |
---|
| 452 | double **A, **Ainv, **Ainv_tp, **Ainv_tp_Ainv, **Boost, **Boost_tp; |
---|
| 453 | double **G[DOF], **M_a, **U, **C_a[DOF], **K, dln_eps[DOF]; |
---|
| 454 | double beta_x, beta_y, sigma_y, a_cst, two_Lc; |
---|
| 455 | |
---|
| 456 | const int n = 2 * DOF; |
---|
| 457 | const int indx[] = { 1, 4, 2, 5, 6, 3 }; |
---|
| 458 | |
---|
| 459 | const double P0 = 1e9 * globval.Energy * q_e / c0; |
---|
| 460 | const double gamma = 1e9 * globval.Energy / m_e; |
---|
| 461 | const double beta = sqrt(1.0 - 1.0 / sqr(gamma)); |
---|
| 462 | |
---|
| 463 | const double N_e = globval.Qb / q_e; |
---|
| 464 | |
---|
| 465 | A = dmatrix(1, n, 1, n); |
---|
| 466 | Ainv = dmatrix(1, n, 1, n); |
---|
| 467 | Ainv_tp = dmatrix(1, n, 1, n); |
---|
| 468 | Ainv_tp_Ainv = dmatrix(1, n, 1, n); |
---|
| 469 | Boost = dmatrix(1, n, 1, n); |
---|
| 470 | Boost_tp = dmatrix(1, n, 1, n); |
---|
| 471 | U = dmatrix(1, n, 1, n); |
---|
| 472 | M_a = dmatrix(1, n, 1, n); |
---|
| 473 | for (i = 0; i < DOF; i++) |
---|
| 474 | G[i] = dmatrix(1, n, 1, n); |
---|
| 475 | for (i = 0; i < DOF; i++) |
---|
| 476 | C_a[i] = dmatrix(1, DOF, 1, DOF); |
---|
| 477 | C_ = dmatrix(1, DOF, 1, DOF); |
---|
| 478 | K = dmatrix(1, DOF, 1, DOF); |
---|
| 479 | |
---|
| 480 | // Compute invariants (in Floquet space): Sigma^-1 = (A^-1)^tp * A^-1 |
---|
| 481 | |
---|
| 482 | for (i = 1; i <= n; i++) |
---|
| 483 | for (j = 1; j <= n; j++) |
---|
| 484 | A[i][j] = A_tps[i - 1][j - 1]; |
---|
| 485 | |
---|
| 486 | dmcopy(A, n, n, U); |
---|
| 487 | dinverse(U, n, Ainv); |
---|
| 488 | dmtranspose(Ainv, n, n, Ainv_tp); |
---|
| 489 | dmmult(Ainv_tp, n, n, Ainv, n, n, Ainv_tp_Ainv); |
---|
| 490 | |
---|
| 491 | for (i = 0; i < DOF; i++) |
---|
| 492 | for (j = 1; j <= n; j++) |
---|
| 493 | for (k = 1; k <= n; k++) |
---|
| 494 | G[i][indx[j - 1]][indx[k - 1]] = Ainv[2 * i + 1][j] |
---|
| 495 | * Ainv[2 * i + 1][k] + Ainv[2 * i + 2][j] * Ainv[2 |
---|
| 496 | * i + 2][k]; |
---|
| 497 | |
---|
| 498 | dmadd(G[0], n, n, G[1], U); |
---|
| 499 | dmadd(U, n, n, G[2], U); |
---|
| 500 | |
---|
| 501 | if (trace) { |
---|
| 502 | printf("\n"); |
---|
| 503 | dmdump(stdout, "G_1:", G[0], n, n, "%11.3e"); |
---|
| 504 | dmdump(stdout, "G_2:", G[1], n, n, "%11.3e"); |
---|
| 505 | dmdump(stdout, "G_3:", G[2], n, n, "%11.3e"); |
---|
| 506 | dmdump(stdout, "Ainv_tp*Ainv:", Ainv_tp_Ainv, n, n, "%11.3e"); |
---|
| 507 | dmdump(stdout, "Sum_a G_a", U, n, n, "%11.3e"); |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | /* Transform from the co-moving to COM frame: |
---|
| 511 | |
---|
| 512 | [ 1 0 0 0 0 0 ] |
---|
| 513 | [ 0 1 0 0 0 0 ] |
---|
| 514 | [ 0 0 1/gamma 0 0 0 ] |
---|
| 515 | [ 0 0 0 1/P0 0 0 ] |
---|
| 516 | [ 0 0 0 0 1/P0 0 ] |
---|
| 517 | [ 0 0 0 0 0 gamma/P0 ] */ |
---|
| 518 | |
---|
| 519 | identity_mat(n, Boost); |
---|
| 520 | Boost[3][3] /= gamma; |
---|
| 521 | Boost[4][4] /= P0; |
---|
| 522 | Boost[5][5] /= P0; |
---|
| 523 | Boost[6][6] *= gamma / P0; |
---|
| 524 | |
---|
| 525 | dmtranspose(Boost, n, n, Boost_tp); |
---|
| 526 | |
---|
| 527 | zero_mat(DOF, C_); |
---|
| 528 | for (i = 0; i < DOF; i++) { |
---|
| 529 | dmmult(Boost_tp, n, n, G[i], n, n, U); |
---|
| 530 | dmmult(U, n, n, Boost, n, n, M_a); |
---|
| 531 | |
---|
| 532 | if (trace) |
---|
| 533 | dmdump(stdout, "M_a:", M_a, n, n, "%11.3e"); |
---|
| 534 | |
---|
| 535 | // Extract the C_a matrices from the momentum components of M_a |
---|
| 536 | for (j = 1; j <= DOF; j++) |
---|
| 537 | for (k = 1; k <= DOF; k++) |
---|
| 538 | C_a[i][j][k] = M_a[DOF + j][DOF + k]; |
---|
| 539 | |
---|
| 540 | dmsmy(C_a[i], DOF, DOF, sqr(P0), C_a[i]); |
---|
| 541 | |
---|
| 542 | if (globval.eps[i] != 0.0) |
---|
| 543 | dmsmy(C_a[i], DOF, DOF, 1.0 / globval.eps[i], U); |
---|
| 544 | else { |
---|
| 545 | cout << "*** do_IBS: zero emittance for plane " << i + 1 |
---|
| 546 | << endl; |
---|
| 547 | exit(1); |
---|
| 548 | } |
---|
| 549 | |
---|
| 550 | dmadd(C_, DOF, DOF, U, C_); |
---|
| 551 | } |
---|
| 552 | |
---|
| 553 | if (trace) { |
---|
| 554 | dmdump(stdout, "C_1:", C_a[0], DOF, DOF, "%11.3e"); |
---|
| 555 | dmdump(stdout, "C_2:", C_a[1], DOF, DOF, "%11.3e"); |
---|
| 556 | dmdump(stdout, "C_3:", C_a[2], DOF, DOF, "%11.3e"); |
---|
| 557 | dmdump(stdout, "C:", C_, DOF, DOF, "%11.3e"); |
---|
| 558 | } |
---|
| 559 | |
---|
| 560 | for (i = 1; i <= DOF; i++) |
---|
| 561 | for (j = 1; j <= DOF; j++) { |
---|
| 562 | // upper bound is infinity |
---|
| 563 | i_ = i; |
---|
| 564 | j_ = j; |
---|
| 565 | K[i][j] = qromb(K_fun, 0.0, 1e8); |
---|
| 566 | } |
---|
| 567 | |
---|
| 568 | if (trace) |
---|
| 569 | dmdump(stdout, "K:", K, DOF, DOF, "%11.3e"); |
---|
| 570 | |
---|
| 571 | // Compute the Coulomb logarithm |
---|
| 572 | beta_x = C_a[0][1][1]; |
---|
| 573 | beta_y = C_a[1][2][2]; |
---|
| 574 | sigma_y = sqrt(beta_y * globval.eps[Y_]); |
---|
| 575 | two_Lc = log(sqr(gamma * globval.eps[X_] * sigma_y / (r_e * beta_x))); |
---|
| 576 | |
---|
| 577 | if (trace) |
---|
| 578 | printf("2(log) = %11.3e\n", two_Lc); |
---|
| 579 | |
---|
| 580 | // include time dilatation and scaling of C matrix |
---|
| 581 | a_cst = L * two_Lc * N_e * sqr(r_e) * c0 / (64.0 * cube(M_PI * beta) |
---|
| 582 | * pow(gamma, 4) * globval.eps[X_] * globval.eps[Y_] |
---|
| 583 | * globval.eps[Z_]); |
---|
| 584 | |
---|
| 585 | // dSigma_ab/dt = a_cst * K_ab, e.g. d<delta^2>/dt = a_cst K_33. |
---|
| 586 | // deps is obtained from C. |
---|
| 587 | |
---|
| 588 | if (trace) |
---|
| 589 | printf("dln_eps:"); |
---|
| 590 | for (i = 0; i < DOF; i++) { |
---|
| 591 | dmmult(C_a[i], DOF, DOF, K, DOF, DOF, U); |
---|
| 592 | // Dt = L/c0 |
---|
| 593 | // dln_eps[i] = a_cst/globval.eps[i]*trace_mat(DOF, U); |
---|
| 594 | globval.D_IBS[i] += a_cst * trace_mat(DOF, U); |
---|
| 595 | if (trace) |
---|
| 596 | printf("%11.3e", dln_eps[i]); |
---|
| 597 | } |
---|
| 598 | if (trace) |
---|
| 599 | printf("\n"); |
---|
| 600 | |
---|
| 601 | free_dmatrix(A, 1, n, 1, n); |
---|
| 602 | free_dmatrix(Ainv, 1, n, 1, n); |
---|
| 603 | free_dmatrix(Ainv_tp, 1, n, 1, n); |
---|
| 604 | free_dmatrix(Ainv_tp_Ainv, 1, n, 1, n); |
---|
| 605 | free_dmatrix(Boost, 1, n, 1, n); |
---|
| 606 | free_dmatrix(Boost_tp, 1, n, 1, n); |
---|
| 607 | free_dmatrix(U, 1, n, 1, n); |
---|
| 608 | free_dmatrix(M_a, 1, n, 1, n); |
---|
| 609 | for (i = 0; i < DOF; i++) |
---|
| 610 | free_dmatrix(G[i], 1, n, 1, n); |
---|
| 611 | for (i = 0; i < DOF; i++) |
---|
| 612 | free_dmatrix(C_a[i], 1, DOF, 1, DOF); |
---|
| 613 | free_dmatrix(C_, 1, DOF, 1, DOF); |
---|
| 614 | free_dmatrix(K, 1, DOF, 1, DOF); |
---|
| 615 | } |
---|
| 616 | |
---|
| 617 | static inline void diff_mat(const tps &B2, const tps &ds, const tps &ps0, |
---|
| 618 | ss_vect<tps> &x) { |
---|
| 619 | } |
---|
| 620 | |
---|
| 621 | }; |
---|
| 622 | |
---|
| 623 | template<typename T> |
---|
| 624 | void get_B2(const double h_ref, const T B[], const ss_vect<T> &xp, T &B2_perp, |
---|
| 625 | T &B2_par) { |
---|
| 626 | // compute |B|^2_perpendicular and |B|^2_parallel |
---|
| 627 | T xn, e[3]; |
---|
| 628 | |
---|
| 629 | xn = 1.0 / sqrt(sqr(1.0 + xp[x_] * h_ref) + sqr(xp[px_]) + sqr(xp[py_])); |
---|
| 630 | e[X_] = xp[px_] * xn; |
---|
| 631 | e[Y_] = xp[py_] * xn; |
---|
| 632 | e[Z_] = (1e0 + xp[x_] * h_ref) * xn; |
---|
| 633 | |
---|
| 634 | // left-handed coordinate system |
---|
| 635 | B2_perp = sqr(B[Y_] * e[Z_] - B[Z_] * e[Y_]) + sqr(B[X_] * e[Y_] - B[Y_] |
---|
| 636 | * e[X_]) + sqr(B[Z_] * e[X_] - B[X_] * e[Z_]); |
---|
| 637 | |
---|
| 638 | // B2_par = sqr(B[X_]*e[X_]+B[Y_]*e[Y_]+B[Z_]*e[Z_]); |
---|
| 639 | } |
---|
| 640 | |
---|
| 641 | template<typename T> |
---|
| 642 | void radiate(ss_vect<T> &x, const double L, const double h_ref, const T B[]) { |
---|
| 643 | T ps0, ps1, ds, B2_perp = 0.0, B2_par = 0.0; |
---|
| 644 | ss_vect<T> xp; |
---|
| 645 | |
---|
| 646 | // large ring: conservation of x' and y' |
---|
| 647 | xp = x; |
---|
| 648 | ps0 = get_p_s(x); |
---|
| 649 | xp[px_] /= ps0; |
---|
| 650 | xp[py_] /= ps0; |
---|
| 651 | |
---|
| 652 | // H = -p_s => ds = H*L |
---|
| 653 | ds = (1.0 + xp[x_] * h_ref + (sqr(xp[px_]) + sqr(xp[py_])) / 2.0) * L; |
---|
| 654 | get_B2(h_ref, B, xp, B2_perp, B2_par); |
---|
| 655 | |
---|
| 656 | if (globval.radiation) { |
---|
| 657 | x[delta_] -= cl_rad * sqr(ps0) * B2_perp * ds; |
---|
| 658 | ps1 = get_p_s(x); |
---|
| 659 | x[px_] = xp[px_] * ps1; |
---|
| 660 | x[py_] = xp[py_] * ps1; |
---|
| 661 | } |
---|
| 662 | |
---|
| 663 | if (globval.emittance) |
---|
| 664 | is_tps<T>::emittance(B2_perp, ds, ps0, xp); |
---|
| 665 | } |
---|
| 666 | |
---|
| 667 | /******************************************************************** |
---|
| 668 | static double get_psi(double irho, double phi, double gap) { |
---|
| 669 | |
---|
| 670 | Purpose: |
---|
| 671 | Correction for magnet gap (longitudinal fringe field) |
---|
| 672 | |
---|
[11] | 673 | Input: |
---|
| 674 | irho: h = 1/rho [1/m] |
---|
| 675 | phi: dipole edge (entrance or exit) angle |
---|
| 676 | gap: full gap between poles |
---|
[3] | 677 | |
---|
| 678 | 2 |
---|
| 679 | K1*gap*h*(1 + sin phi) |
---|
| 680 | psi = ----------------------- * (1 - K2*h*gap*tan phi) |
---|
| 681 | cos phi |
---|
| 682 | |
---|
| 683 | K1 is usually 1/2 |
---|
| 684 | K2 is zero here |
---|
| 685 | For the type of Tanh-like fringe field |
---|
| 686 | |
---|
| 687 | |
---|
| 688 | 18/06/2012 Jianfeng Zhang@LAL |
---|
| 689 | |
---|
| 690 | Fix the bug to calculate psi for the sector magnets which has phi=0 |
---|
| 691 | |
---|
| 692 | *********************************************************************/ |
---|
| 693 | static double get_psi(double irho, double phi, double gap) { |
---|
| 694 | |
---|
| 695 | double psi; |
---|
| 696 | |
---|
| 697 | const double k1 = 0.5, k2 = 0.0; |
---|
| 698 | |
---|
| 699 | // if (phi == 0.0) //NOT hard edge, but sector magnets!!! |
---|
| 700 | // psi = 0.0; |
---|
| 701 | // else |
---|
| 702 | psi = k1 * gap * irho * (1.0 + sqr(sin(dtor(phi)))) / cos(dtor(phi)) |
---|
| 703 | * (1.0 - k2 * gap * irho * tan(dtor(phi))); |
---|
| 704 | |
---|
| 705 | return psi; |
---|
| 706 | } |
---|
| 707 | |
---|
[11] | 708 | /***************************************************************** |
---|
| 709 | * Exact map of the sector dipole |
---|
| 710 | * |
---|
| 711 | * Forest's beam dynamics P.360 Eqn. (12.18) |
---|
| 712 | * Also see Laurent's thesis: P219-220. |
---|
| 713 | * |
---|
| 714 | * Written by Jianfeng Zhang @ LAL, 01/08/2013 |
---|
| 715 | * |
---|
| 716 | * Test version........................................ |
---|
| 717 | * |
---|
| 718 | * ****************************************************************/ |
---|
| 719 | template<typename T> |
---|
| 720 | void dipole_kick(double L, double h_ref, double h_bend, ss_vect<T> &x){ |
---|
| 721 | |
---|
| 722 | ss_vect<T> x0; |
---|
| 723 | T ps, u, dpxf; |
---|
| 724 | |
---|
| 725 | x0 = x; |
---|
| 726 | ps = get_p_s(x0); |
---|
| 727 | |
---|
| 728 | u = sqrt(sqr(1+x0[delta_])-sqr(x0[py_])); |
---|
| 729 | |
---|
| 730 | x[py_] = x0[py_]; |
---|
| 731 | |
---|
| 732 | x[delta_] =x0[delta_]; |
---|
| 733 | |
---|
| 734 | x[px_]= x0[px_]*cos(L*h_ref) + (ps -h_bend*(1/h_ref+x0[x_]))*sin(L*h_ref); |
---|
| 735 | |
---|
| 736 | dpxf = h_ref*sqrt((1+x[delta_])*(1+x[delta_]) - x[px_]*x[px_] - x[py_]*x[py_])-h_bend*(1+h_ref*x[x_]); |
---|
| 737 | |
---|
| 738 | |
---|
| 739 | x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u)); |
---|
| 740 | |
---|
| 741 | x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x[delta_])-sqr(x[px_])-sqr(x[py_]))-dpxf - h_bend); |
---|
| 742 | |
---|
| 743 | x[ct_] += (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u)); |
---|
| 744 | |
---|
| 745 | } |
---|
| 746 | |
---|
| 747 | /***************************************************************** |
---|
| 748 | * pure multipole kick |
---|
| 749 | * |
---|
| 750 | * Forest's beam dynamics P.360 Eqn. (12.18) |
---|
| 751 | * Also see Laurent's thesis: P219-220. |
---|
| 752 | * |
---|
| 753 | * Written by Jianfeng Zhang @ LAL, 01/08/2013 |
---|
| 754 | * |
---|
| 755 | * Test version........................................ |
---|
| 756 | * |
---|
| 757 | * ****************************************************************/ |
---|
| 758 | template<typename T> |
---|
| 759 | void multipole_kick(int Order, double MB[], double L, double h_bend, ss_vect<T> &x){ |
---|
| 760 | |
---|
| 761 | int j=0; |
---|
| 762 | ss_vect<T> x0; |
---|
| 763 | |
---|
| 764 | T ByoBrho = 0, BxoBrho = 0, ByoBrho1=0 ; |
---|
| 765 | |
---|
| 766 | if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) { |
---|
| 767 | x0 = x; |
---|
| 768 | /* compute the total magnetic field Bx and By with Horner's rule */ |
---|
| 769 | ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e) |
---|
| 770 | BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e) |
---|
| 771 | for (j = Order - 1; j >= 1; j--) { |
---|
| 772 | ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax]; |
---|
| 773 | BxoBrho = x0[y_] * ByoBrho + x0[x_] * BxoBrho + MB[HOMmax - j]; |
---|
| 774 | ByoBrho = ByoBrho1; |
---|
| 775 | } |
---|
| 776 | } |
---|
| 777 | |
---|
| 778 | x[px_] -= L * (ByoBrho + h_bend*1); |
---|
| 779 | //vertical kick due to the magnets |
---|
| 780 | x[py_] += L * BxoBrho; |
---|
| 781 | } |
---|
| 782 | |
---|
[3] | 783 | /****************************************************************************/ |
---|
| 784 | /* template<typename T> |
---|
| 785 | void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref,ss_vect<T> &x) |
---|
| 786 | |
---|
| 787 | Purpose: |
---|
| 788 | Refer to Tracy 2 manual. |
---|
| 789 | Calculate multipole kick. |
---|
[11] | 790 | , |
---|
[3] | 791 | (1) The exact Hamiltonian is |
---|
| 792 | |
---|
| 793 | H = A + B where A and B are the kick part defined by |
---|
| 794 | |
---|
| 795 | (exact form of the H of the drift) |
---|
[9] | 796 | A(x,y,-l,px,py,dP) = -sqrt( (1 + delta)^2 - px^2 - py^2 ) + delta |
---|
[3] | 797 | |
---|
| 798 | (exact form of the H of the kick from the magnet) |
---|
[9] | 799 | B(x,y,-l,px,py,dP) = -h*x*sqrt( (1 + delta)^2 - px^2 - py^2 ) + h x + h^2*x^2/2 + int(Re(By+iBx)/Brho) |
---|
[3] | 800 | |
---|
| 801 | so this is the exact Hamitonian for small ring. |
---|
| 802 | |
---|
| 803 | |
---|
| 804 | The kick is given by |
---|
| 805 | e |
---|
[9] | 806 | Dp_x = L* (h_ref*sqrt( (1 + delta)^2 - (px^2 + py^2) ) - h_bend - x*h_bend*h_ref - ---- B_y ) |
---|
[3] | 807 | p_0 |
---|
| 808 | Dp_y = L* e/p_0 * B_x |
---|
| 809 | |
---|
[9] | 810 | x = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * p_x |
---|
[3] | 811 | |
---|
[9] | 812 | y = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * p_y |
---|
[3] | 813 | |
---|
[9] | 814 | CT = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * (1+delta) |
---|
[3] | 815 | |
---|
| 816 | where |
---|
| 817 | e 1 |
---|
| 818 | --- = ----- |
---|
| 819 | p_0 B rho |
---|
| 820 | ==== |
---|
| 821 | \ |
---|
| 822 | (B_y + iB_x) = B rho > (ia_n + b_n ) (x + iy)^n-1 |
---|
| 823 | / |
---|
| 824 | ==== |
---|
| 825 | |
---|
| 826 | (2) The approximate Hamiltonian is |
---|
| 827 | |
---|
| 828 | H = A + B where A and B are defined by |
---|
| 829 | 2 2 |
---|
| 830 | px + py |
---|
| 831 | A(x,y,-l,px,py,dP) = --------- (approximate form of the H of the drift) |
---|
| 832 | 2(1+delta) |
---|
| 833 | 2 2 |
---|
| 834 | B(x,y,-l,px,py,dP) = -h*x*delta + 0.5 h x + int(Re(By+iBx)/Brho) (approximate form of the H of the kick from the magnet) |
---|
| 835 | |
---|
| 836 | so this is the appproximation for large ring |
---|
| 837 | the chromatic term is missing hx*A |
---|
| 838 | |
---|
| 839 | |
---|
| 840 | The kick is given by |
---|
| 841 | |
---|
| 842 | e L L delta L x e L |
---|
| 843 | Dp_x = --- B_y + ------- - ----- , Dp_y = --- B_x |
---|
| 844 | p_0 rho rho^2 p_0 |
---|
| 845 | |
---|
| 846 | //second order of the Hamiltonian motion |
---|
| 847 | x = h_ref * x_0 / (1+delta) * p_x |
---|
| 848 | |
---|
| 849 | y = h_ref * x_0 / (1+delta) * p_y |
---|
| 850 | |
---|
| 851 | CT = h_ref * x_0 / (1+delta) * (p_x^2 + p_y^2)/2/(1+delta) |
---|
| 852 | |
---|
| 853 | where |
---|
| 854 | e 1 |
---|
| 855 | --- = ----- |
---|
| 856 | p_0 B rho |
---|
| 857 | ==== |
---|
| 858 | \ |
---|
| 859 | (B_y + iB_x) = B rho > (ia_n + b_n ) (x + iy)^n-1 |
---|
| 860 | / |
---|
| 861 | ==== |
---|
| 862 | |
---|
| 863 | = B rho [(b_n*x - a_n*y) + i*(bn*y + an*x) + b_(n-1) + i*a_(n-1)] * (x+iy)^n-2 |
---|
| 864 | =...... |
---|
| 865 | |
---|
| 866 | Input: |
---|
| 867 | Order maximum non zero multipole component |
---|
| 868 | MB array of an and bn, magnetic field components |
---|
| 869 | L multipole length |
---|
[11] | 870 | h_bend B_dipole/Brho is the dipole field normalized by the reference momentum; one of the example is the combined dipole |
---|
[3] | 871 | h_ref 1/rho [m^-1] is the curvature of the reference orbit in the dipoles which is used in curvilinear coordinate, |
---|
| 872 | x initial coordinates vector |
---|
| 873 | |
---|
| 874 | Output: |
---|
| 875 | x final coordinates vector |
---|
| 876 | |
---|
| 877 | Return: |
---|
| 878 | none |
---|
| 879 | |
---|
| 880 | Global variables: |
---|
| 881 | none |
---|
| 882 | |
---|
| 883 | Specific functions: |
---|
| 884 | none |
---|
| 885 | |
---|
| 886 | Comments: |
---|
| 887 | 06/11/2012 Jianfeng Nadolski @ LAL |
---|
| 888 | Fix the bug in the kick map of the exact Hamiltonian. |
---|
| 889 | Add the second order of the approximate Hamiltonian. |
---|
| 890 | |
---|
| 891 | **************************************************************************/ |
---|
| 892 | template<typename T> |
---|
| 893 | void thin_kick(int Order, double MB[], double L, double h_bend, double h_ref, |
---|
| 894 | ss_vect<T> &x) { |
---|
| 895 | int j; |
---|
| 896 | T BxoBrho, ByoBrho, ByoBrho1, B[3]; |
---|
| 897 | T sqrtpx, dpx, ps2new, psnew; |
---|
| 898 | T n1, n2; |
---|
| 899 | ss_vect<T> x0, cod; |
---|
| 900 | T u=0.0; |
---|
| 901 | |
---|
| 902 | if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) { |
---|
| 903 | x0 = x; |
---|
| 904 | /* compute the total magnetic field Bx and By with Horner's rule */ |
---|
| 905 | ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e) |
---|
| 906 | BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e) |
---|
| 907 | for (j = Order - 1; j >= 1; j--) { |
---|
| 908 | ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax]; |
---|
| 909 | BxoBrho = x0[y_] * ByoBrho + x0[x_] * BxoBrho + MB[HOMmax - j]; |
---|
| 910 | ByoBrho = ByoBrho1; |
---|
| 911 | } |
---|
| 912 | |
---|
| 913 | |
---|
| 914 | if (globval.radiation || globval.emittance) { |
---|
| 915 | B[X_] = BxoBrho; |
---|
| 916 | B[Y_] = ByoBrho + h_bend; |
---|
| 917 | B[Z_] = 0.0; |
---|
| 918 | radiate(x, L, h_ref, B); |
---|
| 919 | } |
---|
| 920 | |
---|
| 921 | //curvilinear cocodinates, from the dipole components |
---|
| 922 | if (h_ref != 0.0) { |
---|
| 923 | |
---|
| 924 | //exact Hamiltonian |
---|
| 925 | if (globval.H_exact) { |
---|
| 926 | |
---|
| 927 | if(get_p_s(x0)==0) |
---|
| 928 | return; |
---|
| 929 | |
---|
[11] | 930 | //dipole kick map from the exact Hamiltonian; not symplectic |
---|
| 931 | // u = L * h_ref * x0[x_] /get_p_s(x0); |
---|
| 932 | // x[x_] += u * x0[px_]; |
---|
| 933 | // x[y_] += u * x0[py_]; |
---|
| 934 | // x[px_]-= L * (-h_ref * get_p_s(x0) + h_bend + h_ref * h_bend |
---|
| 935 | // * x0[x_] + ByoBrho); |
---|
| 936 | // x[ct_] += u * (1.0+x0[delta_]); |
---|
| 937 | |
---|
| 938 | //kick map due to the dipole vector potential |
---|
| 939 | // x[px_] -= L*(h_bend + h_ref*h_bend*x0[x_]+ByoBrho); |
---|
[3] | 940 | |
---|
| 941 | // sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_])); |
---|
| 942 | // // |
---|
| 943 | // x[px_] = x0[px_]*cos(L*h_ref) + (get_p_s(x0)-h_bend*(h_ref+x0[x_]))*sin(L*h_ref)-L*ByoBrho; |
---|
| 944 | // |
---|
| 945 | // dpx = -x0[px_]*h_ref*sin(L*h_ref)+(get_p_s(x0)-h_bend*(h_ref+x0[x_]))*h_ref*cos(L*h_ref); |
---|
| 946 | // |
---|
| 947 | // x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[py_]))-dpx - h_bend); |
---|
| 948 | // |
---|
| 949 | // x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(sqrtpx*x0[px_])-asin(sqrtpx*x[px_])); |
---|
| 950 | // |
---|
| 951 | // x[ct_]+= (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(sqrtpx*x0[px_])-asin(sqrtpx*x[px_])); |
---|
| 952 | |
---|
| 953 | //map of a wedge |
---|
| 954 | // x[px_] = x0[px_]*cos(L*h_bend) + (get_p_s(x0) - h_ref*x0[x_])*sin(L*h_bend) -L*ByoBrho; |
---|
| 955 | // |
---|
| 956 | // sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_])); |
---|
| 957 | // n1 = x0[px_]*sqrtpx; |
---|
| 958 | // n2 = x[px_]*sqrtpx; |
---|
| 959 | // |
---|
| 960 | // x[x_] = x0[x_]*cos(L*h_bend) + (x0[x_]*x0[px_]*sin(2*L*h_bend)+sqr(sin(L*h_bend))*(2*x0[x_]*get_p_s(x0)-h_ref*sqr(x0[x_]))) |
---|
| 961 | // /(sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[y_]))+get_p_s(x0)*cos(L*h_bend)-x0[px_]*sin(L*h_bend)); |
---|
| 962 | // |
---|
| 963 | // x[y_] += x0[py_]*L*h_bend*h_ref + x0[py_]*h_ref*(asin(n1) - asin(n2)); |
---|
| 964 | // x[ct_]+= (1+x0[delta_])*L*h_bend*h_ref + (1+x0[delta_])*h_ref*(asin(n1) - asin(n2)); |
---|
| 965 | // |
---|
[9] | 966 | }//approximate Hamiltonian |
---|
[10] | 967 | else {//first order map |
---|
[9] | 968 | x[px_] -= L * (ByoBrho + (h_bend - h_ref) / 2.0 + h_ref * h_bend |
---|
[3] | 969 | * x0[x_] - h_ref * x0[delta_]); |
---|
[9] | 970 | x[ct_] += L * h_ref * x0[x_]; |
---|
| 971 | |
---|
[11] | 972 | // second order, cubic term of H, H3; has problem when do tracking...works for THOMX!!!!! |
---|
| 973 | // not symplectic???? |
---|
| 974 | if(0){ |
---|
| 975 | u = L * h_ref * x0[x_] /(1.0+x0[delta_]); |
---|
[9] | 976 | x[x_] += u * x0[px_]; |
---|
| 977 | x[y_] += u * x0[py_]; |
---|
[19] | 978 | // x[px_]+= 0.5*u/x0[x_]*(x0[px_]*x0[px_]+x0[py_]*x0[py_]); |
---|
[9] | 979 | x[ct_] += u*(sqr(x0[px_])+sqr(x0[py_]))/(2.0*(1.0+x0[delta_])); |
---|
[10] | 980 | |
---|
[11] | 981 | // cout << "test......2..... u = " << u<< endl; |
---|
| 982 | } |
---|
| 983 | } |
---|
[3] | 984 | |
---|
| 985 | }//from other straight magnets |
---|
| 986 | else |
---|
[11] | 987 | x[px_] -= L * (ByoBrho + h_bend*1); //h_bend should be trigger on for combined dipoles |
---|
[3] | 988 | |
---|
| 989 | //vertical kick due to the magnets |
---|
| 990 | x[py_] += L * BxoBrho; |
---|
| 991 | |
---|
| 992 | /* |
---|
| 993 | if(h_ref == 0 && h_bend !=0){ |
---|
| 994 | sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_])); |
---|
| 995 | ps2new = sqr(1+x0[delta_]) - sqr(x[px_]) -sqr(x0[py_]); |
---|
| 996 | psnew = sqrt(ps2new); |
---|
| 997 | |
---|
| 998 | x[x_] += 1/h_bend*(psnew - get_p_s(x0)); |
---|
| 999 | x[y_] += x0[py_]/h_bend*(asin(x0[px_]*sqrtpx)-asin(x[px_]*sqrtpx)); |
---|
| 1000 | x[ct_]+= (1+x0[delta_])/h_bend*(asin(x0[px_]*sqrtpx)-asin(x[px_]*sqrtpx)); |
---|
| 1001 | }*/ |
---|
| 1002 | |
---|
| 1003 | |
---|
| 1004 | } |
---|
| 1005 | } |
---|
| 1006 | |
---|
| 1007 | /****************************************************************************/ |
---|
| 1008 | /* template<typename T> |
---|
| 1009 | static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x) |
---|
| 1010 | |
---|
| 1011 | Purpose: |
---|
| 1012 | Compute edge focusing for a dipole |
---|
| 1013 | There is no radiation coming from the edge |
---|
[11] | 1014 | The standard formula used is : (SLAC-75, K. Brown's first order of fringe field but with the correction of 1/(1+delta) in R43) |
---|
[3] | 1015 | |
---|
| 1016 | px = px0 + irho tan(phi) *x0 |
---|
| 1017 | |
---|
| 1018 | |
---|
| 1019 | irho |
---|
| 1020 | pz = pz0 - ------ tan(phi - psi) *z0 |
---|
| 1021 | 1 + dP |
---|
| 1022 | |
---|
| 1023 | for psi definition see its function |
---|
| 1024 | |
---|
| 1025 | Input: |
---|
| 1026 | irho = inverse of curvature radius (rho = 5.36 m for SOLEIL) |
---|
| 1027 | phi = entrance/exit angle of the dipole edge,usually half the |
---|
| 1028 | curvature angle of a dipole |
---|
| 1029 | gap = gap of the dipole for longitudinal fringing field (see psi) |
---|
| 1030 | x = input particle coordinates |
---|
| 1031 | Entrance: bool flag, true, then calculate edge focus and fringe field at the entrance of the dipole |
---|
| 1032 | false,...........................................at the exit fo the dipole |
---|
| 1033 | |
---|
| 1034 | Output: |
---|
| 1035 | x output particle coordinates |
---|
| 1036 | |
---|
| 1037 | Return: |
---|
| 1038 | none |
---|
| 1039 | |
---|
| 1040 | Global variables: |
---|
| 1041 | none |
---|
| 1042 | |
---|
| 1043 | Specific functions: |
---|
| 1044 | none |
---|
| 1045 | |
---|
| 1046 | Comments: |
---|
| 1047 | 05/07/10 Add energy dependence part irho replaced by irho/(1.0+x[delta_]) |
---|
| 1048 | now the chromaticity attribution of the dipole edge is used |
---|
| 1049 | by a simple 1.0+x[delta_], but not a complicated Hamiltonian |
---|
| 1050 | expansion. |
---|
| 1051 | Now chromaticities in Tracy II and Tracy III are the same. |
---|
| 1052 | Modification based on Tracy II soleil version. |
---|
| 1053 | |
---|
| 1054 | April 2011, No energy dependence for H-plane according to SSC-141 note (E. Forest) |
---|
[11] | 1055 | Agreement better with MAD8 (LNLS), SOLEIL small effect (.1). Modified by Laurent Nadolski @ soleil |
---|
| 1056 | Comments: Jianfeng Zhang @ LAL, 05/06/2013 |
---|
| 1057 | This model only works for Soleil lattice, but not works for ThomX and SuperB DR lattice. |
---|
[3] | 1058 | |
---|
| 1059 | 18/06/2012 Jianfeng Zhang @ LAL |
---|
[11] | 1060 | Add the contribution of x' to the x[py_] for small ring like Thom-X, using Alexandre Loulergue's correction. |
---|
| 1061 | This models works for ThomX, Soleil, and SuperB DR lattices; and |
---|
| 1062 | close to the Forest's model and K. Brown's models. |
---|
[3] | 1063 | |
---|
[24] | 1064 | 07/2013 Jianfeng Zhang @ LAL |
---|
| 1065 | Fix the bug in x[py], change |
---|
| 1066 | x0[y]/ (1.0 + x0[delta_]*1) to x0[y_]/ (1.0 + x0[delta_]*0), |
---|
| 1067 | now the vertical chromaticity is not 0.5 as in ELEGANT. |
---|
| 1068 | |
---|
[3] | 1069 | ****************************************************************************/ |
---|
| 1070 | template<typename T> |
---|
| 1071 | static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x, bool Entrance) { |
---|
| 1072 | |
---|
[24] | 1073 | ss_vect<T> x0; |
---|
| 1074 | x0 = x; |
---|
[3] | 1075 | |
---|
| 1076 | if (true) { |
---|
| 1077 | |
---|
[11] | 1078 | // cout << "Dipole edge effect WITH Alex's correction: " << endl; |
---|
[3] | 1079 | |
---|
| 1080 | //with the contribution from the entrance angle of the particle at the edge |
---|
| 1081 | // Original model, written by L. Nadolski, this is a reasonal model!!!! |
---|
| 1082 | // warning: => diverging Taylor map (see SSC-141) |
---|
| 1083 | // x[px_] += irho * tan(dtor(phi)) * x[x_]; |
---|
| 1084 | // x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_] |
---|
| 1085 | // / (1.0 + x[delta_]); |
---|
| 1086 | |
---|
| 1087 | |
---|
| 1088 | // K.Brown 2rd order transfer matrix of the dipole fringe field for sector dipoles (phi=0) |
---|
| 1089 | // compared the fringe field model used in Tracy |
---|
| 1090 | // results: K. Brown's model is wrong for small ring like Thom-X!!!! |
---|
[11] | 1091 | // The original model of tracy is correct for tracking!!!!!!!!!! |
---|
[3] | 1092 | /* |
---|
| 1093 | double psi; |
---|
| 1094 | psi=gap*0.5*irho*(1+sin(dtor(phi))*sin(dtor(phi)))/cos(dtor(phi)); |
---|
| 1095 | |
---|
| 1096 | if(Entrance){ |
---|
| 1097 | x[x_] += irho/2*x[y_]*x[y_]; |
---|
| 1098 | x[px_] += 0.5*x[x_]*x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5*x[y_]*x[y_]; |
---|
| 1099 | x[py_] -= -irho*tan(-dtor(psi))*x[y_]+1*x[x_]*x[y_]+irho*x[px_]*x[y_]+ (irho*psi/cos(-dtor(psi))/cos(-dtor(psi))) * x[y_] |
---|
| 1100 | / (1.0 + x[delta_]); |
---|
| 1101 | } |
---|
| 1102 | else |
---|
| 1103 | { |
---|
| 1104 | x[x_] -= irho/2*x[y_]*x[y_]; |
---|
| 1105 | x[px_] += 0.5*x[x_]*x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5*x[y_]*x[y_]; |
---|
| 1106 | x[py_] -= -irho*tan(-dtor(psi))*x[y_]+1*x[x_]*x[y_]-irho*x[px_]*x[y_]+ (irho*psi/cos(-dtor(psi))/cos(-dtor(psi))) * x[y_] |
---|
| 1107 | / (1.0 + x[delta_]); |
---|
| 1108 | } |
---|
| 1109 | */ |
---|
| 1110 | |
---|
[11] | 1111 | //add the contribution of the entrance momentum of the particle, from Alex Loulergue. |
---|
[3] | 1112 | if(Entrance){ |
---|
[24] | 1113 | x[px_] += irho * tan(dtor(phi) ) * x0[x_]; |
---|
| 1114 | x[py_] -= irho * tan(dtor(phi) + x0[px_]/(1+x0[delta_]*1) - get_psi(irho, phi, gap)/(1+x0[delta_]*0)) * x0[y_] |
---|
| 1115 | / (1.0 + x0[delta_]*0); |
---|
[3] | 1116 | }else{ |
---|
[24] | 1117 | x[px_] += irho * tan(dtor(phi)) * x0[x_]; |
---|
| 1118 | x[py_] -= irho * tan(dtor(phi) - x0[px_]/(1+x0[delta_]*1) - get_psi(irho, phi, gap)/(1+x0[delta_]*0)) * x0[y_] |
---|
| 1119 | / (1.0 + x0[delta_]*0); |
---|
[3] | 1120 | } |
---|
| 1121 | |
---|
| 1122 | } else {//original one |
---|
| 1123 | |
---|
[11] | 1124 | //cout << "Dipole edge effect WITHOUT Alex's correction: " << endl; |
---|
[3] | 1125 | x[px_] += irho * tan(dtor(phi)) * x[x_]; |
---|
[11] | 1126 | x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]/ (1.0 + x[delta_]*1); |
---|
[3] | 1127 | } |
---|
| 1128 | } |
---|
| 1129 | |
---|
| 1130 | |
---|
| 1131 | /****************************************************************** |
---|
[11] | 1132 | template<typename T> |
---|
| 1133 | static void BendCurvature(double irho, double H, ss_vect<T> &x) |
---|
| 1134 | |
---|
| 1135 | The leanding order term T211, T233, and T413 on page 117 & page 118 |
---|
| 1136 | of the dipole edge effect in SLAC-75, using the symplectic |
---|
| 1137 | form in |
---|
[3] | 1138 | E. Forest: "Fringe field in MAD, Part II: Bend curvature in MAD-X for the |
---|
| 1139 | mudule PTC". |
---|
| 1140 | P.10 Eq. (5) and (42). |
---|
| 1141 | |
---|
[11] | 1142 | For an rectangular dipole, the bending curvature term is zero, since 1/R = 0. |
---|
| 1143 | |
---|
| 1144 | irho: curvature of the central orbit inside the dipole |
---|
| 1145 | H: curvature of the entrance/exit pole face of the dipole |
---|
| 1146 | |
---|
| 1147 | Comments: |
---|
| 1148 | Written by Jianfeng Zhang @ LAL, 01/10/2012. |
---|
[3] | 1149 | ******************************************************************/ |
---|
| 1150 | template<typename T> |
---|
[11] | 1151 | static void BendCurvature(double irho, double H, ss_vect<T> &x) { |
---|
[3] | 1152 | |
---|
[11] | 1153 | if(trace) |
---|
| 1154 | cout << "Forest's bend curvature...."<<endl; |
---|
| 1155 | |
---|
[3] | 1156 | T pm2,u, coeff1, coeff2,coeff3; |
---|
| 1157 | |
---|
[11] | 1158 | ss_vect<T> x0; |
---|
[3] | 1159 | |
---|
[11] | 1160 | x0 = x; |
---|
[3] | 1161 | |
---|
[11] | 1162 | pm2 = sqr(1+x0[delta_]) - sqr(x0[px_]); |
---|
| 1163 | u = irho*H/2/pm2; |
---|
[3] | 1164 | |
---|
[11] | 1165 | coeff1 = u * 2*x0[px_]*(1+x0[delta_])/pm2; // |
---|
| 1166 | coeff2 = u * (1+x0[delta_]); // |
---|
| 1167 | coeff3 = u * (pm2 - 2*(1+x0[delta_]))/pm2; |
---|
[3] | 1168 | |
---|
| 1169 | |
---|
[11] | 1170 | x[x_] += x0[x_]/(1-coeff1*x0[y_]*x0[y_]); |
---|
| 1171 | x[px_] -= coeff2*x0[y_]*x0[y_] - irho*H/2*sqr(x0[x_]); |
---|
| 1172 | x[py_] -= 2*coeff2*x[x_]*x0[y_]; |
---|
| 1173 | x[ct_] -= coeff3*x[x_]*x0[y_]*x0[y_]; |
---|
[3] | 1174 | } |
---|
| 1175 | |
---|
| 1176 | |
---|
| 1177 | |
---|
| 1178 | |
---|
| 1179 | |
---|
| 1180 | /**************************************************************** |
---|
| 1181 | template<typename T> |
---|
| 1182 | void p_rot(double phi, ss_vect<T> &x) |
---|
| 1183 | |
---|
| 1184 | Purpose: |
---|
| 1185 | Rotate from the beam coorinate to the dipole face. |
---|
| 1186 | This is a pure geometric operation, no physical meaning. |
---|
| 1187 | |
---|
| 1188 | See P.307 Eq. (10.26) in |
---|
| 1189 | E. Forest "Beam dynamics: A new attitude and framework". |
---|
| 1190 | |
---|
| 1191 | Input: |
---|
| 1192 | phi: entrance or exit angle of the dipole edge |
---|
| 1193 | x: initial cooridinate of the particle |
---|
| 1194 | |
---|
| 1195 | Output: |
---|
| 1196 | |
---|
| 1197 | Return: |
---|
| 1198 | |
---|
| 1199 | ****************************************************************/ |
---|
| 1200 | template<typename T> |
---|
| 1201 | void p_rot(double phi, ss_vect<T> &x) { |
---|
| 1202 | T c, s, ps, p; |
---|
| 1203 | ss_vect<T> x1; |
---|
| 1204 | |
---|
| 1205 | c = cos(dtor(phi)); |
---|
| 1206 | s = sin(dtor(phi)); |
---|
| 1207 | x1 = x; |
---|
| 1208 | ps = get_p_s(x); |
---|
| 1209 | p = c * ps - s * x1[px_]; |
---|
| 1210 | x[x_] = x1[x_] * ps / p; |
---|
| 1211 | x[px_] = s * ps + c * x1[px_]; |
---|
| 1212 | x[y_] += x1[x_] * x1[py_] * s / p; |
---|
| 1213 | x[ct_] += (1.0 + x1[delta_]) * x1[x_] * s / p; |
---|
| 1214 | } |
---|
| 1215 | |
---|
| 1216 | |
---|
| 1217 | /*************************************************************** |
---|
| 1218 | template<typename T> |
---|
| 1219 | void bend_fringe(double hb, ss_vect<T> &x) |
---|
| 1220 | |
---|
| 1221 | Purpose: |
---|
| 1222 | the effect of longitudinal fringe field using exact Hamitonian. |
---|
| 1223 | This is a hard effect model, the fringe field and edge focus |
---|
| 1224 | happens at the entrance/exit point of the dipole pole face. |
---|
| 1225 | |
---|
| 1226 | include only the first order of irho. |
---|
| 1227 | Input: |
---|
| 1228 | |
---|
| 1229 | Output: |
---|
| 1230 | |
---|
| 1231 | Return: |
---|
| 1232 | |
---|
| 1233 | ****************************************************************/ |
---|
| 1234 | template<typename T> |
---|
| 1235 | void bend_fringe(double hb, ss_vect<T> &x) { |
---|
| 1236 | T coeff, u, ps, ps2, ps3; |
---|
| 1237 | ss_vect<T> x1; |
---|
| 1238 | |
---|
| 1239 | coeff = -hb / 2.0; |
---|
| 1240 | x1 = x; |
---|
| 1241 | ps = get_p_s(x); |
---|
| 1242 | ps2 = sqr(ps); |
---|
| 1243 | ps3 = ps * ps2; |
---|
| 1244 | u = 1.0 + 4.0 * coeff * x1[px_] * x1[y_] * x1[py_] / ps3; |
---|
| 1245 | if (u >= 0.0) { |
---|
| 1246 | x[y_] = 2.0 * x1[y_] / (1.0 + sqrt(u)); |
---|
| 1247 | x[x_] = x1[x_] - coeff * sqr(x[y_]) * (ps2 + sqr(x1[px_])) / ps3; |
---|
| 1248 | x[py_] = x1[py_] + 2.0 * coeff * x1[px_] * x[y_] / ps; |
---|
| 1249 | x[ct_] = x1[ct_] - coeff * x1[px_] * sqr(x[y_]) * (1.0 + x1[delta_]) |
---|
| 1250 | / ps3; |
---|
| 1251 | } else {//give the error message and set the unstable x value during the tracking |
---|
| 1252 | //this x value will be discard as an unstable value during the tracking |
---|
| 1253 | printf("bend_fringe: *** Speed of light exceeded!\n"); |
---|
| 1254 | x[y_] = 10; |
---|
| 1255 | x[x_] = 10; |
---|
| 1256 | x[py_] =10; |
---|
| 1257 | x[ct_] = 10; |
---|
| 1258 | //exit_(1); |
---|
| 1259 | } |
---|
| 1260 | } |
---|
| 1261 | |
---|
| 1262 | /*************************************************************** |
---|
| 1263 | template<typename T> |
---|
| 1264 | void bend_fringe(double hb, double gap, ss_vect<T> &x) |
---|
| 1265 | |
---|
| 1266 | Purpose: |
---|
| 1267 | the effect of longitudinal fringe field using exact Hamitonian. |
---|
| 1268 | This is a hard effect model, the fringe field and edge focus |
---|
| 1269 | happens at the entrance/exit point of the dipole pole face. |
---|
| 1270 | |
---|
| 1271 | include the second order of irho. |
---|
| 1272 | |
---|
| 1273 | See E. Forest and et al. |
---|
| 1274 | "Fringe field in MAD Part I: Second order Fringe in MAD-X for the module PTC", |
---|
| 1275 | P. 8-9, Eq. (35) and (34). |
---|
| 1276 | |
---|
| 1277 | Input: |
---|
| 1278 | |
---|
| 1279 | Output: |
---|
| 1280 | |
---|
| 1281 | Return: |
---|
| 1282 | |
---|
| 1283 | ****************************************************************/ |
---|
| 1284 | template<typename T> |
---|
| 1285 | void bend_fringe(double hb, double gap, ss_vect<T> &x) { |
---|
[10] | 1286 | |
---|
| 1287 | bool track=false; |
---|
| 1288 | |
---|
| 1289 | if(track) cout << "bend_fringe(): Forest's model"<<endl; |
---|
| 1290 | |
---|
[3] | 1291 | T coeff1,coeff2,coeff3,coeff4; |
---|
| 1292 | T ps, ps2, ps3, ps5; |
---|
| 1293 | |
---|
| 1294 | ss_vect<T> x1; |
---|
| 1295 | |
---|
| 1296 | gap = gap*0.5; |
---|
| 1297 | |
---|
| 1298 | x1 = x; |
---|
| 1299 | ps = get_p_s(x); |
---|
| 1300 | |
---|
| 1301 | if(ps==0) |
---|
| 1302 | return; |
---|
| 1303 | |
---|
| 1304 | ps2 = sqr(ps); |
---|
| 1305 | ps3 = ps * ps2; |
---|
| 1306 | ps5 = ps2*ps3; |
---|
| 1307 | |
---|
| 1308 | coeff1 = 1.0 - 4.0*sqr(hb)*gap*x1[py_]*x1[y_]/ps3; |
---|
| 1309 | coeff2 = sqr(hb)*gap*2*x1[px_]*(2*sqr(x1[px_])-sqr(1.0+x1[delta_]))/ps5 + |
---|
| 1310 | hb*(ps/(sqr(1+x1[delta_])-sqr(x1[px_])) + 2*sqr(x1[px_])*ps/sqr(sqr(1+x1[delta_])-sqr(x1[px_]))); |
---|
| 1311 | |
---|
| 1312 | coeff3 = hb*(x1[px_]/ps)/(1+sqr(x1[py_])/ps2) - sqr(hb)*gap*((ps2+sqr(x1[px_]))/ps3+sqr(x1[px_])/ps2*(ps2+x1[py_])/ps3); |
---|
| 1313 | coeff4 = -sqr(hb)*gap*4*sqr(x1[px_])*(1+x1[delta_])/ps5; |
---|
| 1314 | |
---|
| 1315 | if (coeff1 >= 0.0) { |
---|
| 1316 | x[y_] = 2.0 * x1[y_] / (1.0 + sqrt(coeff1)); |
---|
| 1317 | } else { |
---|
| 1318 | //give the error message and set the unstable x value during the tracking |
---|
| 1319 | //this x value will be discard as an unstable value during the tracking |
---|
| 1320 | printf("bend_fringe: *** Speed of light exceeded!\n"); |
---|
| 1321 | x[y_] = 10; |
---|
| 1322 | x[x_] = 10; |
---|
| 1323 | x[py_] =10; |
---|
| 1324 | x[ct_] = 10; |
---|
| 1325 | // exit_(1); |
---|
| 1326 | } |
---|
| 1327 | |
---|
| 1328 | x[x_] += 0.5*coeff2*sqr(x[y_]); |
---|
| 1329 | x[py_]-= coeff3*x[y_]; |
---|
| 1330 | x[ct_] -= 0.5*coeff4*sqr(x[y_]); |
---|
| 1331 | } |
---|
| 1332 | |
---|
| 1333 | /**************************************************************************** |
---|
| 1334 | * template<typename T> |
---|
| 1335 | void quad_fringe(double b2, ss_vect<T> &x) |
---|
| 1336 | |
---|
| 1337 | Purpose: |
---|
| 1338 | Compute edge focusing for a quadrupole |
---|
| 1339 | There is no radiation coming from the edge |
---|
| 1340 | |
---|
| 1341 | The standard formula used is using more general form with exponential |
---|
| 1342 | form. If keep up to the 4-th order nonlinear terms, the formula with goes to the |
---|
| 1343 | following: |
---|
| 1344 | (see E. Forest's book (Beam Dynamics: A New Attitude and Framework), p390): |
---|
| 1345 | b2 |
---|
| 1346 | x = x0 +/- ---------- (x0^3 + 3*z0^2*x0) |
---|
| 1347 | 12(1 + dP) |
---|
| 1348 | |
---|
| 1349 | b2 |
---|
| 1350 | px = px0 +/- ---------- (2*x0*y0*py0 - x0^2*px0 - y0^2*py0) |
---|
| 1351 | 4(1 + dP) |
---|
| 1352 | |
---|
| 1353 | b2 |
---|
| 1354 | y = y0 -/+ ---------- (y0^3 + 3*x0^2*y0) |
---|
| 1355 | 12(1 + dP) |
---|
| 1356 | |
---|
| 1357 | b2 |
---|
| 1358 | py = py0 -/+ ---------- (2*x0*y0*px0 - y0^2*py0 - x0^2*py0) |
---|
| 1359 | 4(1 + dP) |
---|
| 1360 | |
---|
| 1361 | dP = dP0; |
---|
| 1362 | |
---|
| 1363 | |
---|
| 1364 | b2 |
---|
| 1365 | tau = tau0 -/+ ----------- (y0^3*px - x0^3*py + 3*x0^2*y*py - 3*y0^2*x0*px) |
---|
| 1366 | 12(1 + dP)^2 |
---|
| 1367 | |
---|
| 1368 | Input: |
---|
| 1369 | b2 = quadrupole strength |
---|
| 1370 | x = input particle coordinates |
---|
| 1371 | |
---|
| 1372 | Output: |
---|
| 1373 | x output particle coordinates |
---|
| 1374 | |
---|
| 1375 | Return: |
---|
| 1376 | none |
---|
| 1377 | |
---|
| 1378 | Global variables: |
---|
| 1379 | none |
---|
| 1380 | |
---|
| 1381 | Specific functions: |
---|
| 1382 | none |
---|
| 1383 | |
---|
| 1384 | Comments: |
---|
| 1385 | Now in Tracy III, no definition "entrance" and "exit", when called in Mpole_pass, |
---|
| 1386 | first call with M --> PB[quad+HOMmax], then |
---|
| 1387 | call with -M --> PB[quad+HOMmax] |
---|
| 1388 | |
---|
| 1389 | ****************************************************************************/ |
---|
| 1390 | template<typename T> |
---|
| 1391 | void quad_fringe(double b2, ss_vect<T> &x) { |
---|
| 1392 | T u, ps; |
---|
| 1393 | |
---|
| 1394 | u = b2 / (12.0 * (1.0 + x[delta_])); |
---|
| 1395 | ps = u / (1.0 + x[delta_]); |
---|
| 1396 | |
---|
| 1397 | x[py_] /= 1.0 - 3.0 * u * sqr(x[y_]); |
---|
| 1398 | x[y_] -= u * cube(x[y_]); |
---|
| 1399 | |
---|
| 1400 | if (globval.Cavity_on) |
---|
| 1401 | x[ct_] -= ps * cube(x[y_]) * x[py_]; //-y^3*py |
---|
| 1402 | |
---|
| 1403 | x[px_] /= 1.0 + 3.0 * u * sqr(x[x_]); //+x^2 |
---|
| 1404 | |
---|
| 1405 | if (globval.Cavity_on) |
---|
| 1406 | x[ct_] += ps * cube(x[x_]) * x[px_]; //+x^3*px |
---|
| 1407 | |
---|
| 1408 | x[x_] += u * cube(x[x_]); //+x^3 |
---|
| 1409 | u = u * 3.0; |
---|
| 1410 | ps = ps * 3.0; |
---|
| 1411 | x[y_] = exp(-u * sqr(x[x_])) * x[y_]; //+x^2*y |
---|
| 1412 | x[py_] = exp(u * sqr(x[x_])) * x[py_]; //+x^2*py |
---|
| 1413 | x[px_] += 2.0 * u * x[x_] * x[y_] * x[py_]; //+2*x*y*py |
---|
| 1414 | |
---|
| 1415 | if (globval.Cavity_on) |
---|
| 1416 | x[ct_] -= ps * sqr(x[x_]) * x[y_] * x[py_]; // -3*x^2*y*py |
---|
| 1417 | |
---|
| 1418 | x[x_] = exp(u * sqr(x[y_])) * x[x_]; //+x*y^2 |
---|
| 1419 | x[px_] = exp(-u * sqr(x[y_])) * x[px_]; // -x^2*px-y^2*px |
---|
| 1420 | x[py_] -= 2.0 * u * x[y_] * x[x_] * x[px_]; //-2*x*y*px |
---|
| 1421 | |
---|
| 1422 | if (globval.Cavity_on) |
---|
| 1423 | x[ct_] += ps * sqr(x[y_]) * x[x_] * x[px_]; // +3*y^2*x*px |
---|
| 1424 | } |
---|
| 1425 | |
---|
| 1426 | /****************************************************************************/ |
---|
| 1427 | /* void Mpole_Pass(CellType &Cell, ss_vect<T> &x) |
---|
| 1428 | |
---|
| 1429 | Purpose: |
---|
| 1430 | multipole pass,for dipole, quadrupole,sextupole,decupole,etc |
---|
| 1431 | Using DA method. |
---|
| 1432 | |
---|
| 1433 | Input: |
---|
| 1434 | |
---|
| 1435 | Output: |
---|
| 1436 | |
---|
| 1437 | |
---|
| 1438 | Return: |
---|
| 1439 | none |
---|
| 1440 | |
---|
| 1441 | Global variables: |
---|
| 1442 | none |
---|
| 1443 | |
---|
| 1444 | Specific functions: |
---|
| 1445 | none |
---|
| 1446 | |
---|
| 1447 | Comments: |
---|
| 1448 | none |
---|
| 1449 | |
---|
| 1450 | ****************************************************************************/ |
---|
| 1451 | |
---|
| 1452 | template<typename T> |
---|
| 1453 | void Mpole_Pass(CellType &Cell, ss_vect<T> &x) { |
---|
| 1454 | int seg = 0; |
---|
| 1455 | double k = 0.0, dL = 0.0, dL1 = 0.0, dL2 = 0.0; |
---|
| 1456 | double dkL1 = 0.0, dkL2 = 0.0, h_ref = 0.0; |
---|
| 1457 | elemtype *elemp; |
---|
| 1458 | MpoleType *M; |
---|
| 1459 | |
---|
| 1460 | ss_vect<T> x1, x2; |
---|
| 1461 | |
---|
| 1462 | |
---|
| 1463 | elemp = &Cell.Elem; |
---|
| 1464 | M = elemp->M; |
---|
| 1465 | |
---|
| 1466 | /* Global -> Local */ |
---|
| 1467 | GtoL(x, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1); |
---|
| 1468 | |
---|
| 1469 | //entrance of the magnet |
---|
| 1470 | if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) { /* fringe fields */ |
---|
| 1471 | |
---|
| 1472 | if (globval.quad_fringe && (M->PB[Quad + HOMmax] != 0.0) && (M->quadFF1 |
---|
| 1473 | == 1)) { |
---|
| 1474 | quad_fringe(M->quadFFscaling * M->PB[Quad + HOMmax], x); |
---|
| 1475 | } |
---|
| 1476 | |
---|
| 1477 | //fringe field and edge focusing of dipoles |
---|
| 1478 | if (M->Pirho != 0.0 && M->dipEdge_effect1 == 1){ |
---|
[11] | 1479 | if (!globval.H_exact) { //big ring */ modified linear term T43 on page 117 & 118 in SLAC-75. |
---|
| 1480 | EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x,true); |
---|
| 1481 | } else {//small and big rings; Forest's model |
---|
| 1482 | if(M->PH1!=0){ |
---|
| 1483 | // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0. |
---|
| 1484 | BendCurvature(M->Pirho, M->PH1, x); |
---|
| 1485 | } |
---|
| 1486 | p_rot(M->PTx1, x); //rotate from cartesian cooridnate to curlinear curved beam coordinate; |
---|
[3] | 1487 | // since the map of a sector dipole is used in Tracy 3 |
---|
| 1488 | bend_fringe(M->Pirho, M->Pgap, x); //fringe field |
---|
| 1489 | } |
---|
| 1490 | } |
---|
| 1491 | } |
---|
| 1492 | |
---|
| 1493 | if (M->Pthick == thick) { |
---|
| 1494 | if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates |
---|
| 1495 | // if (M->n_design == Dip || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates |
---|
| 1496 | h_ref = M->Pirho; |
---|
| 1497 | dL = elemp->PL / M->PN; |
---|
| 1498 | } else {// Cartesian coordinates |
---|
| 1499 | h_ref = 0.0; |
---|
| 1500 | if (M->Pirho == 0.0) |
---|
| 1501 | dL = elemp->PL / M->PN; |
---|
| 1502 | else |
---|
| 1503 | dL = 1.0 / M->Pirho * (sin(dtor(M->PTx1)) + sin(dtor(M->PTx2))) //straight length of the dipole |
---|
| 1504 | / M->PN; |
---|
| 1505 | } |
---|
| 1506 | } |
---|
| 1507 | |
---|
| 1508 | switch (M->Pmethod) { |
---|
| 1509 | |
---|
| 1510 | case Meth_Linear: |
---|
| 1511 | |
---|
| 1512 | case Meth_First: |
---|
| 1513 | if (M->Pthick == thick) { |
---|
| 1514 | /* First Linear */ |
---|
| 1515 | // LinTrans(5L, M->AU55, x); |
---|
| 1516 | k = M->PB[Quad + HOMmax]; |
---|
| 1517 | /* retrieve normal quad component already in AU55 */ |
---|
| 1518 | M->PB[Quad + HOMmax] = 0.0; |
---|
| 1519 | /* Kick w/o quad component */ |
---|
| 1520 | thin_kick(M->Porder, M->PB, elemp->PL, 0.0, 0.0, x); |
---|
| 1521 | /* restore quad component */ |
---|
| 1522 | M->PB[Quad + HOMmax] = k; |
---|
| 1523 | /* Second Linear */ |
---|
| 1524 | // LinTrans(5L, M->AD55, x); |
---|
| 1525 | } else |
---|
| 1526 | /* thin kick */ |
---|
| 1527 | thin_kick(M->Porder, M->PB, 1.0, 0.0, 0.0, x); |
---|
| 1528 | break; |
---|
| 1529 | |
---|
| 1530 | case Meth_Second: |
---|
[11] | 1531 | |
---|
| 1532 | // cout << "Mpole_Pass: Meth_Second not supported" << endl; |
---|
| 1533 | // exit_(0); |
---|
| 1534 | // break; |
---|
[3] | 1535 | |
---|
[11] | 1536 | //specially for the test of the dipole map,see Forest's beam dynamics, p361. eqn.(12.19) |
---|
| 1537 | |
---|
| 1538 | dL1 = 0.5*dL; |
---|
| 1539 | dL2 = dL; |
---|
| 1540 | |
---|
| 1541 | for (seg = 1; seg <= M->PN; seg++) { |
---|
| 1542 | dipole_kick(dL1,M->Pirho,h_ref,x); |
---|
| 1543 | multipole_kick(M->Porder, M->PB, dkL1, M->Pirho, x); |
---|
| 1544 | dipole_kick(dL1,M->Pirho,h_ref,x); |
---|
| 1545 | } |
---|
| 1546 | |
---|
| 1547 | |
---|
| 1548 | |
---|
[3] | 1549 | case Meth_Fourth: //forth order symplectic integrator |
---|
| 1550 | if (M->Pthick == thick) { |
---|
| 1551 | dL1 = c_1 * dL; |
---|
| 1552 | dL2 = c_2 * dL; |
---|
| 1553 | dkL1 = d_1 * dL; |
---|
| 1554 | dkL2 = d_2 * dL; |
---|
| 1555 | |
---|
| 1556 | dcurly_H = 0.0; |
---|
| 1557 | dI4 = 0.0; |
---|
| 1558 | for (seg = 1; seg <= M->PN; seg++) { |
---|
| 1559 | if (globval.emittance && (!globval.Cavity_on) && (M->Pirho |
---|
| 1560 | != 0.0)) { |
---|
| 1561 | dcurly_H += is_tps<tps>::get_curly_H(x); |
---|
| 1562 | dI4 += is_tps<tps>::get_dI4(x); |
---|
| 1563 | } |
---|
| 1564 | |
---|
[11] | 1565 | Drift(dL1, x); |
---|
[3] | 1566 | thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x); |
---|
| 1567 | Drift(dL2, x); |
---|
[11] | 1568 | thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x); |
---|
[3] | 1569 | |
---|
| 1570 | if (globval.emittance && (!globval.Cavity_on) && (M->Pirho |
---|
| 1571 | != 0.0)) { |
---|
| 1572 | dcurly_H += 4.0 * is_tps<tps>::get_curly_H(x); |
---|
| 1573 | dI4 += 4.0 * is_tps<tps>::get_dI4(x); |
---|
| 1574 | } |
---|
[11] | 1575 | |
---|
| 1576 | |
---|
[3] | 1577 | Drift(dL2, x); |
---|
[11] | 1578 | thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x); |
---|
| 1579 | Drift(dL1, x); |
---|
[3] | 1580 | |
---|
| 1581 | if (globval.emittance && (!globval.Cavity_on) && (M->Pirho |
---|
| 1582 | != 0.0)) { |
---|
| 1583 | dcurly_H += is_tps<tps>::get_curly_H(x); |
---|
| 1584 | dI4 += is_tps<tps>::get_dI4(x); |
---|
| 1585 | } |
---|
| 1586 | } |
---|
| 1587 | |
---|
| 1588 | if (globval.emittance && (!globval.Cavity_on) && (M->Pirho != 0)) { |
---|
| 1589 | dcurly_H /= 6.0 * M->PN; |
---|
| 1590 | dI4 *= M->Pirho * (sqr(M->Pirho) + 2.0 |
---|
| 1591 | * M->PBpar[Quad + HOMmax]) / (6.0 * M->PN); |
---|
| 1592 | |
---|
| 1593 | I2 += elemp->PL * sqr(M->Pirho); |
---|
| 1594 | I4 += elemp->PL * dI4; |
---|
| 1595 | I5 += elemp->PL * fabs(cube(M->Pirho)) * dcurly_H; |
---|
| 1596 | } |
---|
| 1597 | } else |
---|
| 1598 | thin_kick(M->Porder, M->PB, 1.0, 0.0, 0.0, x); |
---|
| 1599 | |
---|
| 1600 | break; |
---|
| 1601 | } |
---|
| 1602 | |
---|
| 1603 | //exit of the magnets |
---|
| 1604 | if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) { |
---|
| 1605 | |
---|
| 1606 | /* dipole fringe fields */ |
---|
| 1607 | if (M->Pirho != 0.0 && M->dipEdge_effect2 == 1){ |
---|
[11] | 1608 | if (!globval.H_exact) { //big ring, linear model correction |
---|
| 1609 | EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x,false); |
---|
| 1610 | }else{ |
---|
| 1611 | bend_fringe(-M->Pirho, M->Pgap,x); //Forest's fringe field of the dipole |
---|
| 1612 | p_rot(M->PTx2, x); //rotate back to the cartesian cooridinate |
---|
| 1613 | if(M->PH2!=0){ |
---|
| 1614 | // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0. |
---|
| 1615 | BendCurvature(M->Pirho, M->PH2, x); |
---|
[3] | 1616 | } |
---|
| 1617 | } |
---|
| 1618 | } |
---|
| 1619 | //quadrupole fringe field |
---|
| 1620 | if (globval.quad_fringe && (M->PB[Quad + HOMmax] != 0.0) && (M->quadFF2 |
---|
| 1621 | == 1)) |
---|
| 1622 | quad_fringe(-M->quadFFscaling * M->PB[Quad + HOMmax], x); |
---|
| 1623 | } |
---|
| 1624 | |
---|
| 1625 | /* Local -> Global */ |
---|
| 1626 | LtoG(x, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1); |
---|
| 1627 | } |
---|
| 1628 | |
---|
| 1629 | template<typename T> |
---|
| 1630 | void Marker_Pass(CellType &Cell, ss_vect<T> &X) { |
---|
| 1631 | elemtype *elemp; |
---|
| 1632 | |
---|
| 1633 | elemp = &Cell.Elem; |
---|
| 1634 | /* Global -> Local */ |
---|
| 1635 | GtoL(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 1636 | /* Local -> Global */ |
---|
| 1637 | LtoG(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 1638 | } |
---|
| 1639 | |
---|
| 1640 | /**************************************************************************** |
---|
| 1641 | * void Cav_Pass(CellType *Cell, double *X) |
---|
| 1642 | |
---|
| 1643 | Purpose: |
---|
| 1644 | Tracking through a cavity |
---|
| 1645 | |
---|
| 1646 | Input: |
---|
| 1647 | Cell cavity element to track through |
---|
| 1648 | X input coordinates |
---|
| 1649 | |
---|
| 1650 | Output: |
---|
| 1651 | X output coordinates |
---|
| 1652 | |
---|
| 1653 | Return: |
---|
| 1654 | none |
---|
| 1655 | |
---|
| 1656 | Global variables: |
---|
| 1657 | none |
---|
| 1658 | |
---|
| 1659 | Specific functions: |
---|
| 1660 | none |
---|
| 1661 | |
---|
| 1662 | Comments: |
---|
| 1663 | none |
---|
| 1664 | |
---|
| 1665 | ****************************************************************************/ |
---|
| 1666 | template<typename T> |
---|
| 1667 | void Cav_Pass(CellType &Cell, ss_vect<T> &X) { |
---|
| 1668 | elemtype *elemp; |
---|
| 1669 | CavityType *C; |
---|
| 1670 | T delta; |
---|
| 1671 | |
---|
| 1672 | elemp = &Cell.Elem; |
---|
| 1673 | C = elemp->C; |
---|
| 1674 | if (globval.Cavity_on && C->Pvolt != 0.0) { |
---|
| 1675 | delta = -C->Pvolt / (globval.Energy * 1e9) * sin(2.0 * M_PI * C->Pfreq |
---|
| 1676 | / c0 * X[ct_] + C->phi); |
---|
| 1677 | X[delta_] += delta; |
---|
| 1678 | |
---|
| 1679 | if (globval.radiation) |
---|
| 1680 | globval.dE -= is_double<T>::cst(delta); |
---|
| 1681 | |
---|
| 1682 | if (globval.pathlength) |
---|
| 1683 | X[ct_] -= C->Ph / C->Pfreq * c0; |
---|
| 1684 | } |
---|
| 1685 | } |
---|
| 1686 | |
---|
| 1687 | template<typename T> |
---|
| 1688 | inline void get_Axy(const WigglerType *W, const double z, ss_vect<T> &x, |
---|
| 1689 | T AxoBrho[], T AyoBrho[]) |
---|
| 1690 | |
---|
| 1691 | { |
---|
| 1692 | int i; |
---|
| 1693 | double ky, kz_n; |
---|
| 1694 | T cx, cz, sx, sz, chy, shy; |
---|
| 1695 | |
---|
| 1696 | for (i = 0; i <= 3; ++i) { |
---|
| 1697 | AxoBrho[i] = 0.0; |
---|
| 1698 | AyoBrho[i] = 0.0; |
---|
| 1699 | } |
---|
| 1700 | |
---|
| 1701 | for (i = 0; i < W->n_harm; i++) { |
---|
| 1702 | kz_n = W->harm[i] * 2.0 * M_PI / W->lambda; |
---|
| 1703 | ky = sqrt(sqr(W->kxV[i]) + sqr(kz_n)); |
---|
| 1704 | cx = cos(W->kxV[i] * x[x_]); |
---|
| 1705 | sx = sin(W->kxV[i] * x[x_]); |
---|
| 1706 | chy = cosh(ky * x[y_]); |
---|
| 1707 | shy = sinh(ky * x[y_]); |
---|
| 1708 | sz = sin(kz_n * z); |
---|
| 1709 | |
---|
| 1710 | AxoBrho[0] += W->BoBrhoV[i] / kz_n * cx * chy * sz; |
---|
| 1711 | AyoBrho[0] += W->BoBrhoV[i] * W->kxV[i] / (ky * kz_n) * sx * shy * sz; |
---|
| 1712 | |
---|
| 1713 | // derivatives with respect to x |
---|
| 1714 | AxoBrho[1] -= W->BoBrhoV[i] * W->kxV[i] / kz_n * sx * chy * sz; |
---|
| 1715 | AyoBrho[1] += W->BoBrhoV[i] * sqr(W->kxV[i]) / (ky * kz_n) * cx * shy |
---|
| 1716 | * sz; |
---|
| 1717 | |
---|
| 1718 | // derivatives with respect to y |
---|
| 1719 | AxoBrho[2] += W->BoBrhoV[i] * ky / kz_n * cx * shy * sz; |
---|
| 1720 | AyoBrho[2] += W->BoBrhoV[i] * W->kxV[i] / kz_n * sx * chy * sz; |
---|
| 1721 | |
---|
| 1722 | if (globval.radiation) { |
---|
| 1723 | cz = cos(kz_n * z); |
---|
| 1724 | // derivatives with respect to z |
---|
| 1725 | AxoBrho[3] += W->BoBrhoV[i] * cx * chy * cz; |
---|
| 1726 | AyoBrho[3] += W->BoBrhoV[i] * W->kxV[i] / ky * sx * shy * cz; |
---|
| 1727 | } |
---|
| 1728 | } |
---|
| 1729 | } |
---|
| 1730 | |
---|
| 1731 | /* |
---|
| 1732 | template<typename T> |
---|
| 1733 | inline void get_Axy_map(const FieldMapType *FM, const double z, |
---|
| 1734 | const ss_vect<T> &x, T AxoBrho[], T AyoBrho[]) |
---|
| 1735 | { |
---|
| 1736 | float y, ax0, ax1, ax2, ay0, ay1, ay2; |
---|
| 1737 | |
---|
| 1738 | const float dy = 1e-3, dz = 1e-3; |
---|
| 1739 | |
---|
| 1740 | y = is_double<T>::cst(x[y_]); |
---|
| 1741 | |
---|
| 1742 | if ((z < FM->s_pos[1]) || (z > FM->s_pos[FM->n_s])) { |
---|
| 1743 | cout << scientific << setprecision(3) |
---|
| 1744 | << "get_Axy_map: s out of range " << z << endl; |
---|
| 1745 | exit_(1); |
---|
| 1746 | } |
---|
| 1747 | |
---|
| 1748 | if ((y < FM->y_pos[1]) || (y > FM->y_pos[FM->m_y])) { |
---|
| 1749 | cout << scientific << setprecision(3) |
---|
| 1750 | << "get_Axy_map: y out of range " << y << endl; |
---|
| 1751 | exit_(1); |
---|
| 1752 | } |
---|
| 1753 | |
---|
| 1754 | splin2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->AxoBrho2, FM->m_y, FM->n_s, |
---|
| 1755 | y, z, &ax1); |
---|
| 1756 | AxoBrho[0] = FM->scl*ax1; |
---|
| 1757 | |
---|
| 1758 | splin2(FM->y_pos, FM->s_pos, FM->AyoBrho, FM->AyoBrho2, FM->m_y, FM->n_s, |
---|
| 1759 | y, z, &ay1); |
---|
| 1760 | AyoBrho[0] = FM->scl*ay1; |
---|
| 1761 | |
---|
| 1762 | // derivatives with respect to x |
---|
| 1763 | AxoBrho[1] = FM->scl*0.0; AyoBrho[1] = FM->scl*0.0; |
---|
| 1764 | |
---|
| 1765 | // derivatives with respect to y |
---|
| 1766 | splin2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->AxoBrho2, FM->m_y, FM->n_s, |
---|
| 1767 | y+dy, z, &ax2); |
---|
| 1768 | splin2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->AxoBrho2, FM->m_y, FM->n_s, |
---|
| 1769 | y-dy, z, &ax1); |
---|
| 1770 | splin2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->AxoBrho2, FM->m_y, FM->n_s, |
---|
| 1771 | y, z, &ax0); |
---|
| 1772 | AxoBrho[2] = |
---|
| 1773 | (ax2-ax1)/(2.0*dy) + (ax2+ax1-2.0*ax0)/sqr(dy)*is_tps<T>::set_prm(y_+1); |
---|
| 1774 | AxoBrho[2] *= FM->scl; |
---|
| 1775 | |
---|
| 1776 | splin2(FM->y_pos, FM->s_pos, FM->AyoBrho, FM->AyoBrho2, FM->m_y, FM->n_s, |
---|
| 1777 | y+dy, z, &ay2); |
---|
| 1778 | splin2(FM->y_pos, FM->s_pos, FM->AyoBrho, FM->AyoBrho2, FM->m_y, FM->n_s, |
---|
| 1779 | y-dy, z, &ay1); |
---|
| 1780 | splin2(FM->y_pos, FM->s_pos, FM->AyoBrho, FM->AyoBrho2, FM->m_y, FM->n_s, |
---|
| 1781 | y, z, &ay0); |
---|
| 1782 | AyoBrho[2] = |
---|
| 1783 | (ay2-ay1)/(2.0*dy) + (ay2+ay1-2.0*ay0)/sqr(dy)*is_tps<T>::set_prm(y_+1); |
---|
| 1784 | AyoBrho[2] *= FM->scl; |
---|
| 1785 | |
---|
| 1786 | if (globval.radiation) { |
---|
| 1787 | // derivatives with respect to z |
---|
| 1788 | splin2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->AxoBrho2, FM->m_y, FM->n_s, |
---|
| 1789 | y, z+dz, &ax2); |
---|
| 1790 | splin2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->AxoBrho2, FM->m_y, FM->n_s, |
---|
| 1791 | y, z-dz, &ax1); |
---|
| 1792 | AxoBrho[3] = (ax2-ax1)/(2.0*dz); AxoBrho[3] *= FM->scl; |
---|
| 1793 | |
---|
| 1794 | splin2(FM->y_pos, FM->s_pos, FM->AyoBrho, FM->AyoBrho2, FM->m_y, FM->n_s, |
---|
| 1795 | y, z+dz, &ay2); |
---|
| 1796 | splin2(FM->y_pos, FM->s_pos, FM->AyoBrho, FM->AyoBrho2, FM->m_y, FM->n_s, |
---|
| 1797 | y, z-dz, &ay1); |
---|
| 1798 | AyoBrho[3] = (ay2-ay1)/(2.0*dz); AyoBrho[3] *= FM->scl; |
---|
| 1799 | if (false) |
---|
| 1800 | cout << fixed << setprecision(5) |
---|
| 1801 | << setw(8) << z << setw(9) << is_double<T>::cst(AxoBrho[3]) << endl; |
---|
| 1802 | } |
---|
| 1803 | } |
---|
| 1804 | */ |
---|
| 1805 | |
---|
| 1806 | template<typename T> |
---|
| 1807 | void Wiggler_pass_EF(const elemtype &elem, ss_vect<T> &x) { |
---|
| 1808 | // First order symplectic integrator for wiggler using expanded Hamiltonian |
---|
| 1809 | |
---|
| 1810 | int i, nstep = 0; |
---|
| 1811 | double h, z; |
---|
| 1812 | T AxoBrho[4], AyoBrho[4], psi, hodp, a12, a21, a22, det; |
---|
| 1813 | T d1, d2, a11, c11, c12, c21, c22, x2, B[3]; |
---|
| 1814 | |
---|
| 1815 | switch (elem.Pkind) { |
---|
| 1816 | case Wigl: |
---|
| 1817 | nstep = elem.W->PN; |
---|
| 1818 | break; |
---|
| 1819 | case FieldMap: |
---|
| 1820 | nstep = elem.FM->n_step; |
---|
| 1821 | break; |
---|
| 1822 | default: |
---|
| 1823 | cout << "Wiggler_pass_EF: unknown element type" << endl; |
---|
| 1824 | exit_(1); |
---|
| 1825 | break; |
---|
| 1826 | } |
---|
| 1827 | |
---|
| 1828 | h = elem.PL / nstep; |
---|
| 1829 | z = 0.0; |
---|
| 1830 | for (i = 1; i <= nstep; ++i) { |
---|
| 1831 | switch (elem.Pkind) { |
---|
| 1832 | case Wigl: |
---|
| 1833 | get_Axy(elem.W, z, x, AxoBrho, AyoBrho); |
---|
| 1834 | break; |
---|
| 1835 | case FieldMap: |
---|
| 1836 | // get_Axy_map(elem.FM, z, x, AxoBrho, AyoBrho); |
---|
| 1837 | break; |
---|
| 1838 | default: |
---|
| 1839 | cout << "Wiggler_pass_EF: unknown element type" << endl; |
---|
| 1840 | exit_(1); |
---|
| 1841 | break; |
---|
| 1842 | } |
---|
| 1843 | |
---|
| 1844 | psi = 1.0 + x[delta_]; |
---|
| 1845 | hodp = h / psi; |
---|
| 1846 | a11 = hodp * AxoBrho[1]; |
---|
| 1847 | a12 = hodp * AyoBrho[1]; |
---|
| 1848 | a21 = hodp * AxoBrho[2]; |
---|
| 1849 | a22 = hodp * AyoBrho[2]; |
---|
| 1850 | det = 1.0 - a11 - a22 + a11 * a22 - a12 * a21; |
---|
| 1851 | d1 = hodp * AxoBrho[0] * AxoBrho[1]; |
---|
| 1852 | d2 = hodp * AxoBrho[0] * AxoBrho[2]; |
---|
| 1853 | c11 = (1.0 - a22) / det; |
---|
| 1854 | c12 = a12 / det; |
---|
| 1855 | c21 = a21 / det; |
---|
| 1856 | c22 = (1.0 - a11) / det; |
---|
| 1857 | x2 = c11 * (x[px_] - d1) + c12 * (x[py_] - d2); |
---|
| 1858 | |
---|
| 1859 | x[py_] = c21 * (x[px_] - d1) + c22 * (x[py_] - d2); |
---|
| 1860 | x[px_] = x2; |
---|
| 1861 | x[x_] += hodp * (x[px_] - AxoBrho[0]); |
---|
| 1862 | x[y_] += hodp * x[py_]; |
---|
| 1863 | x[ct_] += h * (sqr((x[px_] - AxoBrho[0]) / psi) + sqr((x[py_] |
---|
| 1864 | - AyoBrho[0]) / psi)) / 2.0; |
---|
| 1865 | |
---|
| 1866 | if (false) |
---|
| 1867 | cout << scientific << setprecision(3) << setw(8) << z << setw(11) |
---|
| 1868 | << is_double<T>::cst(x[x_]) << setw(11) |
---|
| 1869 | << is_double<T>::cst(x[px_]) << setw(11) |
---|
| 1870 | << is_double<T>::cst(x[y_]) << setw(11) |
---|
| 1871 | << is_double<T>::cst(x[py_]) << endl; |
---|
| 1872 | |
---|
| 1873 | if (globval.pathlength) |
---|
| 1874 | x[ct_] += h; |
---|
| 1875 | |
---|
| 1876 | if (globval.radiation || globval.emittance) { |
---|
| 1877 | B[X_] = -AyoBrho[3]; |
---|
| 1878 | B[Y_] = AxoBrho[3]; |
---|
| 1879 | B[Z_] = AyoBrho[1] - AxoBrho[2]; |
---|
| 1880 | radiate(x, h, 0.0, B); |
---|
| 1881 | } |
---|
| 1882 | |
---|
| 1883 | z += h; |
---|
| 1884 | } |
---|
| 1885 | } |
---|
| 1886 | |
---|
| 1887 | template<typename T> |
---|
| 1888 | inline void get_Axy2(const double z, const double kxV, const double kxH, |
---|
| 1889 | const double kz, const double BoBrhoV, const double BoBrhoH, |
---|
| 1890 | const double phi, ss_vect<T> &x, T AxoBrho[], T AyoBrho[]) { |
---|
| 1891 | int i; |
---|
| 1892 | T cx, sx, cz1, cz2, sz1, sz2, chy, shy, kyH, kyV, chx, shx, cy, sy; |
---|
| 1893 | |
---|
| 1894 | for (i = 0; i <= 3; ++i) { |
---|
| 1895 | AxoBrho[i] = 0.0; |
---|
| 1896 | AyoBrho[i] = 0.0; |
---|
| 1897 | } |
---|
| 1898 | |
---|
| 1899 | kyV = sqrt(sqr(kz) + sqr(kxV)); |
---|
| 1900 | kyH = sqrt(sqr(kz) + sqr(kxH)); |
---|
| 1901 | cx = cos(kxV * x[x_]); |
---|
| 1902 | sx = sin(kxV * x[x_]); |
---|
| 1903 | cy = cos(kxH * x[y_]); |
---|
| 1904 | sy = sin(kxH * x[y_]); |
---|
| 1905 | chx = cosh(kyH * x[x_]); |
---|
| 1906 | shx = sinh(kyH * x[x_]); |
---|
| 1907 | chy = cosh(kyV * x[y_]); |
---|
| 1908 | shy = sinh(kyV * x[y_]); |
---|
| 1909 | sz1 = sin(kz * z); |
---|
| 1910 | sz2 = sin(kz * z + phi); |
---|
| 1911 | |
---|
| 1912 | AxoBrho[0] += BoBrhoV / kz * cx * chy * sz1; |
---|
| 1913 | AxoBrho[0] -= BoBrhoH * kxH / (kyH * kz) * shx * sy * sz2; |
---|
| 1914 | AyoBrho[0] += BoBrhoV * kxV / (kyV * kz) * sx * shy * sz1; |
---|
| 1915 | AyoBrho[0] -= BoBrhoH / kz * chx * cy * sz2; |
---|
| 1916 | |
---|
| 1917 | /* derivatives with respect to x */ |
---|
| 1918 | AxoBrho[1] -= BoBrhoV * kxV / kz * sx * chy * sz1; |
---|
| 1919 | AxoBrho[1] -= BoBrhoH * kxH / kz * chx * sy * sz2; |
---|
| 1920 | AyoBrho[1] += BoBrhoV * sqr(kxV) / (kyV * kz) * cx * shy * sz1; |
---|
| 1921 | AyoBrho[1] -= BoBrhoH * kyH / kz * shx * cy * sz2; |
---|
| 1922 | |
---|
| 1923 | /* derivatives with respect to y */ |
---|
| 1924 | AxoBrho[2] += BoBrhoV * kyV / kz * cx * shy * sz1; |
---|
| 1925 | AxoBrho[2] -= BoBrhoH * sqr(kxH) / (kyH * kz) * shx * cy * sz2; |
---|
| 1926 | AyoBrho[2] += BoBrhoV * kxV / kz * sx * chy * sz1; |
---|
| 1927 | AyoBrho[2] += BoBrhoH * kxH / kz * chx * sy * sz2; |
---|
| 1928 | |
---|
| 1929 | if (globval.radiation) { |
---|
| 1930 | cz1 = cos(kz * z); |
---|
| 1931 | cz2 = cos(kz * z + phi); |
---|
| 1932 | /* derivatives with respect to z */ |
---|
| 1933 | AxoBrho[3] += BoBrhoV * cx * chy * cz1; |
---|
| 1934 | AxoBrho[3] -= BoBrhoH * kxH / kyH * shx * sy * cz2; |
---|
| 1935 | AyoBrho[3] += BoBrhoV * kxV / kyV * sx * shy * cz1; |
---|
| 1936 | AyoBrho[3] -= BoBrhoH * chx * cy * cz2; |
---|
| 1937 | } |
---|
| 1938 | } |
---|
| 1939 | |
---|
| 1940 | template<typename T> |
---|
| 1941 | void Wiggler_pass_EF2(int nstep, double L, double kxV, double kxH, double kz, |
---|
| 1942 | double BoBrhoV, double BoBrhoH, double phi, ss_vect<T> &x) { |
---|
| 1943 | // First order symplectic integrator for wiggler using expanded Hamiltonian |
---|
| 1944 | |
---|
| 1945 | int i; |
---|
| 1946 | double h, z; |
---|
| 1947 | T hodp, B[3], px1, px2, px3, py1, py2, AxoBrho[4], AyoBrho[4], psi; |
---|
| 1948 | T px = 0.0, py = 0.0; |
---|
| 1949 | |
---|
| 1950 | h = L / nstep; |
---|
| 1951 | z = 0.0; |
---|
| 1952 | for (i = 1; i <= nstep; ++i) { |
---|
| 1953 | get_Axy2(z, kxV, kxH, kz, BoBrhoV, BoBrhoH, phi, x, AxoBrho, AyoBrho); |
---|
| 1954 | |
---|
| 1955 | psi = 1.0 + x[delta_]; |
---|
| 1956 | hodp = h / psi; |
---|
| 1957 | |
---|
| 1958 | px1 = (x[px_] - (AxoBrho[0] * AxoBrho[1] + AyoBrho[0] * AyoBrho[1]) |
---|
| 1959 | * hodp) * (1 - AyoBrho[2] * hodp); |
---|
| 1960 | px2 = (x[py_] - (AxoBrho[0] * AxoBrho[2] + AyoBrho[0] * AyoBrho[2]) |
---|
| 1961 | * hodp) * AyoBrho[1] * hodp; |
---|
| 1962 | px3 = (1 - AxoBrho[1] * hodp) * (1 - AyoBrho[2] * hodp) - AxoBrho[2] |
---|
| 1963 | * AyoBrho[1] * hodp * hodp; |
---|
| 1964 | |
---|
| 1965 | py1 = (x[py_] - (AxoBrho[0] * AxoBrho[2] + AyoBrho[0] * AyoBrho[2]) |
---|
| 1966 | * hodp) * (1 - AxoBrho[1] * hodp); |
---|
| 1967 | py2 = (x[px_] - (AxoBrho[0] * AxoBrho[1] + AyoBrho[0] * AyoBrho[1]) |
---|
| 1968 | * hodp) * AxoBrho[2] * hodp; |
---|
| 1969 | |
---|
| 1970 | py = (py1 + py2) / px3; |
---|
| 1971 | px = (px1 + px2) / px3; |
---|
| 1972 | x[x_] += hodp * (px - AxoBrho[0]); |
---|
| 1973 | x[y_] += hodp * (py - AyoBrho[0]); |
---|
| 1974 | x[ct_] += h * (sqr((px - AxoBrho[0]) / psi) + sqr((py - AyoBrho[0]) |
---|
| 1975 | / psi)) / 2.0; |
---|
| 1976 | |
---|
| 1977 | if (globval.pathlength) |
---|
| 1978 | x[ct_] += h; |
---|
| 1979 | |
---|
| 1980 | if (globval.radiation || globval.emittance) { |
---|
| 1981 | B[X_] = -AyoBrho[3]; |
---|
| 1982 | B[Y_] = AxoBrho[3]; |
---|
| 1983 | B[Z_] = AyoBrho[1] - AxoBrho[2]; |
---|
| 1984 | radiate(x, h, 0.0, B); |
---|
| 1985 | } |
---|
| 1986 | |
---|
| 1987 | z += h; |
---|
| 1988 | } |
---|
| 1989 | |
---|
| 1990 | x[px_] = px; |
---|
| 1991 | x[py_] = py; |
---|
| 1992 | } |
---|
| 1993 | |
---|
| 1994 | template<typename T> |
---|
| 1995 | inline void get_Axy_EF3(const WigglerType *W, const double z, |
---|
| 1996 | const ss_vect<T> &x, T &AoBrho, T dAoBrho[], T &dp, const bool hor) { |
---|
| 1997 | int i; |
---|
| 1998 | double ky, kz_n; |
---|
| 1999 | T cx, sx, sz, chy, shy, cz; |
---|
| 2000 | |
---|
| 2001 | AoBrho = 0.0; |
---|
| 2002 | dp = 0.0; |
---|
| 2003 | |
---|
| 2004 | for (i = 0; i < 3; i++) |
---|
| 2005 | dAoBrho[i] = 0.0; |
---|
| 2006 | |
---|
| 2007 | for (i = 0; i < W->n_harm; i++) { |
---|
| 2008 | kz_n = W->harm[i] * 2.0 * M_PI / W->lambda; |
---|
| 2009 | ky = sqrt(sqr(W->kxV[i]) + sqr(kz_n)); |
---|
| 2010 | |
---|
| 2011 | cx = cos(W->kxV[i] * x[x_]); |
---|
| 2012 | sx = sin(W->kxV[i] * x[x_]); |
---|
| 2013 | chy = cosh(ky * x[y_]); |
---|
| 2014 | shy = sinh(ky * x[y_]); |
---|
| 2015 | sz = sin(kz_n * z); |
---|
| 2016 | |
---|
| 2017 | if (hor) { |
---|
| 2018 | // A_x/Brho |
---|
| 2019 | AoBrho += W->BoBrhoV[i] / kz_n * cx * chy * sz; |
---|
| 2020 | |
---|
| 2021 | if (globval.radiation) { |
---|
| 2022 | cz = cos(kz_n * z); |
---|
| 2023 | dAoBrho[X_] -= W->BoBrhoV[i] * W->kxV[i] / kz_n * sx * chy * sz; |
---|
| 2024 | dAoBrho[Y_] += W->BoBrhoV[i] * ky / kz_n * cx * shy * sz; |
---|
| 2025 | dAoBrho[Z_] += W->BoBrhoV[i] * cx * chy * cz; |
---|
| 2026 | } |
---|
| 2027 | |
---|
| 2028 | // dp_y |
---|
| 2029 | if (W->kxV[i] == 0.0) |
---|
| 2030 | dp += W->BoBrhoV[i] / kz_n * ky * x[x_] * shy * sz; |
---|
| 2031 | else |
---|
| 2032 | dp += W->BoBrhoV[i] / (W->kxV[i] * kz_n) * ky * sx * shy * sz; |
---|
| 2033 | } else { |
---|
| 2034 | // A_y/Brho |
---|
| 2035 | AoBrho += W->BoBrhoV[i] * W->kxV[i] / (ky * kz_n) * sx * shy * sz; |
---|
| 2036 | |
---|
| 2037 | if (globval.radiation) { |
---|
| 2038 | cz = cos(kz_n * z); |
---|
| 2039 | dAoBrho[X_] += W->BoBrhoV[i] * sqr(W->kxV[i]) / (ky * kz_n) |
---|
| 2040 | * cx * shy * sz; |
---|
| 2041 | dAoBrho[Y_] += W->BoBrhoV[i] * W->kxV[i] / kz_n * sx * chy * sz; |
---|
| 2042 | dAoBrho[Z_] += W->BoBrhoV[i] * W->kxV[i] / ky * sx * shy * cz; |
---|
| 2043 | } |
---|
| 2044 | |
---|
| 2045 | // dp_x |
---|
| 2046 | dp += W->BoBrhoV[i] / kz_n * sqr(W->kxV[i] / ky) * cx * chy * sz; |
---|
| 2047 | } |
---|
| 2048 | } |
---|
| 2049 | } |
---|
| 2050 | |
---|
| 2051 | template<typename T> |
---|
| 2052 | void Wiggler_pass_EF3(const elemtype &elem, ss_vect<T> &x) { |
---|
| 2053 | /* Second order symplectic integrator for insertion devices based on: |
---|
| 2054 | |
---|
| 2055 | E. Forest, et al "Explicit Symplectic Integrator for s-dependent |
---|
| 2056 | Static Magnetic Field" */ |
---|
| 2057 | |
---|
| 2058 | int i; |
---|
| 2059 | double h, z; |
---|
| 2060 | T hd, AxoBrho, AyoBrho, dAxoBrho[3], dAyoBrho[3], dpy, dpx, B[3]; |
---|
| 2061 | |
---|
| 2062 | h = elem.PL / elem.W->PN; |
---|
| 2063 | z = 0.0; |
---|
| 2064 | |
---|
| 2065 | for (i = 1; i <= elem.W->PN; i++) { |
---|
| 2066 | hd = h / (1.0 + x[delta_]); |
---|
| 2067 | |
---|
| 2068 | // 1: half step in z |
---|
| 2069 | z += 0.5 * h; |
---|
| 2070 | |
---|
| 2071 | // 2: half drift in y |
---|
| 2072 | get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2073 | |
---|
| 2074 | x[px_] -= dpx; |
---|
| 2075 | x[py_] -= AyoBrho; |
---|
| 2076 | x[y_] += 0.5 * hd * x[py_]; |
---|
| 2077 | x[ct_] += sqr(0.5) * hd * sqr(x[py_]) / (1.0 + x[delta_]); |
---|
| 2078 | |
---|
| 2079 | get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2080 | |
---|
| 2081 | x[px_] += dpx; |
---|
| 2082 | x[py_] += AyoBrho; |
---|
| 2083 | |
---|
| 2084 | // 3: full drift in x |
---|
| 2085 | get_Axy_EF3(elem.W, z, x, AxoBrho, dAxoBrho, dpy, true); |
---|
| 2086 | |
---|
| 2087 | x[px_] -= AxoBrho; |
---|
| 2088 | x[py_] -= dpy; |
---|
| 2089 | x[x_] += hd * x[px_]; |
---|
| 2090 | x[ct_] += 0.5 * hd * sqr(x[px_]) / (1.0 + x[delta_]); |
---|
| 2091 | |
---|
| 2092 | if (globval.pathlength) |
---|
| 2093 | x[ct_] += h; |
---|
| 2094 | |
---|
| 2095 | get_Axy_EF3(elem.W, z, x, AxoBrho, dAxoBrho, dpy, true); |
---|
| 2096 | |
---|
| 2097 | x[px_] += AxoBrho; |
---|
| 2098 | x[py_] += dpy; |
---|
| 2099 | |
---|
| 2100 | // 4: a half drift in y |
---|
| 2101 | get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2102 | |
---|
| 2103 | x[px_] -= dpx; |
---|
| 2104 | x[py_] -= AyoBrho; |
---|
| 2105 | x[y_] += 0.5 * hd * x[py_]; |
---|
| 2106 | x[ct_] += sqr(0.5) * hd * sqr(x[py_]) / (1.0 + x[delta_]); |
---|
| 2107 | |
---|
| 2108 | get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2109 | |
---|
| 2110 | x[px_] += dpx; |
---|
| 2111 | x[py_] += AyoBrho; |
---|
| 2112 | |
---|
| 2113 | // 5: half step in z |
---|
| 2114 | z += 0.5 * h; |
---|
| 2115 | |
---|
| 2116 | if (globval.radiation || globval.emittance) { |
---|
| 2117 | get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2118 | get_Axy_EF3(elem.W, z, x, AxoBrho, dAxoBrho, dpy, true); |
---|
| 2119 | B[X_] = -dAyoBrho[Z_]; |
---|
| 2120 | B[Y_] = dAxoBrho[Z_]; |
---|
| 2121 | B[Z_] = dAyoBrho[X_] - dAxoBrho[Y_]; |
---|
| 2122 | radiate(x, h, 0.0, B); |
---|
| 2123 | } |
---|
| 2124 | } |
---|
| 2125 | } |
---|
| 2126 | |
---|
| 2127 | template<typename T> |
---|
| 2128 | void Wiggler_Pass(CellType &Cell, ss_vect<T> &X) { |
---|
| 2129 | int seg; |
---|
| 2130 | double L, L1, L2, K1, K2; |
---|
| 2131 | elemtype *elemp; |
---|
| 2132 | WigglerType *W; |
---|
| 2133 | ss_vect<T> X1; |
---|
| 2134 | |
---|
| 2135 | elemp = &Cell.Elem; |
---|
| 2136 | W = elemp->W; |
---|
| 2137 | /* Global -> Local */ |
---|
| 2138 | GtoL(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 2139 | switch (W->Pmethod) { |
---|
| 2140 | |
---|
| 2141 | case Meth_Linear: |
---|
| 2142 | // LinTrans(5L, W->W55, X); |
---|
| 2143 | cout << "Wiggler_Pass: Meth_Linear not supported" << endl; |
---|
| 2144 | exit_(1); |
---|
| 2145 | break; |
---|
| 2146 | |
---|
| 2147 | case Meth_First: |
---|
| 2148 | if ((W->BoBrhoV[0] != 0.0) || (W->BoBrhoH[0] != 0.0)) { |
---|
| 2149 | if (!globval.EPU) |
---|
| 2150 | Wiggler_pass_EF(Cell.Elem, X); |
---|
| 2151 | else { |
---|
| 2152 | Wiggler_pass_EF2(W->PN, elemp->PL, W->kxV[0], W->kxH[0], 2.0 |
---|
| 2153 | * M_PI / W->lambda, W->BoBrhoV[0], W->BoBrhoH[0], |
---|
| 2154 | W->phi[0], X); |
---|
| 2155 | } |
---|
| 2156 | } else |
---|
| 2157 | // drift if field = 0 |
---|
| 2158 | Drift(elemp->PL, X); |
---|
| 2159 | break; |
---|
| 2160 | |
---|
| 2161 | case Meth_Second: |
---|
| 2162 | if ((W->BoBrhoV[0] != 0.0) || (W->BoBrhoH[0] != 0.0)) { |
---|
| 2163 | Wiggler_pass_EF3(Cell.Elem, X); |
---|
| 2164 | } else |
---|
| 2165 | // drift if field = 0 |
---|
| 2166 | Drift(elemp->PL, X); |
---|
| 2167 | break; |
---|
| 2168 | |
---|
| 2169 | case Meth_Fourth: /* 4-th order integrator */ |
---|
| 2170 | L = elemp->PL / W->PN; |
---|
| 2171 | L1 = c_1 * L; |
---|
| 2172 | L2 = c_2 * L; |
---|
| 2173 | K1 = d_1 * L; |
---|
| 2174 | K2 = d_2 * L; |
---|
| 2175 | for (seg = 1; seg <= W->PN; seg++) { |
---|
| 2176 | Drift(L1, X); |
---|
| 2177 | X1 = X; |
---|
| 2178 | thin_kick(W->Porder, W->PBW, K1, 0.0, 0.0, X1); |
---|
| 2179 | X[py_] = X1[py_]; |
---|
| 2180 | Drift(L2, X); |
---|
| 2181 | X1 = X; |
---|
| 2182 | thin_kick(W->Porder, W->PBW, K2, 0.0, 0.0, X1); |
---|
| 2183 | X[py_] = X1[py_]; |
---|
| 2184 | Drift(L2, X); |
---|
| 2185 | X1 = X; |
---|
| 2186 | thin_kick(W->Porder, W->PBW, K1, 0.0, 0.0, X1); |
---|
| 2187 | X[py_] = X1[py_]; |
---|
| 2188 | Drift(L1, X); |
---|
| 2189 | } |
---|
| 2190 | break; |
---|
| 2191 | } |
---|
| 2192 | /* Local -> Global */ |
---|
| 2193 | LtoG(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 2194 | } |
---|
| 2195 | |
---|
| 2196 | #undef eps |
---|
| 2197 | #undef kx |
---|
| 2198 | |
---|
| 2199 | template<typename T> |
---|
| 2200 | void FieldMap_Pass(CellType &Cell, ss_vect<T> &X) { |
---|
| 2201 | |
---|
| 2202 | GtoL(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 2203 | |
---|
| 2204 | Wiggler_pass_EF(Cell.Elem, X); |
---|
| 2205 | |
---|
| 2206 | LtoG(X, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 2207 | } |
---|
| 2208 | |
---|
| 2209 | /******************************************************************** |
---|
| 2210 | void Insertion_Pass(CellType &Cell, ss_vect<T> &x) |
---|
| 2211 | |
---|
| 2212 | Purpose: |
---|
| 2213 | Track vector x through an insertion |
---|
| 2214 | If radiation or cavity on insertion is like a drift |
---|
| 2215 | |
---|
| 2216 | Input: |
---|
| 2217 | Cell element to track through |
---|
| 2218 | x initial coordinates vector |
---|
| 2219 | |
---|
| 2220 | Output: |
---|
| 2221 | x final coordinates vector |
---|
| 2222 | |
---|
| 2223 | Return: |
---|
| 2224 | none |
---|
| 2225 | |
---|
| 2226 | Global variables: |
---|
| 2227 | none |
---|
| 2228 | |
---|
| 2229 | Specific functions: |
---|
| 2230 | LinearInterpolation2 |
---|
| 2231 | Drft |
---|
| 2232 | CopyVec |
---|
| 2233 | |
---|
| 2234 | Comments: |
---|
| 2235 | Outside of interpolation table simulated by putting 1 in x[4] |
---|
| 2236 | 01/07/03 6D tracking activated |
---|
| 2237 | 10/01/05 First order kick part added |
---|
| 2238 | *******************************************************************/ |
---|
| 2239 | |
---|
| 2240 | template<typename T> |
---|
| 2241 | void Insertion_Pass(CellType &Cell, ss_vect<T> &x) { |
---|
| 2242 | |
---|
| 2243 | elemtype *elemp; |
---|
| 2244 | double LN = 0.0; |
---|
| 2245 | T tx1, tz1; /* thetax and thetaz retrieved from |
---|
| 2246 | interpolation routine First order kick*/ |
---|
| 2247 | T tx2, tz2; /* thetax and thetaz retrieved from |
---|
| 2248 | interpolation routine Second order Kick */ |
---|
| 2249 | T d; |
---|
| 2250 | double alpha0 = 0.0; // 1/ brh0 |
---|
| 2251 | double alpha02 = 0.0; // alpha square |
---|
| 2252 | int Nslice = 0; |
---|
| 2253 | int i = 0; |
---|
| 2254 | bool outoftable = false; |
---|
| 2255 | |
---|
| 2256 | elemp = &Cell.Elem; |
---|
| 2257 | Nslice = elemp->ID->PN; |
---|
| 2258 | alpha0 = c0 / globval.Energy * 1E-9 * elemp->ID->scaling1; |
---|
| 2259 | alpha02 = (c0 / globval.Energy * 1E-9) * (c0 / globval.Energy * 1E-9) |
---|
| 2260 | * elemp->ID->scaling2; |
---|
| 2261 | |
---|
| 2262 | // /* Global -> Local */ |
---|
| 2263 | // GtoL(X, Cell->dS, Cell->dT, 0.0, 0.0, 0.0); |
---|
| 2264 | |
---|
| 2265 | // (Nslice+1) drifts, n slice kicks |
---|
| 2266 | LN = elemp->PL / (Nslice + 1); |
---|
| 2267 | Drift(LN, x); |
---|
| 2268 | |
---|
| 2269 | for (i = 1; i <= Nslice; i++) { |
---|
| 2270 | // First order kick map |
---|
| 2271 | if (elemp->ID->firstorder) { |
---|
| 2272 | if (!elemp->ID->linear) { |
---|
| 2273 | //cout << "InsertionPass: Spline\n"; |
---|
| 2274 | SplineInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable, |
---|
| 2275 | 1); |
---|
| 2276 | } else{ |
---|
| 2277 | LinearInterpolation2(x[x_], x[y_], tx1, tz1, Cell, outoftable, |
---|
| 2278 | 1); |
---|
| 2279 | } |
---|
| 2280 | if (outoftable) { |
---|
| 2281 | x[x_] = 1e30; |
---|
| 2282 | return; |
---|
| 2283 | } |
---|
| 2284 | |
---|
| 2285 | d = alpha0 / Nslice; |
---|
| 2286 | x[px_] += d * tx1; |
---|
| 2287 | x[py_] += d * tz1; |
---|
| 2288 | } |
---|
| 2289 | |
---|
| 2290 | // Second order kick map |
---|
| 2291 | if (elemp->ID->secondorder) { |
---|
| 2292 | if (!elemp->ID->linear) { |
---|
| 2293 | //cout << "InsertionPass: Spline\n"; |
---|
| 2294 | SplineInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable, |
---|
| 2295 | 2); |
---|
| 2296 | } else { |
---|
| 2297 | //cout << "InsertionPass: Linear\n"; |
---|
| 2298 | LinearInterpolation2(x[x_], x[y_], tx2, tz2, Cell, outoftable, |
---|
| 2299 | 2); |
---|
| 2300 | } |
---|
| 2301 | if (outoftable) { |
---|
| 2302 | x[x_] = 1e30; |
---|
| 2303 | return; |
---|
| 2304 | } |
---|
| 2305 | |
---|
| 2306 | d = alpha02 / Nslice / (1.0 + x[delta_]); |
---|
| 2307 | x[px_] += d * tx2; |
---|
| 2308 | x[py_] += d * tz2; |
---|
| 2309 | } |
---|
| 2310 | Drift(LN, x); |
---|
| 2311 | } |
---|
| 2312 | // CopyVec(6L, x, Cell->BeamPos); |
---|
| 2313 | |
---|
| 2314 | // /* Local -> Global */ |
---|
| 2315 | // LtoG(X, Cell->dS, Cell->dT, 0.0, 0.0, 0.0); |
---|
| 2316 | } |
---|
| 2317 | |
---|
| 2318 | template<typename T> |
---|
| 2319 | void sol_pass(const elemtype &elem, ss_vect<T> &x) { |
---|
| 2320 | int i; |
---|
| 2321 | double h, z; |
---|
| 2322 | T hd, AxoBrho, AyoBrho, dAxoBrho[3], dAyoBrho[3], dpy, dpx, B[3]; |
---|
| 2323 | |
---|
| 2324 | h = elem.PL / elem.Sol->N; |
---|
| 2325 | z = 0.0; |
---|
| 2326 | |
---|
| 2327 | for (i = 1; i <= elem.Sol->N; i++) { |
---|
| 2328 | hd = h / (1.0 + x[delta_]); |
---|
| 2329 | |
---|
| 2330 | // 1: half step in z |
---|
| 2331 | z += 0.5 * h; |
---|
| 2332 | |
---|
| 2333 | // 2: half drift in y |
---|
| 2334 | AyoBrho = elem.Sol->BoBrho * x[x_] / 2.0; |
---|
| 2335 | dpx = elem.Sol->BoBrho * x[y_] / 2.0; |
---|
| 2336 | // get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2337 | |
---|
| 2338 | x[px_] -= dpx; |
---|
| 2339 | x[py_] -= AyoBrho; |
---|
| 2340 | x[y_] += 0.5 * hd * x[py_]; |
---|
| 2341 | x[ct_] += sqr(0.5) * hd * sqr(x[py_]) / (1.0 + x[delta_]); |
---|
| 2342 | |
---|
| 2343 | AyoBrho = elem.Sol->BoBrho * x[x_] / 2.0; |
---|
| 2344 | dpx = elem.Sol->BoBrho * x[y_] / 2.0; |
---|
| 2345 | // get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2346 | |
---|
| 2347 | x[px_] += dpx; |
---|
| 2348 | x[py_] += AyoBrho; |
---|
| 2349 | |
---|
| 2350 | // 3: full drift in x |
---|
| 2351 | AxoBrho = -elem.Sol->BoBrho * x[y_] / 2.0; |
---|
| 2352 | dpy = -elem.Sol->BoBrho * x[x_] / 2.0; |
---|
| 2353 | // get_Axy_EF3(elem.W, z, x, AxoBrho, dAxoBrho, dpy, true); |
---|
| 2354 | |
---|
| 2355 | x[px_] -= AxoBrho; |
---|
| 2356 | x[py_] -= dpy; |
---|
| 2357 | x[x_] += hd * x[px_]; |
---|
| 2358 | x[ct_] += 0.5 * hd * sqr(x[px_]) / (1.0 + x[delta_]); |
---|
| 2359 | |
---|
| 2360 | if (globval.pathlength) |
---|
| 2361 | x[ct_] += h; |
---|
| 2362 | |
---|
| 2363 | AxoBrho = -elem.Sol->BoBrho * x[y_] / 2.0; |
---|
| 2364 | dpy = -elem.Sol->BoBrho * x[x_] / 2.0; |
---|
| 2365 | // get_Axy_EF3(elem.W, z, x, AxoBrho, dAxoBrho, dpy, true); |
---|
| 2366 | |
---|
| 2367 | x[px_] += AxoBrho; |
---|
| 2368 | x[py_] += dpy; |
---|
| 2369 | |
---|
| 2370 | // 4: a half drift in y |
---|
| 2371 | AyoBrho = elem.Sol->BoBrho * x[x_] / 2.0; |
---|
| 2372 | dpx = elem.Sol->BoBrho * x[y_] / 2.0; |
---|
| 2373 | // get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2374 | |
---|
| 2375 | x[px_] -= dpx; |
---|
| 2376 | x[py_] -= AyoBrho; |
---|
| 2377 | x[y_] += 0.5 * hd * x[py_]; |
---|
| 2378 | x[ct_] += sqr(0.5) * hd * sqr(x[py_]) / (1.0 + x[delta_]); |
---|
| 2379 | |
---|
| 2380 | AyoBrho = elem.Sol->BoBrho * x[x_] / 2.0; |
---|
| 2381 | dpx = elem.Sol->BoBrho * x[y_] / 2.0; |
---|
| 2382 | // get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2383 | |
---|
| 2384 | x[px_] += dpx; |
---|
| 2385 | x[py_] += AyoBrho; |
---|
| 2386 | |
---|
| 2387 | // 5: half step in z |
---|
| 2388 | z += 0.5 * h; |
---|
| 2389 | |
---|
| 2390 | if (globval.radiation || globval.emittance) { |
---|
| 2391 | dAxoBrho[X_] = 0.0; |
---|
| 2392 | dAxoBrho[Y_] = -elem.Sol->BoBrho / 2.0; |
---|
| 2393 | dAxoBrho[Z_] = 0.0; |
---|
| 2394 | dAyoBrho[X_] = elem.Sol->BoBrho / 2.0; |
---|
| 2395 | dAyoBrho[Y_] = 0.0; |
---|
| 2396 | dAyoBrho[Z_] = 0.0; |
---|
| 2397 | // get_Axy_EF3(elem.W, z, x, AyoBrho, dAyoBrho, dpx, false); |
---|
| 2398 | // get_Axy_EF3(elem.W, z, x, AxoBrho, dAxoBrho, dpy, true); |
---|
| 2399 | B[X_] = -dAyoBrho[Z_]; |
---|
| 2400 | B[Y_] = dAxoBrho[Z_]; |
---|
| 2401 | B[Z_] = dAyoBrho[X_] - dAxoBrho[Y_]; |
---|
| 2402 | radiate(x, h, 0.0, B); |
---|
| 2403 | } |
---|
| 2404 | } |
---|
| 2405 | } |
---|
| 2406 | |
---|
| 2407 | template<typename T> |
---|
| 2408 | void Solenoid_Pass(CellType &Cell, ss_vect<T> &ps) { |
---|
| 2409 | |
---|
| 2410 | GtoL(ps, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 2411 | |
---|
| 2412 | sol_pass(Cell.Elem, ps); |
---|
| 2413 | |
---|
| 2414 | LtoG(ps, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 2415 | } |
---|
| 2416 | |
---|
| 2417 | // Matrix model |
---|
| 2418 | |
---|
| 2419 | void GtoL_M(Matrix &X, Vector2 &T) { |
---|
| 2420 | Matrix R; |
---|
| 2421 | |
---|
| 2422 | /* Rotate */ |
---|
| 2423 | R[0][0] = T[0]; R[0][1] = 0.0; R[0][2] = T[1]; R[0][3] = 0.0; |
---|
| 2424 | R[1][0] = 0.0; R[1][1] = T[0]; R[1][2] = 0.0; R[1][3] = T[1]; |
---|
| 2425 | R[2][0] = -T[1]; R[2][1] = 0.0; R[2][2] = T[0]; R[2][3] = 0.0; |
---|
| 2426 | R[3][0] = 0.0; R[3][1] = -T[1]; R[3][2] = 0.0; R[3][3] = T[0]; |
---|
| 2427 | MulLMat(4L, R, X); |
---|
| 2428 | } |
---|
| 2429 | |
---|
| 2430 | void LtoG_M(Matrix &X, Vector2 &T) { |
---|
| 2431 | Matrix R; |
---|
| 2432 | |
---|
| 2433 | /* Rotate */ |
---|
| 2434 | R[0][0] = T[0]; R[0][1] = 0.0; R[0][2] = -T[1]; R[0][3] = 0.0; |
---|
| 2435 | R[1][0] = 0.0; R[1][1] = T[0]; R[1][2] = 0.0; R[1][3] = -T[1]; |
---|
| 2436 | R[2][0] = T[1]; R[2][1] = 0.0; R[2][2] = T[0]; R[2][3] = 0.0; |
---|
| 2437 | R[3][0] = 0.0; R[3][1] = T[1]; R[3][2] = 0.0; R[3][3] = T[0]; |
---|
| 2438 | MulLMat(4, R, X); |
---|
| 2439 | } |
---|
| 2440 | |
---|
| 2441 | /**************************************************************************** |
---|
| 2442 | * void Drift_Pass_M(CellType *Cell, double *xref, vector *x) |
---|
| 2443 | |
---|
| 2444 | Purpose: |
---|
| 2445 | matrix propagation through a drift |
---|
| 2446 | x = D55*x |
---|
| 2447 | xref= drift(xref) |
---|
| 2448 | |
---|
| 2449 | Input: |
---|
| 2450 | xref vector |
---|
| 2451 | x matrix |
---|
| 2452 | |
---|
| 2453 | Output: |
---|
| 2454 | xref |
---|
| 2455 | x |
---|
| 2456 | |
---|
| 2457 | Return: |
---|
| 2458 | none |
---|
| 2459 | |
---|
| 2460 | Global variables: |
---|
| 2461 | none |
---|
| 2462 | |
---|
| 2463 | Specific functions: |
---|
| 2464 | MulLMat, Drft |
---|
| 2465 | |
---|
| 2466 | Comments: |
---|
| 2467 | none |
---|
| 2468 | |
---|
| 2469 | ****************************************************************************/ |
---|
| 2470 | |
---|
| 2471 | void Drift_Pass_M(CellType &Cell, Vector &xref, Matrix &X) { |
---|
| 2472 | |
---|
| 2473 | MulLMat(5, Cell.Elem.D->D55, X); |
---|
| 2474 | Drift(Cell.Elem.PL, xref); |
---|
| 2475 | } |
---|
| 2476 | |
---|
| 2477 | /****************************************************************************/ |
---|
| 2478 | /* void thinkick_M(long Order, double *MB, double L, double irho, pthicktype pthick, |
---|
| 2479 | double *xref, vector *x) |
---|
| 2480 | |
---|
| 2481 | Purpose: |
---|
| 2482 | |
---|
| 2483 | |
---|
| 2484 | Input: |
---|
| 2485 | none |
---|
| 2486 | |
---|
| 2487 | Output: |
---|
| 2488 | none |
---|
| 2489 | |
---|
| 2490 | Return: |
---|
| 2491 | none |
---|
| 2492 | |
---|
| 2493 | Global variables: |
---|
| 2494 | none |
---|
| 2495 | |
---|
| 2496 | Specific functions: |
---|
| 2497 | thinkick |
---|
| 2498 | |
---|
| 2499 | Comments: |
---|
| 2500 | none |
---|
| 2501 | |
---|
| 2502 | ****************************************************************************/ |
---|
| 2503 | void thin_kick_M(int Order, double MB[], double L, double irho, Vector &xref, |
---|
| 2504 | Matrix &x) { |
---|
| 2505 | int i; |
---|
| 2506 | mpolArray MMB; |
---|
| 2507 | Vector z; |
---|
| 2508 | Matrix Mk; |
---|
| 2509 | |
---|
| 2510 | if (2 > Order || Order > HOMmax) |
---|
| 2511 | return; |
---|
| 2512 | for (i = 2; i <= Order; i++) { |
---|
| 2513 | MMB[i + HOMmax - 1] = (i - 1) * MB[i + HOMmax]; |
---|
| 2514 | MMB[HOMmax - i + 1] = (i - 1) * MB[HOMmax - i]; |
---|
| 2515 | } |
---|
| 2516 | z[0] = xref[0]; |
---|
| 2517 | z[1] = 0.0; |
---|
| 2518 | z[2] = xref[2]; |
---|
| 2519 | z[3] = 0.0; |
---|
| 2520 | z[4] = 0.0; |
---|
| 2521 | z[5] = 0.0; |
---|
| 2522 | thin_kick(Order - 1, MMB, L, 0.0, 0.0, z); |
---|
| 2523 | z[1] -= L * sqr(irho); |
---|
| 2524 | UnitMat(5L, Mk); |
---|
| 2525 | Mk[1][0] = z[1]; |
---|
| 2526 | Mk[1][2] = z[3]; |
---|
| 2527 | Mk[3][0] = z[3]; |
---|
| 2528 | Mk[3][2] = -z[1]; |
---|
| 2529 | MulLMat(5L, Mk, x); |
---|
| 2530 | } |
---|
| 2531 | |
---|
| 2532 | static void make3by3(Matrix &A, double a11, double a12, double a13, double a21, |
---|
| 2533 | double a22, double a23, double a31, double a32, double a33) { |
---|
| 2534 | /* |
---|
| 2535 | Set the 3x3 matrix A to: |
---|
| 2536 | (a11 a12 a13) |
---|
| 2537 | A = (a21 a22 a23) |
---|
| 2538 | (a31 a32 a33) |
---|
| 2539 | */ |
---|
| 2540 | |
---|
| 2541 | UnitMat(ss_dim, A); /* set matrix to unit 3x3 matrix */ |
---|
| 2542 | A[0][0] = a11; A[0][1] = a12; A[0][2] = a13; |
---|
| 2543 | A[1][0] = a21; A[1][1] = a22; A[1][2] = a23; |
---|
| 2544 | A[2][0] = a31; A[2][1] = a32; A[2][2] = a33; |
---|
| 2545 | } |
---|
| 2546 | |
---|
| 2547 | /**************************************************************************** |
---|
| 2548 | * void make4by5(vector *A, double a11, double a12, double a15, |
---|
| 2549 | double a21, double a22, double a25, double a33, |
---|
| 2550 | double a34, double a35, double a43, double a44, |
---|
| 2551 | double a45) |
---|
| 2552 | Purpose: |
---|
| 2553 | Constructor for matrix A |
---|
| 2554 | All the matrix terms are explicitly given as input |
---|
| 2555 | |
---|
| 2556 | Input: |
---|
| 2557 | A matrix to initialize |
---|
| 2558 | aij matrix terms |
---|
| 2559 | |
---|
| 2560 | Output: |
---|
| 2561 | A initialized matrix |
---|
| 2562 | |
---|
| 2563 | Return: |
---|
| 2564 | none |
---|
| 2565 | |
---|
| 2566 | Global variables: |
---|
| 2567 | none |
---|
| 2568 | |
---|
| 2569 | Specific functions: |
---|
| 2570 | UnitMat |
---|
| 2571 | |
---|
| 2572 | Comments: |
---|
| 2573 | none |
---|
| 2574 | |
---|
| 2575 | ****************************************************************************/ |
---|
| 2576 | |
---|
| 2577 | static void make4by5(Matrix &A, double a11, double a12, double a15, double a21, |
---|
| 2578 | double a22, double a25, double a33, double a34, double a35, double a43, |
---|
| 2579 | double a44, double a45) { |
---|
| 2580 | UnitMat(ss_dim, A); /* Initializes A to identity matrix */ |
---|
| 2581 | A[0][0] = a11; A[0][1] = a12; A[0][4] = a15; |
---|
| 2582 | A[1][0] = a21; A[1][1] = a22; A[1][4] = a25; |
---|
| 2583 | A[2][2] = a33; A[2][3] = a34; A[2][4] = a35; |
---|
| 2584 | A[3][2] = a43; A[3][3] = a44; A[3][4] = a45; |
---|
| 2585 | } |
---|
| 2586 | |
---|
| 2587 | |
---|
| 2588 | static void mergeto4by5(Matrix &A, Matrix &AH, Matrix &AV) { |
---|
| 2589 | /* |
---|
| 2590 | merges two 3x3 matrices AH (H-plane) and AV (V-plane) into one |
---|
| 2591 | big 4x5 matrix |
---|
| 2592 | |
---|
| 2593 | (ah11 ah12 0 0 ah13) |
---|
| 2594 | (ah21 ah22 0 0 ah23) |
---|
| 2595 | A= ( 0 0 av11 av12 av13) |
---|
| 2596 | ( 0 0 av21 av22 ah13) |
---|
| 2597 | ( 0 0 0 0 1) |
---|
| 2598 | */ |
---|
| 2599 | int i, j; |
---|
| 2600 | |
---|
| 2601 | UnitMat(ss_dim, A); |
---|
| 2602 | for (i = 1; i <= 2; i++) { |
---|
| 2603 | A[i - 1][4] = AH[i - 1][2]; |
---|
| 2604 | A[i + 1][4] = AV[i - 1][2]; |
---|
| 2605 | for (j = 1; j <= 2; j++) { |
---|
| 2606 | A[i - 1][j - 1] = AH[i - 1][j - 1]; |
---|
| 2607 | A[i + 1][j + 1] = AV[i - 1][j - 1]; |
---|
| 2608 | } |
---|
| 2609 | } |
---|
| 2610 | } |
---|
| 2611 | |
---|
| 2612 | void Drift_SetMatrix(int Fnum1, int Knum1) { |
---|
| 2613 | /* |
---|
| 2614 | Make transport matrix for drift from |
---|
| 2615 | familiy Fnum1 and Kid number Knum |
---|
| 2616 | L = L / (1 + dP) |
---|
| 2617 | |
---|
| 2618 | ( 1 L 0 0 0) |
---|
| 2619 | D55 = ( 0 1 0 0 0) |
---|
| 2620 | ( 0 0 1 L 0) |
---|
| 2621 | ( 0 0 0 1 0) |
---|
| 2622 | */ |
---|
| 2623 | |
---|
| 2624 | double L; |
---|
| 2625 | CellType *cellp; |
---|
| 2626 | elemtype *elemp; |
---|
| 2627 | DriftType *D; |
---|
| 2628 | |
---|
| 2629 | if (ElemFam[Fnum1 - 1].nKid <= 0) |
---|
| 2630 | return; |
---|
| 2631 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 2632 | elemp = &cellp->Elem; |
---|
| 2633 | D = elemp->D; |
---|
| 2634 | L = elemp->PL/(1.0+globval.dPparticle); /* L = L / (1 + dP) */ |
---|
| 2635 | make4by5(D->D55, |
---|
| 2636 | 1.0, L, 0.0, 0.0, 1.0, 0.0, |
---|
| 2637 | 1.0, L, 0.0, 0.0, 1.0, 0.0); |
---|
| 2638 | } |
---|
| 2639 | |
---|
| 2640 | static void driftmat(Matrix &ah, double L) { |
---|
| 2641 | L /= 1 + globval.dPparticle; |
---|
| 2642 | make4by5(ah, |
---|
| 2643 | 1.0, L, 0.0, 0.0, 1.0, 0.0, |
---|
| 2644 | 1.0, L, 0.0, 0.0, 1.0, 0.0); |
---|
| 2645 | } |
---|
| 2646 | |
---|
| 2647 | static void quadmat(Matrix &ahv, double L, double k) { |
---|
| 2648 | /* |
---|
| 2649 | creates the avh matrix for a quadrupole |
---|
| 2650 | where av and ah are the horizontal and vertical |
---|
| 2651 | focusing or defocusing matrices |
---|
| 2652 | |
---|
| 2653 | |
---|
| 2654 | 1/2 1/2 |
---|
| 2655 | cos(L* (|K|) ) sin(L* (|K|) ) |
---|
| 2656 | c = --------------- s = --------------- |
---|
| 2657 | 1/2 1/2 |
---|
| 2658 | (1 + Dp) (1 + Dp) |
---|
| 2659 | 1/2 |
---|
| 2660 | sk = (|K|(1+dP)) |
---|
| 2661 | |
---|
| 2662 | - if k > 0 |
---|
| 2663 | H plane V plane |
---|
| 2664 | |
---|
| 2665 | ( c s/sk 0 ) ( ch sh/k 0 ) |
---|
| 2666 | ah = ( sk*s c 0 ) av = ( sk*sh ch 0 ) |
---|
| 2667 | ( 0 0 1 ) ( 0 0 1 ) |
---|
| 2668 | |
---|
| 2669 | |
---|
| 2670 | ( ah11 ah12 0 0 ah13 ) |
---|
| 2671 | ( ah21 ah22 0 0 ah23 ) |
---|
| 2672 | avh = ( 0 0 av11 av12 av13 ) |
---|
| 2673 | ( 0 0 av21 av22 ah13 ) |
---|
| 2674 | ( 0 0 0 0 1 ) */ |
---|
| 2675 | |
---|
| 2676 | double t, sk, sk0, s, c; |
---|
| 2677 | Matrix a, ah, av; |
---|
| 2678 | |
---|
| 2679 | if (k == 0.0) { |
---|
| 2680 | /* pure drift focusing */ |
---|
| 2681 | driftmat(ahv, L); |
---|
| 2682 | return; |
---|
| 2683 | } |
---|
| 2684 | sk0 = sqrt(fabs(k)); |
---|
| 2685 | t = L * sk0 / sqrt(1.0 + globval.dPparticle); |
---|
| 2686 | c = cos(t); |
---|
| 2687 | s = sin(t); |
---|
| 2688 | sk = sk0 * sqrt(1.0 + globval.dPparticle); |
---|
| 2689 | make3by3(a, c, s / sk, 0.0, -sk * s, c, 0.0, 0.0, 0.0, 1.0); |
---|
| 2690 | if (k > 0.0) |
---|
| 2691 | CopyMat(3L, a, ah); |
---|
| 2692 | else |
---|
| 2693 | CopyMat(3L, a, av); |
---|
| 2694 | c = cosh(t); |
---|
| 2695 | s = sinh(t); |
---|
| 2696 | sk = sk0 * sqrt(1.0 + globval.dPparticle); |
---|
| 2697 | make3by3(a, c, s / sk, 0.0, sk * s, c, 0.0, 0.0, 0.0, 1.0); |
---|
| 2698 | if (k > 0.0) |
---|
| 2699 | CopyMat(3L, a, av); |
---|
| 2700 | else |
---|
| 2701 | CopyMat(3L, a, ah); |
---|
| 2702 | mergeto4by5(ahv, ah, av); |
---|
| 2703 | } |
---|
| 2704 | |
---|
| 2705 | /****************************************************************************/ |
---|
| 2706 | /* static void bendmat(vector *M, double L, double irho, double phi1, |
---|
| 2707 | double phi2, double gap, double k) |
---|
| 2708 | |
---|
| 2709 | Purpose: called by Mpole_Setmatrix |
---|
| 2710 | |
---|
| 2711 | For a quadrupole see quadmat routine for explanation |
---|
| 2712 | |
---|
| 2713 | For a dipole |
---|
| 2714 | |
---|
| 2715 | (1 0 0) |
---|
| 2716 | Edge(theta) = (h*tan(theta) 1 0) |
---|
| 2717 | (0 0 1) |
---|
| 2718 | |
---|
| 2719 | (1 0 0) |
---|
| 2720 | Edge(theta) = (-h*tan(theta-psi) 1 0) |
---|
| 2721 | (0 0 1) |
---|
| 2722 | |
---|
| 2723 | 2 |
---|
| 2724 | K1*gap*h*(1 + sin phi) |
---|
| 2725 | psi = -----------------------, K1 = 1/2 |
---|
| 2726 | cos phi |
---|
| 2727 | |
---|
| 2728 | Input: |
---|
| 2729 | L : length [m] |
---|
| 2730 | irho: 1/rho [1/m] |
---|
| 2731 | phi1: entrance edge angle [degres] |
---|
| 2732 | phi2: exit edge angle [degres] |
---|
| 2733 | K : gradient = n/Rho |
---|
| 2734 | |
---|
| 2735 | Output: |
---|
| 2736 | M transfer matrix |
---|
| 2737 | |
---|
| 2738 | Return: |
---|
| 2739 | none |
---|
| 2740 | |
---|
| 2741 | Global variables: |
---|
| 2742 | none |
---|
| 2743 | |
---|
| 2744 | Specific functions: |
---|
| 2745 | quadmat, make3by3 |
---|
| 2746 | UnitMat, MulRMat, psi |
---|
| 2747 | |
---|
| 2748 | Comments: |
---|
| 2749 | none |
---|
| 2750 | |
---|
| 2751 | ****************************************************************************/ |
---|
| 2752 | static void bendmat(Matrix &M, double L, double irho, double phi1, double phi2, |
---|
| 2753 | double gap, double k) { |
---|
| 2754 | /* called by Mpole_Setmatrix |
---|
| 2755 | |
---|
| 2756 | For a quadrupole see quadmat routine for explanation |
---|
| 2757 | |
---|
| 2758 | For a dipole |
---|
| 2759 | |
---|
| 2760 | (1 0 0) |
---|
| 2761 | Edge(theta) = (h*tan(theta) 1 0) |
---|
| 2762 | (0 0 1) |
---|
| 2763 | |
---|
| 2764 | (1 0 0) |
---|
| 2765 | Edge(theta) = (-h*tan(theta-psi) 1 0) |
---|
| 2766 | (0 0 1) |
---|
| 2767 | |
---|
| 2768 | 2 |
---|
| 2769 | K1*gap*h*(1 + sin phi) |
---|
| 2770 | psi = -----------------------, K1 = 1/2 |
---|
| 2771 | cos phi */ |
---|
| 2772 | |
---|
| 2773 | double r, s, c, sk, p, fk, afk; |
---|
| 2774 | Matrix edge, ah, av; |
---|
| 2775 | double coef, scoef; |
---|
| 2776 | |
---|
| 2777 | if (irho == 0.0) { |
---|
| 2778 | /* Quadrupole matrix */ |
---|
| 2779 | quadmat(M, L, k); |
---|
| 2780 | return; |
---|
| 2781 | } |
---|
| 2782 | |
---|
| 2783 | /* For multipole w/ irho !=0 eg dipole */ |
---|
| 2784 | coef = 1.0 + globval.dPparticle; |
---|
| 2785 | scoef = sqrt(coef); |
---|
| 2786 | r = L * irho / scoef; |
---|
| 2787 | |
---|
| 2788 | if (k == 0.0) { |
---|
| 2789 | /* simple vertical dipole magnet */ |
---|
| 2790 | /* simple vertical dipole magnet */ |
---|
| 2791 | /* |
---|
| 2792 | H-plane |
---|
| 2793 | 2 2 2 |
---|
| 2794 | px + py 2 x |
---|
| 2795 | H = -------- - h*x*dP + h --- |
---|
| 2796 | 2*(1+dP) 2 |
---|
| 2797 | |
---|
| 2798 | 2 2 |
---|
| 2799 | dx h h*dP |
---|
| 2800 | --2 + ---- x = ----- |
---|
| 2801 | ds 1+dP 1+dP |
---|
| 2802 | |
---|
| 2803 | |
---|
| 2804 | let be u = Lh/sqrt(1+dP) then the transfert matrix becomes: |
---|
| 2805 | |
---|
| 2806 | ( sin(u) 1- cos(u) ) |
---|
| 2807 | ( cos(u) -------------- ----------------- ) |
---|
| 2808 | ( h sqrt(1+dP) h ) |
---|
| 2809 | ( -sin(u)*sqrt(1+dP)*h cos(u) sin(u)*sqrt(1+dP) ) |
---|
| 2810 | ( 0 0 1 ) |
---|
| 2811 | |
---|
| 2812 | */ |
---|
| 2813 | c = cos(r); |
---|
| 2814 | s = sin(r); |
---|
| 2815 | make3by3(ah, c, s / (irho * scoef), (1.0 - c) / irho, |
---|
| 2816 | -s * scoef * irho, c, s * scoef, 0.0, 0.0, 1.0); |
---|
| 2817 | /* |
---|
| 2818 | V-plane: it is just a drift |
---|
| 2819 | ( L ) |
---|
| 2820 | ( 1 ---- 0) |
---|
| 2821 | ( 1+dP ) |
---|
| 2822 | ( 0 1 0) |
---|
| 2823 | ( 0 0 1) |
---|
| 2824 | */ |
---|
| 2825 | make3by3(av, 1.0, L / coef, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0); |
---|
| 2826 | } else { |
---|
| 2827 | /* gradient bend, k= n/rho^2 */ |
---|
| 2828 | /* |
---|
| 2829 | K = -k -h*h |
---|
| 2830 | p = L*sqrt(|K|)/sqrt(1+dP) |
---|
| 2831 | */ |
---|
| 2832 | fk = -k - irho * irho; |
---|
| 2833 | afk = fabs(fk); |
---|
| 2834 | sk = sqrt(afk); |
---|
| 2835 | p = L * sk / scoef; |
---|
| 2836 | if (fk < 0.0) { |
---|
| 2837 | /* |
---|
| 2838 | H-plane |
---|
| 2839 | 2 2 2 |
---|
| 2840 | px + py 2 x |
---|
| 2841 | H = -------- - h*x*dP + (k+h ) --- |
---|
| 2842 | 2*(1+dP) 2 |
---|
| 2843 | |
---|
| 2844 | 2 2 |
---|
| 2845 | dx k+h h*dP |
---|
| 2846 | --2 + ---- x = ----- |
---|
| 2847 | ds 1+dP 1+dP |
---|
| 2848 | |
---|
| 2849 | |
---|
| 2850 | let be u = Lsqrt(|h*h+b2|)/sqrt(1+dP) |
---|
| 2851 | then the transfert matrix becomes: |
---|
| 2852 | |
---|
| 2853 | ( sin(u) 1- cos(u) ) |
---|
| 2854 | ( cos(u) -------------- ----------------- ) |
---|
| 2855 | ( sk*sqrt(1+dP) |k+h*h|*h ) |
---|
| 2856 | ( -sin(u)*sqrt(1+dP)*sk cos(u) h*sin(u)*sqrt(1+dP)/sk) |
---|
| 2857 | ( 0 0 1 ) |
---|
| 2858 | |
---|
| 2859 | */ |
---|
| 2860 | c = cos(p); |
---|
| 2861 | s = sin(p); |
---|
| 2862 | make3by3(ah, c, s / sk / scoef, irho * (1.0 - c) / (coef * afk), |
---|
| 2863 | -scoef * sk * s, c, scoef * irho / sk * s, 0.0, 0.0, 1.0); |
---|
| 2864 | sk = sqrt(fabs(k)); |
---|
| 2865 | p = L * sk / scoef; |
---|
| 2866 | c = cosh(p); |
---|
| 2867 | s = sinh(p); |
---|
| 2868 | make3by3(av, c, s / sk / scoef, 0.0, sk * s * scoef, c, 0.0, 0.0, |
---|
| 2869 | 0.0, 1.0); |
---|
| 2870 | } else { |
---|
| 2871 | /* vertically focusing */ |
---|
| 2872 | c = cosh(p); |
---|
| 2873 | s = sinh(p); |
---|
| 2874 | make3by3(ah, c, s / sk / scoef, (c - 1.0) * irho / afk, scoef * s |
---|
| 2875 | * sk, c, scoef * s * irho / sk, 0.0, 0.0, 1.0); |
---|
| 2876 | sk = sqrt(fabs(k)); |
---|
| 2877 | p = L * sk / scoef; |
---|
| 2878 | c = cos(p); |
---|
| 2879 | s = sin(p); |
---|
| 2880 | make3by3(av, c, s / sk / scoef, 0.0, -sk * s * scoef, c, 0.0, 0.0, |
---|
| 2881 | 0.0, 1.0); |
---|
| 2882 | } |
---|
| 2883 | } |
---|
| 2884 | /* Edge focusing, no effect due to gap between AU and AD */ |
---|
| 2885 | |
---|
| 2886 | /* |
---|
| 2887 | (1 0 0) |
---|
| 2888 | Edge(theta) = (h*tan(theta) 1 0) |
---|
| 2889 | (0 0 1) |
---|
| 2890 | |
---|
| 2891 | (1 0 0) |
---|
| 2892 | Edge(theta) = (-h*tan(theta-psi) 1 0) |
---|
| 2893 | (0 0 1) |
---|
| 2894 | |
---|
| 2895 | 2 |
---|
| 2896 | K1*gap*h*(1 + sin phi) |
---|
| 2897 | psi = -----------------------, K1 = 1/2 |
---|
| 2898 | cos phi |
---|
| 2899 | |
---|
| 2900 | */ |
---|
| 2901 | if (phi1 != 0.0 || gap > 0.0) { |
---|
| 2902 | UnitMat(3L, edge); |
---|
| 2903 | edge[1][0] = irho * tan(dtor(phi1)); |
---|
| 2904 | MulRMat(3L, ah, edge); /* ah <- ah*edge */ |
---|
| 2905 | if (true) |
---|
| 2906 | edge[1][0] = -irho * tan(dtor(phi1) - get_psi(irho, phi1, gap)) |
---|
| 2907 | / coef; |
---|
| 2908 | else |
---|
| 2909 | edge[1][0] = -irho * tan(dtor(phi1) - get_psi(irho, phi1, gap)); |
---|
| 2910 | MulRMat(3L, av, edge); /* av <- av*edge */ |
---|
| 2911 | } else if (phi2 != 0.0 || gap < 0.0) { |
---|
| 2912 | UnitMat(3L, edge); |
---|
| 2913 | edge[1][0] = irho * tan(dtor(phi2)); |
---|
| 2914 | MulLMat(3L, edge, ah); /* av <- edge*av */ |
---|
| 2915 | if (true) |
---|
| 2916 | edge[1][0] = -irho * tan(dtor(phi2) |
---|
| 2917 | - get_psi(irho, phi2, fabs(gap))) / coef; |
---|
| 2918 | else |
---|
| 2919 | edge[1][0] = -irho * tan(dtor(phi2) |
---|
| 2920 | - get_psi(irho, phi2, fabs(gap))); |
---|
| 2921 | MulLMat(3L, edge, av); /* av <- edge*av */ |
---|
| 2922 | } |
---|
| 2923 | mergeto4by5(M, ah, av); |
---|
| 2924 | } |
---|
| 2925 | /* ************************************************************** |
---|
| 2926 | * Compute transfert matrix for a quadrupole magnet |
---|
| 2927 | the transfert matrix A is plitted into two part |
---|
| 2928 | A = AD55xAU55 |
---|
| 2929 | |
---|
| 2930 | where AD55 is the downstream transfert matrix |
---|
| 2931 | corresponding to a half magnet w/ an exit angle |
---|
| 2932 | and no entrance angle. |
---|
| 2933 | The linear frindge field is taken into account |
---|
| 2934 | |
---|
| 2935 | where AU55 is the upstream transfert matrix |
---|
| 2936 | corresponding to a half magnet w/ an entrance |
---|
| 2937 | angle and no exit angle. |
---|
[11] | 2938 | The linear fringe field is taken into account |
---|
[3] | 2939 | **************************************************************/ |
---|
| 2940 | void Mpole_Setmatrix(int Fnum1, int Knum1, double K) { |
---|
| 2941 | |
---|
| 2942 | |
---|
| 2943 | CellType *cellp; |
---|
| 2944 | elemtype *elemp; |
---|
| 2945 | MpoleType *M; |
---|
| 2946 | |
---|
| 2947 | if (ElemFam[Fnum1 - 1].nKid <= 0) { |
---|
| 2948 | printf("Mpole_Setmatrix: no kids in famile %d\n", Fnum1); |
---|
| 2949 | return; |
---|
| 2950 | } |
---|
| 2951 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 2952 | elemp = &cellp->Elem; |
---|
| 2953 | M = elemp->M; |
---|
| 2954 | |
---|
| 2955 | bendmat(M->AU55, elemp->PL / 2.0, M->Pirho, M->PTx1, 0.0, M->Pgap, K); |
---|
| 2956 | bendmat(M->AD55, elemp->PL / 2.0, M->Pirho, 0.0, M->PTx2, -M->Pgap, K); |
---|
| 2957 | } |
---|
| 2958 | |
---|
| 2959 | /****************************************************************************/ |
---|
| 2960 | /* void Wiggler_Setmatrix(long Fnum1, long Knum1, double L, double lambda, double k0, double kx) |
---|
| 2961 | |
---|
| 2962 | Purpose: |
---|
| 2963 | |
---|
| 2964 | |
---|
| 2965 | Input: |
---|
| 2966 | none |
---|
| 2967 | |
---|
| 2968 | Output: |
---|
| 2969 | none |
---|
| 2970 | |
---|
| 2971 | Return: |
---|
| 2972 | none |
---|
| 2973 | |
---|
| 2974 | Global variables: |
---|
| 2975 | none |
---|
| 2976 | |
---|
| 2977 | Specific functions: |
---|
| 2978 | none |
---|
| 2979 | |
---|
| 2980 | Comments: |
---|
| 2981 | none |
---|
| 2982 | |
---|
| 2983 | ****************************************************************************/ |
---|
| 2984 | void Wiggler_Setmatrix(int Fnum1, int Knum1, double L, double kx, double kz, |
---|
| 2985 | double k0) { |
---|
| 2986 | double t, s, c, k, ky, LL; |
---|
| 2987 | Matrix ah, av; |
---|
| 2988 | double TEMP; |
---|
| 2989 | WigglerType *W; |
---|
| 2990 | |
---|
| 2991 | LL = L / (1.0 + globval.dPparticle); |
---|
| 2992 | if (kx == 0e0) |
---|
| 2993 | make3by3(ah, 1.0, LL, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0); |
---|
| 2994 | else { |
---|
| 2995 | TEMP = kx / kz; |
---|
| 2996 | k = sqrt(TEMP * TEMP * fabs(k0)); |
---|
| 2997 | t = LL * k; |
---|
| 2998 | c = cosh(t); |
---|
| 2999 | s = sinh(t); |
---|
| 3000 | make3by3(ah, c, s / k, 0.0, k * s, c, 0.0, 0.0, 0.0, 1.0); |
---|
| 3001 | } |
---|
| 3002 | if (k0 == 0e0) |
---|
| 3003 | make3by3(av, 1.0, LL, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0); |
---|
| 3004 | else { |
---|
| 3005 | ky = sqrt(kx * kx + kz * kz); |
---|
| 3006 | TEMP = ky / kz; |
---|
| 3007 | k = sqrt(TEMP * TEMP * fabs(k0)); |
---|
| 3008 | t = LL * k; |
---|
| 3009 | c = cos(t); |
---|
| 3010 | s = sin(t); |
---|
| 3011 | make3by3(av, c, s / k, 0.0, -k * s, c, 0.0, 0.0, 0.0, 1.0); |
---|
| 3012 | } |
---|
| 3013 | W = Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem.W; |
---|
| 3014 | mergeto4by5(W->W55, ah, av); |
---|
| 3015 | } |
---|
| 3016 | |
---|
| 3017 | /****************************************************************************/ |
---|
| 3018 | /* void Mpole_Pass_M(CellType *Cell, double *xref, vector *x) |
---|
| 3019 | |
---|
| 3020 | Purpose: |
---|
| 3021 | |
---|
| 3022 | |
---|
| 3023 | Input: |
---|
| 3024 | none |
---|
| 3025 | |
---|
| 3026 | Output: |
---|
| 3027 | none |
---|
| 3028 | |
---|
| 3029 | Return: |
---|
| 3030 | none |
---|
| 3031 | |
---|
| 3032 | Global variables: |
---|
| 3033 | none |
---|
| 3034 | |
---|
| 3035 | Specific functions: |
---|
| 3036 | none |
---|
| 3037 | |
---|
| 3038 | Comments: |
---|
| 3039 | BUG does nothing if quad .... !!!! Meth_first not see |
---|
| 3040 | maybe just put a globval.MatMeth |
---|
| 3041 | 30/03/03 : changed Meth_First by Meth_Fourth |
---|
| 3042 | |
---|
| 3043 | ****************************************************************************/ |
---|
| 3044 | void Mpole_Pass_M(CellType &Cell, Vector &xref, Matrix &x) { |
---|
| 3045 | double k; |
---|
| 3046 | elemtype *elemp; |
---|
| 3047 | MpoleType *M; |
---|
| 3048 | |
---|
| 3049 | elemp = &Cell.Elem; |
---|
| 3050 | M = elemp->M; |
---|
| 3051 | /* Global -> Local */ |
---|
| 3052 | GtoL_M(x, Cell.dT); |
---|
| 3053 | GtoL(xref, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1); |
---|
| 3054 | |
---|
| 3055 | switch (M->Pmethod) { |
---|
| 3056 | |
---|
| 3057 | case Meth_Linear: |
---|
| 3058 | |
---|
| 3059 | case Meth_Fourth: /* Nothing */ |
---|
| 3060 | // Laurent |
---|
| 3061 | // case Meth_First: /* Nothing */ |
---|
| 3062 | /* Tracy integrator */ |
---|
| 3063 | if (M->Pthick == thick) { |
---|
| 3064 | /* thick element */ |
---|
| 3065 | /* First Linear */ |
---|
| 3066 | MulLMat(5, M->AU55, x); |
---|
| 3067 | LinTrans(5, M->AU55, xref); |
---|
| 3068 | k = M->PB[Quad + HOMmax]; |
---|
| 3069 | M->PB[Quad + HOMmax] = 0.0; |
---|
| 3070 | /* Kick */ |
---|
| 3071 | thin_kick_M(M->Porder, M->PB, elemp->PL, 0.0, xref, x); |
---|
| 3072 | thin_kick(M->Porder, M->PB, elemp->PL, 0.0, 0.0, xref); |
---|
| 3073 | M->PB[Quad + HOMmax] = k; |
---|
| 3074 | /* Second Linear */ |
---|
| 3075 | MulLMat(5L, M->AD55, x); |
---|
| 3076 | LinTrans(5L, M->AD55, xref); |
---|
| 3077 | } else { |
---|
| 3078 | /* thin kick */ |
---|
| 3079 | thin_kick_M(M->Porder, M->PB, 1.0, 0.0, xref, x); |
---|
| 3080 | thin_kick(M->Porder, M->PB, 1.0, 0.0, 0.0, xref); |
---|
| 3081 | } |
---|
| 3082 | break; |
---|
| 3083 | } |
---|
| 3084 | |
---|
| 3085 | /* Local -> Global */ |
---|
| 3086 | LtoG_M(x, Cell.dT); |
---|
| 3087 | LtoG(xref, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1); |
---|
| 3088 | } |
---|
| 3089 | |
---|
| 3090 | /****************************************************************************/ |
---|
| 3091 | /* void Wiggler_Pass_M(CellType *Cell, double *xref, vector *x) |
---|
| 3092 | |
---|
| 3093 | Purpose: |
---|
| 3094 | |
---|
| 3095 | |
---|
| 3096 | Input: |
---|
| 3097 | none |
---|
| 3098 | |
---|
| 3099 | Output: |
---|
| 3100 | none |
---|
| 3101 | |
---|
| 3102 | Return: |
---|
| 3103 | none |
---|
| 3104 | |
---|
| 3105 | Global variables: |
---|
| 3106 | none |
---|
| 3107 | |
---|
| 3108 | Specific functions: |
---|
| 3109 | none |
---|
| 3110 | |
---|
| 3111 | Comments: |
---|
| 3112 | none |
---|
| 3113 | |
---|
| 3114 | ****************************************************************************/ |
---|
| 3115 | void Wiggler_Pass_M(CellType &Cell, Vector &xref, Matrix &x) { |
---|
| 3116 | elemtype *elemp; |
---|
| 3117 | WigglerType *W; |
---|
| 3118 | |
---|
| 3119 | elemp = &Cell.Elem; |
---|
| 3120 | W = elemp->W; |
---|
| 3121 | |
---|
| 3122 | /* Global -> Local */ |
---|
| 3123 | GtoL_M(x, Cell.dT); |
---|
| 3124 | GtoL(xref, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 3125 | |
---|
| 3126 | switch (W->Pmethod) { |
---|
| 3127 | case Meth_Linear: /* Nothing */ |
---|
| 3128 | /* Tracy integrator */ |
---|
| 3129 | MulLMat(5, W->W55, x); |
---|
| 3130 | LinTrans(5L, W->W55, xref); |
---|
| 3131 | break; |
---|
| 3132 | } |
---|
| 3133 | |
---|
| 3134 | /* Local -> Global */ |
---|
| 3135 | LtoG_M(x, Cell.dT); |
---|
| 3136 | LtoG(xref, Cell.dS, Cell.dT, 0.0, 0.0, 0.0); |
---|
| 3137 | } |
---|
| 3138 | |
---|
| 3139 | void Insertion_SetMatrix(int Fnum1, int Knum1) { |
---|
| 3140 | /* void Insertion_SetMatrix(int Fnum1, int Knum1) |
---|
| 3141 | |
---|
| 3142 | Purpose: called by Insertion_Init |
---|
| 3143 | Constructs the linear matrices |
---|
| 3144 | K55 kick matrix for one slice |
---|
| 3145 | D55 drift matrix for one slice |
---|
| 3146 | KD55 full linear transport matrix |
---|
| 3147 | |
---|
| 3148 | Input: |
---|
| 3149 | Fnum1 Family number |
---|
| 3150 | Knum1 Kid number |
---|
| 3151 | |
---|
| 3152 | Output: |
---|
| 3153 | none |
---|
| 3154 | |
---|
| 3155 | Return: |
---|
| 3156 | none |
---|
| 3157 | |
---|
| 3158 | Global variables: |
---|
| 3159 | globval |
---|
| 3160 | |
---|
| 3161 | Specific functions: |
---|
| 3162 | LinearInterpDeriv2 |
---|
| 3163 | |
---|
| 3164 | Comments: |
---|
| 3165 | 04/07/03 derivative interpolation around closed orbit |
---|
| 3166 | 10/01/05 First order kick added |
---|
| 3167 | |
---|
| 3168 | Need to be checked energy dependence and so on. */ |
---|
| 3169 | |
---|
| 3170 | int i = 0; |
---|
| 3171 | double L = 0.0; |
---|
| 3172 | CellType *cellp; |
---|
| 3173 | elemtype *elemp; |
---|
| 3174 | InsertionType *ID; |
---|
| 3175 | double alpha0 = 0.0, alpha02 = 0.0; |
---|
| 3176 | double DTHXDX = 0.0, DTHXDZ = 0.0, DTHZDX = 0.0, DTHZDZ = 0.0; |
---|
| 3177 | int Nslice = 0; |
---|
| 3178 | |
---|
| 3179 | if (ElemFam[Fnum1 - 1].nKid <= 0) |
---|
| 3180 | return; |
---|
| 3181 | |
---|
| 3182 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 3183 | elemp = &cellp->Elem; |
---|
| 3184 | ID = elemp->ID; |
---|
| 3185 | Nslice = ID->PN; |
---|
| 3186 | alpha0 = c0 / globval.Energy * 1E-9 * ID->scaling1; |
---|
| 3187 | alpha02 = (c0 / globval.Energy * 1E-9) * c0 / globval.Energy * 1E-9 / (1.0 |
---|
| 3188 | + globval.dPparticle) * ID->scaling2; |
---|
| 3189 | |
---|
| 3190 | UnitMat(6L, ID->D55); |
---|
| 3191 | UnitMat(6L, ID->K55); |
---|
| 3192 | UnitMat(6L, ID->KD55); |
---|
| 3193 | |
---|
| 3194 | // if (globval.radiation == false && globval.Cavity_on == false) |
---|
| 3195 | // { |
---|
| 3196 | /* (Nslice + 1) Drifts for Nslice Kicks */ |
---|
| 3197 | |
---|
| 3198 | /* Drift Matrix */ |
---|
| 3199 | L = elemp->PL / (Nslice + 1) / (1.0 + globval.dPparticle); |
---|
| 3200 | make4by5(ID->D55, 1.0, L, 0.0, 0.0, 1.0, 0.0, 1.0, L, 0.0, 0.0, 1.0, 0.0); |
---|
| 3201 | |
---|
| 3202 | /* First order Kick */ |
---|
| 3203 | if (ID->firstorder) { |
---|
| 3204 | /* quadrupole Kick matrix linearized around closed orbit */ |
---|
| 3205 | if (!ID->linear) { |
---|
| 3206 | // SplineInterpDeriv3(cellp->BeamPos[0], cellp->BeamPos[2], |
---|
| 3207 | // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp); |
---|
| 3208 | } else { |
---|
| 3209 | // LinearInterpDeriv2(cellp->BeamPos[0], cellp->BeamPos[2], |
---|
| 3210 | // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp, 1); |
---|
| 3211 | } |
---|
| 3212 | ID->K55[1][0] = ID->K55[1][0] + alpha0 * DTHXDX / Nslice; |
---|
| 3213 | ID->K55[1][2] = ID->K55[1][2] + alpha0 * DTHXDZ / Nslice; |
---|
| 3214 | ID->K55[3][0] = ID->K55[3][0] + alpha0 * DTHZDX / Nslice; |
---|
| 3215 | ID->K55[3][2] = ID->K55[3][2] + alpha0 * DTHZDZ / Nslice; |
---|
| 3216 | } |
---|
| 3217 | |
---|
| 3218 | /* Second order Kick */ |
---|
| 3219 | if (ID->secondorder) { |
---|
| 3220 | /* quadrupole Kick matrix linearized around closed orbit */ |
---|
| 3221 | if (!ID->linear) { |
---|
| 3222 | // SplineInterpDeriv3(cellp->BeamPos[0], cellp->BeamPos[2], |
---|
| 3223 | // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp); |
---|
| 3224 | } else { |
---|
| 3225 | // LinearInterpDeriv2(cellp->BeamPos[0], cellp->BeamPos[2], |
---|
| 3226 | // &DTHXDX, &DTHXDZ, &DTHZDX, &DTHZDZ, cellp, 2); |
---|
| 3227 | } |
---|
| 3228 | ID->K55[1][0] = ID->K55[1][0] + alpha02 * DTHXDX / Nslice; |
---|
| 3229 | ID->K55[1][2] = ID->K55[1][2] + alpha02 * DTHXDZ / Nslice; |
---|
| 3230 | ID->K55[3][0] = ID->K55[3][0] + alpha02 * DTHZDX / Nslice; |
---|
| 3231 | ID->K55[3][2] = ID->K55[3][2] + alpha02 * DTHZDZ / Nslice; |
---|
| 3232 | } |
---|
| 3233 | |
---|
| 3234 | MulLMat(6L, ID->D55, ID->KD55); |
---|
| 3235 | |
---|
| 3236 | for (i = 1; i <= Nslice; i++) { |
---|
| 3237 | MulLMat(6L, ID->K55, ID->KD55); |
---|
| 3238 | MulLMat(6L, ID->D55, ID->KD55); |
---|
| 3239 | } |
---|
| 3240 | |
---|
| 3241 | // } |
---|
| 3242 | // else |
---|
| 3243 | // { |
---|
| 3244 | // L = elemp->PL/(1.0 + globval.dPparticle); /* L = L/(1 + dP) */ |
---|
| 3245 | // make4by5(ID->KD55, |
---|
| 3246 | // 1.0, L, 0.0, 0.0, 1.0, 0.0, |
---|
| 3247 | // 1.0, L, 0.0, 0.0, 1.0, 0.0); |
---|
| 3248 | // } |
---|
| 3249 | } |
---|
| 3250 | |
---|
| 3251 | /********************************************************************* |
---|
| 3252 | Purpose: called by Elem_Pass_M |
---|
| 3253 | matrix propagation through a insertion kick-driftlike matrix |
---|
| 3254 | x = KD55*x |
---|
| 3255 | xref= insertion(xref) |
---|
| 3256 | |
---|
| 3257 | Input: |
---|
| 3258 | xref vector |
---|
| 3259 | x matrix |
---|
| 3260 | |
---|
| 3261 | Output: |
---|
| 3262 | xref |
---|
| 3263 | x |
---|
| 3264 | |
---|
| 3265 | Return: |
---|
| 3266 | none |
---|
| 3267 | |
---|
| 3268 | Global variables: |
---|
| 3269 | none |
---|
| 3270 | |
---|
| 3271 | Specific functions: |
---|
| 3272 | MulLMat, Drft |
---|
| 3273 | |
---|
| 3274 | Comments: |
---|
| 3275 | 01/07/03 6D tracking activated |
---|
| 3276 | ******************************************************************************/ |
---|
| 3277 | void Insertion_Pass_M(CellType &Cell, Vector &xref, Matrix &M) { |
---|
| 3278 | |
---|
| 3279 | |
---|
| 3280 | elemtype *elemp; |
---|
| 3281 | |
---|
| 3282 | elemp = &Cell.Elem; |
---|
| 3283 | |
---|
| 3284 | /* Global -> Local */ |
---|
| 3285 | // GtoL_M(x, Cell->dT); |
---|
| 3286 | // GtoL(xref, Cell->dS, Cell->dT, 0.0, 0.0, 0.0); |
---|
| 3287 | // if (globval.radiation == false && globval.Cavity_on == false) |
---|
| 3288 | // { |
---|
| 3289 | MulLMat(5, elemp->ID->KD55, M); /* M<-KD55*M */ |
---|
| 3290 | LinTrans(5, elemp->ID->KD55, xref); |
---|
| 3291 | // } |
---|
| 3292 | // else |
---|
| 3293 | // { |
---|
| 3294 | // MulLMat(5L, elemp->ID->D55, M); /* X<-D55*X */ |
---|
| 3295 | // Drft(elemp->PL, elemp->PL/(1.0+xref[4]), xref); |
---|
| 3296 | // } |
---|
| 3297 | |
---|
| 3298 | /* Local -> Global */ |
---|
| 3299 | // LtoG_M(x, Cell->dT); |
---|
| 3300 | // LtoG(xref, Cell->dS, Cell->dT, 0.0, 0.0, 0.0); |
---|
| 3301 | } |
---|
| 3302 | |
---|
| 3303 | /****************************************************************************/ |
---|
| 3304 | /* void getelem(long i, CellType *cellrec) |
---|
| 3305 | |
---|
| 3306 | Purpose: |
---|
| 3307 | assign all the information of i-th element from array Cell[i] to pointer cellrec |
---|
| 3308 | |
---|
| 3309 | Input: |
---|
| 3310 | |
---|
| 3311 | Output: |
---|
| 3312 | none |
---|
| 3313 | |
---|
| 3314 | Return: |
---|
| 3315 | none |
---|
| 3316 | |
---|
| 3317 | Global variables: |
---|
| 3318 | none |
---|
| 3319 | |
---|
| 3320 | Specific functions: |
---|
| 3321 | none |
---|
| 3322 | |
---|
| 3323 | Comments: |
---|
| 3324 | |
---|
| 3325 | ****************************************************************************/ |
---|
| 3326 | void getelem(long i, CellType *cellrec) { |
---|
| 3327 | *cellrec = Cell[i]; |
---|
| 3328 | } |
---|
| 3329 | /****************************************************************************/ |
---|
| 3330 | /* void putelem(long i, CellType *cellrec) |
---|
| 3331 | |
---|
| 3332 | Purpose: |
---|
| 3333 | assign all the information of pointer cellrec to i-th element to array Cell[i] |
---|
| 3334 | |
---|
| 3335 | Input: |
---|
| 3336 | |
---|
| 3337 | Output: |
---|
| 3338 | none |
---|
| 3339 | |
---|
| 3340 | Return: |
---|
| 3341 | none |
---|
| 3342 | |
---|
| 3343 | Global variables: |
---|
| 3344 | none |
---|
| 3345 | |
---|
| 3346 | Specific functions: |
---|
| 3347 | none |
---|
| 3348 | |
---|
| 3349 | Comments: |
---|
| 3350 | |
---|
| 3351 | ****************************************************************************/ |
---|
| 3352 | void putelem(long i, CellType *cellrec) { |
---|
| 3353 | Cell[i] = *cellrec; |
---|
| 3354 | } |
---|
| 3355 | |
---|
| 3356 | /****************************************************************************/ |
---|
| 3357 | /* int GetnKid(const int Fnum1) |
---|
| 3358 | |
---|
| 3359 | Purpose: |
---|
| 3360 | return the number of kid in the family |
---|
| 3361 | |
---|
| 3362 | Input: |
---|
| 3363 | Fnum1 : family number |
---|
| 3364 | |
---|
| 3365 | |
---|
| 3366 | Output: |
---|
| 3367 | none |
---|
| 3368 | |
---|
| 3369 | Return: |
---|
| 3370 | none |
---|
| 3371 | |
---|
| 3372 | Global variables: |
---|
| 3373 | none |
---|
| 3374 | |
---|
| 3375 | Specific functions: |
---|
| 3376 | none |
---|
| 3377 | |
---|
| 3378 | Comments: |
---|
| 3379 | |
---|
| 3380 | ****************************************************************************/ |
---|
| 3381 | int GetnKid(const int Fnum1) { |
---|
| 3382 | return (ElemFam[Fnum1 - 1].nKid); |
---|
| 3383 | } |
---|
| 3384 | |
---|
| 3385 | /****************************************************************************/ |
---|
| 3386 | /* long Elem_GetPos(const int Fnum1, const int Knum1) |
---|
| 3387 | |
---|
| 3388 | Purpose: |
---|
| 3389 | get the element index in the lattice |
---|
| 3390 | |
---|
| 3391 | |
---|
| 3392 | Input: |
---|
| 3393 | Fnum1 : family number of the element |
---|
| 3394 | Knum1 : kid number of the element in Fnum1 |
---|
| 3395 | |
---|
| 3396 | Output: |
---|
| 3397 | none |
---|
| 3398 | |
---|
| 3399 | Return: |
---|
| 3400 | loc : element index in the lattice |
---|
| 3401 | |
---|
| 3402 | Global variables: |
---|
| 3403 | none |
---|
| 3404 | |
---|
| 3405 | Specific functions: |
---|
| 3406 | none |
---|
| 3407 | |
---|
| 3408 | Comments: |
---|
| 3409 | example: |
---|
| 3410 | long FORLIM = GetnKid(ElemIndex("CH")); // get number of CH |
---|
| 3411 | // search element position for Family CH |
---|
| 3412 | for(k=1;k<FORLIM;k++){ |
---|
| 3413 | fprintf(stdout, "elem %d is at position %ld \n", k, Elem_GetPos(ElemIndex("CH"), k)); |
---|
| 3414 | } |
---|
| 3415 | |
---|
| 3416 | |
---|
| 3417 | 21/12/2011 Jianfeng Zhang@ soleil |
---|
| 3418 | Add warning message: when call Elem_GetPos(), the kid index knum1 start from 1 !!!!! |
---|
| 3419 | |
---|
| 3420 | ****************************************************************************/ |
---|
| 3421 | long Elem_GetPos(const int Fnum1, const int Knum1) { |
---|
| 3422 | long int loc; |
---|
| 3423 | |
---|
| 3424 | if(Knum1 < 1){ |
---|
| 3425 | cout << "Elem_GetPos: kid index of the family starts from 1 !!!" << endl; |
---|
| 3426 | cout << "Element: " << ElemFam[Fnum1 - 1].ElemF.PName << "with Fnum: " <<Fnum1<<" Knum: "<<Knum1<<endl; |
---|
| 3427 | exit_(1); |
---|
| 3428 | } |
---|
| 3429 | if (ElemFam[Fnum1 - 1].nKid != 0) |
---|
| 3430 | loc = ElemFam[Fnum1 - 1].KidList[Knum1 - 1]; |
---|
| 3431 | else { |
---|
| 3432 | loc = -1; |
---|
| 3433 | printf("Elem_GetPos: there are no kids in family %d (%s)\n", Fnum1, |
---|
| 3434 | ElemFam[Fnum1 - 1].ElemF.PName); |
---|
[11] | 3435 | exit_(1); |
---|
[3] | 3436 | } |
---|
| 3437 | |
---|
| 3438 | return loc; |
---|
| 3439 | } |
---|
| 3440 | |
---|
| 3441 | static double thirdroot(double a) { |
---|
| 3442 | /* By substitution method */ |
---|
| 3443 | int i; |
---|
| 3444 | double x; |
---|
| 3445 | |
---|
| 3446 | x = 1.0; |
---|
| 3447 | i = 0; |
---|
| 3448 | do { |
---|
| 3449 | i++; |
---|
| 3450 | x = (x + a) / (x * x + 1e0); |
---|
| 3451 | } while (i != 250); |
---|
| 3452 | return x; |
---|
| 3453 | } |
---|
| 3454 | |
---|
| 3455 | void SI_init(void) { |
---|
| 3456 | /* c_1 = 1/(2*(2-2^(1/3))), c_2 = (1-2^(1/3))/(2*(2-2^(1/3))) |
---|
| 3457 | d_1 = 1/(2-2^(1/3)), d_2 = -2^(1/3)/(2-2^(1/3)) */ |
---|
| 3458 | |
---|
| 3459 | double C_gamma, C_u; |
---|
| 3460 | |
---|
| 3461 | c_1 = 1e0 / (2e0 * (2e0 - thirdroot(2e0))); |
---|
| 3462 | c_2 = 0.5e0 - c_1; |
---|
| 3463 | d_1 = 2e0 * c_1; |
---|
| 3464 | d_2 = 1e0 - 2e0 * d_1; |
---|
| 3465 | |
---|
[11] | 3466 | |
---|
| 3467 | |
---|
| 3468 | |
---|
| 3469 | |
---|
[3] | 3470 | // classical radiation |
---|
| 3471 | // C_gamma = 8.846056192e-05; |
---|
| 3472 | // C_gamma = 4 * pi * r_e [m] / ( 3 * (m_e [GeV/c^2] * c^2)^3 ) |
---|
| 3473 | C_gamma = 4.0 * M_PI * r_e / (3.0 * cube(1e-9 * m_e)); |
---|
| 3474 | // P_gamma = e^2 c^3 / 2 / pi * C_gamma (E [GeV])^2 (B [T])^2 |
---|
| 3475 | // p_s = P_s/P, E = P*c, B/(Brho) = p/e |
---|
| 3476 | cl_rad = C_gamma * cube(globval.Energy) / (2.0 * M_PI); |
---|
| 3477 | |
---|
| 3478 | // eletron rest mass [GeV]: slightly off??? |
---|
| 3479 | // m_e_ = 0.5110034e-03; |
---|
| 3480 | // h_bar times c [GeV m] |
---|
| 3481 | // hbar_t_c = 1.9732858e-16; |
---|
| 3482 | // quantum fluctuations |
---|
| 3483 | C_u = 55.0 / (24.0 * sqrt(3.0)); |
---|
| 3484 | q_fluct = 3.0 * C_u * C_gamma * 1e-9 * h_bar * c0 / (4.0 * M_PI * cube(1e-9 |
---|
| 3485 | * m_e)) * pow(globval.Energy, 5.0); |
---|
| 3486 | } |
---|
| 3487 | |
---|
| 3488 | static void Mpole_Print(FILE *f, int Fnum1) { |
---|
| 3489 | elemtype *elemp; |
---|
| 3490 | MpoleType *M; |
---|
| 3491 | |
---|
| 3492 | elemp = &ElemFam[Fnum1 - 1].ElemF; |
---|
| 3493 | M = elemp->M; |
---|
| 3494 | fprintf(f, "Element[%3d ] \n", Fnum1); |
---|
| 3495 | fprintf(f, " Name: %.*s, Kind: mpole, L=% .8E\n", SymbolLength, |
---|
| 3496 | elemp->PName, elemp->PL); |
---|
| 3497 | fprintf(f, " Method: %d, N=%4d\n", M->Pmethod, M->PN); |
---|
| 3498 | } |
---|
| 3499 | |
---|
| 3500 | /**************************************************************************** |
---|
| 3501 | * void Drift_Print(FILE **f, long Fnum1) |
---|
| 3502 | |
---|
| 3503 | Purpose: called by Elem_Print |
---|
| 3504 | Print out drift characteristics in a file |
---|
| 3505 | Name, kind, length, method, slice number |
---|
| 3506 | |
---|
| 3507 | Input: |
---|
| 3508 | Fnum1 Family number |
---|
| 3509 | f pointer on file id |
---|
| 3510 | |
---|
| 3511 | Output: |
---|
| 3512 | none |
---|
| 3513 | |
---|
| 3514 | Return: |
---|
| 3515 | none |
---|
| 3516 | |
---|
| 3517 | Global variables: |
---|
| 3518 | none |
---|
| 3519 | |
---|
| 3520 | Specific functions: |
---|
| 3521 | none |
---|
| 3522 | |
---|
| 3523 | Comments: |
---|
| 3524 | none |
---|
| 3525 | |
---|
| 3526 | ****************************************************************************/ |
---|
| 3527 | static void Drift_Print(FILE *f, int Fnum1) { |
---|
| 3528 | ElemFamType *elemfamp; |
---|
| 3529 | elemtype *elemp; |
---|
| 3530 | |
---|
| 3531 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 3532 | elemp = &elemfamp->ElemF; |
---|
| 3533 | fprintf(f, "Element[%3d ] \n", Fnum1); |
---|
| 3534 | fprintf(f, " Name: %.*s, Kind: drift, L=% .8E\n", SymbolLength, |
---|
| 3535 | elemp->PName, elemp->PL); |
---|
| 3536 | fprintf(f, " nKid:%3d\n\n", elemfamp->nKid); |
---|
| 3537 | } |
---|
| 3538 | |
---|
| 3539 | /**************************************************************************** |
---|
| 3540 | * void Wiggler_Print(FILE **f, long Fnum1) |
---|
| 3541 | |
---|
| 3542 | Purpose: called by Elem_Print |
---|
| 3543 | Print out drift characteristics in a file |
---|
| 3544 | Name, kind, length |
---|
| 3545 | |
---|
| 3546 | Input: |
---|
| 3547 | Fnum1 Family number |
---|
| 3548 | f pointer on file id |
---|
| 3549 | |
---|
| 3550 | Output: |
---|
| 3551 | none |
---|
| 3552 | |
---|
| 3553 | Return: |
---|
| 3554 | none |
---|
| 3555 | |
---|
| 3556 | Global variables: |
---|
| 3557 | none |
---|
| 3558 | |
---|
| 3559 | Specific functions: |
---|
| 3560 | none |
---|
| 3561 | |
---|
| 3562 | Comments: |
---|
| 3563 | none |
---|
| 3564 | |
---|
| 3565 | ****************************************************************************/ |
---|
| 3566 | static void Wiggler_Print(FILE *f, int Fnum1) { |
---|
| 3567 | elemtype *elemp; |
---|
| 3568 | |
---|
| 3569 | elemp = &ElemFam[Fnum1 - 1].ElemF; |
---|
| 3570 | fprintf(f, "Element[%3d ] \n", Fnum1); |
---|
| 3571 | fprintf(f, " Name: %.*s, Kind: wiggler, L=% .8E\n\n", NameLength, |
---|
| 3572 | elemp->PName, elemp->PL); |
---|
| 3573 | } |
---|
| 3574 | |
---|
| 3575 | /**************************************************************************** |
---|
| 3576 | * void Insertion_Print(FILE **f, long Fnum1) |
---|
| 3577 | |
---|
| 3578 | Purpose: called by Elem_Print |
---|
| 3579 | Print out drift characteristics in a file |
---|
| 3580 | Name, kind, length |
---|
| 3581 | |
---|
| 3582 | Input: |
---|
| 3583 | Fnum1 Family number |
---|
| 3584 | f pointer on file id |
---|
| 3585 | |
---|
| 3586 | Output: |
---|
| 3587 | none |
---|
| 3588 | |
---|
| 3589 | Return: |
---|
| 3590 | none |
---|
| 3591 | |
---|
| 3592 | Global variables: |
---|
| 3593 | none |
---|
| 3594 | |
---|
| 3595 | Specific functions: |
---|
| 3596 | none |
---|
| 3597 | |
---|
| 3598 | Comments: |
---|
| 3599 | none |
---|
| 3600 | |
---|
| 3601 | ****************************************************************************/ |
---|
| 3602 | static void Insertion_Print(FILE *f, int Fnum1) { |
---|
| 3603 | elemtype *elemp; |
---|
| 3604 | |
---|
| 3605 | elemp = &ElemFam[Fnum1 - 1].ElemF; |
---|
| 3606 | fprintf(f, "Element[%3d ] \n", Fnum1); |
---|
| 3607 | fprintf(f, " Name: %.*s, Kind: wiggler, L=% .8E\n\n", SymbolLength, |
---|
| 3608 | elemp->PName, elemp->PL); |
---|
| 3609 | } |
---|
| 3610 | |
---|
| 3611 | /**************************************************************************** |
---|
| 3612 | * void Insertion_SetMatrix(long Fnum1, long Knum1) |
---|
| 3613 | |
---|
| 3614 | Purpose: called by Insertion_Init |
---|
| 3615 | Constructs the linear matrices |
---|
| 3616 | K55 kick matrix for one slice |
---|
| 3617 | D55 drift matrix for one slice |
---|
| 3618 | KD55 full linear transport matrix |
---|
| 3619 | |
---|
| 3620 | Input: |
---|
| 3621 | Fnum1 Family number |
---|
| 3622 | Knum1 Kid number |
---|
| 3623 | |
---|
| 3624 | Output: |
---|
| 3625 | none |
---|
| 3626 | |
---|
| 3627 | Return: |
---|
| 3628 | none |
---|
| 3629 | |
---|
| 3630 | Global variables: |
---|
| 3631 | globval |
---|
| 3632 | |
---|
| 3633 | Specific functions: |
---|
| 3634 | LinearInterpDeriv2 |
---|
| 3635 | |
---|
| 3636 | Comments: |
---|
| 3637 | 04/07/03 derivative interpolation around closed orbit |
---|
| 3638 | 10/01/05 First order kick added |
---|
| 3639 | |
---|
| 3640 | Need to be checked energy dependence and so on. |
---|
| 3641 | ****************************************************************************/ |
---|
| 3642 | |
---|
| 3643 | void Elem_Print(FILE *f, int Fnum1) { |
---|
| 3644 | int i; |
---|
| 3645 | |
---|
| 3646 | if (Fnum1 == 0) { |
---|
| 3647 | // print all elements |
---|
| 3648 | for (i = 1; i <= globval.Elem_nFam; i++) |
---|
| 3649 | Elem_Print(f, i); |
---|
| 3650 | return; |
---|
| 3651 | } |
---|
| 3652 | |
---|
| 3653 | switch (ElemFam[Fnum1 - 1].ElemF.Pkind) { |
---|
| 3654 | case drift: |
---|
| 3655 | Drift_Print(f, Fnum1); |
---|
| 3656 | break; |
---|
| 3657 | |
---|
| 3658 | case Mpole: |
---|
| 3659 | Mpole_Print(f, Fnum1); |
---|
| 3660 | break; |
---|
| 3661 | case Wigl: |
---|
| 3662 | Wiggler_Print(f, Fnum1); |
---|
| 3663 | break; |
---|
| 3664 | case FieldMap: |
---|
| 3665 | break; |
---|
| 3666 | case Insertion: |
---|
| 3667 | Insertion_Print(f, Fnum1); |
---|
| 3668 | break; |
---|
| 3669 | case Cavity: |
---|
| 3670 | break; |
---|
| 3671 | case marker: |
---|
| 3672 | break; |
---|
| 3673 | case Spreader: |
---|
| 3674 | break; |
---|
| 3675 | case Recombiner: |
---|
| 3676 | break; |
---|
| 3677 | case Solenoid: |
---|
| 3678 | break; |
---|
| 3679 | case undef: |
---|
| 3680 | break; |
---|
| 3681 | } |
---|
| 3682 | } |
---|
| 3683 | |
---|
| 3684 | double Mpole_GetPB(int Fnum1, int Knum1, int Order); |
---|
| 3685 | |
---|
| 3686 | /**************************************************************************** |
---|
| 3687 | * double Elem_GetKval(long Fnum1, long Knum1, long Order) |
---|
| 3688 | |
---|
| 3689 | Purpose: |
---|
| 3690 | Get K values |
---|
| 3691 | |
---|
| 3692 | Input: |
---|
| 3693 | Fnum1 Famility number |
---|
| 3694 | Knum1 Kids number |
---|
| 3695 | Order mutipole component 1 for dipole, 2 for quadrupole) |
---|
| 3696 | |
---|
| 3697 | Output: |
---|
| 3698 | none |
---|
| 3699 | |
---|
| 3700 | Return: |
---|
| 3701 | 0.0 if drift |
---|
| 3702 | integrated field if multipole |
---|
| 3703 | |
---|
| 3704 | |
---|
| 3705 | Global variables: |
---|
| 3706 | ElemFam |
---|
| 3707 | |
---|
| 3708 | Specific functions: |
---|
| 3709 | Mpole_GetPB |
---|
| 3710 | |
---|
| 3711 | Comments: |
---|
| 3712 | 01/02/03 add return = 0 for marker and cavity |
---|
| 3713 | 22/04/03 Insertion added |
---|
| 3714 | |
---|
| 3715 | ****************************************************************************/ |
---|
| 3716 | double Elem_GetKval(int Fnum1, int Knum1, int Order) { |
---|
| 3717 | double Result = 0.0; |
---|
| 3718 | elemtype *elemp; |
---|
| 3719 | |
---|
| 3720 | if (Fnum1 > 0) { |
---|
| 3721 | elemp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem; |
---|
| 3722 | switch (elemp->Pkind) { |
---|
| 3723 | case drift: |
---|
| 3724 | Result = 0.0; |
---|
| 3725 | break; |
---|
| 3726 | case marker: |
---|
| 3727 | Result = 0.0; |
---|
| 3728 | break; |
---|
| 3729 | case Cavity: |
---|
| 3730 | Result = 0.0; |
---|
| 3731 | break; |
---|
| 3732 | case Mpole: /* KL*/ |
---|
| 3733 | if (elemp->M->Pthick == thick) |
---|
| 3734 | Result = elemp->PL * Mpole_GetPB(Fnum1, Knum1, Order); |
---|
| 3735 | else |
---|
| 3736 | Result = Mpole_GetPB(Fnum1, Knum1, Order); |
---|
| 3737 | break; |
---|
| 3738 | case Wigl: |
---|
| 3739 | Result |
---|
| 3740 | = elemp->PL |
---|
| 3741 | * sqrt(2.0 * Cell[ElemFam[Fnum1 - 1].KidList[Knum1 |
---|
| 3742 | - 1]].Elem.W->PBW[Order + HOMmax]); |
---|
| 3743 | break; |
---|
| 3744 | case FieldMap: |
---|
| 3745 | Result = 0.0; |
---|
| 3746 | break; |
---|
| 3747 | case Insertion: |
---|
| 3748 | Result = 0.0; |
---|
| 3749 | break; |
---|
| 3750 | case Spreader: |
---|
| 3751 | Result = 0.0; |
---|
| 3752 | break; |
---|
| 3753 | case Recombiner: |
---|
| 3754 | Result = 0.0; |
---|
| 3755 | break; |
---|
| 3756 | case Solenoid: |
---|
| 3757 | Result = 0.0; |
---|
| 3758 | break; |
---|
| 3759 | case undef: |
---|
| 3760 | break; |
---|
| 3761 | } |
---|
| 3762 | } else |
---|
| 3763 | Result = 0.0; |
---|
| 3764 | |
---|
| 3765 | return Result; |
---|
| 3766 | } |
---|
| 3767 | |
---|
| 3768 | #define n 4 |
---|
| 3769 | void LinsTrans(Matrix &A, Vector &b) { |
---|
| 3770 | int j; |
---|
| 3771 | Vector c; |
---|
| 3772 | |
---|
| 3773 | CopyVec(n, b, c); /* c=b */ |
---|
| 3774 | LinTrans(n, A, c); /* c<-A*c */ |
---|
| 3775 | for (j = 0; j < n; j++) |
---|
| 3776 | c[j] += A[j][n] * b[n] + A[n][j]; |
---|
| 3777 | CopyVec(n, c, b); /* b=c */ |
---|
| 3778 | } |
---|
| 3779 | #undef n |
---|
| 3780 | |
---|
| 3781 | #define n 4 |
---|
| 3782 | void MulLsMat(Matrix &A, Matrix &B) { |
---|
| 3783 | int i, k; |
---|
| 3784 | Matrix C; |
---|
| 3785 | |
---|
| 3786 | CopyMat(n, B, C); /* C<-B */ |
---|
| 3787 | MulLMat(n, A, C); /* C<-A*C */ |
---|
| 3788 | for (i = 0; i < n; i++) { |
---|
| 3789 | C[i][n] = A[i][n]; |
---|
| 3790 | C[n][i] = 0.0; |
---|
| 3791 | for (k = 0; k < n; k++) { |
---|
| 3792 | C[i][n] += A[i][k] * B[k][n]; |
---|
| 3793 | C[n][i] += A[i][k] * B[n][k]; |
---|
| 3794 | } |
---|
| 3795 | } |
---|
| 3796 | C[n][n] = 1.0; |
---|
| 3797 | CopyMat(n + 1, C, B); /* B<-C */ |
---|
| 3798 | } |
---|
| 3799 | #undef n |
---|
| 3800 | |
---|
| 3801 | /**************************************************************************** |
---|
| 3802 | * void Drift_Alloc(elemtype *Elem) |
---|
| 3803 | |
---|
| 3804 | Purpose: |
---|
| 3805 | Dynamic memory allocation for drift element |
---|
| 3806 | |
---|
| 3807 | Input: |
---|
| 3808 | Pointer on element |
---|
| 3809 | |
---|
| 3810 | Output: |
---|
| 3811 | memory space for drift in Elem->UU.D |
---|
| 3812 | |
---|
| 3813 | Return: |
---|
| 3814 | none |
---|
| 3815 | |
---|
| 3816 | Global variables: |
---|
| 3817 | none |
---|
| 3818 | |
---|
| 3819 | Specific functions: |
---|
| 3820 | none |
---|
| 3821 | |
---|
| 3822 | Comments: |
---|
| 3823 | none |
---|
| 3824 | |
---|
| 3825 | ****************************************************************************/ |
---|
| 3826 | void Drift_Alloc(elemtype *Elem) { |
---|
| 3827 | Elem->D = (DriftType *) malloc(sizeof(DriftType)); |
---|
| 3828 | } |
---|
| 3829 | |
---|
| 3830 | void Mpole_Alloc(elemtype *Elem) { |
---|
| 3831 | int j; |
---|
| 3832 | MpoleType *M; |
---|
| 3833 | |
---|
| 3834 | /* Memory allocation */ |
---|
| 3835 | Elem->M = (MpoleType *) malloc(sizeof(MpoleType)); |
---|
| 3836 | M = Elem->M; |
---|
| 3837 | M->Pmethod = Meth_Fourth; |
---|
| 3838 | M->PN = 0; |
---|
| 3839 | /* Displacement errors */ |
---|
| 3840 | for (j = 0; j <= 1; j++) { |
---|
| 3841 | M->PdSsys[j] = 0.0; |
---|
| 3842 | M->PdSrms[j] = 0.0; |
---|
| 3843 | M->PdSrnd[j] = 0.0; |
---|
| 3844 | } |
---|
| 3845 | M->PdTpar = 0.0; /* Roll angle */ |
---|
| 3846 | M->PdTsys = 0.0; /* systematic Roll errors */ |
---|
| 3847 | M->PdTrms = 0.0; /* random Roll errors */ |
---|
| 3848 | M->PdTrnd = 0.0; /* random seed */ |
---|
| 3849 | for (j = -HOMmax; j <= HOMmax; j++) { |
---|
| 3850 | /* Initializes multipoles strengths to zero */ |
---|
| 3851 | M->PB[j + HOMmax] = 0.0; |
---|
| 3852 | M->PBpar[j + HOMmax] = 0.0; |
---|
| 3853 | M->PBsys[j + HOMmax] = 0.0; |
---|
| 3854 | M->PBrms[j + HOMmax] = 0.0; |
---|
| 3855 | M->PBrnd[j + HOMmax] = 0.0; |
---|
| 3856 | } |
---|
| 3857 | M->Porder = 0; |
---|
| 3858 | M->n_design = 0; |
---|
| 3859 | M->Pirho = 0.0; /* inverse of curvature radius */ |
---|
| 3860 | M->PTx1 = 0.0; /* Entrance angle */ |
---|
| 3861 | M->PTx2 = 0.0; /* Exit angle */ |
---|
[11] | 3862 | M->PH1 = 0.0; /* Entrance pole face curvature*/ |
---|
| 3863 | M->PH2 = 0.0; /* Exit pole face curvature */ |
---|
[3] | 3864 | M->Pgap = 0.0; /* Gap for fringe field ??? */ |
---|
| 3865 | |
---|
| 3866 | M->Pc0 = 0.0; |
---|
| 3867 | M->Pc1 = 0.0; |
---|
| 3868 | M->Ps1 = 0.0; |
---|
| 3869 | M->quadFF1 = 0L; |
---|
| 3870 | M->quadFF2 = 0L; |
---|
| 3871 | M->sextFF1 = 0L; |
---|
| 3872 | M->sextFF2 = 0L; |
---|
| 3873 | M->quadFFscaling = 0.0; |
---|
| 3874 | |
---|
| 3875 | } |
---|
| 3876 | |
---|
| 3877 | /**************************************************************************** |
---|
| 3878 | * void Cav_Alloc(elemtype *Elem) |
---|
| 3879 | |
---|
| 3880 | Purpose: |
---|
| 3881 | Constructor for a cavity element |
---|
| 3882 | |
---|
| 3883 | Input: |
---|
| 3884 | none |
---|
| 3885 | |
---|
| 3886 | Output: |
---|
| 3887 | none |
---|
| 3888 | |
---|
| 3889 | Return: |
---|
| 3890 | none |
---|
| 3891 | |
---|
| 3892 | Global variables: |
---|
| 3893 | none |
---|
| 3894 | |
---|
| 3895 | Specific functions: |
---|
| 3896 | none |
---|
| 3897 | |
---|
| 3898 | Comments: |
---|
| 3899 | none |
---|
| 3900 | |
---|
| 3901 | ****************************************************************************/ |
---|
| 3902 | void Cav_Alloc(elemtype *Elem) { |
---|
| 3903 | CavityType *C; |
---|
| 3904 | |
---|
| 3905 | Elem->C = (CavityType *) malloc(sizeof(CavityType)); |
---|
| 3906 | C = Elem->C; |
---|
| 3907 | C->Pvolt = 0.0; |
---|
| 3908 | C->Pfreq = 0.0; |
---|
| 3909 | C->phi = 0.0; |
---|
| 3910 | C->Ph = 0; |
---|
| 3911 | } |
---|
| 3912 | |
---|
| 3913 | void Wiggler_Alloc(elemtype *Elem) { |
---|
| 3914 | int j; |
---|
| 3915 | WigglerType *W; |
---|
| 3916 | |
---|
| 3917 | Elem->W = (WigglerType *) malloc(sizeof(WigglerType)); |
---|
| 3918 | W = Elem->W; |
---|
| 3919 | W->Pmethod = Meth_Linear; |
---|
| 3920 | W->PN = 0; |
---|
| 3921 | for (j = 0; j <= 1; j++) { |
---|
| 3922 | W->PdSsys[j] = 0.0; |
---|
| 3923 | W->PdSrnd[j] = 0.0; |
---|
| 3924 | } |
---|
| 3925 | W->PdTpar = 0.0; |
---|
| 3926 | W->PdTsys = 0.0; |
---|
| 3927 | W->PdTrnd = 0.0; |
---|
| 3928 | W->n_harm = 0; |
---|
| 3929 | for (j = 0; j < n_harm_max; j++) { |
---|
| 3930 | W->BoBrhoV[j] = 0.0; |
---|
| 3931 | W->BoBrhoH[j] = 0.0; |
---|
| 3932 | W->kxV[j] = 0.0; |
---|
| 3933 | W->kxH[j] = 0.0; |
---|
| 3934 | W->lambda = 0.0; |
---|
| 3935 | W->phi[j] = 0.0; |
---|
| 3936 | } |
---|
| 3937 | for (j = 0; j <= HOMmax; j++) |
---|
| 3938 | W->PBW[j + HOMmax] = 0.0; |
---|
| 3939 | W->Porder = 0; |
---|
| 3940 | } |
---|
| 3941 | |
---|
| 3942 | void FieldMap_Alloc(elemtype *Elem, const bool alloc_fm) { |
---|
| 3943 | FieldMapType *FM; |
---|
| 3944 | |
---|
| 3945 | Elem->FM = (FieldMapType *) malloc(sizeof(FieldMapType)); |
---|
| 3946 | FM = Elem->FM; |
---|
| 3947 | FM->n_step = 0; |
---|
| 3948 | FM->n[X_] = 0; |
---|
| 3949 | FM->n[Y_] = 0; |
---|
| 3950 | FM->n[Z_] = 0; |
---|
| 3951 | FM->scl = 1.0; |
---|
| 3952 | |
---|
| 3953 | /* if (alloc_fm) { |
---|
| 3954 | FM->AxoBrho = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3955 | FM->AxoBrho2 = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3956 | FM->AyoBrho = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3957 | FM->AyoBrho2 = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3958 | |
---|
| 3959 | FM->Bx = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3960 | FM->By = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3961 | FM->Bz = matrix(1, m_max_FM, 1, n_max_FM); |
---|
| 3962 | |
---|
| 3963 | FM->x_pos[1] = 1e30; FM->x_pos[FM->m_x] = -1e30; |
---|
| 3964 | FM->y_pos[1] = 1e30; FM->y_pos[FM->m_y] = -1e30; |
---|
| 3965 | FM->s_pos[1] = 1e30; FM->s_pos[FM->n_s] = -1e30; |
---|
| 3966 | }*/ |
---|
| 3967 | |
---|
| 3968 | // free_vector(FM->x_pos, 1, m_max_FM); free_vector(FM->y_pos, 1, m_max_FM); |
---|
| 3969 | // free_vector(FM->s_pos, 1, n_max_FM); |
---|
| 3970 | // free_matrix(FM->AxoBrho, 1, m_max_FM, 1, n_max_FM); |
---|
| 3971 | // free_matrix(FM->AxoBrho2, 1, m_max_FM, 1, n_max_FM); |
---|
| 3972 | |
---|
| 3973 | // free_matrix(Bx, 1, m_max_FM, 1, n_max_FM); |
---|
| 3974 | // free_matrix(By, 1, m_max_FM, 1, n_max_FM); |
---|
| 3975 | // free_matrix(Bz, 1, m_max_FM, 1, n_max_FM); |
---|
| 3976 | } |
---|
| 3977 | |
---|
| 3978 | /**************************************************************************** |
---|
| 3979 | * void Insertion_Alloc(elemtype *Elem, boolean firstflag, boolean secondflag) |
---|
| 3980 | |
---|
| 3981 | Purpose: called by Insertion_Init and Lat_DealElement |
---|
| 3982 | Dynamic memory allocation for a Insertion element and various |
---|
| 3983 | initializations |
---|
| 3984 | |
---|
| 3985 | Input: |
---|
| 3986 | Elem Element for memory allocation |
---|
| 3987 | firstflag true if first order kick map to be loaded |
---|
| 3988 | secondflag true if second order kick map to be loaded |
---|
| 3989 | |
---|
| 3990 | Output: |
---|
| 3991 | none |
---|
| 3992 | |
---|
| 3993 | Return: |
---|
| 3994 | none |
---|
| 3995 | |
---|
| 3996 | Global variables: |
---|
| 3997 | none |
---|
| 3998 | |
---|
| 3999 | Specific functions: |
---|
| 4000 | none |
---|
| 4001 | |
---|
| 4002 | Comments: |
---|
| 4003 | 10/01/05 First order kick part added |
---|
| 4004 | 4 November 2010 Splitting 1st and 2nd order X and Z axes |
---|
| 4005 | |
---|
| 4006 | ****************************************************************************/ |
---|
| 4007 | |
---|
| 4008 | void Insertion_Alloc(elemtype *Elem) { |
---|
| 4009 | int i = 0, j = 0; |
---|
| 4010 | InsertionType *ID; |
---|
| 4011 | |
---|
| 4012 | Elem->ID = (InsertionType *) malloc(sizeof(InsertionType)); |
---|
| 4013 | ID = Elem->ID; |
---|
| 4014 | |
---|
| 4015 | ID->Pmethod = Meth_Linear; |
---|
| 4016 | ID->PN = 0; |
---|
| 4017 | ID->nx1 = 0; |
---|
| 4018 | ID->nz1 = 0; |
---|
| 4019 | ID->nx2 = 0; |
---|
| 4020 | ID->nz2 = 0; |
---|
| 4021 | |
---|
| 4022 | /* Initialization thetax and thetaz to 0*/ |
---|
| 4023 | |
---|
| 4024 | // first order kick map |
---|
| 4025 | if (ID->firstorder) { |
---|
| 4026 | for (i = 0; i < IDZMAX; i++) { |
---|
| 4027 | for (j = 0; j < IDXMAX; j++) { |
---|
| 4028 | ID->thetax1[i][j] = 0.0; |
---|
| 4029 | ID->thetaz1[i][j] = 0.0; |
---|
| 4030 | } |
---|
| 4031 | } |
---|
| 4032 | } |
---|
| 4033 | |
---|
| 4034 | // second order kick map |
---|
| 4035 | if (ID->secondorder) { |
---|
| 4036 | for (i = 0; i < IDZMAX; i++) { |
---|
| 4037 | for (j = 0; j < IDXMAX; j++) { |
---|
| 4038 | ID->thetax2[i][j] = 0.0; |
---|
| 4039 | ID->thetaz2[i][j] = 0.0; |
---|
| 4040 | } |
---|
| 4041 | } |
---|
| 4042 | } |
---|
| 4043 | |
---|
| 4044 | // stuffs for interpolation |
---|
| 4045 | for (j = 0; j < IDXMAX; j++) { |
---|
| 4046 | ID->tabx1[j] = 0.0; |
---|
| 4047 | ID->tabx2[j] = 0.0; |
---|
| 4048 | } |
---|
| 4049 | |
---|
| 4050 | for (j = 0; j < IDZMAX; j++) { |
---|
| 4051 | ID->tabz1[j] = 0.0; |
---|
| 4052 | ID->tabz2[j] = 0.0; |
---|
| 4053 | } |
---|
| 4054 | |
---|
| 4055 | // filenames |
---|
| 4056 | strcpy(ID->fname1, ""); |
---|
| 4057 | strcpy(ID->fname2, ""); |
---|
| 4058 | |
---|
| 4059 | for (j = 0; j <= 1; j++) { |
---|
| 4060 | ID->PdSsys[j] = 0.0; |
---|
| 4061 | ID->PdSrnd[j] = 0.0; |
---|
| 4062 | } |
---|
| 4063 | ID->PdTpar = 0.0; |
---|
| 4064 | ID->PdTsys = 0.0; |
---|
| 4065 | ID->PdTrnd = 0.0; |
---|
| 4066 | ID->Porder = 0; |
---|
| 4067 | } |
---|
| 4068 | |
---|
| 4069 | void Spreader_Alloc(elemtype *Elem) { |
---|
| 4070 | int k; |
---|
| 4071 | |
---|
| 4072 | Elem->Spr = (SpreaderType *) malloc(sizeof(SpreaderType)); |
---|
| 4073 | |
---|
| 4074 | for (k = 0; k < Spreader_max; k++) |
---|
| 4075 | Elem->Spr->Cell_ptrs[k] = NULL; |
---|
| 4076 | } |
---|
| 4077 | |
---|
| 4078 | void Recombiner_Alloc(elemtype *Elem) { |
---|
| 4079 | Elem->Rec = (RecombinerType *) malloc(sizeof(RecombinerType)); |
---|
| 4080 | } |
---|
| 4081 | |
---|
| 4082 | void Solenoid_Alloc(elemtype *Elem) { |
---|
| 4083 | int j; |
---|
| 4084 | SolenoidType *Sol; |
---|
| 4085 | |
---|
| 4086 | Elem->Sol = (SolenoidType *) malloc(sizeof(SolenoidType)); |
---|
| 4087 | Sol = Elem->Sol; |
---|
| 4088 | Sol->N = 0; |
---|
| 4089 | for (j = 0; j <= 1; j++) { |
---|
| 4090 | Sol->PdSsys[j] = 0.0; |
---|
| 4091 | Sol->PdSrms[j] = 0.0; |
---|
| 4092 | Sol->PdSrnd[j] = 0.0; |
---|
| 4093 | } |
---|
| 4094 | Sol->dTpar = 0.0; |
---|
| 4095 | Sol->dTsys = 0.0; |
---|
| 4096 | Sol->dTrnd = 0.0; |
---|
| 4097 | } |
---|
| 4098 | |
---|
| 4099 | /****************************************************************************/ |
---|
| 4100 | /* void Drift_Init(long Fnum1) |
---|
| 4101 | |
---|
| 4102 | Purpose: |
---|
| 4103 | Constructor of a drift element |
---|
| 4104 | see comments in Drift_SetMatrix |
---|
| 4105 | |
---|
| 4106 | Input: |
---|
| 4107 | Fnum1 Family number |
---|
| 4108 | |
---|
| 4109 | Output: |
---|
| 4110 | none |
---|
| 4111 | |
---|
| 4112 | Return: |
---|
| 4113 | none |
---|
| 4114 | |
---|
| 4115 | Global variables: |
---|
| 4116 | ElemFam |
---|
| 4117 | Cell |
---|
| 4118 | |
---|
| 4119 | Specific functions: |
---|
| 4120 | Drift_Alloc |
---|
| 4121 | Drift_SetMatrix |
---|
| 4122 | |
---|
| 4123 | Comments: |
---|
| 4124 | none |
---|
| 4125 | |
---|
| 4126 | ****************************************************************************/ |
---|
| 4127 | void Drift_Init(int Fnum1) { |
---|
| 4128 | int i; |
---|
| 4129 | ElemFamType *elemfamp; |
---|
| 4130 | CellType *cellp; |
---|
| 4131 | |
---|
| 4132 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4133 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4134 | /* Get in Cell kid # i from Family Fnum1 */ |
---|
| 4135 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4136 | /* Dynamic memory allocation for element */ |
---|
| 4137 | Drift_Alloc(&cellp->Elem); |
---|
| 4138 | /* copy low level routine */ |
---|
| 4139 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4140 | /* Set the drift length */ |
---|
| 4141 | cellp->Elem.PL = elemfamp->ElemF.PL; |
---|
| 4142 | /* set the kind of element */ |
---|
| 4143 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4144 | /* set pointer for the D dynamic space */ |
---|
| 4145 | *cellp->Elem.D = *elemfamp->ElemF.D; |
---|
| 4146 | cellp->dT[0] = 1e0; /* cos = 1 */ |
---|
| 4147 | cellp->dT[1] = 0.0; /* sin = 0 */ |
---|
| 4148 | cellp->dS[0] = 0.0; /* no H displacement */ |
---|
| 4149 | cellp->dS[1] = 0.0; /* no V displacement */ |
---|
| 4150 | /* set Drift matrix */ |
---|
| 4151 | Drift_SetMatrix(Fnum1, i); |
---|
| 4152 | } |
---|
| 4153 | } |
---|
| 4154 | |
---|
| 4155 | static int UpdatePorder(elemtype &Elem) { |
---|
| 4156 | int i, order; |
---|
| 4157 | MpoleType *M; |
---|
| 4158 | |
---|
| 4159 | M = Elem.M; |
---|
| 4160 | if (M->Pirho != 0.0) /* non zero curvature => bend */ |
---|
| 4161 | order = 1; |
---|
| 4162 | else |
---|
| 4163 | /* mutipole */ |
---|
| 4164 | order = 0; |
---|
| 4165 | for (i = -HOMmax; i <= HOMmax; i++) |
---|
| 4166 | if (M->PB[i + HOMmax] != 0.0 && abs(i) > order) |
---|
| 4167 | order = abs(i); |
---|
| 4168 | |
---|
| 4169 | return order; |
---|
| 4170 | } |
---|
| 4171 | |
---|
| 4172 | void Mpole_Init(int Fnum1) { |
---|
| 4173 | double x; |
---|
| 4174 | int i; |
---|
| 4175 | ElemFamType *elemfamp; |
---|
| 4176 | CellType *cellp; |
---|
| 4177 | elemtype *elemp; |
---|
| 4178 | |
---|
| 4179 | /* Pointer on element */ |
---|
| 4180 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4181 | memcpy(elemfamp->ElemF.M->PB, elemfamp->ElemF.M->PBpar, sizeof(mpolArray)); |
---|
| 4182 | /* Update the right multipole order */ |
---|
| 4183 | elemfamp->ElemF.M->Porder = UpdatePorder(elemfamp->ElemF); |
---|
| 4184 | /* Quadrupole strength */ |
---|
| 4185 | x = elemfamp->ElemF.M->PBpar[Quad + HOMmax]; |
---|
| 4186 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4187 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4188 | /* Memory allocation and set everything to zero */ |
---|
| 4189 | Mpole_Alloc(&cellp->Elem); |
---|
| 4190 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4191 | /* set length */ |
---|
| 4192 | cellp->Elem.PL = elemfamp->ElemF.PL; |
---|
| 4193 | /* set element kind (Mpole) */ |
---|
| 4194 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4195 | *cellp->Elem.M = *elemfamp->ElemF.M; |
---|
| 4196 | |
---|
| 4197 | elemp = &cellp->Elem; |
---|
| 4198 | /* set entrance and exit angles */ |
---|
| 4199 | cellp->dT[0] = cos(dtor(elemp->M->PdTpar)); |
---|
| 4200 | cellp->dT[1] = sin(dtor(elemp->M->PdTpar)); |
---|
| 4201 | |
---|
| 4202 | /* set displacement to zero */ |
---|
| 4203 | cellp->dS[0] = 0.0; |
---|
| 4204 | cellp->dS[1] = 0.0; |
---|
| 4205 | |
---|
| 4206 | if (elemp->PL != 0.0 || elemp->M->Pirho != 0.0) { |
---|
| 4207 | /* Thick element or radius non zero element */ |
---|
| 4208 | elemp->M->Pthick = pthicktype(thick); |
---|
| 4209 | /* sin(L*irho/2) =sin(theta/2) half the angle */ |
---|
| 4210 | elemp->M->Pc0 = sin(elemp->PL * elemp->M->Pirho / 2e0); |
---|
| 4211 | /* cos roll: sin(theta/2)*cos(dT) */ |
---|
| 4212 | elemp->M->Pc1 = cellp->dT[0] * elemp->M->Pc0; |
---|
| 4213 | /* sin roll: sin(theta/2)*cos(dT) */ |
---|
| 4214 | elemp->M->Ps1 = cellp->dT[1] * elemp->M->Pc0; |
---|
| 4215 | Mpole_Setmatrix(Fnum1, i, x); |
---|
| 4216 | } else |
---|
| 4217 | /* element as thin lens */ |
---|
| 4218 | elemp->M->Pthick = pthicktype(thin); |
---|
| 4219 | } |
---|
| 4220 | } |
---|
| 4221 | |
---|
| 4222 | #define order 2 |
---|
| 4223 | void Wiggler_Init(int Fnum1) { |
---|
| 4224 | int i; |
---|
| 4225 | double x; |
---|
| 4226 | ElemFamType *elemfamp; |
---|
| 4227 | CellType *cellp; |
---|
| 4228 | elemtype *elemp; |
---|
| 4229 | |
---|
| 4230 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4231 | /* ElemF.M^.PB := ElemF.M^.PBpar; */ |
---|
| 4232 | elemfamp->ElemF.W->Porder = order; |
---|
| 4233 | x = elemfamp->ElemF.W->PBW[Quad + HOMmax]; |
---|
| 4234 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4235 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4236 | Wiggler_Alloc(&cellp->Elem); |
---|
| 4237 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4238 | cellp->Elem.PL = elemfamp->ElemF.PL; |
---|
| 4239 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4240 | *cellp->Elem.W = *elemfamp->ElemF.W; |
---|
| 4241 | |
---|
| 4242 | elemp = &cellp->Elem; |
---|
| 4243 | cellp->dT[0] = cos(dtor(elemp->M->PdTpar)); |
---|
| 4244 | cellp->dT[1] = sin(dtor(elemp->M->PdTpar)); |
---|
| 4245 | |
---|
| 4246 | cellp->dS[0] = 0.0; |
---|
| 4247 | cellp->dS[1] = 0.0; |
---|
| 4248 | Wiggler_Setmatrix(Fnum1, i, cellp->Elem.PL, cellp->Elem.W->kxV[0], 2.0 |
---|
| 4249 | * M_PI / cellp->Elem.W->lambda, x); |
---|
| 4250 | } |
---|
| 4251 | } |
---|
| 4252 | #undef order |
---|
| 4253 | |
---|
| 4254 | /* |
---|
| 4255 | void get_Ax(const int m, const int n, float **By, FieldMapType *FM) |
---|
| 4256 | { |
---|
| 4257 | int j, k; |
---|
| 4258 | |
---|
| 4259 | const double Brho = 1e9*globval.Energy/c0; |
---|
| 4260 | |
---|
| 4261 | FM->m_y = m; FM->n_s = n; |
---|
| 4262 | |
---|
| 4263 | for (j = 1; j <= m; j++) { |
---|
| 4264 | FM->AxoBrho[j][1] = 0.0; |
---|
| 4265 | for (k = 2; k <= n; k++) |
---|
| 4266 | FM->AxoBrho[j][k] = |
---|
| 4267 | FM->AxoBrho[j][k-1] + By[j][k]*(FM->s_pos[k]-FM->s_pos[k-1])/Brho; |
---|
| 4268 | } |
---|
| 4269 | |
---|
| 4270 | splie2(FM->y_pos, FM->s_pos, FM->AxoBrho, FM->m_y, FM->n_s, FM->AxoBrho2); |
---|
| 4271 | } |
---|
| 4272 | |
---|
| 4273 | |
---|
| 4274 | void get_Ay(const int m, const int n, float **Bx, FieldMapType *FM) |
---|
| 4275 | { |
---|
| 4276 | int j, k; |
---|
| 4277 | |
---|
| 4278 | const double Brho = 1e9*globval.Energy/c0; |
---|
| 4279 | |
---|
| 4280 | FM->m_x = m; FM->n_s = n; |
---|
| 4281 | |
---|
| 4282 | for (j = 1; j <= m; j++) { |
---|
| 4283 | FM->AyoBrho[j][1] = 0.0; |
---|
| 4284 | for (k = 2; k <= n; k++) |
---|
| 4285 | FM->AyoBrho[j][k] = |
---|
| 4286 | FM->AyoBrho[j][k-1] - Bx[j][k]*(FM->s_pos[k]-FM->s_pos[k-1])/Brho; |
---|
| 4287 | } |
---|
| 4288 | |
---|
| 4289 | splie2(FM->x_pos, FM->s_pos, FM->AyoBrho, FM->m_x, FM->n_s, FM->AyoBrho2); |
---|
| 4290 | } |
---|
| 4291 | */ |
---|
| 4292 | |
---|
| 4293 | void get_B(const char *file_name, FieldMapType *FM) { |
---|
| 4294 | char line[max_str]; |
---|
| 4295 | int i, j, k, l; |
---|
| 4296 | ifstream inf; |
---|
| 4297 | |
---|
| 4298 | printf("\n"); |
---|
| 4299 | printf("reading field map %s\n", file_name); |
---|
| 4300 | |
---|
| 4301 | file_rd(inf, file_name); |
---|
| 4302 | |
---|
| 4303 | inf.getline(line, max_str); |
---|
| 4304 | // read number of steps |
---|
| 4305 | sscanf(line, "%d,%d,%d", &FM->n[X_], &FM->n[Y_], &FM->n[Z_]); |
---|
| 4306 | // skip comment |
---|
| 4307 | inf.getline(line, max_str); |
---|
| 4308 | |
---|
| 4309 | i = 1; |
---|
| 4310 | j = 0; |
---|
| 4311 | k = 1; |
---|
| 4312 | while (inf.getline(line, max_str) != NULL) { |
---|
| 4313 | j++; |
---|
| 4314 | if (j > FM->n[Y_]) { |
---|
| 4315 | k++; |
---|
| 4316 | j = 1; |
---|
| 4317 | } |
---|
| 4318 | if (k > FM->n[Z_]) { |
---|
| 4319 | i++; |
---|
| 4320 | k = 1; |
---|
| 4321 | } |
---|
| 4322 | |
---|
| 4323 | if ((i > i_max_FM) || (j > j_max_FM) || (k > k_max_FM)) { |
---|
| 4324 | printf("get_B: max index exceeded %d %d %d (%d %d %d)\n", i, j, k, |
---|
| 4325 | i_max_FM, j_max_FM, k_max_FM); |
---|
| 4326 | exit_(1); |
---|
| 4327 | } |
---|
| 4328 | |
---|
| 4329 | sscanf(line, "%lf,%lf,%lf,%lf,%lf,%lf", &FM->xyz[X_][i - 1][j - 1][k |
---|
| 4330 | - 1], &FM->xyz[Y_][i - 1][j - 1][k - 1], &FM->xyz[Z_][i - 1][j |
---|
| 4331 | - 1][k - 1], &FM->B[X_][i - 1][j - 1][k - 1], |
---|
| 4332 | &FM->B[Y_][i - 1][j - 1][k - 1], |
---|
| 4333 | &FM->B[Z_][i - 1][j - 1][k - 1]); |
---|
| 4334 | for (l = 0; l < 3; l++) |
---|
| 4335 | FM->xyz[l][i - 1][j - 1][k - 1] *= 1e-3; |
---|
| 4336 | |
---|
| 4337 | } |
---|
| 4338 | |
---|
| 4339 | printf("no of steps: n_x = %d, n_y = %d, n_z = %d\n", FM->n[X_], FM->n[Y_], |
---|
| 4340 | FM->n[Z_]); |
---|
| 4341 | |
---|
| 4342 | // get_Ay(m, n, FM->Bx, FM); prt_Bx(FM); |
---|
| 4343 | } |
---|
| 4344 | |
---|
| 4345 | void FieldMap_Init(int Fnum1) { |
---|
| 4346 | int i; |
---|
| 4347 | ElemFamType *elemfamp; |
---|
| 4348 | CellType *cellp; |
---|
| 4349 | elemtype *elemp; |
---|
| 4350 | |
---|
| 4351 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4352 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4353 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4354 | FieldMap_Alloc(&cellp->Elem, false); |
---|
| 4355 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4356 | cellp->Elem.PL = elemfamp->ElemF.PL; |
---|
| 4357 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4358 | *cellp->Elem.FM = *elemfamp->ElemF.FM; |
---|
| 4359 | |
---|
| 4360 | elemp = &cellp->Elem; |
---|
| 4361 | cellp->dT[0] = 1.0; |
---|
| 4362 | cellp->dT[1] = 0.0; |
---|
| 4363 | cellp->dS[X_] = 0.0; |
---|
| 4364 | cellp->dS[Y_] = 0.0; |
---|
| 4365 | } |
---|
| 4366 | } |
---|
| 4367 | |
---|
| 4368 | /**************************************************************************** |
---|
| 4369 | * void Cav_Init(long Fnum1) |
---|
| 4370 | |
---|
| 4371 | Purpose: called by Cell_Init() |
---|
| 4372 | Constructor for a cavity |
---|
| 4373 | set the RF voltage, frequency read from the lattice file |
---|
| 4374 | |
---|
| 4375 | Input: |
---|
| 4376 | Fnum1 Family number |
---|
| 4377 | |
---|
| 4378 | Output: |
---|
| 4379 | none |
---|
| 4380 | |
---|
| 4381 | Return: |
---|
| 4382 | none |
---|
| 4383 | |
---|
| 4384 | Global variables: |
---|
| 4385 | ElemFam, Cell |
---|
| 4386 | |
---|
| 4387 | Specific functions: |
---|
| 4388 | none |
---|
| 4389 | |
---|
| 4390 | Comments: |
---|
| 4391 | none |
---|
| 4392 | |
---|
| 4393 | ****************************************************************************/ |
---|
| 4394 | void Cav_Init(int Fnum1) { |
---|
| 4395 | int i; |
---|
| 4396 | ElemFamType *elemfamp; |
---|
| 4397 | elemtype *elemp; |
---|
| 4398 | CellType *cellp; |
---|
| 4399 | |
---|
| 4400 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4401 | elemp = &elemfamp->ElemF; |
---|
| 4402 | for (i = 0; i < elemfamp->nKid; i++) { |
---|
| 4403 | cellp = &Cell[elemfamp->KidList[i]]; |
---|
| 4404 | cellp->Elem = elemfamp->ElemF; |
---|
| 4405 | } |
---|
| 4406 | } |
---|
| 4407 | |
---|
| 4408 | void Marker_Init(int Fnum1) { |
---|
| 4409 | int i; |
---|
| 4410 | ElemFamType *elemfamp; |
---|
| 4411 | elemtype *elemp; |
---|
| 4412 | CellType *cellp; |
---|
| 4413 | |
---|
| 4414 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4415 | elemp = &elemfamp->ElemF; |
---|
| 4416 | for (i = 0; i < elemfamp->nKid; i++) { |
---|
| 4417 | cellp = &Cell[elemfamp->KidList[i]]; |
---|
| 4418 | cellp->Elem = elemfamp->ElemF; |
---|
| 4419 | cellp->dT[0] = 1.0; |
---|
| 4420 | cellp->dT[1] = 0.0; |
---|
| 4421 | cellp->dS[0] = 0.0; |
---|
| 4422 | cellp->dS[1] = 0.0; |
---|
| 4423 | } |
---|
| 4424 | } |
---|
| 4425 | |
---|
| 4426 | /**************************************************************************** |
---|
| 4427 | * INSERTION * |
---|
| 4428 | ****************************************************************************/ |
---|
| 4429 | |
---|
| 4430 | /**************************************************************************** |
---|
| 4431 | * void Insertion_Init(long Fnum1) |
---|
| 4432 | |
---|
| 4433 | Purpose: called by Cell_Init |
---|
| 4434 | Initializes the insertion |
---|
| 4435 | Fills in all the parameters read in the RADIA file |
---|
| 4436 | Constructs the linear matrix |
---|
| 4437 | |
---|
| 4438 | Input: |
---|
| 4439 | Fnum1: family number |
---|
| 4440 | |
---|
| 4441 | Output: |
---|
| 4442 | none |
---|
| 4443 | |
---|
| 4444 | Return: |
---|
| 4445 | none |
---|
| 4446 | |
---|
| 4447 | Global variables: |
---|
| 4448 | none |
---|
| 4449 | |
---|
| 4450 | Specific functions: |
---|
| 4451 | none |
---|
| 4452 | |
---|
| 4453 | Comments: |
---|
| 4454 | none |
---|
| 4455 | |
---|
| 4456 | ****************************************************************************/ |
---|
| 4457 | |
---|
| 4458 | void Insertion_Init(int Fnum1) { |
---|
| 4459 | int i; |
---|
| 4460 | ElemFamType *elemfamp; |
---|
| 4461 | CellType *cellp; |
---|
| 4462 | elemtype *elemp; |
---|
| 4463 | |
---|
| 4464 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4465 | // elemfamp->ElemF.ID->Porder = order; |
---|
| 4466 | // x = elemfamp->ElemF.ID->PBW[Quad + HOMmax]; |
---|
| 4467 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4468 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4469 | Insertion_Alloc(&cellp->Elem); |
---|
| 4470 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4471 | cellp->Elem.PL = elemfamp->ElemF.PL; |
---|
| 4472 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4473 | *cellp->Elem.ID = *elemfamp->ElemF.ID; |
---|
| 4474 | |
---|
| 4475 | elemp = &cellp->Elem; |
---|
| 4476 | |
---|
| 4477 | cellp->dT[0] = cos(dtor(elemp->ID->PdTpar)); |
---|
| 4478 | cellp->dT[1] = sin(dtor(elemp->ID->PdTpar)); |
---|
| 4479 | cellp->dS[0] = 0.0; |
---|
| 4480 | cellp->dS[1] = 0.0; |
---|
| 4481 | |
---|
| 4482 | Insertion_SetMatrix(Fnum1, i); |
---|
| 4483 | } |
---|
| 4484 | } |
---|
| 4485 | |
---|
| 4486 | void Spreader_Init(int Fnum1) { |
---|
| 4487 | int i; |
---|
| 4488 | ElemFamType *elemfamp; |
---|
| 4489 | CellType *cellp; |
---|
| 4490 | |
---|
| 4491 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4492 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4493 | /* Get in Cell kid # i from Family Fnum1 */ |
---|
| 4494 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4495 | /* Dynamic memory allocation for element */ |
---|
| 4496 | Spreader_Alloc(&cellp->Elem); |
---|
| 4497 | /* copy low level routine */ |
---|
| 4498 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4499 | /* set the kind of element */ |
---|
| 4500 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4501 | /* set pointer for the dynamic space */ |
---|
| 4502 | *cellp->Elem.Spr = *elemfamp->ElemF.Spr; |
---|
| 4503 | cellp->dT[0] = 1e0; /* cos = 1 */ |
---|
| 4504 | cellp->dT[1] = 0.0; /* sin = 0 */ |
---|
| 4505 | cellp->dS[0] = 0.0; /* no H displacement */ |
---|
| 4506 | cellp->dS[1] = 0.0; /* no V displacement */ |
---|
| 4507 | } |
---|
| 4508 | } |
---|
| 4509 | |
---|
| 4510 | void Recombiner_Init(int Fnum1) { |
---|
| 4511 | int i; |
---|
| 4512 | ElemFamType *elemfamp; |
---|
| 4513 | CellType *cellp; |
---|
| 4514 | |
---|
| 4515 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4516 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4517 | /* Get in Cell kid # i from Family Fnum1 */ |
---|
| 4518 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4519 | /* Dynamic memory allocation for element */ |
---|
| 4520 | Spreader_Alloc(&cellp->Elem); |
---|
| 4521 | /* copy low level routine */ |
---|
| 4522 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4523 | /* set the kind of element */ |
---|
| 4524 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4525 | /* set pointer for the dynamic space */ |
---|
| 4526 | *cellp->Elem.Rec = *elemfamp->ElemF.Rec; |
---|
| 4527 | cellp->dT[0] = 1e0; /* cos = 1 */ |
---|
| 4528 | cellp->dT[1] = 0.0; /* sin = 0 */ |
---|
| 4529 | cellp->dS[0] = 0.0; /* no H displacement */ |
---|
| 4530 | cellp->dS[1] = 0.0; /* no V displacement */ |
---|
| 4531 | } |
---|
| 4532 | } |
---|
| 4533 | |
---|
| 4534 | void Solenoid_Init(int Fnum1) { |
---|
| 4535 | int i; |
---|
| 4536 | ElemFamType *elemfamp; |
---|
| 4537 | CellType *cellp; |
---|
| 4538 | elemtype *elemp; |
---|
| 4539 | |
---|
| 4540 | /* Pointer on element */ |
---|
| 4541 | elemfamp = &ElemFam[Fnum1 - 1]; |
---|
| 4542 | for (i = 1; i <= elemfamp->nKid; i++) { |
---|
| 4543 | cellp = &Cell[elemfamp->KidList[i - 1]]; |
---|
| 4544 | /* Memory allocation and set everything to zero */ |
---|
| 4545 | Solenoid_Alloc(&cellp->Elem); |
---|
| 4546 | memcpy(cellp->Elem.PName, elemfamp->ElemF.PName, sizeof(partsName)); |
---|
| 4547 | /* set length */ |
---|
| 4548 | cellp->Elem.PL = elemfamp->ElemF.PL; |
---|
| 4549 | /* set element kind */ |
---|
| 4550 | cellp->Elem.Pkind = elemfamp->ElemF.Pkind; |
---|
| 4551 | *cellp->Elem.Sol = *elemfamp->ElemF.Sol; |
---|
| 4552 | |
---|
| 4553 | elemp = &cellp->Elem; |
---|
| 4554 | /* set entrance and exit angles */ |
---|
| 4555 | cellp->dT[0] = 1.0; |
---|
| 4556 | cellp->dT[1] = 0.0; |
---|
| 4557 | |
---|
| 4558 | /* set displacement to zero */ |
---|
| 4559 | cellp->dS[0] = 0.0; |
---|
| 4560 | cellp->dS[1] = 0.0; |
---|
| 4561 | } |
---|
| 4562 | } |
---|
| 4563 | |
---|
| 4564 | /************************************************************************************** |
---|
| 4565 | void Mpole_SetPB(int Fnum1, int Knum1, int Order) |
---|
| 4566 | |
---|
| 4567 | Purpose: |
---|
| 4568 | called by Cell_SetdP |
---|
| 4569 | Update full multipole composent as sum of design, systematic |
---|
| 4570 | and random part; and update the maximum order of the multipole |
---|
| 4571 | component p_order. |
---|
| 4572 | The ramdom error is the multiplication of PBrms and PBrnd |
---|
| 4573 | Compute transport matrix if quadrupole (Order=2) |
---|
| 4574 | Set multipole order to Order if multipole (Order >2) |
---|
| 4575 | |
---|
| 4576 | Input: |
---|
| 4577 | Fnum1 family name |
---|
| 4578 | Knum1 kid number |
---|
| 4579 | Order maximum order of the multipole |
---|
| 4580 | |
---|
| 4581 | Output: |
---|
| 4582 | None |
---|
| 4583 | |
---|
| 4584 | Return: |
---|
| 4585 | None |
---|
| 4586 | |
---|
| 4587 | Gloval variables: |
---|
| 4588 | None |
---|
| 4589 | |
---|
| 4590 | Specific functions: |
---|
| 4591 | |
---|
| 4592 | Comments: |
---|
| 4593 | None |
---|
| 4594 | |
---|
| 4595 | **************************************************************************************/ |
---|
| 4596 | void Mpole_SetPB(int Fnum1, int Knum1, int Order) { |
---|
| 4597 | /* */ |
---|
| 4598 | |
---|
| 4599 | CellType *cellp; /* pointer on the Cell */ |
---|
| 4600 | elemtype *elemp; /* pointer on the Elemetype */ |
---|
| 4601 | MpoleType *M;/* Pointer on the Multipole */ |
---|
| 4602 | |
---|
| 4603 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 4604 | elemp = &cellp->Elem; |
---|
| 4605 | M = elemp->M; |
---|
| 4606 | M->PB[Order + HOMmax] = M->PBpar[Order + HOMmax] + M->PBsys[Order + HOMmax] |
---|
| 4607 | + M->PBrms[Order + HOMmax] * M->PBrnd[Order + HOMmax]; |
---|
| 4608 | if (abs(Order) > M->Porder && M->PB[Order + HOMmax] != 0.0) |
---|
| 4609 | M->Porder = abs(Order); |
---|
| 4610 | if (M->Pmethod == Meth_Linear && Order == 2L) |
---|
| 4611 | Mpole_Setmatrix(Fnum1, Knum1, M->PB[Order + HOMmax]); |
---|
| 4612 | cellconcat = false; |
---|
| 4613 | } |
---|
| 4614 | |
---|
| 4615 | /************************************************************************************** |
---|
| 4616 | double Mpole_GetPB(int Fnum1, int Knum1, int Order) |
---|
| 4617 | |
---|
| 4618 | Purpose: |
---|
| 4619 | Return multipole strength (of order Order) for Knum1 element of |
---|
| 4620 | family Fnum1 |
---|
| 4621 | Order = 2 for normal quadrupole, bn components |
---|
| 4622 | = -2 for skew quadrupole an components |
---|
| 4623 | |
---|
| 4624 | Input: |
---|
| 4625 | Fnum1 family name |
---|
| 4626 | Knum1 kid number |
---|
| 4627 | Order order of the multipole |
---|
| 4628 | |
---|
| 4629 | Output: |
---|
| 4630 | None |
---|
| 4631 | |
---|
| 4632 | Return: |
---|
| 4633 | None |
---|
| 4634 | |
---|
| 4635 | Gloval variables: |
---|
| 4636 | None |
---|
| 4637 | |
---|
| 4638 | Specific functions: |
---|
| 4639 | |
---|
| 4640 | Comments: |
---|
| 4641 | None |
---|
| 4642 | |
---|
| 4643 | **************************************************************************************/ |
---|
| 4644 | double Mpole_GetPB(int Fnum1, int Knum1, int Order) { |
---|
| 4645 | |
---|
| 4646 | MpoleType *M; /* Pointer on the multipole */ |
---|
| 4647 | |
---|
| 4648 | M = Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem.M; |
---|
| 4649 | return (M->PB[Order + HOMmax]); |
---|
| 4650 | } |
---|
| 4651 | |
---|
| 4652 | void Mpole_DefPBpar(int Fnum1, int Knum1, int Order, double PBpar) { |
---|
| 4653 | elemtype *elemp; |
---|
| 4654 | MpoleType *M; |
---|
| 4655 | |
---|
| 4656 | elemp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem; |
---|
| 4657 | M = elemp->M; |
---|
| 4658 | |
---|
| 4659 | M->PBpar[Order + HOMmax] = PBpar; |
---|
| 4660 | } |
---|
| 4661 | |
---|
| 4662 | |
---|
| 4663 | void Mpole_DefPBsys(int Fnum1, int Knum1, int Order, double PBsys) { |
---|
| 4664 | /*Fnum1, Knum1, Order : integer*/ |
---|
| 4665 | elemtype *elemp; |
---|
| 4666 | MpoleType *M; |
---|
| 4667 | |
---|
| 4668 | elemp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem; |
---|
| 4669 | M = elemp->M; |
---|
| 4670 | |
---|
| 4671 | M->PBsys[Order + HOMmax] = PBsys; |
---|
| 4672 | } |
---|
| 4673 | /********************************************************************* |
---|
| 4674 | void Mpole_SetdS(int Fnum1, int Knum1) |
---|
| 4675 | |
---|
| 4676 | Purpose: |
---|
| 4677 | Set misalignment error to the element with Fnum and Knum |
---|
| 4678 | |
---|
| 4679 | Input: |
---|
| 4680 | Fnum1 family number |
---|
| 4681 | Knum1 kid number |
---|
| 4682 | |
---|
| 4683 | |
---|
| 4684 | **********************************************************************/ |
---|
| 4685 | void Mpole_SetdS(int Fnum1, int Knum1) { |
---|
| 4686 | int j; |
---|
| 4687 | CellType *cellp; |
---|
| 4688 | elemtype *elemp; |
---|
| 4689 | MpoleType *M; |
---|
| 4690 | |
---|
| 4691 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 4692 | elemp = &cellp->Elem; |
---|
| 4693 | M = elemp->M; |
---|
| 4694 | for (j = 0; j <= 1; j++) |
---|
| 4695 | cellp->dS[j] = M->PdSsys[j] + M->PdSrms[j] * M->PdSrnd[j]; |
---|
| 4696 | cellconcat = false; |
---|
| 4697 | } |
---|
| 4698 | |
---|
| 4699 | /********************************************************************* |
---|
| 4700 | void Mpole_SetdT(int Fnum1, int Knum1) |
---|
| 4701 | |
---|
| 4702 | Purpose: |
---|
| 4703 | Set rotation error to the element with Fnum and Knum |
---|
| 4704 | |
---|
| 4705 | Input: |
---|
| 4706 | Fnum1 family number |
---|
| 4707 | Knum1 kid number |
---|
| 4708 | |
---|
| 4709 | |
---|
| 4710 | **********************************************************************/ |
---|
| 4711 | void Mpole_SetdT(int Fnum1, int Knum1) { |
---|
| 4712 | CellType *cellp; |
---|
| 4713 | elemtype *elemp; |
---|
| 4714 | MpoleType *M; |
---|
| 4715 | |
---|
| 4716 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 4717 | elemp = &cellp->Elem; |
---|
| 4718 | M = elemp->M; |
---|
| 4719 | cellp->dT[0] = cos(dtor(M->PdTpar + M->PdTsys + M->PdTrms * M->PdTrnd)); |
---|
| 4720 | cellp->dT[1] = sin(dtor(M->PdTpar + M->PdTsys + M->PdTrms * M->PdTrnd)); |
---|
| 4721 | /* Calculate simplified p_rots */ |
---|
| 4722 | M->Pc0 = sin(elemp->PL * M->Pirho / 2e0); |
---|
| 4723 | M->Pc1 = cos(dtor(M->PdTpar)) * M->Pc0; |
---|
| 4724 | M->Ps1 = sin(dtor(M->PdTpar)) * M->Pc0; |
---|
| 4725 | cellconcat = false; |
---|
| 4726 | } |
---|
| 4727 | |
---|
| 4728 | /****************************************************************************/ |
---|
| 4729 | /* double Mpole_GetdT(long Fnum1, long Knum1) |
---|
| 4730 | |
---|
| 4731 | Purpose: |
---|
| 4732 | Return total roll angle of the element |
---|
| 4733 | Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem, which is |
---|
| 4734 | a sum of a design ,a systematic error and a random error part. |
---|
| 4735 | |
---|
| 4736 | |
---|
| 4737 | Input: |
---|
| 4738 | none |
---|
| 4739 | |
---|
| 4740 | Output: |
---|
| 4741 | none |
---|
| 4742 | |
---|
| 4743 | Return: |
---|
| 4744 | none |
---|
| 4745 | |
---|
| 4746 | Global variables: |
---|
| 4747 | none |
---|
| 4748 | |
---|
| 4749 | Specific functions: |
---|
| 4750 | none |
---|
| 4751 | |
---|
| 4752 | Comments: |
---|
| 4753 | none |
---|
| 4754 | |
---|
| 4755 | ****************************************************************************/ |
---|
| 4756 | double Mpole_GetdT(int Fnum1, int Knum1) { |
---|
| 4757 | elemtype *elemp; |
---|
| 4758 | MpoleType *M; |
---|
| 4759 | |
---|
| 4760 | elemp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem; |
---|
| 4761 | M = elemp->M; |
---|
| 4762 | |
---|
| 4763 | return (M->PdTpar + M->PdTsys + M->PdTrms * M->PdTrnd); |
---|
| 4764 | } |
---|
| 4765 | |
---|
| 4766 | /**************************************************************************** |
---|
| 4767 | * void Mpole_DefdTpar(long Fnum1, long Knum1, double PdTpar) |
---|
| 4768 | |
---|
| 4769 | Purpose: |
---|
| 4770 | Set design roll angle to {\ttfamily PdTpar} degrees. |
---|
| 4771 | |
---|
| 4772 | Input: |
---|
| 4773 | none |
---|
| 4774 | |
---|
| 4775 | Output: |
---|
| 4776 | none |
---|
| 4777 | |
---|
| 4778 | Return: |
---|
| 4779 | none |
---|
| 4780 | |
---|
| 4781 | Global variables: |
---|
| 4782 | none |
---|
| 4783 | |
---|
| 4784 | Specific functions: |
---|
| 4785 | none |
---|
| 4786 | |
---|
| 4787 | Comments: |
---|
| 4788 | none |
---|
| 4789 | |
---|
| 4790 | ****************************************************************************/ |
---|
| 4791 | void Mpole_DefdTpar(int Fnum1, int Knum1, double PdTpar) { |
---|
| 4792 | elemtype *elemp; |
---|
| 4793 | MpoleType *M; |
---|
| 4794 | |
---|
| 4795 | elemp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem; |
---|
| 4796 | M = elemp->M; |
---|
| 4797 | |
---|
| 4798 | M->PdTpar = PdTpar; |
---|
| 4799 | } |
---|
| 4800 | |
---|
| 4801 | /**************************************************************************** |
---|
| 4802 | * void Mpole_DefdTsys(long Fnum1, long Knum1, double PdTsys) |
---|
| 4803 | |
---|
| 4804 | Purpose: |
---|
| 4805 | Set systematic roll angle error to {\ttfamily PdTsys} degrees. |
---|
| 4806 | |
---|
| 4807 | Input: |
---|
| 4808 | none |
---|
| 4809 | |
---|
| 4810 | Output: |
---|
| 4811 | none |
---|
| 4812 | |
---|
| 4813 | Return: |
---|
| 4814 | none |
---|
| 4815 | |
---|
| 4816 | Global variables: |
---|
| 4817 | none |
---|
| 4818 | |
---|
| 4819 | Specific functions: |
---|
| 4820 | none |
---|
| 4821 | |
---|
| 4822 | Comments: |
---|
| 4823 | none |
---|
| 4824 | |
---|
| 4825 | ****************************************************************************/ |
---|
| 4826 | void Mpole_DefdTsys(int Fnum1, int Knum1, double PdTsys) { |
---|
| 4827 | elemtype *elemp; |
---|
| 4828 | MpoleType *M; |
---|
| 4829 | |
---|
| 4830 | elemp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]].Elem; |
---|
| 4831 | M = elemp->M; |
---|
| 4832 | |
---|
| 4833 | M->PdTsys = PdTsys; |
---|
| 4834 | } |
---|
| 4835 | |
---|
| 4836 | void Wiggler_SetPB(int Fnum1, int Knum1, int Order) { |
---|
| 4837 | CellType *cellp; |
---|
| 4838 | elemtype *elemp; |
---|
| 4839 | WigglerType *W; |
---|
| 4840 | |
---|
| 4841 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 4842 | elemp = &cellp->Elem; |
---|
| 4843 | W = elemp->W; |
---|
| 4844 | if (abs(Order) > W->Porder) |
---|
| 4845 | W->Porder = abs(Order); |
---|
| 4846 | if (W->Pmethod == Meth_Linear && Order == 2) |
---|
| 4847 | Wiggler_Setmatrix(Fnum1, Knum1, elemp->PL, W->kxV[0], 2.0 * M_PI |
---|
| 4848 | / cellp->Elem.W->lambda, W->PBW[Order + HOMmax]); |
---|
| 4849 | cellconcat = false; |
---|
| 4850 | } |
---|
| 4851 | |
---|
| 4852 | void Wiggler_SetdS(int Fnum1, int Knum1) { |
---|
| 4853 | int j; |
---|
| 4854 | CellType *cellp; |
---|
| 4855 | elemtype *elemp; |
---|
| 4856 | WigglerType *W; |
---|
| 4857 | |
---|
| 4858 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 4859 | elemp = &cellp->Elem; |
---|
| 4860 | W = elemp->W; |
---|
| 4861 | for (j = 0; j <= 1; j++) |
---|
| 4862 | cellp->dS[j] = W->PdSsys[j] + W->PdSrms[j] * W->PdSrnd[j]; |
---|
| 4863 | cellconcat = false; |
---|
| 4864 | if (W->Pmethod == Meth_Linear) |
---|
| 4865 | Wiggler_Setmatrix(Fnum1, Knum1, elemp->PL, W->kxV[0], 2.0 * M_PI |
---|
| 4866 | / cellp->Elem.W->lambda, W->PBW[Quad + HOMmax]); |
---|
| 4867 | cellconcat = false; |
---|
| 4868 | } |
---|
| 4869 | |
---|
| 4870 | void Wiggler_SetdT(int Fnum1, int Knum1) { |
---|
| 4871 | CellType *cellp; |
---|
| 4872 | elemtype *elemp; |
---|
| 4873 | WigglerType *W; |
---|
| 4874 | |
---|
| 4875 | cellp = &Cell[ElemFam[Fnum1 - 1].KidList[Knum1 - 1]]; |
---|
| 4876 | elemp = &cellp->Elem; |
---|
| 4877 | W = elemp->W; |
---|
| 4878 | cellp->dT[0] = cos(dtor(W->PdTpar + W->PdTsys + W->PdTrms * W->PdTrnd)); |
---|
| 4879 | cellp->dT[1] = sin(dtor(W->PdTpar + W->PdTsys + W->PdTrms * W->PdTrnd)); |
---|
| 4880 | if (W->Pmethod == Meth_Linear) |
---|
| 4881 | Wiggler_Setmatrix(Fnum1, Knum1, elemp->PL, W->kxV[0], 2.0 * M_PI |
---|
| 4882 | / cellp->Elem.W->lambda, W->PBW[Quad + HOMmax]); |
---|
| 4883 | cellconcat = false; |
---|
| 4884 | } |
---|