Context Navigation

t2elem.cc @ 32

Last change on this file since 32 was 32, checked in by zhangj, 10 years ago
active the transport of the twiss functions and orbits of the transfer line.
Property svn:executable set to ``*
File size: 133.0 KB

Line
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
91
92	For exact hamitonian:
93
94	(1) ps = sqrt((1+delta)^2 - px_^2 - py_^2)
95
96	(2) using the check of TEAPOT to check ps < 0
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");
115	// exit_(1);
116
117	}
118	}
119	return(p_s);
120	}
121
122	/****************************************************************************/
123	/* Drift(double L, ss_vect<T> &x)
124
125	Purpose:
126	Exact Drift pass in the curvilinear coordinate for the dipole
127
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_refx) 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(Lh_bend)/(1-x0[px_]/get_p_s(x)tan(L*h_bend));
168	x[px_]= x0[px_]cos(Lh_bend) + get_p_s(x)sin(Lh_bend);
169	x[y_] = x0[y_]+x0[px_]x0[x_]tan(Lh_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(Lh_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
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_refx) 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
673	Input:
674	irho: h = 1/rho [1/m]
675	phi: dipole edge (entrance or exit) angle
676	gap: full gap between poles
677
678	2
679	K1gaph*(1 + sin phi)
680	psi = ----------------------- * (1 - K2hgap*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
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(Lh_ref) + (ps -h_bend(1/h_ref+x0[x_]))sin(L*h_ref);
735
736	dpxf = h_refsqrt((1+x[delta_])(1+x[delta_]) - x[px_]x[px_] - x[py_]x[py_])-h_bend(1+h_refx[x_]);
737
738
739	x[y_] += h_ref/h_bendx0[py_]L + x0[py_]/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
740
741	x[x_] = 1/(h_bendh_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_bendL + (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
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.
790	,
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)
796	A(x,y,-l,px,py,dP) = -sqrt( (1 + delta)^2 - px^2 - py^2 ) + delta
797
798	(exact form of the H of the kick from the magnet)
799	B(x,y,-l,px,py,dP) = -hxsqrt( (1 + delta)^2 - px^2 - py^2 ) + h x + h^2*x^2/2 + int(Re(By+iBx)/Brho)
800
801	so this is the exact Hamitonian for small ring.
802
803
804	The kick is given by
805	e
806	Dp_x = L* (h_refsqrt( (1 + delta)^2 - (px^2 + py^2) ) - h_bend - xh_bend*h_ref - ---- B_y )
807	p_0
808	Dp_y = L* e/p_0 * B_x
809
810	x = L* hx_ref / sqrt( (1 + delta)^2 - px^2 - py^2) p_x
811
812	y = L* hx_ref / sqrt( (1 + delta)^2 - px^2 - py^2) p_y
813
814	CT = L* hx_ref / sqrt( (1 + delta)^2 - px^2 - py^2) (1+delta)
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) = -hxdelta + 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_nx - a_ny) + i(bny + anx) + b_(n-1) + ia_(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
870	h_bend B_dipole/Brho is the dipole field normalized by the reference momentum; one of the example is the combined dipole
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 Zhang @ 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=0;
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
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_refh_bend*x0[x_]+ByoBrho);
940
941	// sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
942	// //
943	// x[px_] = x0[px_]cos(Lh_ref) + (get_p_s(x0)-h_bend(h_ref+x0[x_]))sin(Lh_ref)-LByoBrho;
944	//
945	// dpx = -x0[px_]h_refsin(Lh_ref)+(get_p_s(x0)-h_bend(h_ref+x0[x_]))h_refcos(L*h_ref);
946	//
947	// x[x_] = 1/(h_bendh_ref)(h_ref*sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[py_]))-dpx - h_bend);
948	//
949	// x[y_] += h_ref/h_bendx0[py_]L + x0[py_]/h_bend(asin(sqrtpxx0[px_])-asin(sqrtpx*x[px_]));
950	//
951	// x[ct_]+= (1+x0[delta_])h_ref/h_bendL + (1+x0[delta_])/h_bend(asin(sqrtpxx0[px_])-asin(sqrtpx*x[px_]));
952
953	//map of a wedge
954	// x[px_] = x0[px_]cos(Lh_bend) + (get_p_s(x0) - h_refx0[x_])sin(Lh_bend) -LByoBrho;
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(Lh_bend) + (x0[x_]x0[px_]sin(2Lh_bend)+sqr(sin(Lh_bend))(2x0[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(Lh_bend)-x0[px_]sin(Lh_bend));
962	//
963	// x[y_] += x0[py_]Lh_bendh_ref + x0[py_]h_ref*(asin(n1) - asin(n2));
964	// x[ct_]+= (1+x0[delta_])Lh_bendh_ref + (1+x0[delta_])h_ref*(asin(n1) - asin(n2));
965	//
966	}//approximate Hamiltonian
967	else {//first order map
968	x[px_] -= L * (ByoBrho + (h_bend - h_ref) / 2.0 + h_ref * h_bend
969	* x0[x_] - h_ref * x0[delta_]);
970	x[ct_] += L * h_ref * x0[x_];
971
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_]);
976	x[x_] += u * x0[px_];
977	x[y_] += u * x0[py_];
978	// x[px_]+= 0.5u/x0[x_](x0[px_]x0[px_]+x0[py_]x0[py_]);
979	x[ct_] += u(sqr(x0[px_])+sqr(x0[py_]))/(2.0(1.0+x0[delta_]));
980
981	// cout << "test......2..... u = " << u<< endl;
982	}
983	}
984
985	}//from other straight magnets
986	else
987	x[px_] -= L * (ByoBrho + h_bend*1); //h_bend should be trigger on for combined dipoles
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
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)
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)
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.
1058
1059	18/06/2012 Jianfeng Zhang @ LAL
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.
1063
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
1069	****************************************************************************/
1070	template<typename T>
1071	static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x, bool Entrance) {
1072
1073	ss_vect<T> x0;
1074	x0 = x;
1075
1076	if (true) {
1077
1078	// cout << "Dipole edge effect WITH Alex's correction: " << endl;
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!!!!
1091	// The original model of tracy is correct for tracking!!!!!!!!!!
1092	/*
1093	double psi;
1094	psi=gap0.5irho(1+sin(dtor(phi))sin(dtor(phi)))/cos(dtor(phi));
1095
1096	if(Entrance){
1097	x[x_] += irho/2x[y_]x[y_];
1098	x[px_] += 0.5x[x_]x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5x[y_]x[y_];
1099	x[py_] -= -irhotan(-dtor(psi))x[y_]+1x[x_]x[y_]+irhox[px_]x[y_]+ (irhopsi/cos(-dtor(psi))/cos(-dtor(psi))) x[y_]
1100	/ (1.0 + x[delta_]);
1101	}
1102	else
1103	{
1104	x[x_] -= irho/2x[y_]x[y_];
1105	x[px_] += 0.5x[x_]x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5x[y_]x[y_];
1106	x[py_] -= -irhotan(-dtor(psi))x[y_]+1x[x_]x[y_]-irhox[px_]x[y_]+ (irhopsi/cos(-dtor(psi))/cos(-dtor(psi))) x[y_]
1107	/ (1.0 + x[delta_]);
1108	}
1109	*/
1110
1111	//add the contribution of the entrance momentum of the particle, from Alex Loulergue.
1112	if(Entrance){
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);
1116	}else{
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);
1120	}
1121
1122	} else {//original one
1123
1124	//cout << "Dipole edge effect WITHOUT Alex's correction: " << endl;
1125	x[px_] += irho * tan(dtor(phi)) * x[x_];
1126	x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]/ (1.0 + x[delta_]*1);
1127	}
1128	}
1129
1130
1131	/******************************************************************
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
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
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.
1149	******************************************************************/
1150	template<typename T>
1151	static void BendCurvature(double irho, double H, ss_vect<T> &x) {
1152
1153	if(trace)
1154	cout << "Forest's bend curvature...."<<endl;
1155
1156	T pm2,u, coeff1, coeff2,coeff3;
1157
1158	ss_vect<T> x0;
1159
1160	x0 = x;
1161
1162	pm2 = sqr(1+x0[delta_]) - sqr(x0[px_]);
1163	u = irho*H/2/pm2;
1164
1165	coeff1 = u * 2x0[px_](1+x0[delta_])/pm2; //
1166	coeff2 = u * (1+x0[delta_]); //
1167	coeff3 = u * (pm2 - 2*(1+x0[delta_]))/pm2;
1168
1169
1170	x[x_] += x0[x_]/(1-coeff1x0[y_]x0[y_]);
1171	x[px_] -= coeff2x0[y_]x0[y_] - irhoH/2sqr(x0[x_]);
1172	x[py_] -= 2coeff2x[x_]*x0[y_];
1173	x[ct_] -= coeff3x[x_]x0[y_]*x0[y_];
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) {
1286
1287	bool track=false;
1288
1289	if(track) cout << "bend_fringe(): Forest's model"<<endl;
1290
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.0sqr(hb)gapx1[py_]x1[y_]/ps3;
1309	coeff2 = sqr(hb)gap2x1[px_](2*sqr(x1[px_])-sqr(1.0+x1[delta_]))/ps5 +
1310	hb(ps/(sqr(1+x1[delta_])-sqr(x1[px_])) + 2sqr(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)gap4sqr(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.5coeff2sqr(x[y_]);
1329	x[py_]-= coeff3*x[y_];
1330	x[ct_] -= 0.5coeff4sqr(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 + 3z0^2x0)
1347	12(1 + dP)
1348
1349	b2
1350	px = px0 +/- ---------- (2x0y0py0 - x0^2px0 - y0^2*py0)
1351	4(1 + dP)
1352
1353	b2
1354	y = y0 -/+ ---------- (y0^3 + 3x0^2y0)
1355	12(1 + dP)
1356
1357	b2
1358	py = py0 -/+ ---------- (2x0y0px0 - y0^2py0 - x0^2*py0)
1359	4(1 + dP)
1360
1361	dP = dP0;
1362
1363
1364	b2
1365	tau = tau0 -/+ ----------- (y0^3px - x0^3py + 3x0^2ypy - 3y0^2x0px)
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_]; //+2xy*py
1414
1415	if (globval.Cavity_on)
1416	x[ct_] -= ps * sqr(x[x_]) * x[y_] * x[py_]; // -3x^2y*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^2px-y^2px
1420	x[py_] -= 2.0 * u * x[y_] * x[x_] * x[px_]; //-2xy*px
1421
1422	if (globval.Cavity_on)
1423	x[ct_] += ps * sqr(x[y_]) * x[x_] * x[px_]; // +3y^2x*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){
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;
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:
1531
1532	// cout << "Mpole_Pass: Meth_Second not supported" << endl;
1533	// exit_(0);
1534	// break;
1535
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
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
1565	Drift(dL1, x);
1566	thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
1567	Drift(dL2, x);
1568	thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
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	}
1575
1576
1577	Drift(dL2, x);
1578	thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
1579	Drift(dL1, x);
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){
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);
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->scl0.0; AyoBrho[1] = FM->scl0.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.0dy) + (ax2+ax1-2.0ax0)/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.0dy) + (ay2+ay1-2.0ay0)/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.0dz); 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.0dz); 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: 6*1 orbit
2451	x: 6*6 (optics) transfer matrix?
2452
2453	Output:
2454	xref: 6*1 orbit
2455	x: 6*6 (optics) transfer matrix
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=0.0;
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 = ( sks c 0 ) av = ( sksh 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=0.0, sk=0.0, sk0=0.0, s=0.0, c=0.0;
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	K1gaph*(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	K1gaph*(1 + sin phi)
2770	psi = -----------------------, K1 = 1/2
2771	cos phi */
2772
2773	double r=0.0, s=0.0, c=0.0, sk=0.0, p=0.0, fk=0.0, afk=0.0;
2774	Matrix edge, ah, av;
2775	double coef=0.0, scoef=0.0;
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 = -------- - hxdP + 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 = -------- - hxdP + (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	( sksqrt(1+dP) \|k+hh\|*h )
2856	( -sin(u)sqrt(1+dP)sk cos(u) hsin(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	K1gaph*(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 <- ahedge /
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 <- avedge /
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 <- edgeav /
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 <- edgeav /
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.
2938	The linear fringe field is taken into account
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<-KD55M /
3290	LinTrans(5, elemp->ID->KD55, xref);
3291	// }
3292	// else
3293	// {
3294	// MulLMat(5L, elemp->ID->D55, M); /* X<-D55X /
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);
3435	exit_(1);
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
3466
3467
3468
3469
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<-Ac /
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<-AC /
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 */
3862	M->PH1 = 0.0; /* Entrance pole face curvature*/
3863	M->PH2 = 0.0; /* Exit pole face curvature */
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=1;
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=0.0;
4174	int i=1;
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(Lirho/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	}

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: TRACY3/trunk/tracy/tracy/src/t2elem.cc @ 32

Download in other formats: