Context Navigation

source: TRACY3/trunk/tracy/tracy/src/t2ring.cc @ 11

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

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

Line
1	/* Tracy-2
2
3	J. Bengtsson, CBP, LBL 1990 - 1994 Pascal version
4	SLS, PSI 1995 - 1997
5	M. Boege SLS, PSI 1998 C translation
6	L. Nadolski SOLEIL 2002 Link to NAFF, Radia field maps
7	J. Bengtsson NSLS-II, BNL 2004 -
8
9	t2ring.c -- Routines for closed beam lines
10
11	*/
12
13
14	#define n 4
15
16
17	void GetNu(Vector2 &nu, Matrix &M)
18	{
19	int i;
20	double sgn, detp, detm, b, c, th, tv, b2mc, x_arg, y_arg;
21	Matrix M1;
22
23	globval.stable = true;
24
25	CopyMat((long)n, M, M1);
26
27	for (i = 0; i < n; i++)
28	M1[i][i] -= 1.0;
29	detp = DetMat((long)n, M1); /* det(M-I) */
30
31	for (i = 0; i < n; i++)
32	M1[i][i] += 2.0;
33	detm = DetMat((long)n, M1); /* det(M+I) */
34
35	for (i = 0; i < n; i++)
36	M1[i][i] -= 1.0; /* restore M */
37
38	b = (detp-detm)/16.0; c = (detp+detm)/8.0 - 1.0;
39	th = (M1[0][0]+M1[1][1])/2.0; tv = (M1[2][2]+M1[3][3])/2.0;
40
41	b2mc = b*b - c;
42
43	if (b2mc < 0.0) {
44	globval.stable = false; nu[0] = -1.0; nu[1] = -1.0;
45	printf("GetNu: unstable\n");
46	return;
47	}
48
49	sgn = (th > tv)? 1.0 : -1.0;
50
51	x_arg = -b + sgn*sqrt(b2mc);
52	if (fabs(x_arg) <= 1.0) {
53	nu[X_] = acos(x_arg)/(2.0*M_PI);
54	if (M1[0][1] < 0.0) nu[X_] = 1.0 - nu[X_];
55	} else {
56	globval.stable = false; nu[X_] = -1.0;
57	printf("GetNu: unstable (horizontal)\n");
58	}
59
60	y_arg = -b - sgn*sqrt(b2mc);
61	if (fabs(y_arg) <= 1.0) {
62	nu[Y_] = acos(y_arg)/(2.0*M_PI);
63	if (M1[2][3] < 0.0) nu[Y_] = 1.0 - nu[Y_];
64	} else {
65	globval.stable = false; nu[Y_] = -1.0;
66	printf("GetNu: unstable (vertical)\n");
67	return;
68	}
69
70	return;
71	}
72	#undef n
73
74
75	/* implementation */
76
77	/****************************************************************************/
78	/* void Cell_GetABGN(Matrix &M_, Vector2 &alpha, Vector2 &beta, Vector2 &gamma, Vector2 &nu)
79
80	Purpose: called by Ring_Twiss_M and Ring_Twiss
81
82	Get Twiss parameters from the transfer matrix M
83
84	[Nx ]
85	M = [ ]
86	[ Ny]
87
88	where, in the case of mid-plane symmetry
89
90	[ cos(mu) + alpha sin(mu) beta sin(mu) ]
91	N = [ ]
92	[ -gamma sin(mu) cos(mu) - alpha sin(mu) ]
93
94	cos(mu) = Trace(N)/2
95	Alpha = (N11-N22)/(2*sin(mu))
96	beta = N12/sin(mu)
97	gamma = -N21/sin(mu)
98
99	Input:
100	M oneturn matrix
101
102	Output:
103	alpha, beta, gamma vectors of the Twiss parameters
104	nu tune vector
105
106	Return:
107	global.stable
108	global value to check the linear matrix is stable or not,
109	true stable,
110	false unstable
111
112	Global variables:
113	none
114
115	Specific functions:
116	getnu
117
118	Comments:
119	none
120
121	****************************************************************************/
122	void Cell_GetABGN(Matrix &M,
123	Vector2 &alpha, Vector2 &beta, Vector2 &gamma, Vector2 &nu)
124	{
125	int i = 0, j = 0;
126	double c = 0.0, s = 0.0;
127
128	globval.stable = true;
129	for (i = 0; i <= 1; i++) {
130	j = (i+1)*2 - 1;
131	/* c=cos(mu) = Tr(N)/2 */
132	c = (M[j-1][j-1]+M[j][j])/2.0;
133	globval.stable = (fabs(c) < 1.0);
134	if (globval.stable) {
135	// s = sin(mu)
136	s = sqrt(1.0-cc)sgn(M[j-1][j]);
137	alpha[i] = (M[j-1][j-1]-M[j][j])/(2.0*s); beta[i] = M[j-1][j]/s;
138	gamma[i] = -(M[j][j-1]/s);
139	// nu[i] = acos(c)/2/pi;
140	GetNu(nu, M);
141	}
142	}
143	}
144
145
146	#define n 4
147	void Cell_Geteta(long i0, long i1, bool ring, double dP)
148	{
149	long int i = 0, lastpos = 0;
150	int j = 0, k = 0;
151	Vector xref;
152	Vector codbuf[Cell_nLocMax+1];
153	CellType *cellp;
154
155	/* cod for the energy dP - globval.dPcommon / 2e0 */
156	if (ring)
157	GetCOD(globval.CODimax, globval.CODeps, dP-globval.dPcommon/2e0, lastpos);
158	else { /* beam mode */
159	CopyVec(n+2, globval.CODvect, xref); xref[4] = dP - globval.dPcommon/2e0;
160	Cell_Pass(i0, i1, xref, lastpos);
161	}
162
163	/* Store chromatic orbit for elements i0 to i1 in codbuf */
164	for (i = i0; i <= i1; i++)
165	CopyVec(n+2, Cell[i].BeamPos, codbuf[i]);
166
167	/* cod for the energy dP - globval.dPcommon / 2e0 */
168	if (ring)
169	GetCOD(globval.CODimax, globval.CODeps, dP+globval.dPcommon/2e0, lastpos);
170	else { /* beam mode */
171	CopyVec(n+2, globval.CODvect, xref); xref[4] = dP + globval.dPcommon/2e0;
172	Cell_Pass(i0, i1, xref, lastpos);
173	}
174
175	/* cod(dP-globval.dPcommon/2e0) - cod(dP-globval.dPcommon/2e0) */
176	/* eta = ----------------------------------------------------------- */
177	/* globval.dPcommon */
178	/* cod'(dP-globval.dPcommon/2e0) - cod'(dP-globval.dPcommon/2e0) */
179	/* eta'= --------------------------------------------------------------- */
180	/* globval.dPcommon */
181
182	for (i = i0; i <= i1; i++) {
183	cellp = &Cell[i];
184	for (j = 1; j <= 2; j++) {
185	k = j*2 - 1;
186	cellp->Eta[j-1] = (cellp->BeamPos[k-1]-codbuf[i][k-1])/globval.dPcommon;
187	cellp->Etap[j-1] = (cellp->BeamPos[k]-codbuf[i][k])/globval.dPcommon;
188	}
189	}
190	}
191
192	#undef n
193
194	/****************************************************************************/
195	/* void getprm(Matrix &Ascr, Vector2 &alpha, Vector2 &beta)
196
197	Purpose:
198	Computes Twiss parameters alpha and beta from the matrix Ascr
199
200	M oneturn matrix (A and M are symplectic)
201
202	( A11 A12) -1 ( cos(mu) sin(mu))
203	A = ( ) M = (A R A) with R = ( )
204	( A21 A22) 2 2 (-sin(mu) cos(mu))
205	-1 x + px
206	eps = 2*J o a J = ------
207	2
208	2 2 2 2 2 2
209	eps = (A22 + A21 ) x + 2(-A22A12-A21A11) x*px + (A11 + A12 ) px
210	gamma alpha beta
211
212	Input:
213	A matrix
214
215	Output:
216	alpha vector
217	beta vector
218
219	Return:
220	none
221
222	Global variables:
223	none
224
225	Specific functions:
226	none
227
228	Comments:
229	none
230
231	****************************************************************************/
232	#define n 4
233	void getprm(Matrix &Ascr, Vector2 &alpha, Vector2 &beta)
234	{
235	int i, j;
236
237	for (i = 1; i <= 2; i++) {
238	j = i*2 - 1;
239	alpha[i-1] = -(Ascr[j-1][j-1]Ascr[j][j-1] + Ascr[j-1][j]Ascr[j][j]);
240	beta[i-1] = sqr(Ascr[j-1][j-1]) + sqr(Ascr[j-1][j]);
241	}
242	}
243
244	/****************************************************************************/
245	/* void Cell_Twiss_M(long i0, long i1, Matrix &Ascr, bool chroma,
246	bool ring, double dP)
247
248	Purpose: called by Ring_Twiss_M
249	Calculate Twiss parameters from element i0 to element i1
250	Method: extended matrix formalism
251
252	Input:
253	i0 first element
254	i1 last element
255	ring true if a ring
256	chroma true if compute chromaticities
257	dP energy offset
258
259	Output:
260	none
261
262	Return:
263	none
264
265	Global variables:
266	none
267
268	Specific functions:
269	getprm,
270	CopyMat, CopyVec
271	Elem_Pass_M, GetAngle
272
273	Comments:
274	17/07/03 use of M_PI instead of pi
275
276	****************************************************************************/
277	void Cell_Twiss_M(long i0, long i1, Matrix &Ascr, bool chroma,
278	bool ring, double dP)
279	{
280	long int i;
281	int j, k;
282	Vector2 nu1, dnu;
283	Vector xref;
284	Matrix Ascr0, Ascr1;
285	CellType *cellp;
286
287	if (dP != globval.dPparticle) Cell_SetdP(dP);
288
289	/* Init */
290	for (j = 0; j <= 1; j++)
291	nu1[j] = 0.0;
292
293	/* get alpha beta for i0 */
294	cellp = &Cell[i0];
295	getprm(Ascr, cellp->Alpha, cellp->Beta);
296	memcpy(cellp->Nu, nu1, sizeof(Vector2));
297	CopyMat(n+1, Ascr, Ascr0); CopyVec(n+2L, globval.CODvect, xref);
298
299	for (i = i0+1; i <= i1; i++){
300	CopyMat(n+1, Ascr0, Ascr1);
301	/* Ascr1=Elem_MAscr0 /
302	/* xref =Elem(xref) */
303	Elem_Pass_M(i, xref, Ascr1);
304
305	cellp = &Cell[i];
306	/* get alpha and beta for element i */
307	getprm(Ascr1, cellp->Alpha, cellp->Beta);
308
309	for (j = 0; j <= 1; j++) {
310	k = (j+1) * 2 - 1;
311	/* get phase advances */
312	dnu[j] = (GetAngle(Ascr1[k-1][k-1], Ascr1[k-1][k]) -
313	GetAngle(Ascr0[k-1][k-1], Ascr0[k-1][k]))/(2.0*M_PI);
314	if ((cellp->Elem.PL >= 0.0) && (dnu[j] < -1e-16))
315	dnu[j] += 1.0;
316	else if ((cellp->Elem.PL < 0.0) && (dnu[j] > 1e-16))
317	dnu[j] -= 1.0;
318	nu1[j] += dnu[j]; cellp->Nu[j] = nu1[j];
319	/* Only correct for bare lattice */
320	cellp->Eta[j] = Ascr1[k-1][4]; cellp->Etap[j] = Ascr1[k][4];
321	}
322	CopyMat(n+1, Ascr1, Ascr0);
323	}
324	if (chroma) Cell_Geteta(i0, i1, ring, dP);
325	}
326	#undef n
327
328	/****************************************************************************/
329	/* void dagetprm(ss_vect<tps> &Ascr, Vector2 &alpha, Vector2 &beta)
330
331	Purpose: called by Cell_Twiss
332	Compute Twiss parameters alpha and beta from the matrix Ascr
333
334	M oneturn matrix (A and M are symplectic)
335
336	( A11 A12 ) -1 ( cos(mu) sin(mu) )
337	A = ( ) M = (A R A) with R = ( )
338	( A21 A22 ) ( -sin(mu) cos(mu) )
339
340	2 2
341	-1 x + px
342	eps = 2*J o a J = ------
343	2
344	2 2 2 2 2 2
345	eps = (A22 + A21 ) x + 2(-A22A12-A21A11) x*px + (A11 + A12 ) px
346	gamma alpha beta
347
348	A = Asrc
349
350	alpha(i-1) = -(A(2i-1,2i-1)A(2i,2i-1) + A(2i-1,2i)A(2i,2i))
351	beta(i-1) = A(2i-1,2i-1)A(2i-1,2i-1) + A(2i-1,2i)A(2i-1,2i)
352
353	Input:
354	Ascr
355
356	Output:
357	alpha alpha function vector
358	beta beta function vector
359
360	Return:
361	none
362
363	Global variables:
364	none
365
366	Specific functions:
367	getmat
368
369	Comments:
370	none
371
372	****************************************************************************/
373	#define n 4
374	void dagetprm(ss_vect<tps> &Ascr, Vector2 &alpha, Vector2 &beta)
375	{
376	int i = 0, j = 0;
377
378	for (i = 1; i <= 2; i++) {
379	j = i*2 - 1;
380	alpha[i-1] = -(getmat(Ascr, j, j)*getmat(Ascr, j+1, j)
381	+ getmat(Ascr, j, j+1)*getmat(Ascr, j+1, j+1));
382	beta[i-1] = sqr(getmat(Ascr, j, j)) + sqr(getmat(Ascr, j, j+1));
383	}
384	}
385
386	/******************************************************************
387	void Cell_Twiss(long i0, long i1, ss_vect<tps> &Ascr, bool chroma, bool ring,
388	double dP)
389
390
391	Purpose:
392	Get the phase, chromaticites, dispersion functions.
393
394	input:
395	i0: start index of the lattice element
396	i1: end index of the lattice element
397	Ascr: normalized one turn transfer matrix
398	chroma: bool flag to calculate chromaticites
399	ring: for ring lattice or not
400	dP: off momentum
401	********************************************************************/
402	void Cell_Twiss(long i0, long i1, ss_vect<tps> &Ascr, bool chroma, bool ring,
403	double dP)
404	{
405	long int i = 0;
406	int j = 0, k = 0;
407	Vector2 nu1, dnu; /* absolute and relative phase advance */
408	ss_vect<tps> Ascr0, Ascr1;
409	CellType *cellp;
410
411	/* initialization */
412	for (j = 0; j <= 1; j++)
413	{
414	nu1[j] = 0.0;
415	dnu[j] = 0.0;
416	}
417
418	if (globval.radiation) globval.dE = 0.0;
419
420	cellp = &Cell[i0];
421	dagetprm(Ascr, cellp->Alpha, cellp->Beta);
422	memcpy(cellp->Nu, nu1, sizeof(Vector2));
423
424	Ascr0 = Ascr;
425	for (j = 0; j <= n+1; j++)
426	Ascr0[j] += tps(globval.CODvect[j]);
427
428	Ascr1 = Ascr0;
429	for (i = i0; i <= i1; i++)
430	{
431	Elem_Pass(i, Ascr1);
432	cellp = &Cell[i];
433	dagetprm(Ascr1, cellp->Alpha, cellp->Beta);
434	for (j = 0; j <= 1; j++) {
435	k = (j+1)*2 - 1;
436	dnu[j] = (GetAngle(getmat(Ascr1, k, k), getmat(Ascr1, k, k+1)) -
437	GetAngle(getmat(Ascr0, k, k), getmat(Ascr0, k, k+1)))/(2.0*M_PI);
438
439	if ((cellp->Elem.PL >= 0.0) && (dnu[j] < -1e-16))
440	dnu[j] += 1.0;
441	else if ((cellp->Elem.PL < 0.0) && (dnu[j] > 1e-16))
442	dnu[j] -= 1.0;
443
444	nu1[j] += dnu[j];
445
446	cellp->Nu[j] = nu1[j];
447	cellp->Eta[j] = getmat(Ascr1, k, 5)*getmat(Ascr1, 6, 6) -
448	getmat(Ascr1, k, 6)*getmat(Ascr1, 6, 5);
449	cellp->Etap[j] = getmat(Ascr1, k+1, 5)*getmat(Ascr1, 6, 6) -
450	getmat(Ascr1, k+1, 6)*getmat(Ascr1, 6, 5);
451	}
452	Ascr0 = Ascr1;
453	}
454
455	if (chroma && !globval.Cavity_on) Cell_Geteta(i0, i1, ring, dP);
456	}
457
458	#undef n
459
460
461	/****************************************************************************/
462	/* void Ring_Getchrom(double dP)
463
464	Purpose: called by Ring_Twiss_M and Ring_Twiss
465	Compute chromaticites of the ring by numerical differentiation
466	and by evaluating each time the closed orbit
467
468	nu(dP+dPlocal)-nu(dP-dPlocal)
469	xi =-----------------------------
470	2 dPlocal
471
472	Input:
473	dP particle energy offset (should be zero cf comments)
474
475	Output:
476	none
477
478	Return:
479	none
480
481	Global variables:
482	status
483	globval
484
485	Specific functions:
486	Cell_GetABGN, Cell_GetCOD
487
488	Comments:
489	03/01/03 Stability test from GETcod routines slightly modified
490	WARNING : this routine does not give the chromaticities for dP != 0
491	but the local slope of the curve xi=f(dP)
492	16/10/03 Modified convergence test: now done for both planes
493	01/09/10 Modify the convergence criteria on relative diff of the chroma value
494	The previous test does not work well for non zero chromaticities
495	Test convergence changed if chroma closed to zero
496	28/06/11 Fix the bug for the lattice with negative momentum compact factor,
497	by changing expo += 1 to expo += log(5)/LOG10.
498	By Jianfeng Zhang @ SOLEIL
499
500	****************************************************************************/
501	#define n 4
502	#define chromeps 1e-6 /* convergence condition for chromaticity computation */
503	#define LOG10 log(10.0)
504	#define ZEROCHRO 0.1
505
506	void Ring_Getchrom(double dP) {
507	long int lastpos = 0;
508	double NORMchro = 1.0;
509	double dPlocal = 0.0, expo = 0.0, Norm_TEMP = 0.0;
510	int j;
511	Vector2 alpha = { 0.0, 0.0 }, beta = { 0.0, 0.0 }, gamma = { 0.0, 0.0 };
512	Vector2 nu = { 0.0, 0.0 }, nu0 = { 0.0, 0.0 }, TEMP={0.0,0.0}, Chrom={0.0,0.0};;
513
514	if (dP != 0.0)
515	fprintf(stdout,
516	"Ring_Getchrom: Warning this is NOT the CHROMA, dP=%e\n", dP);
517
518	/* initialization */
519	globval.Chrom[0] = 1e38;
520	globval.Chrom[1] = 1e38;
521
522	expo = log(globval.dPcommon) / LOG10;
523	do {
524	for (j = 0; j <= 1; j++)
525	Chrom[j] = globval.Chrom[j];
526
527	dPlocal = exp(expo * LOG10);
528	//Get cod for energy dP - globval.dPcommon
529	GetCOD(globval.CODimax, globval.CODeps, dP - dPlocal * 0.5, lastpos);
530
531	if (!status.codflag) {
532	// if no cod
533	fprintf(stdout, "Ring_Getchrom: Lattice is unstable for"
534	" dP-globval.dPcommon=% .5e\n", dP - dPlocal * 0.5);
535	return;
536	}
537
538	//get tunes for energy dP - globval.dPcommon/2 from oneturn matrix
539	Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu0);
540
541	//Get cod for energy dP+globval.dPcommon
542	GetCOD(globval.CODimax, globval.CODeps, dP + dPlocal * 0.5, lastpos);
543
544	if (!status.codflag) { // if no cod
545	fprintf(stdout, "Ring_Getchrom Lattice is unstable for"
546	" dP+globval.dPcommon=% .5e \n", dP + dPlocal * 0.5);
547	return;
548	}
549
550	//get tunes for energy dP+globval.dPcommon/2 from one turn matrix
551	Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
552
553	if (!globval.stable) {
554	printf("Ring_Getchrom: Lattice is unstable\n");
555	}
556
557	//Get chromaticities by numerical differentiation
558	for (j = 0; j <= 1; j++)
559	globval.Chrom[j] = (nu[j] - nu0[j]) / globval.dPcommon;
560
561	/* Get chromaticities by numerical differentiation*/
562	for (j = 0; j <= 1; j++) {
563	globval.Chrom[j] = (nu[j] - nu0[j]) / dPlocal;
564	TEMP[j] = fabs(globval.Chrom[j] - Chrom[j]);
565	}
566
567	Norm_TEMP = sqrt(TEMP[0] * TEMP[0] + TEMP[1] * TEMP[1]);
568	NORMchro = sqrt(globval.Chrom[0] * globval.Chrom[0] + globval.Chrom[1]
569	* globval.Chrom[1]);
570	// if chroma closed to zero, norm is one for avoiding divergence of convergence test
571	if (NORMchro < ZEROCHRO){
572	NORMchro = 1.0;
573	}
574
575	// TEST CHROMA convergence
576	if (trace) {
577	fprintf(
578	stdout,
579	"\nexpo % e xix = % e xiz = % e Norm_TEMP = %e TEMPX %+e TEMPZ %+e\n",
580	expo, Chrom[0], Chrom[1], Norm_TEMP, TEMP[0], TEMP[1]);
581	fprintf(stdout, "expo % e nux = % e nuz = % e dPlocal= %+e\n",
582	expo, nu0[0], nu0[1], dP - 0.5 * dPlocal);
583	fprintf(stdout, "expo % e nux = % e nuz = % e dPlocal= %+e\n",
584	expo, nu[0], nu[1], dP + 0.5 * dPlocal);
585	}
586	expo += log(5)/LOG10;
587	if (trace)
588	fprintf(stdout, "%+e %+.12e %+.12e %+.12e %+.12e\n", dPlocal,
589	globval.Chrom[0], fabs(globval.Chrom[0] - Chrom[0])
590	/ Chrom[0], globval.Chrom[1], fabs(globval.Chrom[1]
591	- Chrom[1]) / Chrom[1]);
592	} while ((Norm_TEMP > chromeps * NORMchro) \|\| (Norm_TEMP == 0e0));
593
594	status.chromflag = true;
595	}
596
597	#undef n
598	#undef chromeps
599	#undef LOG10
600	#undef ZEROCHRO
601
602
603	/****************************************************************************/
604	/* static void Ring_Twiss_M(bool chroma, double dP)
605
606	Purpose: called by Ring_GetTwiss
607	Calcule les parametres Twiss a l'energie dP
608	si chroma=true, les chromaticites sont calculees
609	Methode matricielle
610
611	M : One turn transfer matrix
612
613	[ cx sx . . ]
614	[ ]
615	-1 [-sx cx . . ]
616	M -> A M A = R = [ ]
617	[ . . cy sy]
618	[ ]
619	[ . . -sy cy]
620
621
622
623	Input:
624	bool true if chromaticities and dispersion to compute
625	else false
626	dP particle energy offset
627
628	Output:
629	none
630
631	Return:
632	none
633
634	Global variables:
635	none
636
637	Specific functions:
638	none
639
640	Comments:
641	none
642
643	****************************************************************************/
644	#define n 4
645	void Ring_Twiss_M(bool chroma, double dP)
646	{
647	long int lastpos = 0;
648	int j;
649	Vector2 alpha={0.0,0.0}, beta={0.0,0.0}, gamma={0.0,0.0}, nu={0.0,0.0};
650	Vector eta0;
651	Matrix R, C, Ascr;
652
653	/* Get closed orbit and compute oneturn matrix */
654	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos);
655
656	if (!status.codflag) /* Check if stable */
657	return;
658	/* compute twiss parameter using oneturn matrix */
659	Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
660
661	/* Get eigenvalues and eigenvectors for the one turn transfer matrix */
662	GDiag((long)n, Cell[globval.Cell_nLoc].S, globval.Ascr, globval.Ascrinv, R,
663	globval.OneTurnMat, globval.Omega, globval.Alphac);
664
665	/* Only correct for bare lattice */
666	for (j = 0; j <= n; j++)
667	eta0[j] = globval.OneTurnMat[j][n];
668
669	UnitMat((long)n, C); SubMat((long)n, globval.OneTurnMat, C);
670
671	if (!InvMat((long)n, C))
672	printf("** matrix is singular\n");
673
674	LinTrans((long)n, C, eta0);
675
676	for (j = 0; j <= n; j++) {
677	globval.Ascr[n][j] = 0.0;
678	globval.Ascr[j][n] = eta0[j];
679	}
680
681	CopyMat(n+1, globval.Ascr, Ascr);
682	Cell_Twiss_M(0, globval.Cell_nLoc, Ascr, chroma, true, dP);
683
684	/* Copies tunes into global variable BUG !!!*/
685	memcpy(globval.TotalTune, Cell[globval.Cell_nLoc].Nu, sizeof(Vector2));
686	status.tuneflag = true;
687
688	if (chroma)
689	{ /* compute chromaticities */
690	Ring_Getchrom(dP);
691	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos);
692	}
693	}
694
695	#undef n
696
697
698	/****************************************************************************/
699	/* void Ring_Twiss(bool chroma, double dP)
700
701	Purpose: called by Ring_GetTwiss
702	Computes twiss parameters for an energy offset dP using da method
703	First call GetCOD() to check whether there is a stable COD, if yes,
704	check the linear one turn matrix is stable or not, then put the linear
705	one turn map
706
707	Input:
708	chroma if true computes chromaticities as well
709	dP energy offset
710
711	Output:
712	none
713
714	Return:
715	none
716
717	Global variables:
718	globval
719	status
720
721	Specific functions:
722	none
723
724	Comments:
725
726
727	****************************************************************************/
728	void Ring_Twiss(bool chroma, double dP)
729	{
730	long int lastpos = 0;
731	int n = 0;
732	Vector2 alpha={0.0, 0.0}, beta={0.0, 0.0};
733	Vector2 gamma={0.0, 0.0}, nu={0.0, 0.0};
734	Matrix R;
735	ss_vect<tps> AScr; /* DA map*/
736
737	if (globval.Cavity_on)
738	n = 6L; /* 6D tracking */
739	else /* 4D tracking */
740	n = 4L;
741
742	if(!trace)
743	printf("\n %d D tracking....\n",n);
744
745	/* Gets closed orbit and computes one turn map around it*/
746	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos);
747
748	if (!status.codflag)
749	return;
750
751	// Check if stable
752	Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
753	// Get eigenvalues and eigenvectors for the one turn transfer matrix
754	GDiag(n, Cell[globval.Cell_nLoc].S, globval.Ascr, globval.Ascrinv, R,
755	globval.OneTurnMat, globval.Omega, globval.Alphac);
756	// Puts zeroes in constant part of da map
757	putlinmat(n, globval.Ascr, AScr);
758
759	if (!globval.Cavity_on)
760	{
761	AScr[delta_] = 0.0;
762	AScr[ct_] = 0.0;
763	}
764
765	Cell_Twiss(0, globval.Cell_nLoc, AScr, chroma, true, dP);
766
767	memcpy(globval.TotalTune, Cell[globval.Cell_nLoc].Nu, sizeof(Vector2));
768	status.tuneflag = true;
769
770	if (chroma && !globval.Cavity_on) {
771	Ring_Getchrom(dP);
772	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos);
773	}
774	}
775
776	/****************************************************************************/
777	/* void Ring_GetTwiss(bool chroma, double dP)
778
779	Purpose:
780	Computes Twiss functions w/ or w/o chromaticities
781	for particle of energy offset dP
782	matrix or DA method
783
784	Also get the momentum compact factor based on the definition:
785	Alphac = delta_C/delta*1/C
786
787	Input:
788	Chroma if true, compute chromaticities and dispersion
789	if false, dispersion is set to zero !!!
790	Dp energy offset
791
792	Output:
793	none
794
795	Return:
796	none
797
798	Global variables:
799	globval, trace
800
801	Specific functions:
802	Ring_Twiss_M, Ring_Twiss
803
804	Comments:
805	none
806
807	****************************************************************************/
808	void Ring_GetTwiss(bool chroma, double dP)
809	{
810
811	if (trace) printf("enter ring_gettwiss\n");
812
813	if (globval.MatMeth) /* matrix method */
814	Ring_Twiss_M(chroma, dP);
815	else /* da method */
816	Ring_Twiss(chroma, dP);
817
818	globval.Alphac = globval.OneTurnMat[ct_][delta_]/Cell[globval.Cell_nLoc].S;
819	if (trace) printf("exit ring_gettwiss\n");
820	}
821
822
823	/* Local variables for Ring_Fittune: */
824	struct LOC_Ring_Fittune
825	{
826	jmp_buf _JL999;
827	};
828
829
830	/****************************************************************************/
831	/* void TransTrace(long i0, long i1, Vector2 &alpha, Vector2 &beta, Vector2 &eta,
832	Vector2 &etap, double *codvect)
833
834	Purpose: TransTwiss
835	Compute Twiss fonction for a tranfert line from element i0 to element i1
836
837	Input:
838	alpha alpha fonctions at the line entrance
839	beta beta fonctions at the line entrance
840	eta disperion fonctions at the line entrance
841	etap dipersion derivatives fonctions at the line entrance
842	codvect closed orbit fonctions at the line entrance
843
844	Output:
845	none
846
847	Return:
848	none
849
850	Global variables:
851	globval, trace
852
853	Specific functions:
854	Cell_Twiss_M, Cell_Pass_M
855
856	Comments:
857	28/10/03 phase advances added
858
859	****************************************************************************/
860	void TransTrace(long i0, long i1, Vector2 &alpha, Vector2 &beta, Vector2 &eta,
861	Vector2 &etap, Vector &codvect)
862	{
863	long i, j, lastpos;
864	double sb;
865	Matrix Ascr;
866
867	UnitMat(6L, Ascr);
868	for (i = 1; i <= 2L; i++) {
869	sb = sqrt(beta[i-1]); j = i*2 - 1;
870	Ascr[j-1][j-1] = sb; Ascr[j-1][j] = 0.0;
871	Ascr[j][j-1] = -(alpha[i-1]/sb); Ascr[j][j] = 1/sb;
872	}
873	Ascr[0][4] = eta[0]; Ascr[1][4] = etap[0];
874	Ascr[2][4] = eta[1]; Ascr[3][4] = etap[1];
875
876	for (i = 0; i < ss_dim; i++) {
877	globval.CODvect[i] = codvect[i];
878	}
879
880	// get twiss parameter for each element
881	Cell_Twiss_M(i0, i1, Ascr, false, false, codvect[4]);
882	memcpy(globval.TotalTune, Cell[globval.Cell_nLoc].Nu, sizeof(Vector2));
883	status.tuneflag = true;
884
885	for (i = 0; i < ss_dim; i++) {
886	globval.CODvect[i] = codvect[i];
887	}
888
889	// compute oneturn matrix
890	UnitMat(5L, globval.OneTurnMat);
891	Cell_Pass_M(0, globval.Cell_nLoc, globval.CODvect, globval.OneTurnMat,
892	lastpos);
893	}
894
895	/****************************************************************************/
896	/* void shiftk(long Elnum, double dk, struct LOC_Ring_Fittune *LINK)
897
898	Purpose:
899
900
901	Input:
902	none
903
904	Output:
905	none
906
907	Return:
908	none
909
910	Global variables:
911	none
912
913	Specific functions:
914	none
915
916	Comments:
917	none
918
919	****************************************************************************/
920	#define dP 0.0
921	void shiftk(long Elnum, double dk, struct LOC_Ring_Fittune *LINK)
922	{
923	CellType *cellp;
924	elemtype *elemp;
925	MpoleType *M;
926
927	cellp = &Cell[Elnum]; elemp = &cellp->Elem; M = elemp->M;
928	M->PBpar[Quad+HOMmax] += dk;
929	Mpole_SetPB(cellp->Fnum, cellp->Knum, (long)Quad);
930	}
931
932
933	/****************************************************************************/
934	/* void checkifstable(struct LOC_Ring_Fittune *LINK)
935
936	Purpose:
937
938
939	Input:
940	none
941
942	Output:
943	none
944
945	Return:
946	none
947
948	Global variables:
949	none
950
951	Specific functions:
952	none
953
954	Comments:
955	none
956
957	****************************************************************************/
958	void checkifstable(struct LOC_Ring_Fittune *LINK)
959	{
960	if (!globval.stable) {
961	printf(" lattice is unstable\n");
962	longjmp(LINK->_JL999, 1);
963	}
964	}
965
966
967	/****************************************************************************/
968	/* void Ring_Fittune(Vector2 &nu, double eps, long nq, long qf, long *qd,
969	double dkL, long imax)
970
971	Purpose: called by fittune
972	Fit tunes using two family of quadrupoles
973	Linear method
974
975	Input:
976	nu target tunes
977	eps precision
978	nq number of quad of family qf and qd
979	qf position of qf magnets
980	qd position of qd magnets
981	dKL variation on strengths
982	imax maximum number of iteration
983
984	Output:
985	none
986
987	Return:
988	none
989
990	Global variables:
991	none
992
993	Specific functions:
994	none
995
996	Comments:
997	17 mars 2004: removed labels
998
999	****************************************************************************/
1000	void Ring_Fittune(Vector2 &nu, double eps, iVector2 &nq, long qf[], long qd[],
1001	double dkL, long imax)
1002	{
1003	struct LOC_Ring_Fittune V;
1004
1005	int i, j, k;
1006	Vector2 nu0, nu1;
1007	Vector dkL1, dnu;
1008	Matrix A;
1009	bool prt = false;
1010
1011	if (setjmp(V._JL999)) return;
1012
1013	if (trace)
1014	printf(" Tune fit, nux =%10.5f, nuy =%10.5f, eps =% .3E,"
1015	" imax =%4ld, dkL = % .5E\n", nu[0], nu[1], eps, imax, dkL);
1016	Ring_GetTwiss(false, dP); checkifstable(&V);
1017	memcpy(nu0, globval.TotalTune, sizeof(Vector2));
1018	i = 0;
1019	do {
1020	i++;
1021	/* First vary kf then kd */
1022	for (j = 1; j <= 2L; j++) {
1023	for (k = 0; k < nq[j-1]; k++) {
1024	if (j == 1)
1025	shiftk(qf[k], dkL, &V); // new value for qf
1026	else
1027	shiftk(qd[k], dkL, &V); // new value for qd
1028	}
1029	Ring_GetTwiss(false, dP);
1030	nu1[0] = globval.TotalTune[0]; nu1[1] = globval.TotalTune[1];
1031	checkifstable(&V);
1032	for (k = 0; k <= 1; k++) {
1033	dnu[k] = nu1[k] - (long)nu1[k] - nu0[k] + (long)nu0[k];
1034	if (fabs(dnu[k]) > 0.5) dnu[k] = 1 - fabs(dnu[k]);
1035	A[k][j-1] = dnu[k]/dkL;
1036	}
1037
1038	/* Restore strength */
1039	for (k = 0; k < nq[j-1]; k++) {
1040	if (j == 1)
1041	shiftk(qf[k], -dkL, &V);
1042	else
1043	shiftk(qd[k], -dkL, &V);
1044	}
1045	}
1046
1047	if (!InvMat(2L, A)) {
1048	printf(" A is singular\n");
1049	return;
1050	}
1051
1052	for (j = 0; j <= 1; j++)
1053	dkL1[j] = nu[j] - nu0[j];
1054	LinTrans(2L, A, dkL1);
1055	for (j = 1; j <= 2; j++) {
1056	for (k = 0; k < nq[j-1]; k++) {
1057	if (j == 1)
1058	shiftk(qf[k], dkL1[j-1], &V);
1059	else
1060	shiftk(qd[k], dkL1[j-1], &V);
1061	}
1062	}
1063	Ring_GetTwiss(false, dP); checkifstable(&V);
1064	memcpy(nu0, globval.TotalTune, sizeof(Vector2));
1065	if (trace)
1066	printf(" Nux = %10.6f(%10.6f), Nuy = %10.6f(%10.6f),"
1067	" QFam1L = % .5E, QFam2L = % .5E @%3d\n",
1068	nu0[0], nu1[0], nu0[1], nu1[1],
1069	Elem_GetKval(Cell[qf[0]].Fnum, 1, (long)Quad),
1070	Elem_GetKval(Cell[qd[0]].Fnum, 1, (long)Quad), i);
1071	} while (sqrt(sqr(nu[0]-nu0[0])+sqr(nu[1]-nu0[1])) >= eps && i != imax);
1072	if (prt)
1073	// print new K-value for 1s element of the family
1074	printf("Reached Nux = %10.6f, Nuy = %10.6f,"
1075	" QFam1L = % .5E, QFam2L = % .5E @iteration %3d\n",
1076	nu[0], nu[1],
1077	Elem_GetKval(Cell[qf[0]].Fnum, 1, (long)Quad),
1078	Elem_GetKval(Cell[qd[0]].Fnum, 1, (long)Quad), i);
1079	}
1080	#undef dP
1081
1082
1083	#define dP 0.0
1084	void shiftkp(long Elnum, double dkp)
1085	{
1086	CellType *cellp;
1087	elemtype *elemp;
1088	MpoleType *M;
1089
1090	cellp = &Cell[Elnum]; elemp = &cellp->Elem; M = elemp->M;
1091	M->PBpar[Sext+HOMmax] += dkp;
1092	Mpole_SetPB(cellp->Fnum, cellp->Knum, (long)Sext);
1093	}
1094
1095
1096	void Ring_Fitchrom(Vector2 &ksi, double eps, iVector2 &ns,
1097	long sf[], long sd[], double dkpL, long imax)
1098	{
1099	bool rad;
1100	long int lastpos;
1101	int i, j, k;
1102	Vector2 ksi0;
1103	Vector dkpL1, dksi;
1104	Matrix A;
1105
1106	if (trace)
1107	printf(" Chromaticity fit, ksix =%10.5f, ksiy =%10.5f, eps =% .3E"
1108	", imax =%4ld, dkpL =%10.5f\n", ksi[0], ksi[1], eps, imax, dkpL);
1109
1110	/* Turn off radiation */
1111	rad = globval.radiation; globval.radiation = false;
1112	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos); Ring_Getchrom(dP);
1113	for (j = 0; j <= 1; j++)
1114	ksi0[j] = globval.Chrom[j];
1115	i = 0;
1116	do {
1117	i++;
1118	/* First vary sf then sd */
1119	for (j = 1; j <= 2; j++) {
1120	for (k = 0; k < ns[j-1]; k++) {
1121	if (j == 1)
1122	shiftkp(sf[k], dkpL);
1123	else
1124	shiftkp(sd[k], dkpL);
1125	}
1126	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos); Ring_Getchrom(dP);
1127	for (k = 0; k <= 1; k++) {
1128	dksi[k] = globval.Chrom[k] - ksi0[k];
1129	A[k][j-1] = dksi[k] / dkpL;
1130	}
1131	/* Restore strength */
1132	for (k = 0; k < ns[j-1]; k++) {
1133	if (j == 1)
1134	shiftkp(sf[k], -dkpL);
1135	else
1136	shiftkp(sd[k], -dkpL);
1137	}
1138	}
1139	if (!InvMat(2L, A)) {
1140	printf(" A is singular\n");
1141	goto _L999;
1142	}
1143	for (j = 0; j <= 1; j++)
1144	dkpL1[j] = ksi[j] - ksi0[j];
1145	LinTrans(2L, A, dkpL1);
1146	for (j = 1; j <= 2; j++) {
1147	for (k = 0; k < ns[j-1]; k++) {
1148	if (j == 1)
1149	shiftkp(sf[k], dkpL1[j-1]);
1150	else
1151	shiftkp(sd[k], dkpL1[j-1]);
1152	}
1153	}
1154	GetCOD(globval.CODimax, globval.CODeps, dP, lastpos); Ring_Getchrom(dP);
1155	for (j = 0; j <= 1; j++)
1156	ksi0[j] = globval.Chrom[j];
1157	if (trace)
1158	printf(" ksix =%10.6f, ksiy =%10.6f, SF = % .5E, SD = % .5E @%3d\n",
1159	ksi0[0], ksi0[1], Elem_GetKval(Cell[sf[0]].Fnum, 1, (long)Sext),
1160	Elem_GetKval(Cell[sd[0]].Fnum, 1, (long)Sext), i);
1161	} while (sqrt(sqr(ksi[0]-ksi0[0])+sqr(ksi[1]-ksi0[1])) >= eps && i != imax);
1162	_L999:
1163	/* Restore radiation */
1164	globval.radiation = rad;
1165	}
1166
1167	#undef dP
1168
1169
1170	#define dP 0.0
1171
1172
1173	/* Local variables for Ring_FitDisp: */
1174	struct LOC_Ring_FitDisp
1175	{
1176	jmp_buf _JL999;
1177	};
1178
1179
1180	static void shiftk_(long Elnum, double dk, struct LOC_Ring_FitDisp *LINK)
1181	{
1182	CellType *cellp;
1183	elemtype *elemp;
1184	MpoleType *M;
1185
1186	cellp = &Cell[Elnum];
1187	elemp = &cellp->Elem;
1188	M = elemp->M;
1189	M->PBpar[Quad+HOMmax] += dk;
1190	Mpole_SetPB(cellp->Fnum, cellp->Knum, (long)Quad);
1191	}
1192
1193
1194	void checkifstable_(struct LOC_Ring_FitDisp *LINK)
1195	{
1196	if (!globval.stable) {
1197	printf(" lattice is unstable\n");
1198	longjmp(LINK->_JL999, 1);
1199	}
1200	}
1201
1202
1203	void Ring_FitDisp(long pos, double eta, double eps, long nq, long q[],
1204	double dkL, long imax)
1205	{
1206	/*pos : integer; eta, eps : double;
1207	nq : integer; var q : fitvect;
1208	dkL : double; imax : integer*/
1209	struct LOC_Ring_FitDisp V;
1210
1211	int i, j;
1212	double dkL1, Eta0, deta;
1213	bool rad = false;
1214
1215	if (setjmp(V._JL999)) goto _L999;
1216
1217	if (trace)
1218	printf(" Dispersion fit, etax =%10.5f, eps =% .3E"
1219	", imax =%4ld, dkL =%10.5f\n",
1220	eta, eps, imax, dkL);
1221	/* Turn off radiation */
1222	rad = globval.radiation; globval.radiation = false;
1223	Ring_GetTwiss(true, dP); checkifstable_(&V);
1224	Eta0 = Cell[pos].Eta[0];
1225	i = 0;
1226	while (fabs(eta-Eta0) > eps && i < imax) {
1227	i++;
1228	for (j = 0; j < nq; j++)
1229	shiftk_(q[j], dkL, &V);
1230	Ring_GetTwiss(true, dP); checkifstable_(&V);
1231	deta = Cell[pos].Eta[0] - Eta0;
1232	if (deta == 0.0) {
1233	printf(" deta is 0\n");
1234	goto _L999;
1235	}
1236	dkL1 = (eta-Eta0)*dkL/deta - dkL;
1237	for (j = 0; j < nq; j++)
1238	shiftk_(q[j], dkL1, &V);
1239	Ring_GetTwiss(true, dP); checkifstable_(&V);
1240	Eta0 = Cell[pos].Eta[0];
1241	if (trace)
1242	printf(" Dispersion = % .5E, kL =% .5E @%3d\n",
1243	Eta0, Elem_GetKval(Cell[q[0]].Fnum, 1, (long)Quad), i);
1244	}
1245	_L999:
1246	/* Restore radiation */
1247	globval.radiation = rad;
1248	}
1249	#undef dP
1250
1251	/* End. */

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

Download in other formats: