Context Navigation

t2elem.cc @ 23

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

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

Context Navigation

source: TRACY3/branches/tracy3-3.10.1b/tracy/tracy/src/t2elem.cc @ 23

Download in other formats: