Context Navigation

nsls-ii_lib.cc @ 32

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

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

Context Navigation

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

Download in other formats: