Context Navigation

source: TRACY3/trunk/tracy/tracy/src/nsls-ii_lib.cc @ 11

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

Property svn:executable set to ``*
File size: 149.2 KB

Line
1	/* NSLS-II specific library
2
3	J. Bengtsson NSLS-II, BNL 2004 -
4
5	T. Shaftan, I. Pinayev, Y. Luo, C. Montag, B. Nash
6
7	*/
8	/* Current revision $Revision: 1.18 $
9	On branch $Name: $
10	Latest change by $Author: zhang $
11	*/
12
13	// global params
14
15	bool DA_bare = false, // dynamic aperture
16	freq_map = false, // include frequency map
17	tune_shift = false; // to calculate tune shift with amplitude and energy
18	int n_aper = 15, // no of dynamical aperture points
19	n_track = 512, // no of turns for tracking
20	n_orbit = 3, //number of integration for orbit correction
21	n_scale = 1; //number to scale the random orbit error
22
23	bool bba = false; // beam based alignment
24	int n_lin = 3,
25	SQ_per_scell = 2,
26	BPM_per_scell = 12,
27	HCM_per_scell = 12,
28	VCM_per_scell = 12;
29	double kick = 0.01e-3; // 0.01 mrad kick for trims
30	int n_stat = 1; // number of statistics;when set the random error;
31
32	double VDweight = 1e3, // weight for vertical dispersion
33	HVweight = 1e0, // weight for coupling Htrim vertical BPM
34	VHweight = 1e0; // weight for coupling Vtrim horizontal BPM
35
36
37	double delta_DA_ = 5.0e-2; // Parameters for dynamic aperture
38
39	// Parameters for frequency map
40	// Note NTURN is set to 10000 (2*NTURN for diffusion)) in "naffutils.h".
41	int n_x = 50, n_y = 30, n_dp = 25, n_tr = 2064;
42	double x_max_FMA = 20e-3, y_max_FMA = 6e-3, delta_FMA = 3e-2;
43	//double x_max_FMA = 20e-3, y_max_FMA = 3e-3, delta_FMA = 3e-2;
44
45	const int max_corr = 100;
46	const int max_bpm = 200;
47
48	int h_corr[max_corr], v_corr[max_corr], bpm_loc[max_bpm]; //non-global varible
49
50	int Fnum_Cart, n_iter_Cart;
51
52	double u_Touschek; // argument for Touschek D(ksi)
53	double chi_m; // argument for IBS D(ksi)
54
55	char in_dir[max_str] = "";
56
57	// Computation result files
58	const char beam_envelope_file[] = "beam_envelope.out";
59
60	// Lattice error and correction files
61	const char CodCorLatFileName[] = "codcorlat.out";
62
63	const char SkewMatFileName[] = "skewmat.out";
64	const char eta_y_FileName[] = "eta_y.out";
65	const char deta_y_FileName[] = "deta_y.out";
66
67
68	const int max_elem = Cell_nLocMax;
69
70
71
72	char ae_file[max_str], fe_file[max_str], ap_file[max_str];
73	int N_BPM, N_HCOR, N_VCOR, N_SKEW, N_COUPLE;
74
75	/* ID correction**********/
76	// maximum number of quadrupole families used for ID correction
77	const int N_Fam_max = 15;
78	// ID correction params
79	int N_calls = 1, N_steps = 1, N_Fam, Q_Fam[N_Fam_max];
80	// Weights for ID correction
81	double scl_nux = 1e2, scl_nuy = 1e2, scl_dbetax = 1e0, scl_dbetay = 1e0;
82	double scl_dnux = 0.1, scl_dnuy = 0.1, ID_step = 0.7;
83
84	int n_sext, sexts[max_elem];
85	double beta0_[max_elem][2], nu0s[max_elem][2], nu0_[2], b2[N_Fam_max];
86	double *SkewRespMat, VertCouple, SkewStrengthCorr, eta_y;
87	double b, w, V, U;
88	double disp_wave_y;
89	Vector2 beta_ref;
90
91	char lat_FileName[max_str];
92
93	// ID_corr global variables
94
95	const int n_b2_max = 500; // max no of quad corrector families
96	const int n_b3_max = 1000; // max no of sextupoles
97	const int max_ID_Fams = 25; // max no of ID families
98
99	int Nsext, Nquad, Nconstr, NconstrO, quad_prms[n_b2_max], id_loc;
100	int n_ID_Fams, ID_Fams[max_ID_Fams];
101	double Ss[n_b3_max], Sq[n_b2_max], sb[2][n_b3_max], sNu[2][n_b3_max];
102	double qb[2][n_b2_max], qb0[2][n_b2_max], sNu0[2][n_b3_max];
103	double qNu0[2][n_b2_max], qNu[2][n_b2_max], IDb[2], IDNu[2];
104	double Nu_X, Nu_Y, Nu_X0, Nu_Y0;
105	double *A1, Xsext, Xsext0, b2Ls_, w1, U1, *V1;
106	double Xoct, b4s, **Aoct;
107	Vector2 dnu0, nu_0;
108
109	ss_vect<tps> map;
110	MNF_struct MNF;
111
112
113	// conversion
114
115	void lwr_case(char str[])
116	{
117	int k;
118
119	for (k = 0; k < (int)strlen(str); k++)
120	str[k] = tolower(str[k]);
121	}
122
123
124	void upr_case(char str[])
125	{
126	int k;
127
128	for (k = 0; k < (int)strlen(str); k++)
129	str[k] = toupper(str[k]);
130	}
131
132
133	// only supported by Redhat
134	#if 0
135	// generate backtrace
136	void prt_trace (void)
137	{
138	const int max_entries = 20;
139
140	void *array[max_entries];
141	size_t size;
142	char **strings;
143	size_t i;
144
145	size = backtrace(array, max_entries);
146	strings = backtrace_symbols(array, size);
147
148	printf("prt_trace: obtained %zd stack frames\n", size);
149
150	for (i = 0; i < size; i++)
151	printf ("%s\n", strings[i]);
152
153	free (strings);
154	}
155	#endif
156
157
158	// file I/O
159
160	// C++
161
162	void file_rd(ifstream &inf, const char file_name[])
163	{
164
165	inf.open(file_name, ios::in);
166	if (!inf.is_open()) {
167	printf("File not found: %s\n", file_name);
168	exit_(-1);
169	}
170	}
171
172
173	void file_wr(ofstream &outf, const char file_name[])
174	{
175
176	outf.open(file_name, ios::out);
177	if (!outf.is_open()) {
178	printf("Could not create file: %s\n", file_name);
179	exit_(-1);
180	}
181	}
182
183
184	// C
185	/***********************************************
186	FILE* file_read(const char file_name[])
187
188	Purpose:
189	Open a file, and return the file ID
190	************************************************/
191	FILE* file_read(const char file_name[])
192	{
193	FILE *fp;
194
195	fp = fopen(file_name, "r");
196	if (fp == NULL) {
197	printf("File not found: %s\n", file_name);
198	exit_(-1);
199	} else
200	return(fp);
201	// avoid compiler warning
202	return NULL;
203	}
204
205
206	FILE* file_write(const char file_name[])
207	{
208	FILE *fp;
209
210	fp = fopen(file_name, "w");
211	if (fp == NULL) {
212	printf("Could not create file: %s\n", file_name);
213	exit_(-1);
214	// avoid compiler warning
215	return NULL;
216	} else
217	return(fp);
218	}
219
220	/* Check for closed orbit if not exit the program*/
221	void chk_cod(const bool cod, const char *proc_name)
222	{
223
224	if (!cod) {
225	printf("%s: closed orbit not found\n", proc_name);
226	exit_(1);
227	}
228	}
229
230
231	void no_sxt(void)
232	{
233	int k;
234
235	cout << endl;
236	cout << "zeroing sextupoles" << endl;
237	for (k = 0; k <= globval.Cell_nLoc; k++)
238	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->Porder >= Sext))
239	SetKpar(Cell[k].Fnum, Cell[k].Knum, Sext, 0.0);
240	}
241
242
243	void get_map(void)
244	{
245	long int lastpos;
246
247
248	getcod(0.0, lastpos);
249
250	map.identity(); map += globval.CODvect;
251	Cell_Pass(0, globval.Cell_nLoc, map, lastpos);
252	map -= globval.CODvect;
253	}
254
255
256	#if ORDER > 1
257
258	tps get_h(void)
259	{
260	ss_vect<tps> map1, R;
261
262	// Parallel transport nonlinear kick to start of lattice,
263	// assumes left to right evaluation.
264
265	if (true)
266	// Dragt-Finn factorization
267	return LieFact_DF(Inv(MNF.A0MNF.A1)mapMNF.A0MNF.A1, R)*R;
268	else {
269	// single Lie exponent
270	danot_(1); map1 = map; danot_(no_tps);
271	return LieFact(Inv(MNF.A0MNF.A1)mapInv(map1)MNF.A0*MNF.A1);
272	}
273	}
274
275	#endif
276
277	void get_m2(const ss_vect<tps> &ps, tps m2[])
278	{
279	int i, j, k;
280
281	k = 0;
282	for (i = 0; i < 2*nd_tps; i++)
283	for (j = i; j < 2*nd_tps; j++) {
284	k++; m2[k-1] = ps[i]*ps[j];
285	}
286	}
287
288
289	ss_vect<tps> get_S(const int n_DOF)
290	{
291	int j;
292	ss_vect<tps> S;
293
294	S.zero();
295	for (j = 0; j < n_DOF; j++) {
296	S[2j] = tps(0.0, 2j+2); S[2j+1] = -tps(0.0, 2j+1);
297	}
298
299	return S;
300	}
301
302
303	ss_vect<tps> tp_S(const int n_DOF, const ss_vect<tps> &A)
304	{
305	int j, jj[ss_dim];
306	ss_vect<tps> S;
307
308	for (j = 1; j <= ss_dim; j++)
309	jj[j-1] = (j <= 2*n_DOF)? 1 : 0;
310
311	S = get_S(n_DOF);
312
313	return -SPInv(A, jj)S;
314	}
315
316
317	void get_dnu(const ss_vect<tps> &A, double dnu[])
318	{
319	int k;
320
321	for (k = 0; k <= 1; k++) {
322	dnu[k] = atan2(A[2k][2k+1], A[2k][2k])/(2.0*M_PI);
323	if (dnu[k] < 0.0) dnu[k] += 1.0;
324	}
325	}
326
327
328	ss_vect<tps> get_A_CS(const ss_vect<tps> &A, double dnu[])
329	{
330	int k;
331	double c, s;
332	ss_vect<tps> Id, R;
333
334	Id.identity(); R.identity(); get_dnu(A, dnu);
335	for (k = 0; k <= 1; k++) {
336
337	c = cos(2.0M_PIdnu[k]); s = sin(2.0M_PIdnu[k]);
338	R[2k] = cId[2k] - sId[2k+1]; R[2k+1] = sId[2k] + cId[2k+1];
339	}
340
341	return A*R;
342	}
343
344
345	void get_twoJ(const int n_DOF, const ss_vect<double> &ps,
346	const ss_vect<tps> &A, double twoJ[])
347	{
348	int j;
349	iVector jj;
350	ss_vect<double> z;
351
352	for (j = 0; j < nv_tps; j++)
353	jj[j] = (j < 2*n_DOF)? 1 : 0;
354
355	z = (PInv(A, jj)*ps).cst();
356
357	for (j = 0; j < n_DOF; j++)
358	twoJ[j] = sqr(z[2j]) + sqr(z[2j+1]);
359	}
360
361
362	double get_curly_H(const double alpha_x, const double beta_x,
363	const double eta_x, const double etap_x)
364	{
365	double curly_H, gamma_x;
366
367	gamma_x = (1.0+pow(alpha_x, 2))/beta_x;
368
369	curly_H = gamma_xsqr(eta_x) + 2.0alpha_xeta_xetap_x + beta_x*sqr(etap_x);
370
371	return curly_H;
372	}
373
374
375	double get_eps_x(void)
376	{
377	bool cav, emit;
378	long int lastpos;
379	double eps_x;
380	ss_vect<tps> A;
381
382	/* Note:
383
384	T
385	M J M = J,
386
387	-1 T \| 0 I \| T \| beta -alpha \|
388	A = -J A J, J = \| \|, A A = \| \|
389	\| -I 0 \| \| -alpha gamma \|
390
391	Transform to Floquet Space:
392
393	-1 T
394	A eta = -J A J eta,
395
396	-1 T -1 T T
397	H~ = ( A eta ) A eta = ( J eta ) A A ( J eta )
398
399	*/
400
401	cav = globval.Cavity_on; emit = globval.emittance;
402
403	globval.Cavity_on = false; globval.emittance = false;
404
405	Ring_GetTwiss(false, 0.0);
406
407	putlinmat(6, globval.Ascr, A); A += globval.CODvect;
408
409	globval.emittance = true;
410
411	Cell_Pass(0, globval.Cell_nLoc, A, lastpos);
412
413	eps_x = 1470.0pow(globval.Energy, 2)I5/(I2-I4);
414
415	printf("\n");
416	printf("eps_x = %5.3f nm.rad\n", eps_x);
417	printf("J_x = %5.3f, J_z = %5.3f\n", 1.0-I4/I2, 2.0+I4/I2);
418
419	globval.Cavity_on = cav; globval.emittance = emit;
420
421	return eps_x;
422	}
423
424	/****************************************************************************/
425	/* void GetEmittance(const int Fnum, const bool prt)
426
427	Purpose:
428	get the emittance
429
430	Input:
431	Fnum family index
432	prt bool flag, true = print the calcuate information
433	Output:
434	none
435
436	Return:
437	emittance
438
439	Global variables:
440
441
442	specific functions:
443
444
445	Comments:
446
447
448	****************************************************************************/
449	void GetEmittance(const int Fnum, const bool prt)
450	{
451	bool emit, rad, cav, path;
452	int i, j, h_RF;
453	long int lastpos, loc;
454	double C, theta, V_RF, phi0, alpha_z, beta_z, gamma_z;
455	double sigma_s, sigma_delta;
456	Vector3 nu;
457	Matrix Ascr;
458	ss_vect<tps> Ascr_map;
459
460	// save state
461	rad = globval.radiation; emit = globval.emittance;
462	cav = globval.Cavity_on; path = globval.pathlength;
463
464	C = Cell[globval.Cell_nLoc].S; /* Circumference of the ring*/
465
466	// damped system
467	globval.radiation = true; globval.emittance = true;
468	globval.Cavity_on = true; globval.pathlength = false;
469
470	Ring_GetTwiss(false, 0.0);
471
472	// radiation loss is computed in Cav_Pass
473
474	// Sum_k( alpha_k) = 2 * U_0 / E
475	// printf("\n");
476	// printf("%24.16e %24.16e\n",
477	// globval.dE,
478	// (globval.alpha_rad[X_]+globval.alpha_rad[Y_]
479	// +globval.alpha_rad[Z_])/2.0);
480
481	globval.U0 = 1e9globval.dEglobval.Energy;
482	V_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Pvolt; //RF voltage
483	h_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Ph; // RF cavity harmonic number
484	phi0 = fabs(asin(globval.U0/V_RF));
485	globval.delta_RF =
486	sqrt(-V_RFcos(M_PI-phi0)(2.0-(M_PI-2.0*(M_PI-phi0))
487	tan(M_PI-phi0))/(abs(globval.Alphac)M_PIh_RF1e9*globval.Energy));
488
489	// compute diffusion coeffs. for eigenvectors [sigma_xx, sigma_yy, sigma_zz]
490	putlinmat(6, globval.Ascr, Ascr_map); Ascr_map += globval.CODvect;
491
492	Cell_Pass(0, globval.Cell_nLoc, Ascr_map, lastpos);
493
494	for (i = 0; i < DOF; i++) {
495	// partition numbers
496	globval.J[i] = 2.0(1.0+globval.CODvect[delta_])globval.alpha_rad[i]
497	/globval.dE;
498	// damping times
499	globval.tau[i] = -C/(c0*globval.alpha_rad[i]);
500	// diffusion coeff. and emittance
501	globval.eps[i] = -globval.D_rad[i]/(2.0*globval.alpha_rad[i]);
502	// fractional tunes
503	nu[i] = atan2(globval.wi[i2], globval.wr[i2])/(2.0*M_PI);
504	if (nu[i] < 0.0) nu[i] = 1.0 + nu[i];
505	}
506
507	// Note, J_x + J_y + J_z not exactly 4 (1st order perturbations)
508	// printf("\n");
509	// printf("%24.16e\n", globval.J[X_]+globval.J[Y_]+globval.J[Z_]);
510
511	// undamped system
512	globval.radiation = false; globval.emittance = false;
513
514	Ring_GetTwiss(false, 0.0);
515
516	/* compute the sigmas arround the lattice:
517
518	sigma = A diag[J_1, J_1, J_2, J_2, J_3, J_3] A^T
519
520	*/
521	for (i = 0; i < 6; i++) {
522	Ascr_map[i] = tps(globval.CODvect[i]);
523	for (j = 0; j < 6; j++)
524	Ascr_map[i] += globval.Ascr[i][j]sqrt(globval.eps[j/2])tps(0.0, j+1);
525	}
526	for (loc = 0; loc <= globval.Cell_nLoc; loc++) {
527	Elem_Pass(loc, Ascr_map);
528	// sigma = A x A^tp
529	getlinmat(6, Ascr_map, Cell[loc].sigma); TpMat(6, Cell[loc].sigma);
530	getlinmat(6, Ascr_map, Ascr); MulLMat(6, Ascr, Cell[loc].sigma);
531	}
532
533	theta = atan2(2e0*Cell[0].sigma[x_][y_],
534	(Cell[0].sigma[x_][x_]-Cell[0].sigma[y_][y_]))/2e0;
535
536	// longitudinal alpha and beta
537	alpha_z =
538	-globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
539	- globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
540	beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
541	gamma_z = (1.0+sqr(alpha_z))/beta_z;
542
543	// bunch size
544	sigma_s = sqrt(beta_z*globval.eps[Z_]);
545	sigma_delta = sqrt(gamma_z*globval.eps[Z_]);
546
547	if (prt) {
548	printf("\n");
549	printf("Emittance:\n");
550	printf("\n");
551	printf("Energy loss per turn [keV]: "
552	"U0 = %3.1f\n",
553	1e-3*globval.U0);
554	printf("Synchronous phase [deg]: "
555	"phi0 = 180 - %4.2f\n",
556	phi0*180.0/M_PI);
557	printf("RF bucket height [%%]: "
558	"delta_RF = %4.2f\n", 1e2*globval.delta_RF);
559	printf("\n");
560	printf("Equilibrium emittance [m.rad]: "
561	"eps_x = %9.3e, eps_y = %9.3e, eps_z = %9.3e\n",
562	globval.eps[X_], globval.eps[Y_], globval.eps[Z_]);
563	printf("Bunch length [mm]: "
564	"sigma_s = %5.3f\n", 1e3*sigma_s);
565	printf("Momentum spread: "
566	"sigma_delta = %9.3e\n", sigma_delta);
567	printf("Partition numbers: "
568	"J_x = %5.3f, J_y = %5.3f, J_z = %5.3f\n",
569	globval.J[X_], globval.J[Y_], globval.J[Z_]);
570	printf("Damping times [msec]: "
571	"tau_x = %3.1f, tau_y = %3.1f, tau_z = %3.1f\n",
572	1e3globval.tau[X_], 1e3globval.tau[Y_], 1e3*globval.tau[Z_]);
573	printf("\n");
574	printf("alphac: "
575	"alphac = %8.4e\n", globval.Alphac);
576	printf("\n");
577	printf("Fractional tunes: "
578	"nu_x = %7.5f, nu_y = %7.5f, nu_z = %7.5f\n",
579	nu[X_], nu[Y_], nu[Z_]);
580	printf(" "
581	"1-nu_x = %7.5f, 1-nu_y = %7.5f, 1-nu_z = %7.5f\n",
582	1e0-nu[X_], 1e0-nu[Y_], 1e0-nu[Z_]);
583	printf("\n");
584	printf("sigmas: "
585	"sigma_x = %5.1f microns, sigma_px = %5.1f urad\n",
586	1e6sqrt(Cell[0].sigma[x_][x_]), 1e6sqrt(Cell[0].sigma[px_][px_]));
587	printf(" "
588	"sigma_y = %5.1f microns, sigma_py = %5.1f urad\n",
589	1e6sqrt(Cell[0].sigma[y_][y_]), 1e6sqrt(Cell[0].sigma[py_][py_]));
590	printf(" "
591	"sigma_s = %5.2f mm, sigma_delta = %8.2e\n",
592	1e3*sqrt(Cell[0].sigma[ct_][ct_]),
593	sqrt(Cell[0].sigma[delta_][delta_]));
594
595	printf("\n");
596	printf("Beam ellipse twist [rad]: tw = %5.3f\n", theta);
597	printf(" [deg]: tw = %5.3f\n", theta*180.0/M_PI);
598	}
599
600	// restore state
601	globval.radiation = rad; globval.emittance = emit;
602	globval.Cavity_on = cav; globval.pathlength = path;
603	}
604
605
606	// output
607
608	void prt_lat(const char *fname, const int Fnum, const bool all)
609	{
610	long int i = 0;
611	double I2, I5, code = 0.0;
612	FILE *outf;
613
614	outf = file_write(fname);
615	fprintf(outf, "# name s code"
616	" alphax betax nux etax etapx");
617	fprintf(outf, " alphay betay nuy etay etapy I5\n");
618	fprintf(outf, "# [m]"
619	" [m] [m]");
620	fprintf(outf, " [m] [m]\n");
621	fprintf(outf, "#\n");
622
623	I2 = 0.0; I5 = 0.0;
624	for (i = 0; i <= globval.Cell_nLoc; i++) {
625	if (all \|\| (Cell[i].Fnum == Fnum)) {
626	switch (Cell[i].Elem.Pkind) {
627	case drift:
628	code = 0.0;
629	break;
630	case Mpole:
631	if (Cell[i].Elem.M->Pirho != 0.0)
632	code = 0.5;
633	else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0)
634	code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]);
635	else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0)
636	code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]);
637	else if (Cell[i].Fnum == globval.bpm)
638	code = 2.0;
639	else
640	code = 0.0;
641	break;
642	default:
643	code = 0.0;
644	break;
645	}
646	fprintf(outf, "%4ld %15s %6.2f %4.1f"
647	" %7.3f %6.3f %6.3f %6.3f %6.3f"
648	" %7.3f %6.3f %6.3f %6.3f %6.3f %8.2e\n",
649	i, Cell[i].Elem.PName, Cell[i].S, code,
650	Cell[i].Alpha[X_], Cell[i].Beta[X_], Cell[i].Nu[X_],
651	Cell[i].Eta[X_], Cell[i].Etap[X_],
652	Cell[i].Alpha[Y_], Cell[i].Beta[Y_], Cell[i].Nu[Y_],
653	Cell[i].Eta[Y_], Cell[i].Etap[Y_], I5);
654	}
655	}
656
657	// fprintf(outf, "\n");
658	// fprintf(outf, "# emittance: %5.3f nm.rad\n", get_eps_x());
659
660	fclose(outf);
661	}
662
663
664	void prt_chrom_lat(void)
665	{
666	long int i;
667	double dbeta_ddelta[Cell_nLocMax][2], detax_ddelta[Cell_nLocMax];
668	double ksi[Cell_nLocMax][2];
669	double code = 0.0;
670	FILE *outf;
671
672	printf("\n");
673	printf("prt_chrom_lat: calling Ring_GetTwiss with delta != 0\n");
674	Ring_GetTwiss(true, globval.dPcommon);
675	for (i = 0; i <= globval.Cell_nLoc; i++) {
676	dbeta_ddelta[i][X_] = Cell[i].Beta[X_];
677	dbeta_ddelta[i][Y_] = Cell[i].Beta[Y_];
678	detax_ddelta[i] = Cell[i].Eta[X_];
679	}
680	printf("prt_chrom_lat: calling Ring_GetTwiss with delta != 0\n");
681	Ring_GetTwiss(true, -globval.dPcommon);
682	ksi[0][X_] = 0.0; ksi[0][Y_] = 0.0;
683	for (i = 0; i <= globval.Cell_nLoc; i++) {
684	dbeta_ddelta[i][X_] -= Cell[i].Beta[X_];
685	dbeta_ddelta[i][Y_] -= Cell[i].Beta[Y_];
686	detax_ddelta[i] -= Cell[i].Eta[X_];
687	dbeta_ddelta[i][X_] /= 2.0*globval.dPcommon;
688	dbeta_ddelta[i][Y_] /= 2.0*globval.dPcommon;
689	detax_ddelta[i] /= 2.0*globval.dPcommon;
690	if (i != 0) {
691	ksi[i][X_] = ksi[i-1][X_]; ksi[i][Y_] = ksi[i-1][Y_];
692	}
693	if (Cell[i].Elem.Pkind == Mpole) {
694	ksi[i][X_] -= Cell[i].Elem.M->PBpar[Quad+HOMmax]
695	Cell[i].Elem.PLCell[i].Beta[X_]/(4.0*M_PI);
696	ksi[i][Y_] += Cell[i].Elem.M->PBpar[Quad+HOMmax]
697	Cell[i].Elem.PLCell[i].Beta[Y_]/(4.0*M_PI);
698	}
699	}
700
701	outf = file_write("chromlat.out");
702	fprintf(outf, "# name s code"
703	" bxex sqrt(bxby) dbx/ddex bxdex/dd"
704	" byex dby/ddex by*dex/dd ksix ksiy"
705	" dbx/dd bx/dd dex/dd\n");
706	fprintf(outf, "# [m] [m]"
707	" [m] [m] [m]");
708	fprintf(outf, " [m] [m] [m]\n");
709	fprintf(outf, "#\n");
710	for (i = 0; i <= globval.Cell_nLoc; i++) {
711	switch (Cell[i].Elem.Pkind) {
712	case drift:
713	code = 0.0;
714	break;
715	case Mpole:
716	if (Cell[i].Elem.M->Pirho != 0)
717	code = 0.5;
718	else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0)
719	code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]);
720	else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0)
721	code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]);
722	else if (Cell[i].Fnum == globval.bpm)
723	code = 2.0;
724	else
725	code = 0.0;
726	break;
727	default:
728	code = 0.0;
729	break;
730	}
731	fprintf(outf, "%4ld %15s %6.2f %4.1f"
732	" %6.3f %8.3f %8.3f %8.3f"
733	" %6.3f %8.3f %8.3f %5.2f %5.2f"
734	" %6.3f %6.3f %6.3f\n",
735	i, Cell[i].Elem.PName, Cell[i].S, code,
736	Cell[i].Beta[X_]*Cell[i].Eta[X_],
737	sqrt(Cell[i].Beta[X_]*Cell[i].Beta[Y_]),
738	dbeta_ddelta[i][X_]*Cell[i].Eta[X_],
739	detax_ddelta[i]*Cell[i].Beta[X_],
740	Cell[i].Beta[Y_]*Cell[i].Eta[X_],
741	dbeta_ddelta[i][Y_]*Cell[i].Eta[X_],
742	detax_ddelta[i]*Cell[i].Beta[Y_],
743	ksi[i][X_], ksi[i][Y_],
744	dbeta_ddelta[i][X_], dbeta_ddelta[i][Y_], detax_ddelta[i]);
745	}
746	fclose(outf);
747	}
748
749	/****************************************************************************/
750	/* void prt_cod(const char *file_name, const int Fnum, const bool all)
751
752	Purpose:
753	print the location, twiss parameters, displacement errors and close orbit
754	at the the location of Fnum or all lattice elements
755
756	Input:
757	file_name file to save the data
758	code in the file:
759	0 drfit
760	0.5 dipole
761	-1 defocusing quadrupole
762	1 focusing quadrupole
763	1.5 sextupole
764	2.0 bpm
765	Fnum family index
766	all true, print the cod at all elements
767	false, print the cod at the family elements
768	Output:
769	none
770
771	Return:
772	none
773
774	Global variables:
775
776
777	specific functions:
778
779
780	Comments:
781
782
783	****************************************************************************/
784	void prt_cod(const char *file_name, const int Fnum, const bool all)
785	{
786	long i;
787	double code = 0.0;
788	FILE *outf;
789	long FORLIM;
790	struct tm *newtime;
791
792	outf = file_write(file_name);
793
794	/* Get time and date */
795	newtime = GetTime();
796
797	fprintf(outf,"# TRACY III -- %s -- %s \n",
798	file_name, asctime2(newtime));
799
800	fprintf(outf, "# name s code betax nux betay nuy"
801	" xcod ycod dSx dSy dipx dipy\n");
802	fprintf(outf, "# [m] [m] [m] "
803	" [mm] [mm] [mm] [mm] [mrad] [mrad]\n");
804	fprintf(outf, "#\n");
805
806	FORLIM = globval.Cell_nLoc;
807	for (i = 0L; i <= FORLIM; i++) {
808	if (all \|\| (Cell[i].Fnum == Fnum)) {
809	switch (Cell[i].Elem.Pkind) {
810	case drift:
811	code = 0.0;
812	break;
813	case Mpole:
814	if (Cell[i].Elem.M->Pirho != 0)
815	code = 0.5;
816	else if (Cell[i].Elem.M->PBpar[Quad+HOMmax] != 0)
817	code = sgn(Cell[i].Elem.M->PBpar[Quad+HOMmax]);
818	else if (Cell[i].Elem.M->PBpar[Sext+HOMmax] != 0)
819	code = 1.5*sgn(Cell[i].Elem.M->PBpar[Sext+HOMmax]);
820	else if (Cell[i].Fnum == globval.bpm)
821	code = 2.0;
822	else
823	code = 0.0;
824	break;
825	default:
826	code = 0.0;
827	break;
828	}
829	/* COD is in local coordinates */
830	fprintf(outf, "%4ld %.*s %6.2f %4.1f %6.3f %6.3f %6.3f %6.3f"
831	" %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n",
832	i, SymbolLength, Cell[i].Elem.PName, Cell[i].S, code,
833	Cell[i].Beta[X_], Cell[i].Nu[X_],
834	Cell[i].Beta[Y_], Cell[i].Nu[Y_],
835	1e3Cell[i].BeamPos[x_], 1e3Cell[i].BeamPos[y_],
836	1e3Cell[i].dS[X_], 1e3Cell[i].dS[Y_], //displacement error
837	-1e3*Elem_GetKval(Cell[i].Fnum, Cell[i].Knum, Dip), //H kick angle, for correctors
838	1e3*Elem_GetKval(Cell[i].Fnum, Cell[i].Knum, -Dip)); //kick angle, for correctors
839	}
840	}
841	fclose(outf);
842	}
843
844
845	void prt_beampos(const char *file_name)
846	{
847	int k;
848	ofstream outf;
849
850	file_wr(outf, file_name);
851
852	outf << "# k s name x px y py delta ct" << endl;
853	outf << "#" << endl;
854
855	for (k = 0; k <= globval.Cell_nLoc; k++)
856	outf << scientific << setprecision(5)
857	<< setw(5) << k << setw(11) << Cell[k].Elem.PName
858	<< setw(13) << Cell[k].BeamPos << endl;
859
860	outf.close();
861	}
862
863
864	/****************************************************************
865	void CheckAlignTol(const char *OutputFile)
866
867	Purpose:
868	check aligment errors of individual magnets on giders
869	the dT and roll angle are all printed out
870	****************************************************************/
871	// misalignments
872	void CheckAlignTol(const char *OutputFile)
873	// check aligment errors of individual magnets on giders
874	// the dT and roll angle are all printed out
875	{
876	int i,j;
877	int n_girders;
878	int gs_Fnum, ge_Fnum;
879	int gs_nKid, ge_nKid;
880	int dip_Fnum,dip_nKid;
881	int loc, loc_gs, loc_ge;
882	char * name;
883	double s;
884	double PdSsys[2], PdSrms[2], PdSrnd[2], dS[2], dT[2];
885	fstream fout;
886
887	gs_Fnum = globval.gs; gs_nKid = GetnKid(gs_Fnum);
888	ge_Fnum = globval.ge; ge_nKid = GetnKid(ge_Fnum);
889	if (gs_nKid == ge_nKid)
890	n_girders= gs_nKid;
891	else {
892	cout << " The numbers of GS and GE not same. " << endl;
893	exit (1);
894	}
895
896	fout.open(OutputFile,ios::out);
897	if(!fout) {
898	cout << "error in opening the file " << endl;
899	exit_(0);
900	}
901
902	fout << "Girders, Quads, Sexts: " << endl;
903	for (i = 1; i <= n_girders; i++){
904	fout << i << ":" << endl;
905	loc_gs = Elem_GetPos(gs_Fnum, i); loc_ge = Elem_GetPos(ge_Fnum, i);
906
907	loc = loc_gs;
908	PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
909	PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
910	PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
911	PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
912	PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
913	PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
914	dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
915	dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
916	s = Cell[loc].S; name = Cell[loc].Elem.PName;
917	fout << " " << name << " " << loc << " " << s
918	<< " " << PdSsys[X_] << " " << PdSsys[Y_]
919	<< " " << PdSrms[X_] << " " << PdSrms[Y_]
920	<< " " << PdSrnd[X_] << " " << PdSrnd[Y_]
921	<< " " << Cell[loc].Elem.M->PdTrms << " "
922	<< Cell[loc].Elem.M->PdTrnd << " " << dS[X_] << " " << dS[Y_]
923	<< " " << atan2( dT[1], dT[0] ) << endl;
924
925	for (j = loc_gs+1; j < loc_ge; j++) {
926	if ((Cell[j].Elem.Pkind == Mpole) &&
927	(Cell[j].Elem.M->n_design >= Quad \|\|
928	Cell[j].Elem.M->n_design >= Sext)) {
929	loc = j;
930	PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
931	PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
932	PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
933	PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
934	PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
935	PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
936	dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
937	dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
938	s = Cell[loc].S; name=Cell[loc].Elem.PName;
939	fout << " " << name << " " << loc << " " << s
940	<< " " << PdSsys[X_] << " " << PdSsys[Y_]
941	<< " " << PdSrms[X_] << " " << PdSrms[Y_]
942	<< " " << PdSrnd[X_] << " " << PdSrnd[Y_]
943	<< " " << Cell[loc].Elem.M->PdTrms << " "
944	<< Cell[loc].Elem.M->PdTrnd
945	<< " " << dS[X_] << " " << dS[Y_]
946	<< " " << atan2( dT[1], dT[0] ) << endl;
947	}
948	}
949
950	loc = loc_ge;
951	PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
952	PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
953	PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
954	PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
955	PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
956	PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
957	dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
958	dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
959	s=Cell[loc].S; name=Cell[loc].Elem.PName;
960	fout << " " << name << " " << loc << " " << s
961	<< " " << PdSsys[X_] << " " << PdSsys[Y_]
962	<< " " << PdSrms[X_] << " " << PdSrms[Y_]
963	<< " " << PdSrnd[X_] << " " << PdSrnd[Y_]
964	<< " " << Cell[loc].Elem.M->PdTrms
965	<< " " << Cell[loc].Elem.M->PdTrnd
966	<< " " << dS[X_] << " " << dS[Y_]
967	<< " " << atan2( dT[1], dT[0] ) << endl;
968
969	}
970
971	fout << " " << endl;
972	fout << "Dipoles: " << endl;
973	dip_Fnum = ElemIndex("B1"); dip_nKid = GetnKid(dip_Fnum);
974	for (i = 1; i <= dip_nKid; i++){
975	loc = Elem_GetPos(dip_Fnum, i);
976	PdSsys[X_] = Cell[loc].Elem.M->PdSsys[X_];
977	PdSsys[Y_] = Cell[loc].Elem.M->PdSsys[Y_];
978	PdSrms[X_] = Cell[loc].Elem.M->PdSrms[X_];
979	PdSrms[Y_] = Cell[loc].Elem.M->PdSrms[Y_];
980	PdSrnd[X_] = Cell[loc].Elem.M->PdSrnd[X_];
981	PdSrnd[Y_] = Cell[loc].Elem.M->PdSrnd[Y_];
982	dS[X_] = Cell[loc].dS[X_]; dS[Y_] = Cell[loc].dS[Y_];
983	dT[0] = Cell[loc].dT[0]; dT[1] = Cell[loc].dT[1];
984	s = Cell[loc].S; name = Cell[loc].Elem.PName;
985	fout << " " << name << " " << loc << " " << s
986	<< " " << PdSsys[X_] << " " << PdSsys[Y_]
987	<< " " << PdSrms[X_] << " " << PdSrms[Y_]
988	<< " " << PdSrnd[X_] << " " << PdSrnd[Y_]
989	<< " " << Cell[loc].Elem.M->PdTrms
990	<< " " << Cell[loc].Elem.M->PdTrnd
991	<< " " << dS[X_] << " " << dS[Y_]
992	<< " " << atan2( dT[1], dT[0] ) << endl;
993	}
994
995	fout.close();
996	}
997
998	/*********************************************************************
999	void misalign_rms_elem(const int Fnum, const int Knum,
1000	const double dx_rms, const double dy_rms,
1001	const double dr_rms, const bool new_rnd)
1002
1003	Purpose:
1004	Set random misalignment to the element with Fnum and Knum
1005
1006	Input:
1007	Fnum family number
1008	Knum kid number
1009	dx_rms rms value of the error in x
1010	dy_rms rms value of the error in y
1011	dr_rms rms value of the error in rotation around s
1012	new_rnd flag to turn on/off using the new random seed
1013
1014	**********************************************************************/
1015	void misalign_rms_elem(const int Fnum, const int Knum,
1016	const double dx_rms, const double dy_rms,
1017	const double dr_rms, const bool new_rnd)
1018	{
1019	long int loc;
1020	MpoleType *mp;
1021
1022	const bool normal = true;
1023
1024	loc = Elem_GetPos(Fnum, Knum); mp = Cell[loc].Elem.M;
1025
1026	mp->PdSrms[X_] = dx_rms; mp->PdSrms[Y_] = dy_rms; mp->PdTrms = dr_rms;
1027	if (new_rnd) {
1028	if (normal) {
1029	mp->PdSrnd[X_] = normranf(); mp->PdSrnd[Y_] = normranf();
1030	mp->PdTrnd = normranf();
1031	} else {
1032	mp->PdSrnd[X_] = ranf(); mp->PdSrnd[Y_] = ranf();
1033	mp->PdTrnd = ranf();
1034	}
1035	}
1036
1037	Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum);
1038	}
1039
1040
1041	/*********************************************************************
1042	void misalign_sys_elem(const int Fnum, const int Knum,
1043	const double dx_sys, const double dy_sys,
1044	const double dr_sys)
1045
1046	Purpose:
1047	Set systematic misalignment to the elements in "type"
1048
1049	Input:
1050	Fnum family number
1051	Knum kid number
1052	dx_sys sys value of the error in x
1053	dy_sys sys value of the error in y
1054	dr_sys sys value of the error in rotation around s
1055
1056	**********************************************************************/
1057	void misalign_sys_elem(const int Fnum, const int Knum,
1058	const double dx_sys, const double dy_sys,
1059	const double dr_sys)
1060	{
1061	long int loc;
1062	MpoleType *mp;
1063
1064	loc = Elem_GetPos(Fnum, Knum); mp = Cell[loc].Elem.M;
1065
1066	mp->PdSsys[X_] = dx_sys; mp->PdSsys[Y_] = dy_sys; mp->PdTsys = dr_sys;
1067
1068	Mpole_SetdS(Fnum, Knum); Mpole_SetdT(Fnum, Knum);
1069	}
1070
1071	/*********************************************************************
1072	void misalign_rms_fam(const int Fnum,
1073	const double dx_rms, const double dy_rms,
1074	const double dr_rms, const bool new_rnd)
1075
1076	Purpose:
1077	Set systematic misalignment to the elements on the girders
1078
1079	Input:
1080	Fnum family number
1081	dx_rms rms value of the error in x
1082	dy_rms rms value of the error in y
1083	dr_rms rms value of the error in rotation around s
1084	new_rnd flag to generate random number
1085	**********************************************************************/
1086	void misalign_rms_fam(const int Fnum,
1087	const double dx_rms, const double dy_rms,
1088	const double dr_rms, const bool new_rnd)
1089	{
1090	int i;
1091
1092	for (i = 1; i <= GetnKid(Fnum); i++)
1093	misalign_rms_elem(Fnum, i, dx_rms, dy_rms, dr_rms, new_rnd);
1094	}
1095
1096	void misalign_sys_fam(const int Fnum,
1097	const double dx_sys, const double dy_sys,
1098	const double dr_sys)
1099	{
1100	int i;
1101
1102	for (i = 1; i <= GetnKid(Fnum); i++)
1103	misalign_sys_elem(Fnum, i, dx_sys, dy_sys, dr_sys);
1104	}
1105
1106	/*********************************************************************
1107	void misalign_rms_type(const int type,
1108	const double dx_rms, const double dy_rms,
1109	const double dr_rms, const bool new_rnd)
1110
1111	Purpose:
1112	Set random misalignment to the elements in "type"
1113
1114	Input:
1115	type element type,has the following value
1116	All
1117	dipole
1118	quad
1119	sext
1120	bpm
1121	girder
1122
1123	dx_rms rms value of the error in x
1124	dy_rms rms value of the error in y
1125	dr_rms rms value of the error in rotation around s
1126	new_rnd flag to turn on/off using the new random seed
1127
1128	**********************************************************************/
1129	void misalign_rms_type(const int type,
1130	const double dx_rms, const double dy_rms,
1131	const double dr_rms, const bool new_rnd)
1132	{
1133	long int k;
1134
1135	if ((type >= All) && (type <= HOMmax)) {
1136	for (k = 1; k <= globval.Cell_nLoc; k++) {
1137	if ((Cell[k].Elem.Pkind == Mpole) &&
1138	((type == Cell[k].Elem.M->n_design) \|\|
1139	((type == All) &&
1140	((Cell[k].Fnum != globval.gs) && (Cell[k].Fnum != globval.ge))))) {
1141	// if all: skip girders
1142	misalign_rms_elem(Cell[k].Fnum, Cell[k].Knum,
1143	dx_rms, dy_rms, dr_rms, new_rnd);
1144	}
1145	}
1146	} else {
1147	printf("misalign_rms_type: incorrect type %d\n", type); exit_(1);
1148	}
1149	}
1150
1151	/*********************************************************************
1152	void misalign_sys_type(const int type,
1153	const double dx_sys, const double dy_sys,
1154	const double dr_sys)
1155
1156	Purpose:
1157	Set systematic misalignment to the elements in "type"
1158
1159	Input:
1160	type element type,has the following value
1161	All
1162	dipole
1163	quad
1164	sext
1165	bpm
1166	girder
1167
1168	dx_sys sys value of the error in x
1169	dy_sys sys value of the error in y
1170	dr_sys sys value of the error in rotation around s
1171
1172	**********************************************************************/
1173	void misalign_sys_type(const int type,
1174	const double dx_sys, const double dy_sys,
1175	const double dr_sys)
1176	{
1177	long int k;
1178
1179	if ((type >= All) && (type <= HOMmax)) {
1180	for (k = 1; k <= globval.Cell_nLoc; k++) {
1181	if ((Cell[k].Elem.Pkind == Mpole) &&
1182	((type == Cell[k].Elem.M->n_design) \|\|
1183	((type == All) &&
1184	((Cell[k].Fnum != globval.gs) && (Cell[k].Fnum != globval.ge))))) {
1185	// if all: skip girders
1186	misalign_sys_elem(Cell[k].Fnum, Cell[k].Knum,
1187	dx_sys, dy_sys, dr_sys);
1188	}
1189	}
1190	} else {
1191	printf("misalign_sys_type: incorrect type %d\n", type); exit_(1);
1192	}
1193	}
1194
1195	/*********************************************************************
1196	void misalign_rms_girders(const int gs, const int ge,
1197	const double dx_rms, const double dy_rms,
1198	const double dr_rms, const bool new_rnd)
1199
1200	Purpose:
1201	Set systematic misalignment to the elements on the girders
1202
1203	Input:
1204	gs girder start marker
1205	ge girder end marker
1206	dx_rms rms value of the error in x
1207	dy_rms rms value of the error in y
1208	dr_rms rms value of the error in rotation around s
1209
1210	**********************************************************************/
1211	void misalign_rms_girders(const int gs, const int ge,
1212	const double dx_rms, const double dy_rms,
1213	const double dr_rms, const bool new_rnd)
1214	{
1215	int i, k, n_girders, n_ge, n_gs;
1216	long int loc_gs, loc_ge, j;
1217	double s_gs, s_ge, dx_gs[2], dx_ge[2], s;
1218
1219	n_gs = GetnKid(gs); n_ge = GetnKid(ge);
1220
1221	if (n_gs == n_ge)
1222	n_girders = n_gs;
1223	else {
1224	cout << "set_girders: no of GS != no of GE" << endl;
1225	exit (1);
1226	}
1227
1228	misalign_rms_fam(gs, dx_rms, dy_rms, dr_rms, new_rnd);
1229	misalign_rms_fam(ge, dx_rms, dy_rms, dr_rms, new_rnd);
1230
1231	for (i = 1; i <= n_girders; i++) {
1232	loc_gs = Elem_GetPos(gs, i); loc_ge = Elem_GetPos(ge, i);
1233	s_gs = Cell[loc_gs].S; s_ge = Cell[loc_ge].S;
1234
1235	// roll for a rigid boby
1236	// Note, girders needs to be introduced as gs->ge pairs
1237	Cell[loc_ge].Elem.M->PdTrnd = Cell[loc_gs].Elem.M->PdTrnd;
1238	Mpole_SetdT(ge, i);
1239
1240	for (k = 0; k <= 1; k++) {
1241	dx_gs[k] = Cell[loc_gs].dS[k]; dx_ge[k] = Cell[loc_ge].dS[k];
1242	}
1243
1244	// move elements onto mis-aligned girder
1245	for (j = loc_gs+1; j < loc_ge; j++) {
1246	if ((Cell[j].Elem.Pkind == Mpole) \|\| (Cell[j].Fnum == globval.bpm)) {
1247	s = Cell[j].S;
1248	for (k = 0; k <= 1; k++)
1249	Cell[j].Elem.M->PdSsys[k]
1250	= dx_gs[k] + (dx_ge[k]-dx_gs[k])*(s-s_gs)/(s_ge-s_gs);
1251	Cell[j].Elem.M->PdTsys =
1252	Cell[loc_gs].Elem.M->PdTrms*Cell[loc_gs].Elem.M->PdTrnd;
1253	}
1254	}
1255	}
1256	}
1257
1258	/*********************************************************************
1259	void misalign_sys_girders(const int gs, const int ge,
1260	const double dx_sys, const double dy_sys,
1261	const double dr_sys)
1262
1263	Purpose:
1264	Set systematic misalignment to the elements on the girders
1265
1266	Input:
1267	gs girder start marker
1268	ge girder end marker
1269	dx_sys sys value of the error in x
1270	dy_sys sys value of the error in y
1271	dr_sys sys value of the error in rotation around s
1272
1273	**********************************************************************/
1274	void misalign_sys_girders(const int gs, const int ge,
1275	const double dx_sys, const double dy_sys,
1276	const double dr_sys)
1277	{
1278	int i, k, n_girders, n_ge, n_gs;
1279	long int loc_gs, loc_ge, j;
1280	double s_gs, s_ge, dx_gs[2], dx_ge[2], s;
1281
1282	n_gs = GetnKid(gs); n_ge = GetnKid(ge);
1283
1284	if (n_gs == n_ge)
1285	n_girders = n_gs;
1286	else {
1287	cout << "set_girders: no of GS != no of GE" << endl;
1288	exit (1);
1289	}
1290
1291	misalign_sys_fam(gs, dx_sys, dy_sys, dr_sys);
1292	misalign_sys_fam(ge, dx_sys, dy_sys, dr_sys);
1293
1294	for (i = 1; i <= n_girders; i++) {
1295	loc_gs = Elem_GetPos(gs, i); loc_ge = Elem_GetPos(ge, i);
1296	s_gs = Cell[loc_gs].S; s_ge = Cell[loc_ge].S;
1297
1298	// roll for a rigid boby
1299	// Note, girders needs to be introduced as gs->ge pairs
1300	Cell[loc_ge].Elem.M->PdTrnd = Cell[loc_gs].Elem.M->PdTrnd;
1301	Mpole_SetdT(ge, i);
1302
1303	for (k = 0; k <= 1; k++) {
1304	dx_gs[k] = Cell[loc_gs].dS[k]; dx_ge[k] = Cell[loc_ge].dS[k];
1305	}
1306
1307	// move elements onto mis-aligned girder
1308	for (j = loc_gs+1; j < loc_ge; j++) {
1309	if ((Cell[j].Elem.Pkind == Mpole) \|\| (Cell[j].Fnum == globval.bpm)) {
1310	s = Cell[j].S;
1311	for (k = 0; k <= 1; k++)
1312	Cell[j].Elem.M->PdSsys[k]
1313	= dx_gs[k] + (dx_ge[k]-dx_gs[k])*(s-s_gs)/(s_ge-s_gs);
1314	Cell[j].Elem.M->PdTsys =
1315	Cell[loc_gs].Elem.M->PdTrms*Cell[loc_gs].Elem.M->PdTrnd;
1316	}
1317	}
1318	}
1319	}
1320
1321	/***************************************************************************
1322	void LoadAlignTol(const char *AlignFile, const bool Scale_it,
1323	const double Scale, const bool new_rnd, const int k)
1324
1325	Purpose:
1326	Load misalignment error from an external file, then set the error
1327	to the corresponding elements.
1328
1329	Input:
1330	AlignFile file from which the misalignment error is readed
1331	Scale_it flag to turn on/off the scale of the error
1332	Scale value to scale the error
1333	new_rnd use new random number or not
1334	k
1335
1336	15/10/2013 Jianfeng Zhang @ LAL
1337	Replace fgets() by getline() to fix the bug to read line of
1338	arbritrary length.
1339
1340	****************************************************************************/
1341	void LoadAlignTol(const char *AlignFile, const bool Scale_it,
1342	const double Scale, const bool new_rnd, const int k)
1343	{
1344	char Name[max_str]=" ", type[max_str]=" ";
1345	char *line = NULL;
1346	size_t len = 0;
1347	ssize_t read;
1348	int Fnum = 0, seed_val = 0, LineNum = 0;
1349	double dx = 0.0, dy = 0.0, dr = 0.0; // x and y misalignments [m] and roll error [rad]
1350	double dr_deg = 0.0;
1351	bool rms = false, set_rnd = false;
1352	bool prt = false;
1353	FILE *fp;
1354
1355	fp = file_read(AlignFile);
1356	if(fp == NULL){
1357	printf("LoadAlignTol(): Error! Failure to read file %s !\n",AlignFile);
1358	exit_(1);
1359	}
1360
1361	if (new_rnd)
1362	printf("set alignment errors\n");
1363	else
1364	printf("scale alignment errors: %4.2f\n", Scale);
1365
1366	while ((read=getline(&line, &len, fp)) != -1) {
1367
1368	LineNum++;
1369	if(!prt){
1370	printf("Line # %ld \n",LineNum);
1371	printf("Retrieved line of length %zu : \n",read);
1372	printf("%s \n",line);
1373	}
1374
1375	//check for whether to set seed
1376	if ((strstr(line, "#") == NULL) && line[0] != '\n' &&
1377	line[0]!='\r' && (strcmp(line, "\r\n") != 0)) {
1378	sscanf(line, "%s", Name);
1379	if (strcmp("seed", Name) == 0) {
1380	set_rnd = true;
1381	sscanf(line, "%*s %d", &seed_val);
1382	seed_val += 2*k;
1383	printf("setting random seed to %d\n", seed_val);
1384	iniranf(seed_val);
1385	} else {
1386	sscanf(line,"%*s %s %lf %lf %lf", type, &dx, &dy, &dr);
1387	dr_deg = dr*180.0/M_PI;
1388
1389	if (strcmp(type, "rms") == 0){
1390	rms = true;
1391	}
1392	else if (strcmp(type, "sys") == 0){
1393	rms = false;
1394	}
1395	else {
1396	printf("LoadAlignTol(): Error! Need to specify misalignment error type: <rms> or <sys> of element %s! \n",
1397	Name);
1398	exit_(1);
1399	}
1400
1401	if (rms && !set_rnd) {
1402	printf("LoadAlignTol(): Error! Seed of <rms> errors is not defined \n");
1403	exit_(1);
1404	}
1405
1406	if (Scale_it) {
1407	dx = Scale; dy = Scale; dr *= Scale;
1408	}
1409
1410	if(!prt){
1411	printf("Set <%s> misalignment errors to: [%s] \n",type, Name);
1412	printf("dx = %e, dy = %e, dr = %e \n \n",dx,dy,dr);
1413	}
1414
1415	if (strcmp("all", Name) == 0) {
1416
1417	if(rms)
1418	misalign_rms_type(All, dx, dy, dr_deg, new_rnd);
1419	else
1420	misalign_sys_type(All, dx, dy, dr_deg);
1421	} else if (strcmp("girder", Name) == 0) {
1422
1423	if (rms)
1424	misalign_rms_girders(globval.gs, globval.ge, dx, dy, dr_deg,
1425	new_rnd);
1426	else
1427	misalign_sys_girders(globval.gs, globval.ge, dx, dy, dr_deg);
1428	} else if (strcmp("dipole", Name) == 0) {
1429
1430	if (rms)
1431	misalign_rms_type(Dip, dx, dy, dr_deg, new_rnd);
1432	else
1433	misalign_sys_type(Dip, dx, dy, dr_deg);
1434	} else if (strcmp("quad", Name) == 0) {
1435
1436	if (rms)
1437	misalign_rms_type(Quad, dx, dy, dr_deg, new_rnd);
1438	else
1439	misalign_sys_type(Quad, dx, dy, dr_deg);
1440	} else if (strcmp("sext", Name) == 0) {
1441
1442	if (rms)
1443	misalign_rms_type(Sext, dx, dy, dr_deg, new_rnd);
1444	else
1445	misalign_sys_type(Sext, dx, dy, dr_deg);
1446	} else if (strcmp("bpm", Name) == 0) {
1447
1448	if (rms)
1449	misalign_rms_fam(globval.bpm, dx, dy, dr_deg, new_rnd);
1450	else
1451	misalign_sys_fam(globval.bpm, dx, dy, dr_deg);
1452	} else {
1453	Fnum = ElemIndex(Name);
1454	if(Fnum > 0) {
1455	if (rms)
1456	misalign_rms_fam(Fnum, dx, dy, dr_deg, new_rnd);
1457	else
1458	misalign_sys_fam(Fnum, dx, dy, dr_deg);
1459	} else{
1460	printf("LoadAlignTol(): Error! Undefined element: %s \n", Name);
1461	exit_(1);
1462	}
1463	}
1464	}
1465	} else
1466	continue;
1467	}
1468	free(line);
1469	fclose(fp);
1470	}
1471
1472	/****************************************************************************/
1473	/* void set_aper_elem(const int Fnum, const int Knum,
1474	const double Dxmin, const double Dxmax,
1475	const double Dymin, const double Dymax)
1476
1477	Purpose:
1478	Define vacuum chamber to the element
1479
1480
1481	Input:
1482	none
1483
1484	Output:
1485	none
1486
1487	Return:
1488	none
1489
1490	Global variables:
1491	none
1492
1493	Specific functions:
1494	none
1495
1496	Comments:
1497	Add k=0 for the default start point "start" in tracy.
1498
1499	****************************************************************************/
1500	// apertures
1501
1502	void set_aper_elem(const int Fnum, const int Knum,
1503	const double Dxmin, const double Dxmax,
1504	const double Dymin, const double Dymax)
1505	{
1506	int k;
1507
1508	if(Fnum==0 && Knum==0)
1509	k=0;
1510	else
1511	k = Elem_GetPos(Fnum, Knum);
1512
1513	Cell[k].maxampl[X_][0] = Dxmin;
1514	Cell[k].maxampl[X_][1] = Dxmax;
1515	Cell[k].maxampl[Y_][0] = Dymin;
1516	Cell[k].maxampl[Y_][1] = Dymax;
1517	}
1518
1519	/****************************************************************************/
1520	/* void set_aper_fam(const int Fnum,
1521	const double Dxmin, const double Dxmax,
1522	const double Dymin, const double Dymax)
1523	Purpose:
1524	Define vacuum chamber to the family "Fnum"
1525
1526
1527	Input:
1528	none
1529
1530	Output:
1531	none
1532
1533	Return:
1534	none
1535
1536	Global variables:
1537	none
1538
1539	Specific functions:
1540	none
1541
1542	Comments:
1543	Enumeration of special variables:
1544	enum { All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6 }
1545
1546	****************************************************************************/
1547	void set_aper_fam(const int Fnum,
1548	const double Dxmin, const double Dxmax,
1549	const double Dymin, const double Dymax)
1550	{
1551	int k;
1552
1553	for (k = 1; k <= GetnKid(Fnum); k++)
1554	set_aper_elem(Fnum, k, Dxmin, Dxmax, Dymin, Dymax);
1555	}
1556
1557	/****************************************************************************/
1558	/* void set_aper_type(const int type, const double Dxmin, const double Dxmax,
1559	const double Dymin, const double Dymax)
1560
1561	Purpose:
1562	Define vacuum chamber to the elements with the same "type", either for "all" element, or "quad","sext", etc.
1563
1564
1565	Input:
1566	none
1567
1568	Output:
1569	none
1570
1571	Return:
1572	none
1573
1574	Global variables:
1575	none
1576
1577	Specific functions:
1578	none
1579
1580	Comments:
1581	1) Enumeration of special variables:
1582	enum { All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6 }
1583	2) change the start point from k=1 to k=0, in order to set the aperture for
1584	the default start point "begin" which is not a lattice element.
1585
1586	****************************************************************************/
1587	void set_aper_type(const int type, const double Dxmin, const double Dxmax,
1588	const double Dymin, const double Dymax)
1589	{
1590	long int k;
1591
1592	if (type >= All && type <= HOMmax)
1593	{
1594	// for(k = 1; k <= globval.Cell_nLoc; k++)
1595	for(k = 0; k <= globval.Cell_nLoc; k++)
1596	if (((Cell[k].Elem.Pkind == Mpole) &&
1597	(Cell[k].Elem.M->n_design == type)) \|\| (type == All))
1598	set_aper_elem(Cell[k].Fnum, Cell[k].Knum, Dxmin, Dxmax, Dymin, Dymax);
1599	}
1600	else
1601	printf("set_aper_type: bad design type %d\n", type);
1602	}
1603
1604
1605	/****************************************************************************/
1606	/* void LoadApers(const char *AperFile, const double scl_x, const double scl_y)
1607
1608	Purpose:
1609	Read vacuum chamber from file "AperFile", and scale the values with scl_x and scl_y.
1610	Assign vacuum chamber to all elements, if there are definition of
1611	chamber at quad and sext, then assign the specific values at these elements.
1612
1613
1614	Input:
1615	AperFile: file with chamber definition
1616	scl_x : scale factor of x limitation
1617	scl_y : scale factor of y limitation
1618
1619	Output:
1620	none
1621
1622	Return:
1623	none
1624
1625	Global variables:
1626	none
1627
1628	Specific functions:
1629	none
1630
1631	Comments:
1632	Enumeration of special variables:
1633	enum { All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6 }
1634
1635	10/2010 Comments by Jianfeng Zhang
1636
1637	****************************************************************************/
1638	void LoadApers(const char *AperFile, const double scl_x, const double scl_y)
1639	{
1640	char line[max_str], Name[max_str];
1641	int Fnum;
1642	double dxmin, dxmax, dymin, dymax; // min and max x and apertures
1643	FILE *fp;
1644
1645	bool prt = true;
1646
1647	fp = file_read(AperFile);
1648
1649	printf("\n");
1650	printf("...Load and Set Apertures.\n");
1651
1652	while (fgets(line, max_str, fp) != NULL) {
1653	if (strstr(line, "#") == NULL)
1654	{
1655	sscanf(line,"%s %lf %lf %lf %lf",
1656	Name, &dxmin, &dxmax, &dymin, &dymax);
1657	dxmin *= scl_x;
1658	dxmax *= scl_x;
1659	dymin *= scl_y;
1660	dymax *= scl_y;
1661
1662	if (strcmp("all", Name)==0) {
1663	if(prt)
1664	printf("setting all apertures to"
1665	" dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1666	dxmin, dxmax, dymin, dymax);
1667	set_aper_type(All, dxmin, dxmax, dymin, dymax);
1668	// ini_aper(dxmin, dxmax, dymin, dymax);
1669	}
1670
1671	else if (strcmp("quad", Name)==0) {
1672	if(prt)
1673	printf("setting apertures at all quads to"
1674	" dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1675	dxmin, dxmax, dymin, dymax);
1676	set_aper_type(Quad, dxmin, dxmax, dymin, dymax);
1677	}
1678
1679	else if (strcmp("sext", Name) == 0) {
1680	if(prt)
1681	printf("setting apertures at all sextupoles to"
1682	" dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1683	dxmin, dxmax, dymin, dymax);
1684	set_aper_type(Sext, dxmin, dxmax, dymin, dymax);
1685	}
1686
1687	else {
1688	Fnum = ElemIndex(Name);
1689	if(Fnum > 0) {
1690	if(prt)
1691	printf("setting apertures at all %s to"
1692	" dxmin = %e, dxmax = %e, dymin = %e, dymax = %e\n",
1693	Name, dxmin, dxmax, dymin, dymax);
1694	set_aper_fam(Fnum, dxmin, dxmax, dymin, dymax);
1695	}
1696	else
1697	printf("LoadApers: lattice does not contain element %s\n", Name);
1698	}
1699	}
1700	else
1701	printf("%s", line);
1702	}
1703
1704	fclose(fp);
1705	}
1706
1707
1708	double get_L(const int Fnum, const int Knum)
1709	{
1710	return Cell[Elem_GetPos(Fnum, Knum)].Elem.PL;
1711	}
1712
1713
1714	void set_L(const int Fnum, const int Knum, const double L)
1715	{
1716
1717	Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L;
1718	}
1719
1720
1721	void set_L(const int Fnum, const double L)
1722	{
1723	int k;
1724
1725	for (k = 1; k <= GetnKid(Fnum); k++)
1726	set_L(Fnum, k, L);
1727	}
1728
1729
1730	void set_dL(const int Fnum, const int Knum, const double dL)
1731	{
1732
1733	Cell[Elem_GetPos(Fnum, Knum)].Elem.PL += dL;
1734	}
1735
1736
1737	// multipole components
1738	/*************************************************************
1739	void get_bn_design_elem(const int Fnum, const int Knum,
1740	const int n, double &bn, double &an)
1741
1742	Purpose:
1743	Get the n-th order design value (bn, an) of the element
1744	with family index "Fnum" and the kid number "Knum".
1745
1746	Input:
1747	Fnum: family index
1748	Knum: kid index
1749	n: n-th design order of the element
1750	bn: bn component of the element
1751	an: an component of the element
1752
1753	Output:
1754	None
1755
1756	Return:
1757	bn, an
1758
1759	Comments:
1760
1761	11/2010 comments by jianfeng Zhang
1762	**************************************************************/
1763	void get_bn_design_elem(const int Fnum, const int Knum,
1764	const int n, double &bn, double &an)
1765	{
1766	elemtype elem;
1767
1768	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1769
1770	bn = elem.M->PBpar[HOMmax+n]; an = elem.M->PBpar[HOMmax-n];
1771	}
1772
1773
1774	void get_bnL_design_elem(const int Fnum, const int Knum,
1775	const int n, double &bnL, double &anL)
1776	{
1777	elemtype elem;
1778
1779	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1780
1781	bnL = elem.M->PBpar[HOMmax+n]; anL = elem.M->PBpar[HOMmax-n];
1782
1783	if (elem.PL != 0.0) {
1784	bnL = elem.PL; anL = elem.PL;
1785	}
1786	}
1787
1788
1789	void set_bn_design_elem(const int Fnum, const int Knum,
1790	const int n, const double bn, const double an)
1791	{
1792	elemtype elem;
1793
1794	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1795
1796	elem.M->PBpar[HOMmax+n] = bn; elem.M->PBpar[HOMmax-n] = an;
1797
1798	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1799	}
1800
1801
1802	void set_dbn_design_elem(const int Fnum, const int Knum,
1803	const int n, const double dbn, const double dan)
1804	{
1805	elemtype elem;
1806
1807	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1808
1809	elem.M->PBpar[HOMmax+n] += dbn; elem.M->PBpar[HOMmax-n] += dan;
1810
1811	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1812	}
1813
1814
1815	void set_bn_design_fam(const int Fnum,
1816	const int n, const double bn, const double an)
1817	{
1818	int k;
1819
1820	for (k = 1; k <= GetnKid(Fnum); k++)
1821	set_bn_design_elem(Fnum, k, n, bn, an);
1822	}
1823
1824
1825	void set_dbn_design_fam(const int Fnum,
1826	const int n, const double dbn, const double dan)
1827	{
1828	int k;
1829
1830	for (k = 1; k <= GetnKid(Fnum); k++)
1831	set_dbn_design_elem(Fnum, k, n, dbn, dan);
1832	}
1833
1834
1835	void set_bnL_design_elem(const int Fnum, const int Knum,
1836	const int n, const double bnL, const double anL)
1837	{
1838	elemtype elem;
1839
1840	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1841
1842	if (elem.PL != 0.0) {
1843	elem.M->PBpar[HOMmax+n] = bnL/elem.PL;
1844	elem.M->PBpar[HOMmax-n] = anL/elem.PL;
1845	} else {
1846	// thin kick
1847	elem.M->PBpar[HOMmax+n] = bnL; elem.M->PBpar[HOMmax-n] = anL;
1848	}
1849
1850	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1851	}
1852
1853
1854	void set_dbnL_design_elem(const int Fnum, const int Knum,
1855	const int n, const double dbnL, const double danL)
1856	{
1857	elemtype elem;
1858
1859	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1860
1861	if (elem.PL != 0.0) {
1862	elem.M->PBpar[HOMmax+n] += dbnL/elem.PL;
1863	elem.M->PBpar[HOMmax-n] += danL/elem.PL;
1864	} else {
1865	// thin kick
1866	elem.M->PBpar[HOMmax+n] += dbnL; elem.M->PBpar[HOMmax-n] += danL;
1867	}
1868
1869	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
1870	}
1871
1872
1873	void set_dbnL_design_fam(const int Fnum,
1874	const int n, const double dbnL, const double danL)
1875	{
1876	int k;
1877
1878	for (k = 1; k <= GetnKid(Fnum); k++)
1879	set_bnL_design_elem(Fnum, k, n, dbnL, danL);
1880	}
1881
1882
1883	void set_bnL_design_fam(const int Fnum,
1884	const int n, const double bnL, const double anL)
1885	{
1886	int k;
1887
1888	for (k = 1; k <= GetnKid(Fnum); k++)
1889	set_bnL_design_elem(Fnum, k, n, bnL, anL);
1890	}
1891
1892
1893	void set_bnL_design_type(const int type,
1894	const int n, const double bnL, const double anL)
1895	{
1896	long int k;
1897
1898	if ((type >= Dip) && (type <= HOMmax)) {
1899	for (k = 1; k <= globval.Cell_nLoc; k++)
1900	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
1901	set_bnL_design_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL, anL);
1902	} else
1903	printf("Bad type argument to set_bnL_design_type()\n");
1904	}
1905
1906	/***********************************************************************
1907	void set_bnL_sys_elem(const int Fnum, const int Knum,
1908	const int n, const double bnL, const double anL)
1909
1910	Purpose:
1911	Set field error to the kid of a family, errors are absolute
1912	integrated field strength, replace the previous value!!!
1913
1914
1915	Input:
1916	Fnum family index
1917	Knum kids index
1918	n order of the error
1919	bnL absolute integrated B component for the n-th error
1920	anL absolute integrated A component for the n-th error
1921
1922
1923
1924	Output:
1925	None
1926
1927	Return:
1928	None
1929
1930	Global variables
1931	None
1932
1933	Specific functions:
1934	None
1935
1936	Comments:
1937	None
1938	**********************************************************************/
1939	void set_bnL_sys_elem(const int Fnum, const int Knum,
1940	const int n, const double bnL, const double anL)
1941	{
1942	elemtype elem;
1943
1944	const bool prt = false;
1945
1946	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
1947
1948	if (elem.PL != 0.0) {
1949	elem.M->PBsys[HOMmax+n] = bnL/elem.PL;
1950	elem.M->PBsys[HOMmax-n] = anL/elem.PL;
1951	} else {
1952	// thin kick
1953	elem.M->PBsys[HOMmax+n] = bnL; elem.M->PBsys[HOMmax-n] = anL;
1954	}
1955
1956	Mpole_SetPB(Fnum, Knum, n); //set for Bn component
1957	Mpole_SetPB(Fnum, Knum, -n); //set for An component
1958
1959	if (prt)
1960	printf("set the n=%d component of %s to %e %e\n",
1961	n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
1962	bnL, elem.M->PBsys[HOMmax+n]);
1963	}
1964
1965	/***********************************************************************
1966	void set_bnL_sys_fam(const int Fnum,
1967	const int n, const double bnL, const double anL)
1968
1969	Purpose:
1970	Set field error to all the kids of a family,errors are absolute
1971	integrated field strength
1972
1973
1974
1975	Input:
1976	Fnum family index
1977	n order of the error
1978	bnL absolute integrated B component for the n-th error
1979	anL absolute integrated A component for the n-th error
1980
1981
1982
1983	Output:
1984	None
1985
1986	Return:
1987	None
1988
1989	Global variables
1990	None
1991
1992	Specific functions:
1993	None
1994
1995	Comments:
1996	None
1997	**********************************************************************/
1998	void set_bnL_sys_fam(const int Fnum,
1999	const int n, const double bnL, const double anL)
2000	{
2001	int k;
2002
2003	for (k = 1; k <= GetnKid(Fnum); k++)
2004	set_bnL_sys_elem(Fnum, k, n, bnL, anL);
2005	}
2006
2007	/***********************************************************************
2008	void set_bnL_sys_type(const int type,
2009	const int n, const double bnL, const double anL)
2010
2011	Purpose:
2012	Set field error to the type, errors are absolute
2013	integrated field strength
2014
2015
2016
2017	Input:
2018	type type name
2019	n order of the error
2020	bnL absolute integrated B component for the n-th error
2021	anL absolute integrated A component for the n-th error
2022
2023
2024
2025	Output:
2026	None
2027
2028	Return:
2029	None
2030
2031	Global variables
2032	None
2033
2034	Specific functions:
2035	None
2036
2037	Comments:
2038	Attension:
2039	This code has a bug, in the for loop, when looking for the
2040	quadrupole family, it will also find the skew quadrupole family,
2041	since their Pkind == Mpole & n_design =2.
2042	It means the multipole error is additionally set N times,
2043	and N is the number of skew quadrupole in the lattice.
2044
2045	10/2010 Comments by Jianfeng Zhang
2046
2047	**********************************************************************/
2048	void set_bnL_sys_type(const int type,
2049	const int n, const double bnL, const double anL)
2050	{
2051	long int k;
2052
2053	if (type >= Dip && type <= HOMmax) {
2054	for(k = 1; k <= globval.Cell_nLoc; k++)
2055	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2056	set_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL, anL);
2057	} else
2058	printf("Bad type argument to set_bnL_sys_type()\n");
2059	}
2060
2061
2062	void set_bnL_rms_elem(const int Fnum, const int Knum,
2063	const int n, const double bnL, const double anL,
2064	const bool new_rnd)
2065	{
2066	elemtype elem;
2067
2068	bool normal = true, prt = false;
2069
2070	elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
2071
2072	if (elem.PL != 0.0) {
2073	elem.M->PBrms[HOMmax+n] = bnL/elem.PL;
2074	elem.M->PBrms[HOMmax-n] = anL/elem.PL;
2075	} else {
2076	// thin kick
2077	elem.M->PBrms[HOMmax+n] = bnL; elem.M->PBrms[HOMmax-n] = anL;
2078	}
2079
2080	if(new_rnd){
2081	if (normal) {
2082	elem.M->PBrnd[HOMmax+n] = normranf();
2083	elem.M->PBrnd[HOMmax-n] = normranf();
2084	} else {
2085	elem.M->PBrnd[HOMmax+n] = ranf(); elem.M->PBrnd[HOMmax-n] = ranf();
2086	}
2087	}
2088
2089	if (prt)
2090	printf("set_bnL_rms_elem: Fnum = %d, Knum = %d"
2091	", bnL = %e, anL = %e %e %e\n",
2092	Fnum, Knum, bnL, anL,
2093	elem.M->PBrms[HOMmax+n], elem.M->PBrms[HOMmax-n]);
2094
2095	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
2096	}
2097
2098
2099	void set_bnL_rms_fam(const int Fnum,
2100	const int n, const double bnL, const double anL,
2101	const bool new_rnd)
2102	{
2103	int k;
2104
2105	for (k = 1; k <= GetnKid(Fnum); k++)
2106	set_bnL_rms_elem(Fnum, k, n, bnL, anL, new_rnd);
2107	}
2108
2109
2110	void set_bnL_rms_type(const int type,
2111	const int n, const double bnL, const double anL,
2112	const bool new_rnd)
2113	{
2114	long int k;
2115
2116	if (type >= Dip && type <= HOMmax) {
2117	for(k = 1; k <= globval.Cell_nLoc; k++)
2118	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2119	set_bnL_rms_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL, anL, new_rnd);
2120	} else
2121	printf("Bad type argument to set_bnL_rms_type()\n");
2122	}
2123
2124	/***********************************************************************
2125	void set_bnr_sys_elem(const int Fnum, const int Knum,
2126	const int n, const double bnr, const double anr)
2127
2128	Purpose:
2129	Set field error to the kid of a family, errors are relative to
2130	design values for (Dip, Quad, Sext, ...)
2131
2132	If r0 != 0, call this function
2133
2134	Input:
2135	Fnum family index
2136	Knum kids index
2137	n order of the error
2138	bnL relative integrated B component for the n-th error
2139	anL relative integrated A component for the n-th error
2140
2141
2142
2143	Output:
2144	None
2145
2146	Return:
2147	None
2148
2149	Global variables
2150	None
2151
2152	Specific functions:
2153	None
2154
2155	Comments:
2156	!!!!!!!!!
2157	This function desn't work for dipole, since
2158	mp->PBpar[HOMmax+1] is not assigned a value when dipole is
2159	read in t2lat.cc.
2160
2161	comments by Jianfeng Zhang 10/2010 @ soleil
2162	**********************************************************************/
2163	void set_bnr_sys_elem(const int Fnum, const int Knum,
2164	const int n, const double bnr, const double anr)
2165	{
2166	int nd;
2167	MpoleType *mp;
2168	bool prt = false;
2169
2170	mp = Cell[Elem_GetPos(Fnum, Knum)].Elem.M; nd = mp->n_design;
2171	// errors are relative to design values for (Dip, Quad, Sext, ...)
2172	mp->PBsys[HOMmax+n] = bnr*mp->PBpar[HOMmax+nd];
2173	mp->PBsys[HOMmax-n] = anr*mp->PBpar[HOMmax+nd];
2174	//update the total field strength PB and maximum order p_order
2175	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
2176
2177	if (prt)
2178	printf("set the n=%d component of %s to %e %e %e\n",
2179	n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
2180	bnr, mp->PBpar[HOMmax+nd], mp->PBsys[HOMmax+n]);
2181	}
2182
2183	/***********************************************************************
2184	void set_bnr_sys_fam(const int Fnum,
2185	const int n, const double bnr, const double anr)
2186
2187	Purpose:
2188	Set field error to all the kids of a family, errors are relative
2189	to design values for (Dip, Quad, Sext, ...).
2190
2191	If r0 != 0, call this function
2192
2193	Input:
2194	Fnum family index
2195	n order of the error
2196	bnL relative integrated B component for the n-th error
2197	anL relative integrated A component for the n-th error
2198
2199
2200
2201	Output:
2202	None
2203
2204	Return:
2205	None
2206
2207	Global variables
2208	None
2209
2210	Specific functions:
2211	None
2212
2213	Comments:
2214	None
2215	**********************************************************************/
2216	void set_bnr_sys_fam(const int Fnum,
2217	const int n, const double bnr, const double anr)
2218	{
2219	int k;
2220
2221	for (k = 1; k <= GetnKid(Fnum); k++)
2222	set_bnr_sys_elem(Fnum, k, n, bnr, anr);
2223	}
2224
2225	/***********************************************************************
2226	void set_bnr_sys_type(const int type,
2227	const int n, const double bnr, const double anr)
2228
2229	Purpose:
2230	Set field error to the type, errors are relative
2231	to design values for (Dip, Quad, Sext, ...)
2232
2233	If r0 != 0, call this function
2234
2235	Input:
2236	type type name
2237	n order of the error
2238	bnL relative integrated B component for the n-th error
2239	anL relative integrated A component for the n-th error
2240
2241
2242
2243	Output:
2244	None
2245
2246	Return:
2247	None
2248
2249	Global variables
2250	None
2251
2252	Specific functions:
2253	None
2254
2255	Comments:
2256	None
2257	**********************************************************************/
2258	void set_bnr_sys_type(const int type,
2259	const int n, const double bnr, const double anr)
2260	{
2261	long int k;
2262
2263	if (type >= Dip && type <= HOMmax) {
2264	for(k = 1; k <= globval.Cell_nLoc; k++)
2265	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2266	set_bnr_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, bnr, anr);
2267	} else
2268	printf("Bad type argument to set_bnr_sys_type()\n");
2269	}
2270
2271	/***********************************************************************
2272	void set_bnr_rms_elem(const int Fnum, const int Knum,
2273	const int n, const double bnr, const double anr,
2274	const bool new_rnd)
2275
2276	Purpose:
2277	Update rms field error PBrms and random value of the rms field error PBrnd
2278	to a elment, errors are relative to design values PBpar for (Dip, Quad,
2279	Sext, ...).
2280
2281	If r0 != 0, call this function
2282
2283	Input:
2284	Fnum family index
2285	Knum kid index
2286	n order of the error
2287	bnL relative integrated B component for the n-th error
2288	anL relative integrated A component for the n-th error
2289	new_rnd true, generate random value for the rms field error
2290	false, NOT generate random value for the rms field error
2291
2292	Output:
2293	None
2294
2295	Return:
2296	None
2297
2298	Global variables
2299	None
2300
2301	Specific functions:
2302	None
2303
2304	Comments:
2305	None
2306	**********************************************************************/
2307	void set_bnr_rms_elem(const int Fnum, const int Knum,
2308	const int n, const double bnr, const double anr,
2309	const bool new_rnd)
2310	{
2311	int nd;
2312	MpoleType *mp;
2313
2314	bool normal = true, prt = false;
2315
2316	mp = Cell[Elem_GetPos(Fnum, Knum)].Elem.M; nd = mp->n_design;
2317	// errors are relative to design values for (Dip, Quad, Sext, ...)
2318	mp->PBrms[HOMmax+n] = bnr*mp->PBpar[HOMmax+nd];
2319	mp->PBrms[HOMmax-n] = anr*mp->PBpar[HOMmax+nd];
2320	if(new_rnd){
2321	if (normal) {
2322	mp->PBrnd[HOMmax+n] = normranf(); mp->PBrnd[HOMmax-n] = normranf();
2323	} else {
2324	mp->PBrnd[HOMmax+n] = ranf(); mp->PBrnd[HOMmax-n] = ranf();
2325	}
2326	}
2327
2328	if (prt)
2329	printf("set_bnr_rms_elem: Fnum = %d, Knum = %d"
2330	", bnr = %e, anr = %e %e %e\n",
2331	Fnum, Knum, bnr, anr, mp->PBrms[HOMmax+n], mp->PBrms[HOMmax-n]);
2332
2333	Mpole_SetPB(Fnum, Knum, n); Mpole_SetPB(Fnum, Knum, -n);
2334	}
2335
2336	/***********************************************************************
2337	void set_bnr_rms_fam(const int Fnum,
2338	const int n, const double bnr, const double anr,
2339	const bool new_rnd)
2340
2341	Purpose:
2342	Update rms field error PBrms and random value of the rms field error PBrnd
2343	to all the kid elments in a family, errors are relative to design values PBpar
2344	for (Dip, Quad, Sext, ...).
2345
2346	If r0 != 0, call this function
2347
2348	Input:
2349	Fnum family index
2350	n order of the error
2351	bnr relative integrated B component for the n-th error
2352	anr relative integrated A component for the n-th error
2353	new_rnd true, generate random value for the rms field error
2354	false, NOT generate random value for the rms field error
2355
2356	Output:
2357	None
2358
2359	Return:
2360	None
2361
2362	Global variables
2363	None
2364
2365	Specific functions:
2366	None
2367
2368	Comments:
2369	None
2370	**********************************************************************/
2371	void set_bnr_rms_fam(const int Fnum,
2372	const int n, const double bnr, const double anr,
2373	const bool new_rnd)
2374	{
2375	int k;
2376
2377	for (k = 1; k <= GetnKid(Fnum); k++)
2378	set_bnr_rms_elem(Fnum, k, n, bnr, anr, new_rnd);
2379	}
2380
2381	/***********************************************************************
2382	void set_bnr_rms_type(const int type,
2383	const int n, const double bnr, const double anr,
2384	const bool new_rnd)
2385
2386	Purpose:
2387	Update rms field error PBrms and random value of the rms field error PBrnd
2388	to all the kid elments in a family, errors are relative to design values PBpar
2389	for (Dip, Quad, Sext, ...).
2390
2391	If r0 != 0, call this function
2392
2393	Input:
2394	type type of the elment,maybe 'dipole', 'quadrupole', 'sextupole', 'all'
2395	n order of the error
2396	bnr relative integrated B component for the n-th error
2397	anr relative integrated A component for the n-th error
2398	new_rnd true, generate random value for the rms field error
2399	false, NOT generate random value for the rms field error
2400
2401	Output:
2402	None
2403
2404	Return:
2405	None
2406
2407	Global variables
2408	None
2409
2410	Specific functions:
2411	None
2412
2413	Comments:
2414	None
2415	**********************************************************************/
2416	void set_bnr_rms_type(const int type,
2417	const int n, const double bnr, const double anr,
2418	const bool new_rnd)
2419	{
2420	long int k;
2421
2422	if (type >= Dip && type <= HOMmax) {
2423	for(k = 1; k <= globval.Cell_nLoc; k++)
2424	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == type))
2425	set_bnr_rms_elem(Cell[k].Fnum, Cell[k].Knum, n, bnr, anr, new_rnd);
2426	} else
2427	printf("Bad type argument to set_bnr_rms_type()\n");
2428	}
2429
2430
2431	double get_Wiggler_BoBrho(const int Fnum, const int Knum)
2432	{
2433	return Cell[Elem_GetPos(Fnum, Knum)].Elem.W->BoBrhoV[0];
2434	}
2435
2436	/****************************************************************************/
2437	/* void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
2438
2439	Purpose:
2440	Set the integrated field strength for the specific element in the family
2441	of wiggers.
2442
2443	Input:
2444	Fnum family name
2445	BoBrhoV integrated field strength
2446
2447	Output:
2448	none
2449
2450	Return:
2451	none
2452
2453	Global variables:
2454	none
2455
2456	Specific functions:
2457	none
2458
2459	Comments:
2460	Called by set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
2461
2462	****************************************************************************/
2463	void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
2464	{
2465	Cell[Elem_GetPos(Fnum, Knum)].Elem.W->BoBrhoV[0] = BoBrhoV;
2466	Cell[Elem_GetPos(Fnum, Knum)].Elem.W->PBW[HOMmax+Quad] = -sqr(BoBrhoV)/2.0;
2467	Wiggler_SetPB(Fnum, Knum, Quad);
2468	}
2469
2470	/****************************************************************************/
2471	/* void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
2472
2473	Purpose:
2474	Set the integrated field strength for the family of wiggers
2475
2476	Input:
2477	Fnum family name
2478	BoBrhoV integrated field strength
2479
2480	Output:
2481	none
2482
2483	Return:
2484	none
2485
2486	Global variables:
2487	none
2488
2489	Specific functions:
2490	none
2491
2492	Comments:
2493	none
2494
2495	****************************************************************************/
2496	void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
2497	{
2498	int k;
2499
2500	for (k = 1; k <= GetnKid(Fnum); k++)
2501	set_Wiggler_BoBrho(Fnum, k, BoBrhoV);
2502	}
2503
2504	/****************************************************************************/
2505	/* void set_ID_scl(const int Fnum, const int Knum, const double scl)
2506
2507	Purpose:
2508	Set the scale factor for the spicified element of wigger, insertion devices,
2509	or fieldMap.
2510	Called by set_ID_scl(const int Fnum, const double scl).
2511
2512	Input:
2513	Fnum family name
2514	Knum kid index of Fname
2515	scl scaling factor
2516
2517	Output:
2518	none
2519
2520	Return:
2521	none
2522
2523	Global variables:
2524	none
2525
2526	Specific functions:
2527	none
2528
2529	Comments:
2530	none
2531
2532	****************************************************************************/
2533	void set_ID_scl(const int Fnum, const int Knum, const double scl)
2534	{
2535	int k;
2536	WigglerType* W;
2537
2538	switch (Cell[Elem_GetPos(Fnum, Knum)].Elem.Pkind) {
2539	case Wigl:
2540	// scale the ID field
2541	W = Cell[Elem_GetPos(Fnum, Knum)].Elem.W;
2542	for (k = 0; k < W->n_harm; k++) {
2543	W->BoBrhoH[k] = scl*ElemFam[Fnum-1].ElemF.W->BoBrhoH[k];
2544	W->BoBrhoV[k] = scl*ElemFam[Fnum-1].ElemF.W->BoBrhoV[k];
2545	}
2546	break;
2547	case Insertion:
2548	Cell[Elem_GetPos(Fnum, Knum)].Elem.ID->scaling1 = scl;
2549	Cell[Elem_GetPos(Fnum, Knum)].Elem.ID->scaling2 = scl;
2550	break;
2551	case FieldMap:
2552	Cell[Elem_GetPos(Fnum, Knum)].Elem.FM->scl = scl;
2553	break;
2554	default:
2555	cout << "set_ID_scl: unknown element type" << endl;
2556	exit_(1);
2557	break;
2558	}
2559	}
2560
2561	/****************************************************************************/
2562	/* void set_ID_scl(const int Fnum, const double scl)
2563
2564	Purpose:
2565	Set the scale factor for the specified wigger, insertion devices, or fieldMap.
2566
2567	Input:
2568	Fnum family name
2569	scl scaling factor
2570
2571	Output:
2572	none
2573
2574	Return:
2575	none
2576
2577	Global variables:
2578	none
2579
2580	Specific functions:
2581	none
2582
2583	Comments:
2584	none
2585
2586	****************************************************************************/
2587	void set_ID_scl(const int Fnum, const double scl)
2588	{
2589	int k;
2590
2591	for (k = 1; k <= GetnKid(Fnum); k++)
2592	set_ID_scl(Fnum, k, scl);
2593	}
2594
2595	/***********************************************************************
2596	void SetFieldValues_fam(const int Fnum, const bool rms, const double r0,
2597	const int n, const double Bn, const double An,
2598	const bool new_rnd)
2599
2600	Purpose:
2601	Set field error to all the kids in a family
2602
2603	Input:
2604	Fnum family name
2605	rms flag to set rms field error
2606	r0 radius at which error is measured
2607	n order of the error
2608	Bn integrated B component for the n-th error
2609	An integrated A component for the n-th error
2610	new_rnd bool flag to set new random number
2611
2612
2613	Output:
2614	None
2615
2616	Return:
2617	None
2618
2619	Global variables
2620	None
2621
2622	Specific functions:
2623	None
2624
2625	Comments:
2626	None
2627	**********************************************************************/
2628	void SetFieldValues_fam(const int Fnum, const bool rms, const double r0,
2629	const int n, const double Bn, const double An,
2630	const bool new_rnd)
2631	{
2632	int N;
2633	double bnr, anr;
2634
2635	N = Cell[Elem_GetPos(Fnum, 1)].Elem.M->n_design;
2636	if (r0 == 0.0) {
2637	// input is: (b_n L), (a_n L)
2638	if(rms)
2639	set_bnL_rms_fam(Fnum, n, Bn, An, new_rnd);
2640	else
2641	set_bnL_sys_fam(Fnum, n, Bn, An);
2642	} else {
2643	bnr = Bn/pow(r0, n-N); anr = An/pow(r0, n-N);
2644	if(rms)
2645	set_bnr_rms_fam(Fnum, n, bnr, anr, new_rnd);
2646	else
2647	set_bnr_sys_fam(Fnum, n, bnr, anr);
2648	}
2649	}
2650
2651
2652	/***********************************************************************
2653	void SetFieldValues_type(const int N, const bool rms, const double r0,
2654	const int n, const double Bn, const double An,
2655	const bool new_rnd)
2656
2657	Purpose:
2658	Set field error to the type
2659
2660	Input:
2661	N type name
2662	rms flag to set rms field error
2663	r0 radius at which error is measured
2664	n order of the error
2665	Bn integrated B component for the n-th error
2666	An integrated A component for the n-th error
2667	new_rnd bool flag to set new random number
2668
2669
2670	Output:
2671	None
2672
2673	Return:
2674	None
2675
2676	Global variables
2677	None
2678
2679	Specific functions:
2680	None
2681
2682	Comments:
2683	None
2684	**********************************************************************/
2685	void SetFieldValues_type(const int N, const bool rms, const double r0,
2686	const int n, const double Bn, const double An,
2687	const bool new_rnd)
2688	{
2689	double bnr, anr;
2690
2691	if (r0 == 0.0) {
2692	// input is: (b_n L), (a_n L)
2693	if(rms)
2694	set_bnL_rms_type(N, n, Bn, An, new_rnd);
2695	else
2696	set_bnL_sys_type(N, n, Bn, An);
2697	} else {
2698	bnr = Bn/pow(r0, n-N); anr = An/pow(r0, n-N);
2699	if(rms)
2700	set_bnr_rms_type(N, n, bnr, anr, new_rnd);
2701	else
2702	set_bnr_sys_type(N, n, bnr, anr);
2703	}
2704	}
2705
2706
2707	/***********************************************************************
2708	void SetFieldErrors(const char *name, const bool rms, const double r0,
2709	const int n, const double Bn, const double An,
2710	const bool new_rnd)
2711
2712	Purpose:
2713	Set field error to the type or element
2714
2715	Input:
2716	name type name or element name
2717	rms flag to set rms field error
2718	r0 radius at which error is measured
2719	n order of the error
2720	Bn integrated B component for the n-th error
2721	An integrated A component for the n-th error
2722	new_rnd flag to set new random number for the error
2723
2724
2725	Output:
2726	None
2727
2728	Return:
2729	None
2730
2731	Global variables
2732	None
2733
2734	Specific functions:
2735	None
2736
2737	Comments:
2738	10/2010 Comments by Jianfeng Zhang
2739	**********************************************************************/
2740	void SetFieldErrors(const char *name, const bool rms, const double r0,
2741	const int n, const double Bn, const double An,
2742	const bool new_rnd)
2743	{
2744	int Fnum;
2745
2746	if (strcmp("all", name) == 0) {
2747	printf("all: not yet implemented\n");
2748	} else if (strcmp("dip", name) == 0) {
2749	SetFieldValues_type(Dip, rms, r0, n, Bn, An, new_rnd);
2750	} else if (strcmp("quad", name) == 0) {
2751	SetFieldValues_type(Quad, rms, r0, n, Bn, An, new_rnd);
2752	} else if (strcmp("sext", name) == 0) {
2753	SetFieldValues_type(Sext, rms, r0, n, Bn, An, new_rnd);
2754	} else {
2755	Fnum = ElemIndex(name);
2756	if(Fnum > 0)
2757	SetFieldValues_fam(Fnum, rms, r0, n, Bn, An, new_rnd);
2758	else
2759	printf("SetFieldErrors: undefined element %s\n", name);
2760	}
2761	}
2762
2763
2764
2765	/********************************************************************
2766	char* get_prm(void)
2767
2768	Purpose:
2769	Find the stoke talbe symbol "\t" in the string, and
2770	then saved and return the position of the pointer.
2771
2772	symbol for end of the line:
2773	\n
2774	Comments:
2775	10/2010 By Jianfeng Zhang
2776	Fix the bug for reading the end of line symbol "\n" , "\r",'\r\n'
2777	at different operation system
2778	*********************************************************************/
2779	char* get_prm(void)
2780	{
2781	char *prm;
2782
2783	prm = strtok(NULL, " \t");
2784
2785	//if ((prm == NULL) \|\| (strcmp(prm, "\r\n") == 0)) {
2786	if ((prm == NULL) \|\| (strcmp(prm, "\n") == 0)\|\| (strcmp(prm, "\r") == 0)\|\| (strcmp(prm, "\r\n") == 0)) {
2787	printf("get_prm: incorrect format\n");
2788	exit_(1);
2789	}
2790
2791	return prm;
2792	}
2793
2794	/**************************************************************************
2795	* void LoadFieldErr(const char *FieldErrorFile, const bool Scale_it,
2796	const double Scale, const bool new_rnd)
2797	*
2798	*
2799	*
2800	**************************************************************************/
2801	void LoadFieldErr(const char *FieldErrorFile, const bool Scale_it,
2802	const double Scale, const bool new_rnd)
2803	{
2804	bool rms, set_rnd;
2805	char line[max_str], name[max_str], type[max_str], *prm;
2806	int k, n, seed_val;
2807	double Bn, An, r0;
2808	FILE *inf;
2809
2810	const bool prt = true;
2811
2812	inf = file_read(FieldErrorFile);
2813
2814	set_rnd = false;
2815	printf("\n");
2816	while (fgets(line, max_str, inf) != NULL) {
2817	if (strstr(line, "#") == NULL) {
2818	// check for whether to set new seed
2819	sscanf(line, "%s", name);
2820	if (strcmp("seed", name) == 0) {
2821	set_rnd = true;
2822	sscanf(line, "%*s %d", &seed_val);
2823	// printf("setting random seed to %d\n", seed_val);
2824	// iniranf(seed_val);
2825	printf("\n");
2826	printf("setting new random numbers (%d)\n", seed_val);
2827	} else {
2828	sscanf(line, "%*s %s %lf", type, &r0);
2829	printf("%-4s %3s %7.1le", name, type, r0);
2830	rms = (strcmp("rms", type) == 0)? true : false;
2831	if (rms && !set_rnd) {
2832	printf("LoadFieldErr: seed not defined\n");
2833	exit_(1);
2834	}
2835	// skip first three parameters
2836	strtok(line, " \t");
2837	for (k = 1; k <= 2; k++)
2838	strtok(NULL, " \t");
2839	while (((prm = strtok(NULL, " \t")) != NULL) &&
2840	(strcmp(prm, "\r\n") != 0)) {
2841	sscanf(prm, "%d", &n);
2842	prm = get_prm();
2843	sscanf(prm, "%lf", &Bn);
2844	prm = get_prm();
2845	sscanf(prm, "%lf", &An);
2846	if (Scale_it)
2847	{Bn = Scale; An = Scale;}
2848	if (prt)
2849	printf(" %2d %9.1e %9.1e\n", n, Bn, An);
2850	// convert to normalized multipole components
2851	SetFieldErrors(name, rms, r0, n, Bn, An, true);
2852	}
2853	}
2854	} else
2855	printf("%s", line);
2856	}
2857
2858	fclose(inf);
2859	}
2860
2861
2862	// correction algorithms
2863
2864	// control of orbit
2865
2866	// closed orbit correction by n_orbit iterations
2867	bool CorrectCOD(int n_orbit)
2868	{
2869	bool cod;
2870	int i;
2871	long int lastpos;
2872	Vector2 mean, sigma, max;
2873
2874	cod = getcod(0.0, lastpos);
2875	if (cod) {
2876	codstat(mean, sigma, max, globval.Cell_nLoc, false); // get orbit stats
2877	printf("\n");
2878	printf("Initial RMS orbit (BPMs): x = %7.1e mm, y = %7.1e mm\n",
2879	1e3sigma[X_], 1e3sigma[Y_]);
2880	codstat(mean, sigma, max, globval.Cell_nLoc, true);
2881	printf("Initial RMS orbit (all): x = %7.1e mm, y = %7.1e mm\n",
2882	1e3sigma[X_], 1e3sigma[Y_]);
2883
2884	for (i = 0; i < n_orbit; i++){
2885	lsoc(1, globval.bpm, globval.hcorr, 1); // correct horizontal orbit
2886	lsoc(1, globval.bpm, globval.vcorr, 2); // correct vertical orbit
2887	cod = getcod(0.0, lastpos); // find closed orbit
2888	if (!cod) break;
2889	}
2890
2891	if (cod) {
2892	codstat(mean, sigma, max, globval.Cell_nLoc, false);
2893	printf("Corrected RMS orbit (BPMs): x = %7.1e mm, y = %7.1e mm\n",
2894	1e3sigma[X_], 1e3sigma[Y_]);
2895	codstat(mean, sigma, max, globval.Cell_nLoc, true);
2896	printf("Corrected RMS orbit (all): x = %7.1e mm, y = %7.1e mm\n",
2897	1e3sigma[X_], 1e3sigma[Y_]);
2898	}
2899	}
2900
2901	return cod;
2902	}
2903
2904
2905	void Align_BPMs(const int n)
2906	{
2907	// Align BPMs to adjacent multipoles.
2908
2909	bool aligned;
2910	int i, j, k;
2911	long int loc;
2912
2913	const int n_step = 5;
2914
2915	printf("\n");
2916	for (i = 1; i <= GetnKid(globval.bpm); i++) {
2917	loc = Elem_GetPos(globval.bpm, i);
2918
2919	if ((loc == 1) \|\| (loc == globval.Cell_nLoc)) {
2920	printf("Align_BPMs: BPM at entrance or exit of lattice: %ld\n", loc);
2921	exit_(1);
2922	}
2923
2924	j = 1; aligned = false;
2925	do {
2926	if ((Cell[loc-j].Elem.Pkind == Mpole) &&
2927	(Cell[loc-j].Elem.M->n_design == n)) {
2928	for (k = 0; k <= 1; k++)
2929	Cell[loc].Elem.M->PdSsys[k] = Cell[loc-j].dS[k];
2930	printf("aligned BPM no %1d to %s\n", i, Cell[loc-j].Elem.PName);
2931	aligned = true; break;
2932	} else if ((Cell[loc+j].Elem.Pkind == Mpole) &&
2933	(Cell[loc+j].Elem.M->n_design == n)) {
2934	for (k = 0; k <= 1; k++)
2935	Cell[loc].Elem.M->PdSsys[k] = Cell[loc+j].dS[k];
2936	printf("aligned BPM no %1d to %s\n", i, Cell[loc+j].Elem.PName);
2937	aligned = true; break;
2938	}
2939
2940	j++;
2941	} while (j <= n_step);
2942
2943	if (aligned)
2944	Mpole_SetdS(globval.bpm, i);
2945	else
2946	printf("Align_BPMs: no multipole adjacent to BPM no %d\n", i);
2947	}
2948	}
2949
2950	/* *****************************************************
2951	void get_bare()
2952
2953	Purpose:
2954	store values of the optics function at the sextupoles
2955	*/
2956	void get_bare()
2957	{
2958	long int j, k;
2959
2960	n_sext = 0;
2961	for (j = 0; j <= globval.Cell_nLoc; j++) {
2962	if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design >= Sext)) {
2963	n_sext++; sexts[n_sext-1] = j;
2964	for (k = 0; k <= 1; k++) {
2965	beta0_[n_sext-1][k] = Cell[j].Beta[k];
2966	nu0s[n_sext-1][k] = Cell[j].Nu[k];
2967	}
2968	}
2969	}
2970
2971	// beta-functions for normalization
2972	beta_ref[X_] = Cell[globval.Cell_nLoc].Beta[X_];
2973	beta_ref[Y_] = Cell[globval.Cell_nLoc].Beta[Y_];
2974
2975	nu0_[X_] = globval.TotalTune[X_]; nu0_[Y_] = globval.TotalTune[Y_];
2976	}
2977
2978
2979	void get_dbeta_dnu(double m_dbeta[], double s_dbeta[],
2980	double m_dnu[], double s_dnu[])
2981	{
2982	int k;
2983	long int j, ind;
2984	double dbeta, dnu;
2985
2986	Ring_GetTwiss(false, 0.0);
2987
2988	for (k = 0; k <= 1; k++) {
2989	m_dbeta[k] = 0.0; s_dbeta[k] = 0.0; m_dnu[k] = 0.0; s_dnu[k] = 0.0;
2990	}
2991
2992	for (j = 0; j < n_sext; j++) {
2993	ind = sexts[j];
2994	for (k = 0; k <= 1; k++) {
2995	dbeta = (Cell[ind].Beta[k]-beta0_[j][k])/beta0_[j][k];
2996	m_dbeta[k] += dbeta; s_dbeta[k] += sqr(dbeta);
2997	dnu = Cell[ind].Nu[k] - nu0s[j][k];
2998	m_dnu[k] += dnu; s_dnu[k] += sqr(dnu);
2999	}
3000	}
3001
3002	for (k = 0; k <= 1; k++) {
3003	m_dbeta[k] /= n_sext; m_dnu[k] /= n_sext;
3004	s_dbeta[k] = sqrt((s_dbeta[k]-n_sext*sqr(m_dbeta[k]))/(n_sext-1));
3005	s_dnu[k] = sqrt((s_dnu[k]-n_sext*sqr(m_dnu[k]))/(n_sext-1));
3006	}
3007	}
3008
3009
3010	bool CorrectCOD_N(const char *ae_file, const int n_orbit,
3011	const int n, const int k)
3012	{
3013	bool cod = false;
3014	int i;
3015	double m_dbeta[2], s_dbeta[2], m_dnu[2], s_dnu[2];
3016
3017	// Clear trim setpoints
3018	set_bnL_design_fam(globval.hcorr, Dip, 0.0, 0.0);
3019	set_bnL_design_fam(globval.vcorr, Dip, 0.0, 0.0);
3020
3021	// load misalignments
3022	LoadAlignTol(ae_file, true, 1.0, true, k);
3023	for (i = 1; i <= n; i++) {
3024	// Scale the rms values
3025	LoadAlignTol(ae_file, true, (double)i/(double)n, false, k);
3026
3027	if (bba) {
3028	// Beam based alignment
3029	Align_BPMs(Quad);
3030	}
3031
3032	cod = CorrectCOD(n_orbit);
3033
3034	if (!cod) break;
3035
3036	get_dbeta_dnu(m_dbeta, s_dbeta, m_dnu, s_dnu);
3037	printf("\n");
3038	printf("RMS dbeta_x/beta_x = %4.2f%%, dbeta_y/beta_y = %4.2f%%\n",
3039	1e2s_dbeta[X_], 1e2s_dbeta[Y_]);
3040	printf("RMS dnu_x = %7.5f, dnu_y = %7.5f\n",
3041	s_dnu[X_], s_dnu[Y_]);
3042	}
3043
3044	return cod;
3045	}
3046
3047
3048	// Control of vertical beam size
3049
3050	void FindSQ_SVDmat(double SkewRespMat, double U,
3051	double *V, double w, int N_COUPLE, int N_SKEW)
3052	{
3053	int i, j;
3054
3055	const double cut = 1e-10; // cut value for SVD singular values
3056
3057	for (i = 1; i <= N_COUPLE; i++)
3058	for (j = 1; j <= N_SKEW; j++)
3059	U[i][j] = SkewRespMat[i][j];
3060
3061	// prepare matrices for SVD
3062	dsvdcmp(U, N_COUPLE, N_SKEW, w, V);
3063
3064	// zero singular values
3065	printf("\n");
3066	printf("singular values:\n");
3067	printf("\n");
3068	for (i = 1; i <= N_SKEW; i++) {
3069	printf("%11.3e", w[i]);
3070	if (w[i] < cut ) {
3071	w[i] = 0.0;
3072	printf(" (zeroed)");
3073	if (i % 8 == 0) printf("\n");
3074	}
3075	}
3076	if (i % 8 != 0) printf("\n");
3077	}
3078
3079
3080	void FindMatrix(double **SkewRespMat, const double deta_y_max)
3081	{
3082	// Ring_GetTwiss(true, 0.0) should be called in advance
3083	int i, j, k;
3084	long int loc;
3085	double nuX, nuY, alpha, eta_y_max;
3086	double *etaSQ;
3087	double betaSQ, nuSQ, betaBPM, nuBPM;
3088	double betaHC, nuHC, betaVC, nuVC;
3089	FILE SkewMatFile, fp;
3090
3091	const int Xi = 1, Yi = 2;
3092	const double pi = M_PI, twopi = 2.0*M_PI;
3093
3094
3095	etaSQ = dvector(1, N_SKEW); betaSQ = dmatrix(1, N_SKEW, 1, 2);
3096	nuSQ = dmatrix(1, N_SKEW, 1, 2);
3097	betaBPM = dmatrix(1, N_BPM, 1, 2); nuBPM = dmatrix(1, N_BPM, 1, 2);
3098	betaHC = dmatrix(1, N_HCOR, 1, 2); nuHC = dmatrix(1, N_HCOR, 1, 2);
3099	betaVC = dmatrix(1, N_VCOR, 1, 2); nuVC = dmatrix(1, N_VCOR, 1, 2);
3100
3101	nuX = globval.TotalTune[X_]; nuY = globval.TotalTune[Y_];
3102
3103	for (i = 1; i <= N_SKEW; i++) {
3104	loc = Elem_GetPos(globval.qt, i);
3105	etaSQ[i] = Cell[loc].Eta[X_];
3106	betaSQ[i][Xi] = Cell[loc].Beta[X_]; betaSQ[i][Yi] = Cell[loc].Beta[Y_];
3107	nuSQ[i][Xi] = Cell[loc].Nu[X_]; nuSQ[i][Yi] = Cell[loc].Nu[Y_];
3108	} // for i=1..N_SKEW
3109
3110	for (i = 1; i <= N_BPM; i++) {
3111	betaBPM[i][Xi] = Cell[bpm_loc[i-1]].Beta[X_];
3112	betaBPM[i][Yi] = Cell[bpm_loc[i-1]].Beta[Y_];
3113	nuBPM[i][Xi] = Cell[bpm_loc[i-1]].Nu[X_];
3114	nuBPM[i][Yi] = Cell[bpm_loc[i-1]].Nu[Y_];
3115	} // for i=1..N_BPM
3116
3117	for (i = 1; i <= N_HCOR; i++) {
3118	betaHC[i][Xi] = Cell[h_corr[i-1]].Beta[X_];
3119	betaHC[i][Yi] = Cell[h_corr[i-1]].Beta[Y_];
3120	nuHC[i][Xi] = Cell[h_corr[i-1]].Nu[X_];
3121	nuHC[i][Yi] = Cell[h_corr[i-1]].Nu[Y_];
3122	} // for i=1..N_HCOR
3123
3124	for (i = 1; i <= N_VCOR; i++) {
3125	betaVC[i][Xi] = Cell[v_corr[i-1]].Beta[X_];
3126	betaVC[i][Yi] = Cell[v_corr[i-1]].Beta[Y_];
3127	nuVC[i][Xi] = Cell[v_corr[i-1]].Nu[X_];
3128	nuVC[i][Yi] = Cell[v_corr[i-1]].Nu[Y_];
3129	} // for i=1..N_VCOR
3130
3131
3132	for (i = 1; i <= N_SKEW; i++) {
3133	// looking for term for vertical dispersion
3134	alpha = etaSQ[i];
3135	// printf("For skew quad %3d kick is %9.2e\n",i,alpha);
3136	for (j = 1; j <= N_BPM; j++) {
3137	SkewRespMat[j][i] = VDweight0.5alphasqrt(betaSQ[i][Yi]betaBPM[j][Yi])
3138	cos(twopifabs(nuSQ[i][Yi]-nuBPM[j][Yi])-pinuY)/sin(pinuY);
3139	} // for (j=1; j<=N_BPM; j++)
3140
3141	// loking for coupling of horizontal trim to vertical BPM
3142	for (k = 1; k <= N_HCOR; k++) {
3143	// find v-kick by i-th skew quad due to the k-th h-trim
3144	alpha = 0.5sqrt(betaSQ[i][Xi]betaHC[k][Xi])*
3145	cos(twopifabs(nuSQ[i][Xi]-nuHC[k][Xi])-pinuX)/sin(pi*nuX);
3146	// find vertical orbit due to the kick
3147	for (j = 1; j <= N_BPM; j++)
3148	SkewRespMat[N_BPM+(k-1)*N_HCOR+j][i] =
3149	HVweight0.5alphasqrt(betaSQ[i][Yi]betaBPM[j][Yi])*
3150	cos(twopifabs(nuSQ[i][Yi]-nuBPM[j][Yi])-pinuY)/sin(pi*nuY);
3151	} //for (k=1; k<=N_HCOR; k++)
3152
3153	//loking for coupling of vertical trim to horizontal BPM
3154	for (k = 1; k <= N_VCOR; k++) {
3155	// find h-kick by i-th skew quad due to the k-th v-trim
3156	alpha = 0.5sqrt(betaSQ[i][Yi]betaVC[k][Yi])*
3157	cos(twopifabs(nuSQ[i][Yi]-nuVC[k][Yi])-pinuY)/sin(pi*nuY);
3158	// find horizontal orbit due to the kick
3159	for (j = 1; j <= N_BPM; j++)
3160	SkewRespMat[N_BPM+N_BPMN_HCOR+(k-1)N_VCOR+j][i] =
3161	VHweight0.5alphasqrt(betaSQ[i][Xi]betaBPM[j][Xi])*
3162	cos(twopifabs(nuSQ[i][Xi]-nuBPM[j][Xi])-pinuX)/sin(pi*nuX);
3163	} //for (k=1; k<=N_VCOR; k++)
3164	} // for i=1..N_SKEW
3165
3166
3167	SkewMatFile = file_write(SkewMatFileName);
3168	for (i = 1; i <= N_SKEW; i++) {
3169	for (j = 1; j <= N_COUPLE; j++)
3170	fprintf(SkewMatFile, "%9.2e ", SkewRespMat[j][i]);
3171	fprintf(SkewMatFile, "\n");
3172	}
3173	fclose(SkewMatFile);
3174
3175	fp = file_write(deta_y_FileName);
3176	eta_y_max = 0.0;
3177	for (j = 1; j <= N_BPM; j++) {
3178	eta_y[j] = 0.0;
3179	for (i = 1; i <= N_SKEW; i++)
3180	if (i % SQ_per_scell == 0) {
3181	eta_y[j] += 0.5etaSQ[i]sqrt(betaSQ[i][Yi]*betaBPM[j][Yi])
3182	cos(twopifabs(nuSQ[i][Yi]-nuBPM[j][Yi])-pinuY)/sin(pinuY);
3183	}
3184	eta_y_max = max(fabs(eta_y[j]), eta_y_max);
3185	}
3186	for (j = 1; j <= N_BPM; j++)
3187	eta_y[j] /= eta_y_max;
3188
3189	for (j = 1; j <= N_BPM; j++) {
3190	eta_y[j] *= deta_y_max;
3191	fprintf(fp, "%6.3f %10.3e\n", Cell[bpm_loc[j-1]].S, 1e3*eta_y[j]);
3192	}
3193	fclose(fp);
3194
3195	free_dvector(etaSQ, 1, N_SKEW); free_dmatrix(betaSQ, 1, N_SKEW, 1, 2);
3196	free_dmatrix(nuSQ, 1, N_SKEW, 1, 2);
3197	free_dmatrix(betaBPM, 1, N_BPM, 1, 2); free_dmatrix(nuBPM, 1, N_BPM, 1, 2);
3198	free_dmatrix(betaHC, 1, N_HCOR, 1, 2); free_dmatrix(nuHC, 1, N_HCOR, 1, 2);
3199	free_dmatrix(betaVC, 1, N_VCOR, 1, 2); free_dmatrix(nuVC, 1, N_VCOR, 1, 2);
3200	} // FindMatrix
3201
3202
3203	void ini_skew_cor(const double deta_y_max)
3204	{
3205	int k;
3206
3207	// No of skew quads, BPMs, and correctors
3208	N_SKEW = GetnKid(globval.qt);
3209
3210	N_BPM = 0;
3211	for (k = 1; k <= GetnKid(globval.bpm); k++) {
3212	N_BPM++;
3213
3214	if (N_BPM > max_bpm) {
3215	printf("ini_skew_cor: max no of BPMs exceeded %d (%d)\n",
3216	N_BPM, max_bpm);
3217	exit_(1);
3218	}
3219
3220	bpm_loc[N_BPM-1] = Elem_GetPos(globval.bpm, k);
3221	}
3222
3223	N_HCOR = 0;
3224	h_corr[N_HCOR++] = Elem_GetPos(globval.hcorr, 1);
3225	h_corr[N_HCOR++] = Elem_GetPos(globval.hcorr, GetnKid(globval.hcorr)/3);
3226	h_corr[N_HCOR++] = Elem_GetPos(globval.hcorr, 2*GetnKid(globval.hcorr)/3);
3227
3228	N_VCOR = 0;
3229	v_corr[N_VCOR++] = Elem_GetPos(globval.vcorr, 1);
3230	v_corr[N_VCOR++] = Elem_GetPos(globval.vcorr, GetnKid(globval.vcorr)/3);
3231	v_corr[N_VCOR++] = Elem_GetPos(globval.vcorr, 2*GetnKid(globval.vcorr)/3);
3232
3233	N_COUPLE = N_BPM*(1+N_HCOR+N_VCOR);
3234
3235	SkewRespMat = dmatrix(1, N_COUPLE, 1, N_SKEW);
3236	VertCouple = dvector(1, N_COUPLE);
3237	SkewStrengthCorr = dvector(1, N_SKEW);
3238	b = dvector(1, N_COUPLE); w = dvector(1, N_SKEW);
3239	V = dmatrix(1, N_SKEW, 1, N_SKEW); U = dmatrix(1, N_COUPLE, 1, N_SKEW);
3240	eta_y = dvector(1, N_BPM);
3241
3242	printf("\n");
3243	printf("Number of trims: horizontal = %d, vertical = %d\n",
3244	N_HCOR, N_VCOR);
3245	printf("Number of BPMs: %6d\n", N_BPM);
3246	printf("Number of skew quads: %6d\n", N_SKEW);
3247	printf("Number of elements in skew vector: %6d\n", N_COUPLE);
3248
3249	// find matrix
3250	Ring_GetTwiss(true, 0.0);
3251
3252	printf("\n");
3253	printf("Looking for response matrix\n");
3254	FindMatrix(SkewRespMat, deta_y_max);
3255
3256	printf("Looking for SVD matrices\n");
3257	FindSQ_SVDmat(SkewRespMat, U, V, w, N_COUPLE, N_SKEW);
3258	}
3259
3260
3261	void FindCoupVector(double *VertCouple)
3262	{
3263	bool cod;
3264	long i,j;
3265	long lastpos;
3266	double orbitP, orbitN;
3267
3268	orbitP = dvector(1, N_BPM); orbitN = dvector(1, N_BPM);
3269
3270	// Find vertical dispersion
3271	Cell_Geteta(0, globval.Cell_nLoc, true, 0e0);
3272
3273	for (i = 1; i <= N_BPM; i++)
3274	VertCouple[i] = VDweight*Cell[bpm_loc[i-1]].Eta[Y_];
3275	// Finished finding vertical dispersion
3276
3277	// Find off diagonal terms for horizontal trims
3278	for (j = 1; j <= N_HCOR; j++) {
3279	// positive kick: "+Dip" for horizontal
3280	SetdKLpar(Cell[h_corr[j-1]].Fnum, Cell[h_corr[j-1]].Knum, +Dip, kick);
3281	cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3282	for (i = 1; i <= N_BPM; i++)
3283	orbitP[i] = Cell[bpm_loc[i-1]].BeamPos[y_];
3284
3285	//negative kick: "+Dip" for horizontal
3286	SetdKLpar(Cell[h_corr[j-1]].Fnum, Cell[h_corr[j-1]].Knum, +Dip, -2*kick);
3287	cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3288	for (i = 1; i <= N_BPM; i++)
3289	orbitN[i] = Cell[bpm_loc[i-1]].BeamPos[y_];
3290
3291	// restore trim valueL: "+Dip" for horizontal
3292	SetdKLpar(Cell[h_corr[j-1]].Fnum, Cell[h_corr[j-1]].Knum, +Dip, kick);
3293
3294	for (i = 1; i <= N_BPM; i++)
3295	VertCouple[N_BPM+(j-1)*N_HCOR+i] =
3296	HVweight(orbitN[i]-orbitP[i])0.5/kick; // sign reversal
3297	} // hcorr cycle
3298
3299
3300	// Find off diagonal terms for vertical trims
3301	for (j = 1; j <= N_VCOR; j++){
3302	// positive kick: "-Dip" for vertical
3303	SetdKLpar(Cell[v_corr[j-1]].Fnum, Cell[v_corr[j-1]].Knum, -Dip, kick);
3304	cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3305	for (i = 1; i <= N_BPM; i++)
3306	orbitP[i] = Cell[bpm_loc[i-1]].BeamPos[x_];
3307
3308	// negative kick: "-Dip" for vertical
3309	SetdKLpar(Cell[v_corr[j-1]].Fnum, Cell[v_corr[j-1]].Knum, -Dip, -2*kick);
3310	cod = getcod(0.0, lastpos); chk_cod(cod, "FindCoupVector");
3311	for (i = 1; i <= N_BPM; i++)
3312	orbitN[i] = Cell[bpm_loc[i-1]].BeamPos[x_];
3313
3314	// restore corrector: "-Dip" for vertical
3315	SetdKLpar(Cell[v_corr[j-1]].Fnum, Cell[v_corr[j-1]].Knum, -Dip, kick);
3316
3317	for (i = 1; i <= N_BPM; i++)
3318	VertCouple[N_BPM+N_BPMN_HCOR+(j-1)N_VCOR+i] =
3319	VHweight(orbitP[i]-orbitN[i])0.5/kick;
3320	} // vcorr cycle
3321
3322	free_dvector(orbitP, 1, N_BPM); free_dvector(orbitN, 1, N_BPM);
3323	} // FindCoupVector
3324
3325
3326	void SkewStat(double VertCouple[])
3327	{
3328	int i;
3329	double max, rms, sk;
3330
3331	// statistics for skew quadrupoles
3332	max = 0.0; rms = 0.0;
3333	for(i = 1; i <= N_SKEW; i++) {
3334	sk = GetKLpar(globval.qt, i, -Quad);
3335	if (fabs(sk) > max) max = fabs(sk);
3336	rms += sqr(sk);
3337	}
3338	rms = sqrt(rms/N_SKEW);
3339	printf("Rms skew strength: %8.2e+/-%8.2e\n", max, rms);
3340
3341	// statistics for vertical dispersion function
3342	max = 0.0; rms = 0.0;
3343	for(i = 1; i <= N_BPM; i++) {
3344	if (fabs(VertCouple[i]) > max) max = fabs(VertCouple[i]/VDweight);
3345	rms += sqr(VertCouple[i]/VDweight);
3346	}
3347	rms = sqrt(rms/N_BPM);
3348	printf("Max vertical dispersion: %8.2e+/-%8.2e mm\n", 1e3max, 1e3rms);
3349
3350	// statistics for off diagonal terms of response matrix (trims->bpms)
3351	max = 0.0; rms = 0.0;
3352	for(i = N_BPM+1; i <= N_BPM*(1+N_HCOR); i++) {
3353	if (fabs(VertCouple[i]) > max) max = fabs(VertCouple[i]/HVweight);
3354	rms += sqr(VertCouple[i]/HVweight);
3355	}
3356	rms = sqrt(rms/(N_HCOR*N_BPM));
3357	printf("Max horizontal coupling: %8.2e+/-%8.2e mm/mrad\n", max, rms);
3358
3359	max = 0.0; rms = 0.0;
3360	for(i = N_BPM*(1+N_HCOR)+1; i <= N_COUPLE; i++) {
3361	if (fabs(VertCouple[i]) > max) max = fabs(VertCouple[i]/VHweight);
3362	rms += sqr(VertCouple[i]/VHweight);
3363	}
3364	rms = sqrt(rms/(N_VCOR*N_BPM));
3365	printf("Max vertical coupling: %8.2e+/-%8.2e mm/mrad\n", max, rms);
3366	}
3367
3368
3369	void corr_eps_y(void)
3370	{
3371	int i, j;
3372	FILE *outf;
3373
3374	// Clear skew quad setpoints
3375	set_bnL_design_fam(globval.qt, Quad, 0.0, 0.0);
3376
3377	// Find coupling vector
3378	printf("\n");
3379	printf("Looking for coupling error\n");
3380	FindCoupVector(VertCouple);
3381
3382	//Find and print coupling statistics
3383	printf("\n");
3384	printf("Before correction\n");
3385	SkewStat(VertCouple);
3386
3387	// Coupling Correction
3388	printf("\n");
3389	for (i = 1; i <= n_lin; i++) {
3390	printf("Looking for correction\n");
3391
3392	//Find Correcting Settings to skew quadrupoles
3393	for (j = 1; j <= N_BPM; j++)
3394	b[j] = VDweight*eta_y[j] - VertCouple[j];
3395
3396	for (j = N_BPM+1; j <= N_COUPLE; j++)
3397	b[j] = -VertCouple[j];
3398
3399	dsvbksb(U, w, V, N_COUPLE, N_SKEW, b, SkewStrengthCorr);
3400
3401	printf("Applying correction\n");
3402	// Add correction
3403	for (j = 1; j <= N_SKEW; j++)
3404	SetdKLpar(globval.qt, j, -Quad, SkewStrengthCorr[j]);
3405
3406	printf("\n");
3407	printf("Looking for coupling error\n");
3408	// Find coupling vector
3409	FindCoupVector(VertCouple);
3410
3411	printf("\n");
3412	printf("After run %d of correction\n", i);
3413	// Find and print coupling statistics
3414	SkewStat(VertCouple);
3415	} // End of coupling Correction
3416
3417	outf = file_write(eta_y_FileName);
3418	for (i = 0; i <= globval.Cell_nLoc; i++)
3419	fprintf(outf, "%4d %7.3f %s %6.3f %10.3e %10.3e\n",
3420	i, Cell[i].S, Cell[i].Elem.PName,
3421	Cell[i].Nu[Y_], 1e3Cell[i].Eta[Y_], 1e3Cell[i].Etap[Y_]);
3422	fclose(outf);
3423
3424	FindCoupVector(VertCouple);
3425	}
3426
3427
3428	/*******************************************************************
3429	void get_IDs(void)
3430
3431	Purpose:
3432	Get the ID.
3433
3434	Input:
3435	none
3436
3437	Output:
3438	none
3439
3440	Comments:
3441	Print out the scaling 2 of Insertion.
3442
3443	********************************************************************/
3444	// Control of IDs
3445
3446	void get_IDs(void)
3447	{
3448	int k;
3449
3450	printf("\n");
3451	n_ID_Fams = 0;
3452	for (k = 0; k < globval.Elem_nFam; k++)
3453	switch (ElemFam[k].ElemF.Pkind) {
3454	case Wigl:
3455	printf("found ID family: %s %12.5e\n",
3456	ElemFam[k].ElemF.PName, ElemFam[k].ElemF.W->BoBrhoV[0]);
3457	n_ID_Fams++; ID_Fams[n_ID_Fams-1] = k + 1;
3458	break;
3459	case Insertion:
3460	printf("found ID family: %s scaling1 = %12.5e, scaling2 = %12.5e\n",
3461	ElemFam[k].ElemF.PName, ElemFam[k].ElemF.ID->scaling1, ElemFam[k].ElemF.ID->scaling2);
3462	n_ID_Fams++; ID_Fams[n_ID_Fams-1] = k + 1;
3463	break;
3464	case FieldMap:
3465	printf("found ID family: %s %12.5e\n",
3466	ElemFam[k].ElemF.PName, ElemFam[k].ElemF.FM->scl);
3467	n_ID_Fams++; ID_Fams[n_ID_Fams-1] = k + 1;
3468	break;
3469	default:
3470	break;
3471	}
3472	}
3473
3474	/****************************************************************************/
3475	/* void set_IDs(const double scl)
3476
3477	Purpose:
3478	Set the scale factor for all the wigger, insertion devices, or fieldMap.
3479
3480	Input:
3481	scl scaling factor
3482
3483	Output:
3484	none
3485
3486	Return:
3487	none
3488
3489	Global variables:
3490	none
3491
3492	Specific functions:
3493	none
3494
3495	Comments:
3496	See also:
3497	set_ID_scl(const int Fnum, const double scl)
3498	set_ID_scl(const int Fnum, const int knum, const double scl)
3499
3500	****************************************************************************/
3501	void set_IDs(const double scl)
3502	{
3503	int k;
3504
3505	printf("\n");
3506	for (k = 0; k < n_ID_Fams; k++) {
3507	switch (ElemFam[ID_Fams[k]-1].ElemF.Pkind) {
3508	case Wigl:
3509	printf("setting ID family: %s %12.5e\n",
3510	ElemFam[ID_Fams[k]-1].ElemF.PName,
3511	scl*ElemFam[ID_Fams[k]-1].ElemF.W->BoBrhoV[0]);
3512
3513	set_Wiggler_BoBrho(ID_Fams[k],
3514	scl*ElemFam[ID_Fams[k]-1].ElemF.W->BoBrhoV[0]);
3515	break;
3516	case Insertion:
3517	printf("setting ID family: %s %12.5e\n",
3518	ElemFam[ID_Fams[k]-1].ElemF.PName, scl);
3519
3520	set_ID_scl(ID_Fams[k], scl);
3521	break;
3522	case FieldMap:
3523	printf("setting ID family: %s %12.5e\n",
3524	ElemFam[ID_Fams[k]-1].ElemF.PName, scl);
3525
3526	set_ID_scl(ID_Fams[k], scl);
3527	break;
3528	default:
3529	cout << "set_IDs: unknown element type" << endl;
3530	exit_(1);
3531	break;
3532	}
3533	}
3534	}
3535
3536
3537	void reset_quads(void)
3538	{
3539	int k;
3540
3541	if (N_Fam > N_Fam_max) {
3542	printf("reset_quads: N_Fam > N_Fam_max: %d (%d)\n", N_Fam, N_Fam_max);
3543	exit_(0);
3544	}
3545
3546	for (k = 0; k < N_Fam; k++) {
3547	// Note, actual values can differ from the original values
3548	/* printf("setting quad family: %s %12.5e\n",
3549	ElemFam[Q_Fam[k]-1].ElemF.PName,
3550	ElemFam[Q_Fam[k]-1].ElemF.M->PBpar[HOMmax+Quad]);
3551
3552	set_bn_design_fam(Q_Fam[k], Quad,
3553	ElemFam[Q_Fam[k]-1].ElemF.M->PBpar[HOMmax+Quad], 0.0);*/
3554
3555	printf("setting quad family: %s %12.5e\n",
3556	ElemFam[Q_Fam[k]-1].ElemF.PName, b2[k]);
3557
3558	set_bn_design_fam(Q_Fam[k], Quad, b2[k], 0.0);
3559	}
3560	}
3561
3562
3563	void SVD(const int m, const int n, double **M, double beta_nu[],
3564	double b2Ls_[], const double s_cut, const bool first)
3565	{
3566	int i, j;
3567
3568	const bool prt = true;
3569
3570	if (first) {
3571	for (i = 1; i <= m; i++)
3572	for (j = 1; j <= n; j++)
3573	U1[i][j] = M[i][j];
3574
3575	dsvdcmp(U1, m, n, w1, V1);
3576
3577	if (prt) {
3578	printf("\n");
3579	printf("singular values:\n");
3580	printf("\n");
3581	}
3582
3583	for (i = 1; i <= n; i++) {
3584	if (prt) printf("%11.3e", w1[i]);
3585	if (w1[i] < s_cut) {
3586	w1[i] = 0.0;
3587	if (prt) printf(" (zeroed)");
3588	}
3589	if (prt) if (i % 8 == 0) printf("\n");
3590	}
3591	if (prt) if (n % 8 != 0) printf("\n");
3592	}
3593
3594	dsvbksb(U1, w1, V1, m, n, beta_nu, b2Ls_);
3595	}
3596
3597
3598	void quad_config()
3599	{
3600	int i, j;
3601	double an;
3602
3603	if (N_Fam > N_Fam_max) {
3604	printf("quad_config: N_Fam > N_Fam_max: %d (%d)\n", N_Fam, N_Fam_max);
3605	exit_(0);
3606	}
3607
3608	Nquad = 0;
3609	for (i = 0; i < N_Fam; i++) {
3610	for (j = 1; j <= GetnKid(Q_Fam[i]); j++) {
3611	Nquad++;
3612
3613	if (Nquad > n_b2_max) {
3614	printf("quad_config: max no of quadrupoles exceeded %d (%d)\n",
3615	Nquad, n_b2_max);
3616	exit_(1);
3617	}
3618
3619	quad_prms[Nquad-1] = Elem_GetPos(Q_Fam[i], j);
3620
3621	if (j == 1) get_bn_design_elem(Q_Fam[i], j, Quad, b2[i], an);
3622	}
3623	}
3624
3625	printf("\n");
3626	printf("quad_config: Nquad = %d\n", Nquad);
3627	}
3628
3629
3630	bool get_SQ(void)
3631	{
3632	int j, k;
3633	Vector2 alpha3[3], beta3[3], nu3[3], eta3[3], etap3[3];
3634	FILE *outf;
3635
3636	const bool prt = false;
3637
3638	/* Note, IDs are split for evaluation of the driving terms at the center:
3639	id1 1, 2
3640	id2 1, 2
3641	... */
3642
3643	// Get Twiss params, no dispersion
3644	Ring_GetTwiss(false, 0.0);
3645
3646	if (!status.codflag \|\| !globval.stable) return false;
3647
3648	// Get global tunes
3649	Nu_X = globval.TotalTune[X_]; Nu_Y = globval.TotalTune[Y_];
3650
3651	if (prt) {
3652	printf("\n");
3653	printf("nu_x = %8.12f, nu_y = %8.12f\n", Nu_X, Nu_Y);
3654
3655	// Get Twiss params in sext
3656	printf("\n");
3657	printf("listing lattice functions at sextupoles\n");
3658	printf("\n");
3659
3660	outf = file_write("latfunS.out");
3661
3662	fprintf(outf, "s betax nux betay nuy\n");
3663	}
3664
3665	printf("\n");
3666	Nsext = 0;
3667	for (k = 0; k < globval.Cell_nLoc; k++) {
3668	if ((Cell[k].Elem.Pkind == Mpole) && (Cell[k].Elem.M->n_design == Sext)) {
3669	Nsext++;
3670
3671	if (Nsext > n_b3_max) {
3672	printf("get_SQ: max no of sextupoles exceeded %d (%d)\n",
3673	Nsext, n_b3_max);
3674	exit_(1);
3675	}
3676
3677	Ss[Nsext-1] = Cell[k].S;
3678
3679	for (j = 0; j <= 1; j++) {
3680	sb[j][Nsext-1] = Cell[k].Beta[j];
3681	sNu[j][Nsext-1] = Cell[k].Nu[j] - nu_0[j];
3682	}
3683
3684	if (prt) {
3685	printf("%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3686	Cell[k].Elem.PName, Ss[Nsext-1],
3687	sb[X_][Nsext-1], sNu[X_][Nsext-1]-nu_0[X_],
3688	sb[Y_][Nsext-1], sNu[Y_][Nsext-1]-nu_0[Y_]);
3689	fprintf(outf, "%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3690	Cell[k].Elem.PName, Ss[Nsext-1],
3691	sb[X_][Nsext-1], sNu[X_][Nsext-1]-nu_0[X_],
3692	sb[Y_][Nsext-1], sNu[Y_][Nsext-1]-nu_0[Y_]);
3693	}
3694	}
3695	}
3696
3697	if (prt) fclose(outf);
3698
3699	// Number of sexts in the ring
3700	printf("No of sextupoles = %d\n", Nsext);
3701
3702	if (prt) {
3703	// Get Twiss params in quads
3704	printf("\n");
3705	printf("listing lattice functions at quadrupoles\n");
3706	printf("\n");
3707
3708	outf = file_write("latfunQ.out");
3709
3710	fprintf(outf, "s name betax nux betay nuy\n");
3711	}
3712
3713	for (k = 0; k < Nquad; k++) {
3714	Sq[k] = Cell[quad_prms[k]].S;
3715	for (j = 0; j <= 1; j++) {
3716	if (Cell[quad_prms[k]].Elem.M->Pthick == thick) {
3717	get_twiss3(quad_prms[k], alpha3, beta3, nu3, eta3, etap3);
3718	qb[j][k] = beta3[Y_][j]; qNu[j][k] = nu3[Y_][j] - nu_0[j];
3719	} else {
3720	qb[j][k] = Cell[quad_prms[k]].Beta[j];
3721	qNu[j][k] = Cell[quad_prms[k]].Nu[j] - nu_0[j];
3722	}
3723	}
3724
3725	if (prt) {
3726	printf("%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3727	Cell[quad_prms[k]].Elem.PName, Sq[k], qb[X_][k],
3728	qNu[X_][k], qb[Y_][k], qNu[Y_][k]);
3729
3730	fprintf(outf, "%s %6.3f %8.5f %8.5f %8.5f %8.5f\n",
3731	Cell[quad_prms[k]].Elem.PName, Sq[k], qb[X_][k],
3732	qNu[X_][k], qb[Y_][k], qNu[Y_][k]);
3733	}
3734	}
3735
3736	if (prt) fclose(outf);
3737
3738	// Number of quads in the ring
3739	printf("No of quads = %d\n", Nquad);
3740
3741	return true;
3742	}
3743
3744
3745	double Bet(double bq, double nus, double nuq, double NuQ)
3746	{
3747	return bqcos(2.0M_PI(2.0fabs(nus-nuq)-NuQ))/(2.0sin(2.0M_PI*NuQ));
3748	}
3749
3750
3751	double Nus(double bq, double nus, double nuq, double NuQ)
3752	{
3753	double Nu, sgn;
3754
3755	sgn = ((nus-nuq) <= 0)? -1: 1;
3756
3757	Nu = -bqsgn(sin(2.0M_PINuQ)+sin(2.0M_PI(2.0*fabs(nus-nuq)-NuQ)))
3758	/(8.0M_PIsin(2.0M_PINuQ));
3759
3760	return Nu;
3761	}
3762
3763	/*****************************************************************
3764	void A_matrix(void)
3765
3766	Purpose:
3767	solove A in the matrix b = A.x
3768	using SVD.
3769	Used for ID correction.
3770
3771	Comments:
3772	change scl_dbeta to scl_dbetax and scl_dbetay
3773	change scl_dnu to scl_dnux and scl_dnuy
3774	change scl_nu to scl_nux and scl_nuy
3775
3776
3777	See also X_vector ()
3778	******************************************************************/
3779	void A_matrix(void)
3780	{
3781	int k, j;
3782	double BtX, BtY, NuX, NuY;
3783
3784	const bool prt = false;
3785
3786	// Defining Twiss in undisturbed quads
3787	for (k = 0; k < Nquad; k++)
3788	for (j = 0; j <= 1; j++) {
3789	qb0[j][k] = qb[j][k];
3790	qNu0[j][k] = qNu[j][k];
3791	}
3792
3793	// Defining Twiss in undisturbed sexts
3794	for (k = 0; k < Nsext; k++)
3795	for (j = 0; j <= 1; j++)
3796	sNu0[j][k] = sNu[j][k];
3797
3798	// Now creating matrix A in X=A*B2L
3799	for (k = 1; k <= Nsext; k++) {
3800	for (j = 1; j <= Nquad; j++) {
3801	BtX = Bet(qb0[X_][j-1], sNu0[X_][k-1], qNu0[X_][j-1], Nu_X0);
3802	NuX = -Nus(qb0[X_][j-1], sNu0[X_][k-1], qNu0[X_][j-1], Nu_X0);
3803	BtY = -Bet(qb0[Y_][j-1], sNu0[Y_][k-1], qNu0[Y_][j-1], Nu_Y0);
3804	NuY = Nus(qb0[Y_][j-1], sNu0[Y_][k-1], qNu0[Y_][j-1], Nu_Y0);
3805	A1[k][j] = scl_dbetax*BtX;
3806	A1[k+Nsext][j] = scl_dbetay*BtY;
3807	A1[k+2Nsext][j] = scl_dnuxNuX;
3808	A1[k+3Nsext][j] = scl_dnuyNuY;
3809	}
3810	}
3811	// Now adding 2 more constraints for global tunes
3812	for (j = 1; j <= Nquad; j++) {
3813	A1[4Nsext+1][j] = -scl_nuxqb0[X_][j-1]/(4.0*M_PI);
3814	A1[4Nsext+2][j] = scl_nuyqb0[Y_][j-1]/(4.0*M_PI);
3815	}
3816
3817	if (prt) {
3818	printf("\n");
3819	printf("AA:\n");
3820	printf("\n");
3821	for (k = 1; k <= Nconstr; k++) {
3822	for (j = 1; j <= Nquad; j++)
3823	printf(" %10.3e", A1[k][j]);
3824	printf("\n");
3825	}
3826	}
3827	}
3828
3829	/****************************************************
3830	void X_vector(const bool first)
3831
3832	Purpose:
3833	solove x in the matrix b = A.x
3834	using SVD.
3835	Used for ID correction.
3836
3837	Comments:
3838	change scl_dbeta to scl_dbetax and scl_dbetay
3839	change scl_dnu to scl_dnux and scl_dnuy
3840	change scl_nu to scl_nux and scl_nuy
3841
3842	See also:
3843	void A_matrix(void)
3844	*****************************************************/
3845	void X_vector(const bool first)
3846	{
3847	int k;
3848
3849	const bool prt = false;
3850
3851	dnu0[X_] = globval.TotalTune[X_] - Nu_X0;
3852	dnu0[Y_] = globval.TotalTune[Y_] - Nu_Y0;
3853
3854	if (first) {
3855	// Initial fill of X
3856	for (k = 1; k <= Nsext; k++) {
3857	Xsext0[k] = sb[X_][k-1]; Xsext0[k+Nsext] = sb[Y_][k-1];
3858	Xsext0[k+2Nsext] = sNu[X_][k-1]; Xsext0[k+3Nsext] = sNu[Y_][k-1];
3859	}
3860	Xsext0[4Nsext+1] = 0.0; Xsext0[4Nsext+2] = 0.0;
3861	} else {
3862	// Now substracting from X in X=A*B2L
3863	for (k = 1; k <= Nsext; k++) {
3864	Xsext[k] = scl_dbetax*(Xsext0[k]-sb[X_][k-1])/sb[X_][k-1];
3865	Xsext[k+Nsext] = scl_dbetay*(Xsext0[k+Nsext]-sb[Y_][k-1])/sb[Y_][k-1];
3866	Xsext[k+2Nsext] = scl_dnux(Xsext0[k+2*Nsext]-sNu[X_][k-1]+dnu0[X_]/2.0);
3867	Xsext[k+3Nsext] = scl_dnuy(Xsext0[k+3*Nsext]-sNu[Y_][k-1]+dnu0[Y_]/2.0);
3868	}
3869	Xsext[4Nsext+1] = scl_nux(Nu_X0-globval.TotalTune[X_]);
3870	Xsext[4Nsext+2] = scl_nuy(Nu_Y0-globval.TotalTune[Y_]);
3871	}
3872
3873	if (prt) {
3874	printf("\n");
3875	printf("X:\n");
3876	printf("\n");
3877	if (first) {
3878	for (k = 1; k <= Nconstr; k++) {
3879	printf(" %10.3e", Xsext0[k]);
3880	if (k % 10 == 0) printf("\n");
3881	}
3882	if (Nconstr % 10 != 0) printf("\n");
3883	} else {
3884	for (k = 1; k <= Nconstr; k++) {
3885	printf(" %10.3e", Xsext[k]);
3886	if (k % 10 == 0) printf("\n");
3887	}
3888	if (Nconstr % 10 != 0) printf("\n");
3889	}
3890	}
3891	}
3892
3893
3894	/*
3895	void ini_ID_corr(void)
3896
3897	Purpose:
3898	Initial ID, used for insertion Device correction
3899
3900
3901	*/
3902	void ini_ID_corr(void)
3903	{
3904
3905	// store ID families
3906	get_IDs();
3907
3908	// zero ID's
3909	set_IDs(0.0);
3910
3911	// Configuring quads (1 --> C means thin quads located in the middle of 1s)
3912	quad_config();
3913
3914	// Configuring quads (1 --> C means thin quads located in the middle of 1s)
3915	// Read Betas and Nus
3916	get_SQ();
3917	Nconstr = 4*Nsext + 2;
3918
3919	// Note, allocated vectors and matrices are deallocated in ID_corr
3920	Xsext = dvector(1, Nconstr); Xsext0 = dvector(1, Nconstr);
3921	b2Ls_ = dvector(1, Nquad); A1 = dmatrix(1, Nconstr, 1, Nquad);
3922	U1 = dmatrix(1, Nconstr, 1, Nquad); w1 = dvector(1, Nquad);
3923	V1 = dmatrix(1, Nquad, 1, Nquad);
3924
3925	// shift zero point to center of ID
3926	// nu_0[X_] = Cell[id_loc].Nu[X_]; nu_0[Y_] = Cell[id_loc].Nu[Y_];
3927	nu_0[X_] = 0.0; nu_0[Y_] = 0.0;
3928
3929	// Defining undisturbed tunes
3930	Nu_X0 = globval.TotalTune[X_]; Nu_Y0 = globval.TotalTune[Y_];
3931
3932	// Set-up matrix A in X=A*b2Ls_
3933	A_matrix();
3934
3935	// Now fill the X in X=A*b2Ls_
3936	X_vector(true);
3937	}
3938
3939
3940	void W_diag(void)
3941	{
3942	double bxf, byf, nxf, nyf, b2Lsum;
3943	int k;
3944
3945	bxf = 0.0; byf = 0.0; nxf = 0.0; nyf = 0.0;
3946	for (k = 1; k <= Nsext; k++) {
3947	bxf += sqr(Xsext[k]);
3948	byf += sqr(Xsext[k+Nsext]);
3949	nxf += sqr(Xsext[k+2*Nsext]);
3950	nyf += sqr(Xsext[k+3*Nsext]);
3951	}
3952
3953	dnu0[X_] = globval.TotalTune[X_] - Nu_X0;
3954	dnu0[Y_] = globval.TotalTune[Y_] - Nu_Y0;
3955
3956	b2Lsum = 0.0;
3957	for (k = 1; k <= Nquad; k++)
3958	b2Lsum += sqr(b2Ls_[k]);
3959
3960	printf("\n");
3961	printf("Residuals: beta [%%], dnu : \n");
3962	printf("dbeta_x: %6.2f dbeta_y: %6.2f nu_x: %12.6e nu_y: %12.6e\n",
3963	sqrt(bxf)/Nsext1e2, sqrt(byf)/Nsext1e2,
3964	sqrt(nxf)/Nsext, sqrt(nyf)/Nsext);
3965	printf("tune shift: dnu_x = %8.5f, dnu_y = %8.5f\n", dnu0[X_], dnu0[Y_]);
3966	printf("Sum b2Ls_: %12.6e\n", sqrt(b2Lsum)/Nquad);
3967	}
3968
3969
3970
3971
3972
3973
3974	/*
3975	For ID correction, from Mao-Sen Qiu at Taiwan light source.
3976	*/
3977
3978	bool ID_corr0(void)
3979	{
3980	int i, j, k;
3981	double b2, a2;
3982	FILE *outf;
3983
3984	printf("\n");
3985	printf("Before correction!\n");
3986
3987	set_IDs(1.0);
3988	//set_IDs((double) i);
3989
3990	get_SQ(); // Read Betas and Nus
3991	X_vector(false); // Fill in dX in dX=A*db2Ls_
3992	W_diag(); // Get statistics
3993
3994	outf = file_write("ID_before_corr.out");
3995	fprintf(outf, "# dbeta_x/beta_x dbeta_y/beta_y dnu_x dnu_y\n");
3996	fprintf(outf, "# [%%] [%%]\n");
3997	fprintf(outf, "#\n");
3998	for (k = 1; k <= Nsext; k++)
3999	fprintf(outf, "%6.1f %6.2f %6.2f %10.3e %10.3e\n", Ss[k-1], 1e2Xsext[k], 1e2Xsext[k+Nsext], Xsext[k+2Nsext], Xsext[k+3Nsext]);
4000	fclose(outf);
4001
4002	// Allow for repeated calls to ID_corr, allocation is done in ini_ID_corr.
4003	if (true) {
4004	free_dvector(Xsext, 1, Nconstr);
4005	free_dvector(Xsext0, 1, Nconstr);
4006	free_dvector(b2Ls_, 1, Nquad);
4007	free_dmatrix(A1, 1, Nconstr, 1, Nquad);
4008	free_dmatrix(U1, 1, Nconstr, 1, Nquad);
4009	free_dvector(w1, 1, Nquad);
4010	free_dmatrix(V1, 1, Nquad, 1, Nquad);
4011	}
4012
4013	printf("\n");
4014	printf("End of before correction!\n");
4015
4016	return true;
4017	}
4018
4019
4020	/******************************************************************
4021	bool ID_corr(const int N_calls, const int N_steps)
4022
4023	Purpose:
4024	Compensate the field error introduced by insertion device.
4025
4026	******************************************************************/
4027	bool ID_corr(const int N_calls, const int N_steps)
4028	{
4029	int i, j, k;
4030	double b2, a2;
4031	FILE *outf;
4032
4033	printf("\n");
4034	printf("ID matching begins!\n");
4035
4036	outf = file_write("ID_corr.out");
4037	for (i = 1; i <= N_steps; i++) { //This brings ID strength in steps
4038	set_IDs((double)i/(double)N_steps);
4039
4040	get_SQ(); // Read Betas and Nus
4041	X_vector(false); // Fill in dX in dX=A*db2Ls_
4042	W_diag(); // Get statistics
4043	for (j = 1; j <= N_calls; j++) {
4044	SVD(Nconstr, Nquad, A1, Xsext, b2Ls_, 1e0, j == 1);
4045
4046	if ((i == N_steps) && (j == N_calls)) {
4047	fprintf(outf, "b_2:\n");
4048	fprintf(outf, "\n");
4049	}
4050
4051	// add quad strengths (db2Ls_)
4052	for (k = 1; k <= Nquad; k++) {
4053	set_dbnL_design_elem(Cell[quad_prms[k-1]].Fnum,
4054	Cell[quad_prms[k-1]].Knum, Quad,
4055	-ID_step*b2Ls_[k], 0.0);
4056
4057	if ((i == N_steps) && (j == N_calls)) {
4058	get_bn_design_elem(Cell[quad_prms[k-1]].Fnum,
4059	Cell[quad_prms[k-1]].Knum, Quad, b2, a2);
4060	fprintf(outf, "%10s %2d %8.5f\n",
4061	Cell[quad_prms[k-1]].Elem.PName, k, b2);
4062	}
4063	}
4064
4065	printf("\n");
4066	printf("Iteration: %2d\n", j);
4067	if (get_SQ()) {
4068	X_vector(false); // Fill in dX in dX=A*db2Ls_
4069	W_diag(); // Get statistics
4070
4071	printglob();
4072	} else {
4073	printf("ID_corr: correction failed\n");
4074	// restore lattice
4075	set_IDs(0.0); reset_quads();
4076	return false;
4077	}
4078	}
4079	}
4080	fclose(outf);
4081
4082	outf = file_write("ID_corr_res.out");
4083	fprintf(outf, "# dbeta_x/beta_x dbeta_y/beta_y dnu_x dnu_y\n");
4084	fprintf(outf, "# [%%] [%%]\n");
4085	fprintf(outf, "#\n");
4086	for (k = 1; k <= Nsext; k++)
4087	fprintf(outf, "%6.1f %6.2f %6.2f %10.3e %10.3e\n",
4088	Ss[k-1], 1e2Xsext[k], 1e2Xsext[k+Nsext],
4089	Xsext[k+2Nsext], Xsext[k+3Nsext]);
4090	fclose(outf);
4091
4092	// Allow for repeated calls to ID_corr, allocation is done in ini_ID_corr.
4093	if (false) {
4094	free_dvector(Xsext, 1, Nconstr); free_dvector(Xsext0, 1, Nconstr);
4095	free_dvector(b2Ls_, 1, Nquad); free_dmatrix(A1, 1, Nconstr, 1, Nquad);
4096	free_dmatrix(U1, 1, Nconstr, 1, Nquad); free_dvector(w1, 1, Nquad);
4097	free_dmatrix(V1, 1, Nquad, 1, Nquad);
4098	}
4099
4100	printf("\n");
4101	printf("ID matching ends!\n");
4102
4103	return true;
4104	}
4105
4106	/* End ID correction functions *****/
4107
4108
4109	void get_param(const char *param_file)
4110	/*
4111	* Read parameter file in input file xxxx.prm
4112	*
4113	*/
4114	{
4115	char line[max_str], name[max_str], str[max_str], s_prm[max_str];
4116	char lat_file[max_str], flat_file[max_str], IDCq_name[max_str][8];
4117	int k;
4118	double f_prm;
4119	FILE *inf;
4120
4121	const bool prt = true;
4122
4123	if (prt) {
4124	printf("\n");
4125	printf("reading in %s\n", param_file);
4126	}
4127
4128	inf = file_read(param_file);
4129
4130	// read parameter file
4131	// initialization
4132	strcpy(ae_file, ""); strcpy(fe_file, ""); strcpy(ap_file, "");
4133
4134	if (prt) printf("\n");
4135	// loop over all lines of the file
4136	while (fgets(line, max_str, inf) != NULL) {
4137	if (prt) printf("%s", line);
4138	if (strstr(line, "#") == NULL) {
4139	// get initial command token
4140	sscanf(line, "%s", name);
4141	// input files *************************************
4142	if (strcmp("in_dir", name) == 0)
4143	sscanf(line, "%*s %s", in_dir);
4144	else if (strcmp("ae_file", name) == 0){
4145	sscanf(line, "%*s %s", str);
4146	sprintf(ae_file,"%s/%s", in_dir, str);
4147	} else if (strcmp("fe_file", name) == 0) {
4148	sscanf(line, "%*s %s", str);
4149	sprintf(fe_file, "%s/%s", in_dir, str);
4150	} else if (strcmp("ap_file", name) == 0) {
4151	sscanf(line, "%*s %s", str);
4152	sprintf(ap_file, "%s/%s", in_dir, str);
4153	} else if (strcmp("lat_file", name) == 0) {
4154	sscanf(line, "%*s %s", lat_file);
4155	sprintf(lat_FileName, "%s/%s", in_dir, lat_file);
4156	Read_Lattice(lat_FileName);
4157	} else if (strcmp("flat_file", name) == 0) {
4158	sscanf(line, "%*s %s", flat_file);
4159	strcat(flat_file, "flat_file.dat"); rdmfile(flat_file);
4160	// ********< parameters >**************************
4161	} else if (strcmp("s_cut", name) == 0) {
4162	sscanf(line, "%*s %lf", &f_prm);
4163	setrancut(f_prm);
4164	} else if (strcmp("n_stat", name) == 0)
4165	sscanf(line, "%*s %d", &n_stat);
4166	else if (strcmp("n_scale", name) == 0)
4167	sscanf(line, "%*s %d", &n_scale);
4168	else if (strcmp("n_orbit", name) == 0)
4169	sscanf(line, "%*s %d", &n_orbit);
4170	else if (strcmp("bpm_name", name) == 0) {
4171	sscanf(line, "%*s %s", s_prm);
4172	globval.bpm = ElemIndex(s_prm);
4173	} else if (strcmp("h_corr", name) == 0) {
4174	sscanf(line, "%*s %s", s_prm);
4175	globval.hcorr = ElemIndex(s_prm);
4176	} else if (strcmp("v_corr", name) == 0) {
4177	sscanf(line, "%*s %s", s_prm);
4178	globval.vcorr = ElemIndex(s_prm);
4179	} else if (strcmp("gs", name) == 0) {
4180	sscanf(line, "%*s %s", s_prm);
4181	globval.gs = ElemIndex(s_prm);
4182	} else if (strcmp("ge", name) == 0) {
4183	sscanf(line, "%*s %s", s_prm);
4184	globval.ge = ElemIndex(s_prm);
4185	} else if (strcmp("DA_bare", name) == 0) {
4186	sscanf(line, "%*s %s", s_prm);
4187	DA_bare = (strcmp(s_prm, "true") == 0)? true : false;
4188	} else if (strcmp("delta", name) == 0) {
4189	sscanf(line, "%*s %lf", &delta_DA_);
4190	} else if (strcmp("freq_map", name) == 0) {
4191	sscanf(line, "%*s %s %d %d %d %d %lf %lf %lf",
4192	s_prm, &n_x, &n_y, &n_dp, &n_tr,
4193	&x_max_FMA, &y_max_FMA, &delta_FMA);
4194	freq_map = (strcmp(s_prm, "true") == 0)? true : false;
4195	} else if (strcmp("tune_shift", name) == 0) {
4196	sscanf(line, "%*s %s", s_prm);
4197	tune_shift = (strcmp(s_prm, "true") == 0)? true : false;
4198	} else if (strcmp("qt", name) == 0) {
4199	sscanf(line, "%*s %s", s_prm);
4200	globval.qt = ElemIndex(s_prm);
4201	} else if (strcmp("disp_wave_y", name) == 0)
4202	sscanf(line, "%*s %lf", &disp_wave_y);
4203	else if (strcmp("n_lin", name) == 0)
4204	sscanf(line, "%*s %d", &n_lin);
4205	else if (strcmp("VDweight", name) == 0)
4206	sscanf(line, "%*s %lf", &VDweight);
4207	else if (strcmp("HVweight", name) == 0)
4208	sscanf(line, "%*s %lf", &HVweight);
4209	else if (strcmp("VHweight", name) == 0)
4210	sscanf(line, "%*s %lf", &VHweight);
4211	else if (strcmp("N_calls", name) == 0) // ID correction parameters
4212	sscanf(line, "%*s %d", &N_calls);
4213	else if (strcmp("N_steps", name) == 0)
4214	sscanf(line, "%*s %d", &N_steps);
4215	else if (strcmp("N_Fam", name) == 0)
4216	sscanf(line, "%*s %d", &N_Fam);
4217	else if (strcmp("IDCquads", name) == 0) {
4218	sscanf(line, "%*s %s %s %s %s %s %s %s %s",
4219	IDCq_name[0], IDCq_name[1], IDCq_name[2], IDCq_name[3],
4220	IDCq_name[4], IDCq_name[5], IDCq_name[6], IDCq_name[7]);
4221	// } else if (strcmp("scl_nu", name) == 0)
4222	// sscanf(line, "%*s %lf", &scl_nu);
4223	} else {
4224	printf("bad line in %s: %s\n", param_file, line);
4225	exit_(1);
4226	}
4227	} else
4228	printf("%s", line);
4229	}
4230
4231	fclose(inf);
4232
4233	if (N_calls > 0) {
4234	if (N_Fam > N_Fam_max) {
4235	printf("get_param: N_Fam > N_Fam_max: %d (%d)\n",
4236	N_Fam, N_Fam_max);
4237	exit_(0);
4238	}
4239
4240	for (k = 0; k < N_Fam; k++)
4241	Q_Fam[k] = ElemIndex(IDCq_name[k]);
4242	}
4243	}
4244
4245
4246	// control of vertical beam size
4247
4248	// output procedures
4249
4250	// Finding statistics on orbit in the sextupoles and maximal trim settings
4251	void Orb_and_Trim_Stat(void)
4252	{
4253	int i, j, N;
4254	int SextCounter = 0;
4255	int bins[5] = { 0, 0, 0, 0, 0 };
4256	double bin = 40.0e-6; // bin size for stat
4257	double tr; // trim strength
4258
4259	Vector2 Sext_max, Sext_sigma, TrimMax, orb;
4260
4261	Sext_max[X_] = 0.0; Sext_max[Y_] = 0.0;
4262	Sext_sigma[X_] = 0.0; Sext_sigma[Y_] = 0.0;
4263	TrimMax[X_] = 0.0; TrimMax[Y_] = 0.0;
4264	N = globval.Cell_nLoc; SextCounter = 0;
4265	for (i = 0; i <= N; i++) {
4266	if ((Cell[i].Elem.Pkind == Mpole) && (Cell[i].Elem.M->n_design == Sext)) {
4267	SextCounter++;
4268	orb[X_] = Cell[i].BeamPos[x_]; orb[Y_] = Cell[i].BeamPos[y_];
4269	Sext_sigma[X_] += orb[X_]orb[X_]; Sext_sigma[Y_] += orb[Y_]orb[Y_];
4270	if (fabs(orb[X_]) > Sext_max[X_]) Sext_max[X_] = fabs(orb[X_]);
4271	if (fabs(orb[Y_]) > Sext_max[Y_]) Sext_max[Y_] = fabs(orb[Y_]);
4272	j = (int) (sqrt(orb[X_]orb[X_]+orb[Y_]orb[Y_])/bin);
4273	if (j > 4) j = 4;
4274	bins[j]++;
4275	} // sextupole handling
4276
4277	if (Cell[i].Fnum == globval.hcorr) {
4278	tr = Cell[i].Elem.M->PBpar[HOMmax+Dip];
4279	if (fabs(tr) > TrimMax[X_]) TrimMax[X_] = fabs(tr);
4280	} // horizontal trim handling
4281	if (Cell[i].Fnum == globval.vcorr) {
4282	tr = Cell[i].Elem.M->PBpar[HOMmax-Dip];
4283	if (fabs(tr) > TrimMax[Y_]) TrimMax[Y_] = fabs(tr);
4284	} // vertical trim handling
4285	} // looking throught the cells
4286
4287	Sext_sigma[X_] = sqrt(Sext_sigma[X_]/SextCounter);
4288	Sext_sigma[Y_] = sqrt(Sext_sigma[Y_]/SextCounter);
4289
4290	printf("In sextupoles maximal horizontal orbit is:"
4291	" %5.3f mm with sigma %5.3f mm\n",
4292	1e3Sext_max[X_], 1e3Sext_sigma[X_]);
4293	printf("and maximal vertical orbit is: "
4294	" %5.3f mm with sigma %5.3f mm.\n",
4295	1e3Sext_max[Y_], 1e3Sext_sigma[Y_]);
4296
4297	for (i = 0; i < 4; i++) {
4298	printf("There are %3d sextupoles with offset between ", bins[i]);
4299	printf(" %5.3f mm and %5.3f mm\n", ibin1e3, (i+1)bin1e3);
4300	}
4301	printf("There are %3d sextupoles with offset ", bins[4]);
4302	printf("more than %5.3f mm \n", 4e3*bin);
4303	printf("Maximal hcorr is %6.3f mrad, maximal vcorr is %6.3f mrad\n",
4304	1e3TrimMax[X_], 1e3TrimMax[Y_]);
4305	}
4306
4307
4308	void prt_codcor_lat(void)
4309	{
4310	int i;
4311	FILE *CodCorLatFile;
4312
4313
4314	CodCorLatFile = file_write(CodCorLatFileName);
4315
4316	fprintf(CodCorLatFile, "# name s sqrt(BxBy) betaX nuX"
4317	" betaY nuY etaX etaX*betaY nuX-nuY \n");
4318	fprintf(CodCorLatFile, "# [m] [m] [m] [m]"
4319	" [m] [m*m] \n");
4320
4321	for (i = 0; i <= globval.Cell_nLoc; i++){
4322	fprintf(CodCorLatFile, "%4d:", i);
4323
4324	if (i == 0)
4325	fprintf(CodCorLatFile, "%.*s", 6, "begin ");
4326	else
4327	fprintf(CodCorLatFile, "%.*s", 6, Cell[i].Elem.PName);
4328
4329	fprintf(CodCorLatFile, "%7.3f %5.2f %5.2f %7.4f %5.2f %7.4f"
4330	" %6.3f %6.3f %6.3f\n",
4331	Cell[i].S, sqrt(Cell[i].Beta[X_]*Cell[i].Beta[Y_]),
4332	Cell[i].Beta[X_], Cell[i].Nu[X_], Cell[i].Beta[Y_], Cell[i].Nu[Y_],
4333	Cell[i].Eta[X_], Cell[i].Eta[X_]*Cell[i].Beta[Y_],
4334	Cell[i].Nu[X_]-Cell[i].Nu[Y_]);
4335	}
4336	fclose(CodCorLatFile);
4337	}
4338
4339	void prt_beamsizes()
4340	{
4341	int k;
4342	FILE *fp;
4343
4344	fp = file_write(beam_envelope_file);
4345
4346	fprintf(fp,"# k s name s_xx s_pxpx s_xpx s_yy s_pypy s_ypy theta_xy\n");
4347	for(k = 0; k <= globval.Cell_nLoc; k++){
4348	fprintf(fp,"%4d %10s %e %e %e %e %e %e %e %e\n",
4349	k, Cell[k].Elem.PName, Cell[k].S,
4350	Cell[k].sigma[x_][x_], Cell[k].sigma[px_][px_],
4351	Cell[k].sigma[x_][px_],
4352	Cell[k].sigma[y_][y_], Cell[k].sigma[py_][py_],
4353	Cell[k].sigma[y_][py_],
4354	atan2(2e0*Cell[k].sigma[x_][y_],
4355	Cell[k].sigma[x_][x_]-Cell[k].sigma[y_][y_])/2e0*180.0/M_PI);
4356	}
4357
4358	fclose(fp);
4359	}
4360
4361
4362	float f_int_Touschek(const float u)
4363	{
4364	double f;
4365
4366	if (u > 0.0)
4367	f = (1.0/u-log(1.0/u)/2.0-1.0)*exp(-u_Touschek/u);
4368	else
4369	f = 0.0;
4370
4371	return f;
4372	}
4373
4374	/****************************************************************************/
4375	/* double Touschek(const double Qb, const double delta_RF,
4376	const double eps_x, const double eps_y,
4377	const double sigma_delta, const double sigma_s)
4378
4379	Purpose:
4380	Calculate Touschek lifetime
4381	using the normal theoretical formula
4382	the momentum acceptance is RF acceptance delta_RF
4383
4384	Parameters:
4385	Qb electric charge per bunch
4386	delta_RF RF momentum acceptance
4387	eps_x emittance x
4388	eps_y emittance y
4389	sigma_delta longitudinal beam size
4390	sigma_s
4391
4392	Input:
4393	none
4394
4395	Output:
4396
4397
4398	Return:
4399	Touschek lifetime [second]
4400
4401	Global variables:
4402
4403
4404	Specific functions:
4405
4406
4407	Comments:
4408	none
4409
4410	****************************************************************************/
4411	double Touschek(const double Qb, const double delta_RF,
4412	const double eps_x, const double eps_y,
4413	const double sigma_delta, const double sigma_s)
4414	{
4415	long int k;
4416	double tau, sigma_x, sigma_y, sigma_xp, L, curly_H;
4417
4418	const double gamma = 1e9*globval.Energy/m_e, N_e = Qb/q_e;
4419
4420	printf("\n");
4421	printf("Qb = %4.2f nC, delta_RF = %4.2f%%"
4422	", eps_x = %9.3e nm.rad, eps_y = %9.3e nm.rad\n",
4423	1e9Qb, 1e2delta_RF, 1e9eps_x, 1e9eps_y);
4424	printf("sigma_delta = %8.2e, sigma_s = %4.2f mm\n",
4425	sigma_delta, 1e3*sigma_s);
4426
4427	tau = 0.0;
4428	for(k = 1; k <= globval.Cell_nLoc; k++) {
4429	L = Cell[k].S - Cell[k-1].S;
4430
4431	curly_H = get_curly_H(Cell[k].Alpha[X_], Cell[k].Beta[X_],
4432	Cell[k].Eta[X_], Cell[k].Etap[X_]);
4433
4434	// compute estimated beam sizes for given
4435	// hor., ver. emittance, sigma_s, and sigma_delta (x ~ 0)
4436	sigma_x = sqrt(Cell[k].Beta[X_]eps_x+sqr(sigma_deltaCell[k].Eta[X_]));
4437	sigma_y = sqrt(Cell[k].Beta[Y_]*eps_y);
4438	sigma_xp = (eps_x/sigma_x)sqrt(1.0+curly_Hsqr(sigma_delta)/eps_x);
4439
4440	u_Touschek = sqr(delta_RF/(gamma*sigma_xp));
4441
4442	tau += qromb(f_int_Touschek, 0.0, 1.0)/(sigma_xsigma_ysigma_xp)*L;
4443	fflush(stdout);
4444	}
4445
4446	tau = N_esqr(r_e)c0/(8.0M_PIcube(gamma)sigma_s)
4447	/(sqr(delta_RF)*Cell[globval.Cell_nLoc].S);
4448
4449	printf("\n");
4450	printf("Touschek lifetime [hrs]: %10.3e\n", 1.0/(3600.0*tau));
4451	return(1.0/(tau));
4452	}
4453
4454	/****************************************************************************/
4455	/* void mom_aper(double &delta, double delta_RF, const long int k,
4456	const int n_turn, const bool positive)
4457
4458
4459	Purpose:
4460	Using Binary search to determine momentum aperture at location k.
4461	If the particle is stable at the last step of tracking, so the
4462	momentum acceptance is RF momentum accpetance.
4463
4464	Input:
4465	delta particle energy deviation
4466	delta_RF RF momentum acceptance
4467	k element index
4468	n_turn number of tracking turns
4469
4470	Output:
4471	none
4472	Return:
4473	delta searched momentum acceptance
4474
4475	Global variables:
4476	none
4477
4478	specific functions:
4479	Cell_Pass
4480
4481	Comments:
4482	nsls-ii version.
4483	See also:
4484	soleil version MomentumAcceptance( ) in soleillib.cc
4485
4486	****************************************************************************/
4487	void mom_aper(double &delta, double delta_RF, const long int k,
4488	const int n_turn, const bool positive)
4489	{
4490	// Binary search to determine momentum aperture at location k.
4491	int j;
4492	long int lastpos;
4493	double delta_min, delta_max;
4494	Vector x;
4495
4496	const double eps = 1e-4;
4497
4498	delta_min = 0.0;
4499	delta_max = positive ? fabs(delta_RF) : -fabs(delta_RF);
4500
4501	while (fabs(delta_max-delta_min) > eps)
4502	{
4503	delta = (delta_max+delta_min)/2.0; /* binary serach*/
4504
4505	// propagate initial conditions
4506	CopyVec(6, globval.CODvect, x); Cell_Pass(0, k, x, lastpos);
4507	// generate Touschek event
4508	x[delta_] += delta;
4509	// complete one turn
4510	Cell_Pass(k, globval.Cell_nLoc, x, lastpos);
4511
4512	if (lastpos < globval.Cell_nLoc)
4513	// particle lost
4514	delta_max = delta;
4515	else {
4516	// track
4517	for(j = 0; j < n_turn; j++) {
4518	Cell_Pass(0, globval.Cell_nLoc, x, lastpos);
4519
4520	if ((delta_max > delta_RF) \|\| (lastpos < globval.Cell_nLoc)) {
4521	// particle lost
4522	delta_max = delta;
4523	break;
4524	}
4525	}
4526
4527	if ((delta_max <= delta_RF) && (lastpos == globval.Cell_nLoc))
4528	// particle not lost
4529	delta_min = delta;
4530	}
4531	}
4532	}
4533
4534	/****************************************************************************/
4535	/* double Touschek(const double Qb, const double delta_RF,const bool consistent,
4536	const double eps_x, const double eps_y,
4537	const double sigma_delta, double sigma_s,
4538	const int n_turn, const bool aper_on,
4539	double sum_delta[][2], double sum2_delta[][2])
4540
4541	Purpose:
4542	Calculate the Touschek lifetime
4543	using the normal theoretical formula
4544	the momentum acceptance is determined by the RF acceptance delta_RF
4545	and the momentum aperture at each element location
4546	which is tracked over n turns,
4547	the vacuum chamber limitation is read from a outside file
4548
4549	Parameters:
4550	Qb electric charge per bunch
4551	delta_RF RF momentum acceptance
4552	consistent
4553	eps_x emittance x
4554	eps_y emittance y
4555	sigma_delta longitudinal beam size
4556	sigma_s
4557	n_turn number of turns for tracking
4558	aper_on must be TRUE
4559	to indicate vacuum chamber setting is read from a outside file
4560	sum_delta[][2]
4561	sum2_delta[][2]
4562
4563
4564	Input:
4565	none
4566
4567	Output:
4568
4569
4570	Return:
4571	Touschek lifetime [second]
4572
4573	Global variables:
4574
4575
4576	Specific functions:
4577	mom_aper
4578
4579	Comments:
4580	none
4581
4582	****************************************************************************/
4583	double Touschek(const double Qb, const double delta_RF,const bool consistent,
4584	const double eps_x, const double eps_y,
4585	const double sigma_delta, double sigma_s,
4586	const int n_turn, const bool aper_on,
4587	double sum_delta[][2], double sum2_delta[][2])
4588	{
4589	bool cav, aper;
4590	long int k;
4591	double rate, delta_p, delta_m, curly_H0, curly_H1, L;
4592	double sigma_x, sigma_y, sigma_xp;
4593
4594	const bool prt = false;
4595
4596	// const char file_name[] = "Touschek.out";
4597	const double eps = 1e-5, gamma = 1e9*globval.Energy/m_e, N_e = Qb/q_e;
4598
4599	cav = globval.Cavity_on; aper = globval.Aperture_on;
4600
4601	globval.Cavity_on = true;
4602
4603	Ring_GetTwiss(true, 0.0);
4604
4605	globval.Aperture_on = aper_on;
4606
4607	printf("\n");
4608	printf("Qb = %4.2f nC, delta_RF = %4.2f%%"
4609	", eps_x = %9.3e m.rad, eps_y = %9.3e m.rad\n",
4610	1e9Qb, 1e2delta_RF, eps_x, eps_y);
4611	printf("sigma_delta = %8.2e, sigma_s = %4.2f mm\n",
4612	sigma_delta, 1e3*sigma_s);
4613
4614	printf("\n");
4615	printf("Momentum aperture:" );
4616	printf("\n");
4617
4618	delta_p = delta_RF; mom_aper(delta_p, delta_RF, 0, n_turn, true);
4619	delta_m = -delta_RF; mom_aper(delta_m, delta_RF, 0, n_turn, false);
4620	delta_p = min(delta_RF, delta_p); delta_m = max(-delta_RF, delta_m);
4621	sum_delta[0][0] += delta_p; sum_delta[0][1] += delta_m;
4622	sum2_delta[0][0] += sqr(delta_p); sum2_delta[0][1] += sqr(delta_m);
4623
4624	rate = 0.0; curly_H0 = -1e30;
4625	for (k = 1; k <= globval.Cell_nLoc; k++) {
4626	L = Cell[k].S - Cell[k-1].S;
4627
4628	curly_H1 = get_curly_H(Cell[k].Alpha[X_], Cell[k].Beta[X_],
4629	Cell[k].Eta[X_], Cell[k].Etap[X_]);
4630
4631	if (fabs(curly_H0-curly_H1) > eps) {
4632	mom_aper(delta_p, delta_RF, k, n_turn, true);
4633	delta_m = -delta_p; mom_aper(delta_m, delta_RF, k, n_turn, false);
4634	delta_p = min(delta_RF, delta_p); delta_m = max(-delta_RF, delta_m);
4635	printf("%4ld %6.2f %3.2lf%% %3.2lf%%\n",
4636	k, Cell[k].S, 1e2delta_p, 1e2delta_m);
4637	curly_H0 = curly_H1;
4638	}
4639
4640	sum_delta[k][0] += delta_p;
4641	sum_delta[k][1] += delta_m;
4642	sum2_delta[k][0] += sqr(delta_p);
4643	sum2_delta[k][1] += sqr(delta_m);
4644
4645	if (prt)
4646	printf("%4ld %6.2f %3.2lf %3.2lf\n",
4647	k, Cell[k].S, 1e2delta_p, 1e2delta_m);
4648
4649	if (!consistent) {
4650	// compute estimated beam sizes for given
4651	// hor., ver. emittance, sigma_s, and sigma_delta (x ~ 0)
4652	sigma_x = sqrt(Cell[k].Beta[X_]eps_x+sqr(sigma_deltaCell[k].Eta[X_]));
4653	sigma_y = sqrt(Cell[k].Beta[Y_]*eps_y);
4654	sigma_xp = (eps_x/sigma_x)sqrt(1.0+curly_H1sqr(sigma_delta)/eps_x);
4655	} else {
4656	// use self-consistent beam sizes
4657	sigma_x = sqrt(Cell[k].sigma[x_][x_]);
4658	sigma_y = sqrt(Cell[k].sigma[y_][y_]);
4659	sigma_xp = sqrt(Cell[k].sigma[px_][px_]);
4660	sigma_s = sqrt(Cell[k].sigma[ct_][ct_]);
4661	}
4662
4663	u_Touschek = sqr(delta_p/(gamma*sigma_xp));
4664
4665	rate += qromb(f_int_Touschek, 0.0, 1.0)
4666	/(sigma_xsigma_xpsigma_ysqr(delta_p))L;
4667
4668	fflush(stdout);
4669	}
4670
4671	rate = N_esqr(r_e)c0/(8.0M_PIcube(gamma)sigma_s)
4672	/Cell[globval.Cell_nLoc].S;
4673
4674	printf("\n");
4675	printf("Touschek lifetime [hrs]: %4.2f\n", 1.0/(3600.0*rate));
4676
4677	globval.Cavity_on = cav; globval.Aperture_on = aper;
4678
4679	return 1/rate;
4680	}
4681
4682
4683	double f_IBS(const double chi_m)
4684	{
4685	// Interpolated
4686
4687	double f, ln_chi_m;
4688
4689	const double A = 1.0, B = 0.579, C = 0.5;
4690
4691	ln_chi_m = log(chi_m); f = A + Bln_chi_m + Csqr(ln_chi_m);
4692
4693	return f;
4694	}
4695
4696
4697	float f_int_IBS(const float chi)
4698	{
4699	// Integrated
4700
4701	double f;
4702
4703	f = exp(-chi)*log(chi/chi_m)/chi;
4704
4705	return f;
4706	}
4707
4708
4709	void IBS(const double Qb,
4710	const double eps_SR[], double eps[],
4711	const double alpha_z, const double beta_z)
4712	{
4713	/* The equilibrium emittance is given by:
4714
4715	sigma_delta^2 = tau_delta*(D_delta_IBS + D_delta_SR)
4716
4717	D_x = <D_delta*curly_H> =>
4718
4719	eps_x = eps_x_SR + eps_x_IBS, D_x_IBS ~ 1/eps_x
4720
4721	i.e., the 2nd term is essentially constant. From
4722
4723	eps_x^2 = eps_xeps_x_SR + eps_xeps_x_IBS
4724
4725	one obtains the IBS limited emittance (eps_x_SR = 0, eps_x^2=eps_x_IBStau_xD_x_IBS->cst):
4726
4727	eps_x,IBS = sqrt(tau_x(eps_xD_x_IBS)
4728
4729	and by solving for eps_x
4730
4731	eps_x = eps_x_SR/2 + sqrt(eps_x_SR^2/4+eps_x_IBS^2)
4732
4733	*/
4734
4735	long int k;
4736	double D_x, D_delta, b_max, L, gamma_z;
4737	double sigma_x, sigma_xp, sigma_y, sigma_p, sigma_s, sigma_delta;
4738	double incr, C, curly_H, eps_IBS[3], eps_x_tot, sigma_E;
4739	double sigma_s_SR, sigma_delta_SR, sigma_delta_IBS;
4740
4741	const bool integrate = false;
4742	const double gamma = 1e9*globval.Energy/m_e, N_e = Qb/q_e;
4743
4744	// bunch size
4745	gamma_z = (1.0+sqr(alpha_z))/beta_z;
4746	sigma_s_SR = sqrt(beta_z*eps_SR[Z_]);
4747	sigma_delta_SR = sqrt(gamma_z*eps_SR[Z_]);
4748
4749	sigma_s = sqrt(beta_zeps[Z_]); sigma_delta = sqrt(gamma_zeps[Z_]);
4750
4751	printf("\n");
4752	printf("Qb = %4.2f nC\n", 1e9*Qb);
4753	printf("eps_x_SR = %9.3e m.rad, eps_x = %9.3e m.rad\n",
4754	eps_SR[X_], eps[X_]);
4755	printf("eps_y_SR = %9.3e m.rad, eps_y = %9.3e m.rad\n",
4756	eps_SR[Y_], eps[Y_]);
4757	printf("eps_z_SR = %9.3e, eps_z = %9.3e\n",
4758	eps_SR[Z_], eps[Z_]);
4759	printf("alpha_z = %9.3e, beta_z = %9.3e\n",
4760	alpha_z, beta_z);
4761	printf("sigma_s_SR = %9.3e mm, sigma_s = %9.3e mm\n",
4762	1e3sigma_s_SR, 1e3sigma_s);
4763	printf("sigma_delta_SR = %9.3e, sigma_delta = %9.3e\n",
4764	sigma_delta_SR, sigma_delta);
4765
4766	D_delta = 0.0; D_x = 0.0;
4767	for(k = 1; k <= globval.Cell_nLoc; k++) {
4768	L = Cell[k].Elem.PL;
4769
4770	curly_H = get_curly_H(Cell[k].Alpha[X_], Cell[k].Beta[X_],
4771	Cell[k].Eta[X_], Cell[k].Etap[X_]);
4772
4773	// compute estimated beam sizes for given
4774	// hor., ver. emittance, sigma_s, and sigma_delta (x ~ 0)
4775	sigma_x = sqrt(Cell[k].Beta[X_]eps[X_]+sqr(sigma_deltaCell[k].Eta[X_]));
4776	sigma_y = sqrt(Cell[k].Beta[Y_]*eps[Y_]);
4777	sigma_xp = (eps[X_]/sigma_x)sqrt(1.0+curly_Hsqr(sigma_delta)/eps[X_]);
4778
4779	sigma_E = gammam_eq_e*sigma_delta;
4780	sigma_p = gammam_eq_e*sigma_xp;
4781
4782	// b_max = sqrt(2.0M_PI)/pow(N_e/(sigma_xsigma_y*sigma_s), 1.0/3.0);
4783	b_max = 2.0sqrt(M_PI)/pow(N_e/(sigma_xsigma_y*sigma_s), 1.0/3.0);
4784
4785	chi_m = r_esqr(m_eq_e/sigma_E)/b_max;
4786	// chi_m = r_esqr(m_eq_e/sigma_p)/b_max;
4787
4788	if (!integrate)
4789	incr = f_IBS(chi_m)/sigma_y*L;
4790	else
4791	incr = qromb(f_int_IBS, chi_m, 1e5chi_m)/sigma_yL;
4792
4793	D_delta += incr; D_x += incr*curly_H;
4794	}
4795
4796	C = Cell[globval.Cell_nLoc].S;
4797
4798	if (true) {
4799	eps_x_tot = 1.0;
4800	D_delta = N_esqr(r_e)c0/(32.0M_PIcube(gamma)eps_x_totsigma_sC);
4801	D_x = N_esqr(r_e)c0/(32.0M_PIcube(gamma)eps_x_totsigma_sC);
4802	eps_IBS[X_] = sqrt(globval.tau[X_]*D_x);
4803	eps_x_tot = eps_SR[X_](1.0+sqrt(1.0+4.0sqr(eps_IBS[X_]/eps_SR[X_])))/2.0;
4804	D_delta /= eps_x_tot; D_x /= eps_x_tot;
4805	sigma_delta_IBS = sqrt(globval.tau[Z_]*D_delta);
4806	eps[X_] = eps_x_tot;
4807
4808	eps[Y_] = eps_SR[Y_]/eps_SR[X_]*eps[X_];
4809	sigma_delta = sqrt(sqr(sigma_delta_SR)+sqr(sigma_delta_IBS));
4810	eps_IBS[Z_] = sigma_s*sigma_delta;
4811	eps[Z_] = eps_SR[Z_] + eps_IBS[Z_];
4812	} else {
4813	D_delta = N_esqr(r_e)c0/(32.0M_PIcube(gamma)eps[X_]sigma_sC);
4814	D_x = N_esqr(r_e)c0/(32.0M_PIcube(gamma)eps[X_]sigma_sC);
4815	eps_IBS[X_] = globval.tau[X_]*D_x;
4816	eps_IBS[Z_] = globval.tau[Z_]*D_delta;
4817	eps[X_] = eps_SR[X_] + eps_IBS[X_];
4818	eps[Y_] = eps_SR[Y_]/eps_SR[X_]*eps[X_];
4819	sigma_delta_IBS = sqrt(globval.tau[Z_]*D_delta);
4820	sigma_delta = sqrt(sqr(sigma_delta_SR)+sqr(sigma_delta_IBS));
4821	// approx.
4822	eps_IBS[Z_] = sigma_s*sigma_delta;
4823	eps[Z_] = eps_SR[Z_] + eps_IBS[Z_];
4824	}
4825
4826	// bunch size
4827	sigma_s = sqrt(beta_z*eps[Z_]);
4828	sigma_delta = sqrt(gamma_z*eps[Z_]);
4829
4830	printf("\n");
4831	printf("D_x,IBS = %9.3e\n", D_x);
4832	printf("eps_x,IBS = %5.3f nm.rad\n", 1e9*eps_IBS[X_]);
4833	printf("eps_x,tot = %5.3f nm.rad\n", 1e9*eps[X_]);
4834	printf("eps_y,tot = %5.3f nm.rad\n", 1e9*eps[Y_]);
4835	printf("eps_z,tot = %9.3e\n", eps[Z_]);
4836	printf("D_delta,IBS = %9.3e\n", D_delta);
4837	printf("sigma_delta,IBS = %9.3e\n", sigma_delta_IBS);
4838	printf("sigma_delta,tot = %9.3e\n", sigma_delta);
4839	printf("sigma_s,tot = %9.3e\n", sigma_s);
4840	}
4841
4842
4843	void rm_space(char *name)
4844	{
4845	int i, k;
4846
4847	i = 0;
4848	while (name[i] == ' ')
4849	i++;
4850
4851	for (k = i; k <= (int)strlen(name)+1; k++)
4852	name[k-i] = name[k];
4853	}
4854
4855	/****************************************************************
4856	* void get_bn(char file_name[], int n, const bool prt)
4857	*
4858	*
4859	****************************************************************/
4860	void get_bn(char file_name[], int n, const bool prt)
4861	{
4862	char line[max_str], str[max_str], str1[max_str], token, name;
4863	int n_prm, Fnum, Knum, order;
4864	double bnL, bn, C, L;
4865	FILE inf, fp_lat;
4866
4867	inf = file_read(file_name); fp_lat = file_write("get_bn.lat");
4868
4869	// if n = 0: go to last data set
4870	if (n == 0) {
4871	while (fgets(line, max_str, inf) != NULL )
4872	if (strstr(line, "n = ") != NULL) sscanf(line, "n = %d", &n);
4873
4874	fclose(inf); inf = file_read(file_name);
4875	}
4876
4877	if (prt) {
4878	printf("\n");
4879	printf("reading values (n=%d): %s\n", n, file_name);
4880	printf("\n");
4881	}
4882
4883	sprintf(str, "n = %d", n);
4884	do
4885	fgets(line, max_str, inf);
4886	while (strstr(line, str) == NULL);
4887
4888	fprintf(fp_lat, "\n");
4889	n_prm = 0;
4890	while (fgets(line, max_str, inf) != NULL) {
4891	if (strcmp(line, "\n") == 0) break;
4892	n_prm++;
4893	name = strtok(line, "(");
4894	rm_space(name);
4895	strcpy(str, name); Fnum = ElemIndex(str);
4896	strcpy(str1, name); upr_case(str1);
4897	token = strtok(NULL, ")"); sscanf(token, "%d", &Knum);
4898	strtok(NULL, "="); token = strtok(NULL, "\n");
4899	sscanf(token, "%lf %d", &bnL, &order);
4900	if (prt) printf("%6s(%2d) = %10.6f %d\n", name, Knum, bnL, order);
4901	if (Fnum != 0) {
4902	if (order == 0)
4903	SetL(Fnum, bnL);
4904	else
4905	SetbnL(Fnum, order, bnL);
4906
4907	L = GetL(Fnum, 1);
4908	if (Knum == 1)
4909	if (order == 0)
4910	fprintf(fp_lat, "%s: Drift, L = %8.6f;\n", str1, bnL);
4911	else {
4912	bn = (L != 0.0)? bnL/L : bnL;
4913	if (order == Quad)
4914	fprintf(fp_lat, "%s: Quadrupole, L = %8.6f, K = %10.6f, N = Nquad"
4915	", Method = Meth;\n", str1, L, bn);
4916	else if (order == Sext)
4917	fprintf(fp_lat, "%s: Sextupole, L = %8.6f, K = %10.6f"
4918	", N = 1, Method = Meth;\n", str1, L, bn);
4919	else {
4920	fprintf(fp_lat, "%s: Multipole, L = %8.6f"
4921	", N = 1, Method = Meth,\n", str1, L);
4922	fprintf(fp_lat, " HOM = (%d, %13.6f, %3.1f);\n",
4923	order, bn, 0.0);
4924	}
4925	}
4926	} else {
4927	printf("element %s not found\n", name);
4928	exit_(1);
4929	}
4930	}
4931	if (prt) printf("\n");
4932
4933	C = Cell[globval.Cell_nLoc].S; recalc_S();
4934	if (prt)
4935	printf("New Cell Length: %5.3f (%5.3f)\n", Cell[globval.Cell_nLoc].S, C);
4936
4937	fclose(inf); fclose(fp_lat);
4938	}
4939
4940
4941	double get_dynap(const double delta)
4942	{
4943	char str[max_str];
4944	int i;
4945	double x_aper[n_aper], y_aper[n_aper], DA;
4946	FILE *fp;
4947
4948	const int prt = true;
4949
4950	fp = file_write("dynap.out");
4951	dynap(fp, 5e-3, 0.0, 0.1e-3, n_aper, n_track, x_aper, y_aper, false, prt);
4952	fclose(fp);
4953	DA = get_aper(n_aper, x_aper, y_aper);
4954
4955	if (true) {
4956	sprintf(str, "dynap_dp%3.1f.out", 1e2*delta);
4957	fp = file_write(str);
4958	dynap(fp, 5e-3, delta, 0.1e-3, n_aper, n_track,
4959	x_aper, y_aper, false, prt);
4960	fclose(fp);
4961	DA += get_aper(n_aper, x_aper, y_aper);
4962
4963	for (i = 0; i < nv_; i++)
4964	globval.CODvect[i] = 0.0;
4965	sprintf(str, "dynap_dp%3.1f.out", -1e2*delta);
4966	fp = file_write(str);
4967	dynap(fp, 5e-3, -delta, 0.1e-3, n_aper,
4968	n_track, x_aper, y_aper, false, prt);
4969	fclose(fp);
4970	DA += get_aper(n_aper, x_aper, y_aper);
4971	}
4972
4973	return DA/3.0;
4974	}
4975
4976
4977	double get_chi2(long int n, double x[], double y[], long int m, Vector b)
4978	{
4979	/* Compute chi2 for polynomial fit */
4980
4981	int i, j;
4982	double sum, z;
4983
4984	sum = 0.0;
4985	for (i = 0; i < n; i++) {
4986	z = 0.0;
4987	for (j = 0; j <= m; j++)
4988	z += b[j]*pow(x[i], j);
4989	sum += pow(y[i]-z, 2);
4990	}
4991	return(sum/n);
4992	}
4993
4994
4995	void pol_fit(int n, double x[], double y[], int order, Vector &b,
4996	double &sigma)
4997	{
4998	/* Polynomial fit by linear chi-square */
4999
5000	int i, j, k;
5001	Matrix T1;
5002	bool prt = false;
5003
5004	const double sigma_k = 1.0, chi2 = 4.0;
5005
5006	/* initialization */
5007	for (i = 0; i <= order; i++) {
5008	b[i] = 0.0;
5009	for (j = 0; j <= order; j++)
5010	T1[i][j] = 0.0;
5011	}
5012
5013	/* polynomal fiting */
5014	for (i = 0; i < n; i++)
5015	for (j = 0; j <= order; j++) {
5016	b[j] += y[i]*pow(x[i], j)/pow(sigma_k, 2);
5017	for (k = 0; k <= order; k++)
5018	T1[j][k] += pow(x[i], j+k)/pow(sigma_k, 2);
5019	}
5020
5021	if (InvMat(order+1, T1)) {
5022	LinTrans(order+1, T1, b); sigma = get_chi2(n, x, y, order, b);
5023	if(prt){
5024	printf("\n");
5025	printf(" n Coeff.\n");
5026	for (i = 0; i <= order; i++)
5027	printf("%3d %10.3e+/-%8.2e\n", i, b[i], sigmasqrt(chi2T1[i][i]));
5028	}
5029	} else {
5030	printf("pol_fit: Matrix is singular\n");
5031	exit_(1);
5032	}
5033	}
5034
5035
5036	void get_ksi2(const double d_delta)
5037	{
5038	const int n_points = 20, order = 5;
5039
5040	int i, n;
5041	double delta[2n_points+1], nu[2][2n_points+1], sigma;
5042	Vector b;
5043	FILE *fp;
5044
5045	fp = file_write("chrom2.out");
5046	n = 0;
5047	for (i = -n_points; i <= n_points; i++) {
5048	n++; delta[n-1] = i*(double)d_delta/(double)n_points;
5049	Ring_GetTwiss(false, delta[n-1]);
5050	nu[0][n-1] = globval.TotalTune[X_]; nu[1][n-1] = globval.TotalTune[Y_];
5051	fprintf(fp, "%5.2f %8.5f %8.5f\n", 1e2*delta[n-1], nu[0][n-1], nu[1][n-1]);
5052	}
5053	printf("\n");
5054	printf("Horizontal chromaticity:\n");
5055	pol_fit(n, delta, nu[0], order, b, sigma);
5056	printf("\n");
5057	printf("Vertical chromaticity:\n");
5058	pol_fit(n, delta, nu[1], order, b, sigma);
5059	printf("\n");
5060	fclose(fp);
5061	}
5062
5063
5064	bool find_nu(const int n, const double nus[], const double eps, double &nu)
5065	{
5066	bool lost;
5067	int k;
5068
5069	const bool prt = false;
5070
5071	if (prt) printf("\n");
5072	k = 0;
5073	while ((k < n) && (fabs(nus[k]-nu) > eps)) {
5074	if (prt) printf("nu = %7.5f(%7.5f) eps = %7.1e\n", nus[k], nu, eps);
5075	k++;
5076	}
5077
5078	if (k < n) {
5079	if (prt) printf("nu = %7.5f(%7.5f)\n", nus[k], nu);
5080	nu = nus[k]; lost = false;
5081	} else {
5082	if (prt) printf("lost\n");
5083	lost = true;
5084	}
5085
5086	return !lost;
5087	}
5088
5089
5090	bool get_nu(const double Ax, const double Ay, const double delta,
5091	const double eps, double &nu_x, double &nu_y)
5092	{
5093	const int n_turn = 512, n_peaks = 5;
5094
5095	bool lost, ok_x, ok_y;
5096	char str[max_str];
5097	int i;
5098	long int lastpos, lastn, n;
5099	double x[n_turn], px[n_turn], y[n_turn], py[n_turn];
5100	double nu[2][n_peaks], A[2][n_peaks];
5101	Vector x0;
5102	FILE *fp;
5103
5104	const bool prt = false;
5105	const char file_name[] = "track.out";
5106
5107	// complex FFT in Floquet space
5108	x0[x_] = Ax; x0[px_] = 0.0; x0[y_] = Ay; x0[py_] = 0.0;
5109	LinTrans(4, globval.Ascrinv, x0);
5110	track(file_name, x0[x_], x0[px_], x0[y_], x0[py_], delta,
5111	n_turn, lastn, lastpos, 1, 0.0);
5112	if (lastn == n_turn) {
5113	GetTrack(file_name, &n, x, px, y, py);
5114	sin_FFT((int)n, x, px); sin_FFT((int)n, y, py);
5115
5116	if (prt) {
5117	strcpy(str, file_name); strcat(str, ".fft");
5118	fp = file_write(str);
5119	for (i = 0; i < n; i++)
5120	fprintf(fp, "%5.3f %12.5e %8.5f %12.5e %8.5f\n",
5121	(double)i/(double)n, x[i], px[i], y[i], py[i]);
5122	fclose(fp);
5123	}
5124
5125	GetPeaks1(n, x, n_peaks, nu[X_], A[X_]);
5126	GetPeaks1(n, y, n_peaks, nu[Y_], A[Y_]);
5127
5128	ok_x = find_nu(n_peaks, nu[X_], eps, nu_x);
5129	ok_y = find_nu(n_peaks, nu[Y_], eps, nu_y);
5130
5131	lost = !ok_x \|\| !ok_y;
5132	} else
5133	lost = true;
5134
5135	return !lost;
5136	}
5137
5138
5139	// tune shift with amplitude
5140	void dnu_dA(const double Ax_max, const double Ay_max, const double delta)
5141	{
5142	bool ok;
5143	int i;
5144	double nu_x, nu_y, Ax, Ay, Jx, Jy;
5145	Vector ps;
5146	FILE *fp;
5147
5148	const int n_ampl = 25;
5149	const double A_min = 1e-3;
5150	// const double eps0 = 0.04, eps = 0.02;
5151	// const double eps0 = 0.025, eps = 0.02;
5152	const double eps0 = 0.04, eps = 0.015;
5153
5154	Ring_GetTwiss(false, 0.0);
5155
5156	ok = false;
5157	fp = file_write("dnu_dAx.out");
5158	fprintf(fp, "# x[mm] y[mm] J_x J_y nu_x nu_y\n");
5159	fprintf(fp, "#\n");
5160	for (i = -n_ampl; i <= n_ampl; i++) {
5161	Ax = i*Ax_max/n_ampl;
5162	if (Ax == 0.0) Ax = A_min;
5163	Ay = A_min;
5164	ps[x_] = Ax; ps[px_] = 0.0; ps[y_] = Ay; ps[py_] = 0.0; getfloqs(ps);
5165	Jx = (sqr(ps[x_])+sqr(ps[px_]))/2.0; Jy = (sqr(ps[y_])+sqr(ps[py_]))/2.0;
5166	if (!ok) {
5167	nu_x = fract(globval.TotalTune[X_]); nu_y = fract(globval.TotalTune[Y_]);
5168	ok = get_nu(Ax, Ay, delta, eps0, nu_x, nu_y);
5169	} else
5170	ok = get_nu(Ax, Ay, delta, eps, nu_x, nu_y);
5171	if (ok)
5172	fprintf(fp, "%10.3e %10.3e %10.3e %10.3e %7.5f %7.5f\n",
5173	1e3Ax, 1e3Ay, 1e6Jx, 1e6Jy, fract(nu_x), fract(nu_y));
5174	else
5175	fprintf(fp, "# %10.3e %10.3e particle lost\n", 1e3Ax, 1e3Ay);
5176	}
5177	fclose(fp);
5178
5179	ok = false;
5180	fp = file_write("dnu_dAy.out");
5181	fprintf(fp, "# x[mm] y[mm] J_x J_y nu_x nu_y\n");
5182	fprintf(fp, "#\n");
5183	for (i = -n_ampl; i <= n_ampl; i++) {
5184	Ax = A_min; Ay = i*Ay_max/n_ampl;
5185	Jx = pow(Ax, 2)/(2.0*Cell[globval.Cell_nLoc].Beta[X_]);
5186	if (Ay == 0.0) Ay = A_min;
5187	Jy = pow(Ay, 2)/(2.0*Cell[globval.Cell_nLoc].Beta[Y_]);
5188	// if (i == 1) exit_(0);
5189	if (!ok) {
5190	nu_x = fract(globval.TotalTune[X_]); nu_y = fract(globval.TotalTune[Y_]);
5191	ok = get_nu(Ax, Ay, delta, eps0, nu_x, nu_y);
5192	} else
5193	ok = get_nu(Ax, Ay, delta, eps, nu_x, nu_y);
5194	if (ok)
5195	fprintf(fp, "%10.3e %10.3e %10.3e %10.3e %7.5f %7.5f\n",
5196	1e3Ax, 1e3Ay, 1e6Jx, 1e6Jy, fract(nu_x), fract(nu_y));
5197	else
5198	fprintf(fp, "# %10.3e %10.3e particle lost\n", 1e3Ax, 1e3Ay);
5199	}
5200	fclose(fp);
5201	}
5202
5203
5204	bool orb_corr(const int n_orbit)
5205	{
5206	bool cod = false;
5207	int i;
5208	long lastpos;
5209	Vector2 xmean, xsigma, xmax;
5210
5211	printf("\n");
5212
5213	// FitTune(qf, qd, nu_x, nu_y);
5214	// printf("\n");
5215	// printf(" qf = %8.5f qd = %8.5f\n",
5216	// GetKpar(qf, 1, Quad), GetKpar(qd, 1, Quad));
5217
5218	globval.CODvect.zero();
5219	for (i = 1; i <= n_orbit; i++) {
5220	cod = getcod(0.0, lastpos);
5221	if (cod) {
5222	codstat(xmean, xsigma, xmax, globval.Cell_nLoc, false);
5223	printf("\n");
5224	printf("RMS orbit [mm]: (%8.1e+/-%7.1e, %8.1e+/-%7.1e)\n",
5225	1e3xmean[X_], 1e3xsigma[X_], 1e3xmean[Y_], 1e3xsigma[Y_]);
5226	lsoc(1, globval.bpm, globval.hcorr, 1);
5227	lsoc(1, globval.bpm, globval.vcorr, 2);
5228	cod = getcod(0.0, lastpos);
5229	if (cod) {
5230	codstat(xmean, xsigma, xmax, globval.Cell_nLoc, false);
5231	printf("RMS orbit [mm]: (%8.1e+/-%7.1e, %8.1e+/-%7.1e)\n",
5232	1e3xmean[X_], 1e3xsigma[X_], 1e3xmean[Y_], 1e3xsigma[Y_]);
5233	} else
5234	printf("orb_corr: failed\n");
5235	} else
5236	printf("orb_corr: failed\n");
5237	}
5238
5239	prt_cod("orb_corr.out", globval.bpm, true);
5240
5241	return cod;
5242	}
5243
5244
5245	double get_code(CellType &Cell)
5246	{
5247	double code = 0.0;
5248
5249	switch (Cell.Elem.Pkind) {
5250	case drift:
5251	code = 0.0;
5252	break;
5253	case Mpole:
5254	if (Cell.Elem.M->Pirho != 0.0) {
5255	code = 0.5;
5256	} else if (Cell.Elem.M->PBpar[Quad+HOMmax] != 0)
5257	code = sgn(Cell.Elem.M->PBpar[Quad+HOMmax]);
5258	else if (Cell.Elem.M->PBpar[Sext+HOMmax] != 0)
5259	code = 1.5*sgn(Cell.Elem.M->PBpar[Sext+HOMmax]);
5260	else if (Cell.Fnum == globval.bpm)
5261	code = 2.0;
5262	else
5263	code = 0.0;
5264	break;
5265	default:
5266	code = 0.0;
5267	break;
5268	}
5269
5270	return code;
5271	}
5272
5273
5274	void get_alphac(void)
5275	{
5276	CellType Cell;
5277
5278	getelem(globval.Cell_nLoc, &Cell);
5279	globval.Alphac = globval.OneTurnMat[ct_][delta_]/Cell.S;
5280	}
5281
5282	/******************************************************************
5283	void get_alphac2(void)
5284
5285	Purpose:
5286	Calculate momentum compact factor up to third order.
5287	Method:
5288	track the orbit offset c*tau at the first element,
5289	at 11 different energy offset(-10-3 to 10-3), then
5290	use polynomal to fit the momentum compact faction factor
5291	up to 3rd order.
5292	The initial coorinates are (x_co,px_co,y_co,py_co,delta,0).
5293
5294	Input:
5295	none
5296
5297	Output:
5298	none
5299
5300	Return:
5301	alphac
5302
5303	specific functions:
5304	none
5305
5306	Comments:
5307	change the initial coorinates from (0 0 0 0 delta 0) to
5308	(x_co,px_co,y_co,py_co,delta,0). i.e., track
5309	around the off-momentum close orbit but not zero orbit.
5310
5311	*******************************************************************/
5312
5313	void get_alphac2(void)
5314	{
5315	/* Note, do not extract from M[5][4], i.e. around delta
5316	dependent fixed point. */
5317
5318	const int n_points = 5;
5319	const double d_delta = 1e-3;
5320	//const double d_delta = 1e-4;
5321
5322	int i, n;
5323	long int lastpos; // last tracking position when the particle is stable
5324	double delta[2n_points+1], alphac[2n_points+1], sigma;
5325	Vector x, b;
5326	CellType Cell;
5327	bool cod;
5328
5329	globval.pathlength = false;
5330	getelem(globval.Cell_nLoc, &Cell);
5331	n = 0;
5332
5333	for (i = -n_points; i <= n_points; i++) {
5334	n++;
5335	delta[n-1] = i*(double)d_delta/(double)n_points;
5336
5337	// for (j = 0; j < nv_; j++)
5338	// x[j] = 0.0;
5339	// x[delta_] = delta[n-1];
5340
5341	cod = getcod(delta[n-1], lastpos);
5342	x = globval.CODvect;
5343
5344
5345	Cell_Pass(0, globval.Cell_nLoc, x, lastpos);
5346	alphac[n-1] = x[ct_]/Cell.S; // this is not mcf but DT/T
5347	}
5348	pol_fit(n, delta, alphac, 3, b, sigma);
5349	printf("Momentum compaction factor:\n");
5350	printf(" alpha1 = %10.4e alpha2 = %10.4e alpha3 = %10.4e\n",
5351	b[1], b[2], b[3]);
5352	}
5353
5354
5355	float f_bend(float b0L[])
5356	{
5357	long int lastpos;
5358	Vector ps;
5359
5360	const int n_prt = 10;
5361
5362	n_iter_Cart++;
5363
5364	SetbnL_sys(Fnum_Cart, Dip, b0L[1]);
5365
5366	ps.zero();
5367	Cell_Pass(Elem_GetPos(Fnum_Cart, 1)-1, Elem_GetPos(Fnum_Cart, 1),
5368	ps, lastpos);
5369
5370	if (n_iter_Cart % n_prt == 0)
5371	cout << scientific << setprecision(3)
5372	<< setw(4) << n_iter_Cart
5373	<< setw(11) << ps[x_] << setw(11) << ps[px_] << setw(11) << ps[ct_]
5374	<< setw(11) << b0L[1] << endl;
5375
5376	return sqr(ps[x_]) + sqr(ps[px_]);
5377	}
5378
5379
5380	void bend_cal_Fam(const int Fnum)
5381	{
5382	/* Adjusts b1L_sys to zero the orbit for a given gradient. */
5383	const int n_prm = 1;
5384
5385	int iter;
5386	float b0L, *xi, fret;
5387	Vector ps;
5388
5389	const float ftol = 1e-15;
5390
5391	b0L = vector(1, n_prm); xi = matrix(1, n_prm, 1, n_prm);
5392
5393	cout << endl;
5394	cout << "bend_cal: " << ElemFam[Fnum-1].ElemF.PName << ":" << endl;
5395
5396	Fnum_Cart = Fnum; b0L[1] = 0.0; xi[1][1] = 1e-3;
5397
5398	cout << endl;
5399	n_iter_Cart = 0;
5400	powell(b0L, xi, n_prm, ftol, &iter, &fret, f_bend);
5401
5402	free_vector(b0L, 1, n_prm); free_matrix(xi, 1, n_prm, 1, n_prm);
5403	}
5404
5405
5406	void bend_cal(void)
5407	{
5408	long int k;
5409
5410	for (k = 1; k <= globval.Elem_nFam; k++)
5411	if ((ElemFam[k-1].ElemF.Pkind == Mpole) &&
5412	(ElemFam[k-1].ElemF.M->Pirho != 0.0) &&
5413	(ElemFam[k-1].ElemF.M->PBpar[Quad+HOMmax] != 0.0))
5414	if (ElemFam[k-1].nKid > 0) bend_cal_Fam(k);
5415	}
5416
5417
5418	double h_ijklm(const tps &h, const int i, const int j, const int k,
5419	const int l, const int m)
5420	{
5421	int i1;
5422	iVector jj;
5423
5424	for (i1 = 0; i1 < nv_tps; i1++)
5425	jj[i1] = 0;
5426	jj[x_] = i; jj[px_] = j; jj[y_] = k; jj[py_] = l; jj[delta_] = m;
5427	return h[jj];
5428	}
5429
5430
5431	ss_vect<tps> get_A(const double alpha[], const double beta[],
5432	const double eta[], const double etap[])
5433	{
5434	ss_vect<tps> A, Id;
5435
5436	Id.identity();
5437
5438	A.identity();
5439	A[x_] = sqrt(beta[X_])*Id[x_];
5440	A[px_] = -alpha[X_]/sqrt(beta[X_])Id[x_] + 1.0/sqrt(beta[X_])Id[px_];
5441	A[y_] = sqrt(beta[Y_])*Id[y_];
5442	A[py_] = -alpha[Y_]/sqrt(beta[Y_])Id[y_] + 1.0/sqrt(beta[Y_])Id[py_];
5443
5444	A[x_] += eta[X_]Id[delta_]; A[px_] += etap[X_]Id[delta_];
5445
5446	return A;
5447	}
5448
5449
5450	void get_ab(const ss_vect<tps> &A,
5451	double alpha[], double beta[], double eta[], double etap[])
5452	{
5453	int k;
5454	ss_vect<tps> A_Atp;
5455
5456	A_Atp = A*tp_S(2, A);
5457
5458	for (k = 0; k <= 1; k++) {
5459	eta[k] = A[2*k][delta_];
5460	etap[k] = A[2*k+1][delta_];
5461	alpha[k] = -A_Atp[2k][2k+1];
5462	beta[k] = A_Atp[2k][2k];
5463	}
5464	}
5465
5466
5467	void set_tune(const char file_name1[], const char file_name2[], const int n)
5468	{
5469	const int n_b2 = 8;
5470
5471	char line[max_str], names[n_b2][max_str];
5472	int j, k, Fnum;
5473	double b2s[n_b2], nu[2];
5474	ifstream inf1, inf2;
5475	FILE *fp_lat;
5476
5477	file_rd(inf1, file_name1); file_rd(inf2, file_name2);
5478
5479	// skip 1st and 2nd line
5480	inf1.getline(line, max_str); inf2.getline(line, max_str);
5481	inf1.getline(line, max_str); inf2.getline(line, max_str);
5482
5483	inf1.getline(line, max_str);
5484	sscanf(line, "%s %s %*s %s %s %s %s",
5485	names[0], names[1], names[2], names[3]);
5486	inf2.getline(line, max_str);
5487	sscanf(line, "%s %s %*s %s %s %s %s",
5488	names[4], names[5], names[6], names[7]);
5489
5490	printf("\n");
5491	printf("quads:");
5492	for (k = 0; k < n_b2; k++) {
5493	lwr_case(names[k]); rm_space(names[k]);
5494	printf(" %s", names[k]);
5495	}
5496	printf("\n");
5497
5498	// skip 4th line
5499	inf1.getline(line, max_str); inf2.getline(line, max_str);
5500
5501	fp_lat = file_write("set_tune.lat");
5502	while (inf1.getline(line, max_str) != NULL) {
5503	sscanf(line, "%d%lf %lf %lf %lf %lf %lf",
5504	&j, &nu[X_], &nu[Y_], &b2s[0], &b2s[1], &b2s[2], &b2s[3]);
5505	inf2.getline(line, max_str);
5506	sscanf(line, "%d%lf %lf %lf %lf %lf %lf",
5507	&j, &nu[X_], &nu[Y_], &b2s[4], &b2s[5], &b2s[6], &b2s[7]);
5508
5509	if (j == n) {
5510	printf("\n");
5511	printf("%3d %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f %8.5f\n",
5512	n,
5513	b2s[0], b2s[1], b2s[2], b2s[3], b2s[4], b2s[5], b2s[6], b2s[7]);
5514
5515	for (k = 0; k < n_b2; k++) {
5516	Fnum = ElemIndex(names[k]);
5517	set_bn_design_fam(Fnum, Quad, b2s[k], 0.0);
5518
5519	fprintf(fp_lat, "%s: Quadrupole, L = %8.6f, K = %10.6f, N = Nquad"
5520	", Method = Meth;\n",
5521	names[k], ElemFam[Fnum-1].ElemF.PL, b2s[k]);
5522	}
5523	break;
5524	}
5525	}
5526
5527	fclose(fp_lat);
5528	}
5529
5530
5531
5532
5533

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

Download in other formats: