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