source: TRACY3/trunk/tracy/tracy/src/t2cell.cc @ 3

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

Initiale import
  • Property svn:executable set to *
File size: 25.1 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*/
10
11
12long    ntransfmat;
13Matrix  transfmat[maxtransfmat];
14long    kicks[maxtransfmat][maxkicks];
15
16/**************************************************************
17template<typename T>
18inline bool CheckAmpl(const ss_vect<T> &x, const long int loc)
19
20  Purpose:
21     Check whether particle coordinate x are outside of
22      the aperture limitation, then return the bool flag not_lost,
23      if true, then beam is lost.
24     
25  Return:
26     not_lost
27     true,    if beam is not lost
28     false,   if beam is lost
29     
30
31***************************************************************/
32template<typename T>
33inline bool CheckAmpl(const ss_vect<T> &x, const long int loc)
34{
35  bool  not_lost;
36
37  if (globval.Aperture_on)
38    not_lost = is_double<T>::cst(x[x_]) > Cell[loc].maxampl[X_][0] &&
39               is_double<T>::cst(x[x_]) < Cell[loc].maxampl[X_][1] && 
40               is_double<T>::cst(x[y_]) > Cell[loc].maxampl[Y_][0]&&
41               is_double<T>::cst(x[y_]) < Cell[loc].maxampl[Y_][1];
42  else
43    not_lost = is_double<T>::cst(x[x_]) > -max_ampl &&
44               is_double<T>::cst(x[x_]) < max_ampl &&
45               is_double<T>::cst(x[y_]) > -max_ampl&&
46               is_double<T>::cst(x[y_]) < max_ampl;
47
48  if (!not_lost) {
49    if (is_double<T>::cst(x[x_]) < Cell[loc].maxampl[X_][0] ||
50        is_double<T>::cst(x[x_]) > Cell[loc].maxampl[X_][1])
51      status.lossplane = 1;
52    else if (is_double<T>::cst(x[y_]) < Cell[loc].maxampl[Y_][0] ||
53             is_double<T>::cst(x[y_]) > Cell[loc].maxampl[Y_][1])
54      status.lossplane = 2;
55           
56    if (trace)
57      printf("CheckAmpl: Particle lost in plane %d at element:"
58             " %5ld s = %10.5f, x = %12.5e, z= %12.5e\n",
59             status.lossplane, loc, Cell[loc].S,
60             is_double<T>::cst(x[x_]), is_double<T>::cst(x[y_]));
61  }
62
63  return not_lost;
64}
65
66
67template<typename T>
68void Elem_Pass(const long i, ss_vect<T> &x)
69{
70
71  switch (Cell[i].Elem.Pkind) {
72    case drift:
73      Drift_Pass(Cell[i], x);
74      break;
75    case Mpole:
76      Mpole_Pass(Cell[i], x);
77      break;
78    case Wigl:
79      Wiggler_Pass(Cell[i], x);
80      break;
81    case FieldMap:
82      FieldMap_Pass(Cell[i], x);
83      break;
84    case Insertion:
85      Insertion_Pass(Cell[i], x);
86      break;
87    case Cavity:
88      Cav_Pass(Cell[i], x);
89      break;
90    case marker:
91      Marker_Pass(Cell[i], x);
92      break;
93    case Spreader:
94      break;
95    case Recombiner:
96      break;
97    case Solenoid:
98      Solenoid_Pass(Cell[i], x);
99      break;
100    default:
101      printf("Elem_Pass ** undefined type\n");
102      break;
103  }
104
105  if (globval.IBS) {
106//    trace = i == 34;
107    is_tps<T>::do_IBS(Cell[i].Elem.PL, x);
108  }
109
110  is_tps<T>::get_ps(x, Cell[i]);
111}
112
113
114void Elem_Pass_M(const long i, Vector &xref, Matrix &x)
115{
116  /* Purpose:
117       Transport vector xref through matrix x
118       xref = Mi(xref)
119       x    = M*x   M: transport matrix of element i                         */
120
121  switch (Cell[i].Elem.Pkind) {
122    case drift:
123      Drift_Pass_M(Cell[i], xref, x);
124      break;
125    case Mpole:
126      Mpole_Pass_M(Cell[i], xref, x);
127      break;
128    case Wigl:
129      Wiggler_Pass_M(Cell[i], xref, x);
130      break;
131    case Insertion:
132      Insertion_Pass_M(Cell[i], xref, x);
133      break;
134    case Cavity:   /* nothing */
135      break;
136    case marker:   /* nothing */
137      break;
138    default:
139      fprintf(stdout,"Elem_Pass_M: ** undefined type\n");
140      break;
141  }
142
143  Cell[i].BeamPos = xref;
144}
145
146
147void Cell_SetdP(const double dP)
148{
149  int          i, j;
150  ElemFamType  *elemfamp;
151  elemtype     *elemp;
152
153  globval.dPparticle = dP;
154 
155  for (i = 1; i <= globval.Elem_nFam; i++) {
156    elemfamp  = &ElemFam[i-1]; elemp = &elemfamp->ElemF;
157    switch (elemp->Pkind) {
158    case drift:
159      for (j = 1; j <= elemfamp->nKid; j++)
160        Drift_SetMatrix(i, j);
161      break;
162    case Mpole:
163      for (j = 1; j <= elemfamp->nKid; j++)
164        Mpole_SetPB(i, j, 2);
165      break;
166    case Insertion:
167      for (j = 1; j <= elemfamp->nKid; j++)
168        Insertion_SetMatrix(i, j);
169      break;
170    case Wigl:
171      for (j = 1; j <= elemfamp->nKid; j++)
172        Wiggler_SetPB(i, j, 2);
173      break;
174    case FieldMap:
175      break;
176    case Cavity:
177      break;
178    case marker:
179      break;
180    case Spreader:
181      break;
182    case Recombiner:
183      break;
184    case Solenoid:
185      break;
186    default:
187      printf("** Cell_SetdP: undefined type\n");
188      break;
189    }
190  }
191}
192
193/****************************************************************************/
194/* template<typename T>
195void Cell_Pass(const long i0, const long i1, ss_vect<T> &x, long &lastpos)
196
197   Purpose:
198       Called by Cell_GetCOD().
199       First check the initial x at position i0 is stable or not, if yes,
200       then call Elem_Pass() to compute DA map from element position i0 to i1
201       for the particle with initial coordinates x, and check whether the particle is
202       lost at the aperture located at this element, finally return the value of lastpos.
203
204   Input:
205       i0             starting position
206       i1             ending position of the element in the lattice
207       x              initial coordinates
208       lastpos        position of the last element when the particle is stable 
209
210   Output:
211       map     DA map
212       lastpos pointer to last position of the particle
213
214   Return:
215           x     position of the particle
216      lastpos   
217                =position of the last element, particle is not lost
218               !=position of the last element, particle is lost
219   Global variables:
220       globval, status
221
222   Specific functions:
223       CheckAmpl
224       Elem_Pass
225
226   Comments:
227       29/12/02 remove label replaced by return
228
229****************************************************************************/
230template<typename T>
231void Cell_Pass(const long i0, const long i1, ss_vect<T> &x, long &lastpos)
232{
233  long int  i = 0;
234 
235  if (globval.MatMeth && (x[delta_] != globval.dPparticle))
236    Cell_SetdP(is_double<T>::cst(x[delta_]));
237   
238  if (globval.radiation) 
239    globval.dE = 0.0;
240
241  if (globval.emittance) 
242    {
243      I2 = 0.0; 
244      I4 = 0.0; 
245      I5 = 0.0;
246      for (i = 0; i < DOF; i++)
247        globval.D_rad[i] = 0.0;
248    }
249
250  if (globval.IBS)
251    for (i = 0; i < DOF; i++)
252      globval.D_IBS[i] = 0.0;
253   
254  /* tracking and check particle is lost or not*/             
255  if (!CheckAmpl(x, i0))
256    lastpos = i0;
257  else 
258    {
259      lastpos = i1;
260     for (i = i0; i <= i1; i++) 
261    //    for (i = i0+1L; i <= i1; i++)
262      {   
263        Elem_Pass(i, x);
264          if (!CheckAmpl(x, i)) 
265             {  lastpos = i;  break; }
266      }
267    }
268}
269
270/**************************************************************************/
271/* void Cell_Pass(const long i0, const long i1, tps &sigma, long &lastpos)
272
273   Purpose:
274       compute DA map from position i0 to i1 with sigma matrix
275       Called by Cell_GetCOD()
276
277   Input:
278       i0  starting position
279       i1 ending position
280
281   Output:
282       map     DA map
283       lastpos pointer to last position
284
285   Return:
286       none
287
288   Global variables:
289       globval, status
290
291   Specific functions:
292       CheckAmpl
293       Elem_Pass
294
295   Comments:
296       Specific for tracy 3
297
298****************************************************************************/
299void Cell_Pass(const long i0, const long i1, tps &sigma, long &lastpos)
300{
301  const int  n = 9;
302
303  int           i, j, jj[n][nv_tps];
304  ss_vect<tps>  Id, A;
305
306  const double  deps = 1e-20;
307
308  Id.identity();
309
310  map = Id + globval.CODvect; Cell_Pass(0, i0, map, lastpos);
311
312  if (lastpos == i0) {
313    map = Id + map.cst(); Cell_Pass(i0, i1, map, lastpos);
314
315    if (lastpos == i1) {
316      // x_1 = zeta(x_0) => f_1(x) = f_0(zeta^-1(x))
317
318      // deterministic part
319      sigma = sigma*Inv(map-map.cst());
320
321      if (globval.emittance) {
322        // stochastic part
323
324        for (i = 0; i < n; i++)
325          for (j = 0; j < nv_tps; j++)
326            jj[i][j] = 0;
327
328        jj[0][x_]  = 2; jj[1][x_]  = 1; jj[1][px_]    = 1; jj[2][px_]    = 2;
329        jj[3][y_]  = 2; jj[4][y_]  = 1; jj[4][py_]    = 1; jj[5][py_]    = 2;
330        jj[6][ct_] = 2; jj[7][ct_] = 1; jj[7][delta_] = 1; jj[8][delta_] = 2;
331
332        putlinmat(6, globval.Ascr, A); sigma = sigma*A;
333
334        for (i = 0; i < 3; i++) {
335          if (globval.eps[i] > deps) {
336            sigma.pook(jj[3*i], sigma[jj[3*i]]-globval.D_rad[i]/2.0);
337            sigma.pook(jj[3*i+2], sigma[jj[3*i+2]]-globval.D_rad[i]/2.0);
338          }
339        }
340
341        sigma = sigma*Inv(A);
342      }
343    }
344  }
345}
346
347
348void Cell_Pass_M(long i0, long i1, Vector &xref, Matrix &mat, long &lastpos)
349{
350  long i;
351
352  if (xref[4] != globval.dPparticle) Cell_SetdP(xref[4]);
353
354  if (!CheckAmpl(xref, i0))
355    lastpos = i0;
356  else {
357    Cell[i0].BeamPos = xref;
358    lastpos = i1;
359    for (i = i0+1; i <= i1; i++) {
360      Elem_Pass_M(i, xref, mat);
361      if (!CheckAmpl(xref, i)) { lastpos = i; break; }
362//      CopyVec(6, xref, Cell[i0].BeamPos);
363    }
364  }
365}
366
367
368#define n 4
369
370bool linearel(long i)
371{
372  /* Purpose: called by Cell_Concat
373       bool function
374         true if linear element
375           ie element is:
376             straight section
377             dipole w/s index (w/o skew quad)
378             normal quadrupole (w/o skew quad)
379             insertion (focalisation)
380             wiggler (focalisation)
381             RF cavity
382             marker
383         false otherwise                                                     */
384
385  bool     status = false;
386  CellType *cellp;
387  elemtype *elemp;
388
389  cellp = &Cell[i]; elemp = &cellp->Elem;
390
391  switch (elemp->Pkind) {
392  case drift: /* straight section */
393    status = true;
394    break;
395  case Mpole:
396    if (elemp->M->Pthick == thick && elemp->M->Porder <= Quad &&
397        elemp->M->PB[HOMmax-Quad] == 0e0)
398      status = true; /* normal quad */
399    else
400      status = false;
401    break;
402  case Wigl:
403    status = true;
404    break;
405  case FieldMap:
406    status = true;
407    break;
408  case Insertion:
409    status = true;
410    break;
411  case Cavity:
412    status = true;
413    break;
414  case marker:
415    status = true;
416    break;
417  case Spreader:
418    status = false;
419    break;
420  case Recombiner:
421    status = false;
422    break;
423  case Solenoid:
424    status = false;
425    break;
426  case undef:
427    break;
428  }
429 
430  return status;
431}
432
433void GtoL_dP(Matrix &mat, Vector2 &dT)
434{
435  long     k = 0;
436  Vector2  dS0;
437  Vector   x;
438
439  dS0[0] = 0.0; dS0[1] = 0.0;
440
441  for (k = 0; k < n; k++)
442    x[k] = mat[k][n];
443
444  GtoL(x, dS0, dT, 0.0, 0.0, 0.0);
445
446  for (k = 0; k < n; k++)
447    mat[k][n] = x[k];
448}
449
450
451static void LtoG_dP(Matrix &mat, Vector2 &dT)
452{
453  long     k;
454  Vector2  dS0;
455  Vector   x;
456
457  dS0[0] = 0.0; dS0[1] = 0.0;
458
459  for (k = 0; k < n; k++)
460    x[k] = mat[k][n];
461
462  LtoG(x, dS0, dT, 0.0, 0.0, 0.0);
463
464  for (k = 0; k < n; k++)
465    mat[k][n] = x[k];
466}
467
468
469void Cell_Concat(double dP)
470{
471  long      j = 0, i1 = 0;
472  double    PB2 = 0.0;
473  CellType  *cellp;
474  elemtype  *elemp;
475  MpoleType *M;
476
477  if (dP != globval.dPparticle) {
478    Cell_SetdP(dP); cellconcat = false;
479  }
480 
481  if (cellconcat) return;
482
483  if (trace) printf("concatenating\n");
484
485  cellconcat = true; i1 = 0; ntransfmat = 1;
486  UnitMat(n+1, transfmat[ntransfmat-1]);
487  kicks[ntransfmat-1][0] = 0;
488
489  do {
490    while ((linearel(i1+1) && (i1+1) <= globval.Cell_nLoc)) {
491      i1++;
492      cellp  = &Cell[i1]; elemp = &cellp->Elem;
493      switch (elemp->Pkind) {
494        case drift:
495          MulLsMat(elemp->D->D55, transfmat[ntransfmat-1]);
496          break;
497
498        case Mpole:
499          M = elemp->M;
500          GtoL_M(transfmat[ntransfmat-1], cellp->dT);
501          GtoL_dP(transfmat[ntransfmat-1], cellp->dT);
502          GtoL(transfmat[ntransfmat-1][n], cellp->dS, cellp->dT,
503          M->Pc0, M->Pc1, M->Ps1);
504          if (M->Pthick == thick)
505          { /* Drift Kick Drift */
506            MulLsMat(M->AU55, transfmat[ntransfmat-1]);
507            /* Assuming there is no quadrupole kick */
508            PB2 = M->PB[Quad+HOMmax];
509            M->PB[Quad+HOMmax] = 0.0;
510            thin_kick(M->Porder, M->PB, elemp->PL, 0.0, 0.0,
511                      transfmat[ntransfmat-1][n]);
512            M->PB[Quad+HOMmax] = PB2;
513            MulLsMat(M->AD55, transfmat[ntransfmat-1]);
514          } else {
515            /* Dipole kick */
516            thin_kick(M->Porder, M->PB, 1.0, 0.0, 0.0,
517                      transfmat[ntransfmat-1][n]);
518          }
519          LtoG_M(transfmat[ntransfmat-1], cellp->dT);
520          LtoG_dP(transfmat[ntransfmat-1], cellp->dT);
521          LtoG(transfmat[ntransfmat-1][n], cellp->dS, cellp->dT,
522               M->Pc0, M->Pc1, M->Ps1);
523          break;
524
525        case Wigl:
526          MulLsMat(elemp->W->W55, transfmat[ntransfmat-1]);
527          break;
528
529        case Insertion:
530           MulLsMat(elemp->ID->KD55, transfmat[ntransfmat-1]);
531          break;
532
533        case Cavity:   /* nothing */
534          break;
535
536        case marker:   /* nothing */
537          break;
538
539        default:
540          printf("**Cell_Concat: undefined type\n");
541          break;
542      }
543    }
544    j = 0;
545    while (!linearel(i1+j+1) && (i1+j+1) <= globval.Cell_nLoc) {
546      j++;
547      if (j >= maxkicks) {
548        printf("Cell_Concat maxkicks exceeded: %ld (%d)\n", j, maxkicks);
549        exit_(1);
550      }
551      cellp  = &Cell[i1+j]; elemp = &cellp->Elem;
552
553      if (elemp->Pkind != Mpole) continue;
554
555      M = elemp->M;
556
557      GtoL_M(transfmat[ntransfmat-1], cellp->dT);
558      GtoL_dP(transfmat[ntransfmat-1], cellp->dT);
559      GtoL(transfmat[ntransfmat-1][n], cellp->dS, cellp->dT, M->Pc0,
560           M->Pc1, M->Ps1);
561
562      if (M->Pthick == thick)
563        MulLsMat(M->AU55, transfmat[ntransfmat-1]);
564
565      kicks[ntransfmat-1][j-1] = i1 + j; kicks[ntransfmat-1][j] = 0;
566
567      ntransfmat++;
568      if (ntransfmat >= maxtransfmat) {
569        printf("Cell_Concat maxtransfmat exceeded: %ld (%d)\n",
570               ntransfmat, maxtransfmat);
571        exit_(1);
572      }
573      UnitMat(n+1, transfmat[ntransfmat-1]);
574      kicks[ntransfmat-1][0] = 0;
575
576      if (M->Pthick == thick)
577        MulLsMat(M->AD55, transfmat[ntransfmat-1]);
578
579      LtoG_M(transfmat[ntransfmat-1], cellp->dT);
580      LtoG_dP(transfmat[ntransfmat-1], cellp->dT);
581      LtoG(transfmat[ntransfmat-1][n], cellp->dS, cellp->dT, M->Pc0,
582           M->Pc1, M->Ps1);
583    }
584    i1 += j;
585  } while (i1 != globval.Cell_nLoc);
586}
587
588#undef n
589
590
591#define n 4
592void Cell_fPass(ss_vect<double> &x, long &lastpos)
593{
594  /* Purpose:
595       Compute the oneturn matrix and propagates xref through it
596       Nota: f means full                                                    */
597
598  long       i, j;
599  double     PB2;
600  CellType   *cellp;
601  elemtype   *elemp;
602  MpoleType  *M;
603
604
605  if (!CheckAmpl(x, 1))
606    lastpos = 1;
607  else { 
608    lastpos = globval.Cell_nLoc;
609    for (i = 0; i < ntransfmat; i++) {
610      LinsTrans(transfmat[i], x);
611      j = 0;
612      while (kicks[i][j] != 0) {
613        j++;
614        cellp = &Cell[kicks[i][j-1]]; elemp = &cellp->Elem; M = elemp->M;
615        CopyVec(n, x, Cell[kicks[i][j-1]].BeamPos);
616        if (M->Pthick == thick) {
617          PB2 = M->PB[Quad+HOMmax];
618          M->PB[Quad+HOMmax] = 0.0;
619          thin_kick(M->Porder, M->PB, elemp->PL, 0.0, 0.0, x);
620          M->PB[Quad+HOMmax] = PB2;
621        } else
622          thin_kick(M->Porder, M->PB, 1.0, 0.0, 0.0, x);
623
624        if (!CheckAmpl(x, kicks[i][j-1])) {
625          lastpos = kicks[i][j-1]; return;
626        }
627      }
628    }
629  }
630}
631#undef n
632
633
634#define n 4
635void Cell_fPass_M(ss_vect<double> &xref, Matrix &mat, long &lastpos)
636{
637  /* Purpose: called by Cell_MatGetCOD
638       Compute the oneturn matrix and propagates xref through it using
639       matrix formalism
640       Nota: f means full
641       
642   Input:
643       lastpos last position if unstable
644       x  starting vector
645
646   Output:
647       mat oneturnmatrix
648
649   Return:
650       none
651
652   Global variables:
653       transfmat contains transfert matrix for each linear element
654       ntransfmat number of transfer matrices
655       Cell contains all elements
656       
657   Specific functions:
658        CheckAmpl, MulLsMat, LinsTrans
659        thin_kick_M, thin_kick                                               */
660
661  long       i= 0, j = 0;
662  double     PB2 = 0.0;
663  CellType   *cellp;
664  elemtype   *elemp;
665  MpoleType  *M;
666
667  if (!CheckAmpl(xref, 1))
668    lastpos = 1;
669  else {
670    lastpos = globval.Cell_nLoc;
671    for (i = 0; i < ntransfmat; i++) {
672      MulLsMat(transfmat[i], mat); LinsTrans(transfmat[i], xref);
673      j = 0;
674      while (kicks[i][j] != 0) {
675        j++;
676        cellp  = &Cell[kicks[i][j-1]]; elemp = &cellp->Elem; M = elemp->M;
677
678        if (M->Pthick == thick) {
679          PB2 = M->PB[Quad+HOMmax];
680          M->PB[Quad+HOMmax] = 0.0;
681          thin_kick_M(M->Porder, M->PB, elemp->PL, 0.0, xref, mat);
682          thin_kick(M->Porder, M->PB, elemp->PL, 0.0, 0.0, xref);
683          M->PB[Quad+HOMmax] = PB2;
684        } else {
685          thin_kick_M(M->Porder, M->PB, 1.0, 0.0, xref, mat);
686          thin_kick(M->Porder, M->PB, 1.0, 0.0, 0.0, xref);
687        }
688
689        if (!CheckAmpl(xref, kicks[i][j-1])) {
690          lastpos = kicks[i][j-1]; return;
691        }
692      }
693    }
694  }
695}
696#undef n
697
698
699#define n 4
700bool Cell_GetCOD_M(long imax, double eps, double dP, long &lastpos)
701{
702  long    i = 0, j = 0;
703  double  dxabs = 0.0;
704  Vector  x0, x1, dx;
705  Matrix  A;
706
707  if (globval.MatMeth) Cell_Concat(dP);
708
709  CopyVec(n, globval.CODvect, x0);
710  x0[delta_] = dP; x0[ct_] = 0.0; i = 0;
711
712  do {
713    i++;
714    UnitMat(n+2, globval.OneTurnMat); x1 = x0;
715
716    Cell_fPass_M(x1, globval.OneTurnMat, lastpos); /* compute oneturn matrix */
717//     Cell_Pass_M(0, globval.Cell_nLoc, x1, globval.OneTurnMat, lastpos);
718
719    if (lastpos == globval.Cell_nLoc) globval.CODvect = x0;
720
721    CopyVec(n, x0, dx); SubVec(n, x1, dx);
722    CopyMat(n, globval.OneTurnMat, A);
723
724    /* A = A - Id */
725    for (j = 0; j < n; j++)
726      A[j][j]--;
727
728    if (InvMat(n, A)) {
729      if (lastpos == globval.Cell_nLoc) {
730        LinTrans(n, A, dx);
731        for (j = 0; j < n; j++)
732          x0[j] += dx[j];
733      }
734    } else
735      printf("  *** A is singular\n");
736
737    dxabs = xabs(4, dx);
738
739    if (trace) {
740      printf("--- CODLOOP%3ld, Err=% .3E/% .3E\n", i, dxabs, eps);
741      printf("% .5E % .5E\n", x0[0], x0[1]);
742      printf("% .5E % .5E\n", x0[2], x0[3]);
743      printf("% .5E % .5E\n", x0[4], x0[5]);
744    }
745  } while (i < imax && dxabs > eps && lastpos == globval.Cell_nLoc);
746
747  x0 = globval.CODvect; Cell_Pass(0, globval.Cell_nLoc, x0, lastpos);
748
749  if (dxabs <= eps && lastpos == globval.Cell_nLoc)
750    status.codflag = true;
751  else {
752    printf("Cell_GetCOD_M: GetCOD failed\n");
753    status.codflag = false;
754  }
755
756  return status.codflag;
757}
758#undef n
759
760/****************************************************************************/
761/* void Cell_GetCOD(long imax, double eps, double dP, long *lastpos)
762
763   Purpose:  called by Ring_Getchrom, getcod
764              Looks for a chromatic closed orbit at dP  and computes oneturn map
765              Search at precision eps and with a maximum of imax iterations
766              using DA method
767              If can't find the COD, then write beampos.dat, cod.out, flat_file_dbg.dat for debug;
768              and also print information for the particle lost.
769              If find the cod, then get the one turn map and cell pass, and
770              return status.codflag. The searched COD is saved in globval.CODvect
771             
772   Input:
773       imax number of iteration for cod search
774       dP   particle energy offset
775       eps  accuracy for cod search
776
777   Output:
778       
779
780   Return:
781      globval.CODvect   COD for the particle with energy offset dP
782       status.codflag   if true, particle is stable,
783                        if false, particle is unstable.
784
785   Global variables:
786       none
787
788   Specific functions:
789       Cell_Pass
790       
791       
792   Comments:
793       Called by getCOD, Ring_Pass
794
795       27/06/11   Correct the bug for the negative momentum compact factor; if the momentum
796                  compact is negative, then switch the Stable Fixed Point to Unstable Fixed Point.
797                  For Cell_GetCOD_M, the tracking is in 4D, so there is no switch of
798     fixed point for the lattice with negative momentum compact factor. 
799      28/06/11   Correct the one turn map for the lattice with negative momentum compact factor, now the one turn map is
800      tracked around the COD.
801****************************************************************************/
802bool Cell_getCOD(long imax, double eps, double dP, long &lastpos)
803{
804  long             j = 0, n = 0, n_iter = 0, h_RF = 0;
805  double           dxabs = 0.0;
806  iVector          jj;
807  ss_vect<double>  x0, x1, dx;
808  ss_vect<tps>     I, dx0, map;
809 
810  if(globval.Cavity_on)
811    n = 6L;
812  else
813    n = 4L;
814   
815  globval.dPparticle = dP;
816
817  if (n == 6) 
818    {
819    // initial guess is zero for 3 D.O.F.
820      x0[x_] = 0.0; 
821      x0[px_] = 0.0; 
822      x0[y_] = 0.0; 
823      x0[py_] = 0.0;
824      x0[delta_] = 0.0; 
825      x0[ct_] = 0.0;
826    } 
827  else 
828    {
829    // or eta*dP for 2 1/2 D.O.F.
830      x0[x_] = Cell[0].Eta[X_]*dP; 
831      x0[px_] = Cell[0].Etap[X_]*dP;
832      x0[y_] = Cell[0].Eta[Y_]*dP; 
833      x0[py_] = Cell[0].Etap[Y_]*dP;
834      x0[delta_] = dP; 
835      x0[ct_] = 0.0;
836    }
837
838  for (j = 0; j < 2*nd_tps; j++)
839    jj[j] = (j < n)? 1 : 0;
840
841  if (trace) 
842    {
843      cout << endl;
844      cout << "Cell_getCOD:" << endl;
845    }
846   
847  n_iter = 0; 
848  I.identity();
849 
850  do 
851    { //search for COD
852      n_iter++; 
853      map.identity(); 
854      map += x0;
855
856      Cell_Pass(0, globval.Cell_nLoc, map, lastpos); 
857
858      if (lastpos == globval.Cell_nLoc) {
859        x1 = map.cst(); 
860        dx = x0 - x1; 
861        dx0 = PInv(map-I-x1, jj)*dx;
862        dxabs = xabs(n, dx); 
863        x0 += dx0.cst();
864        } 
865      else{
866        prt_beampos("beampos.dat");
867        prt_cod("cod.out", globval.bpm, true);
868        prtmfile("flat_file_dbg.dat");
869//      prt_trace();
870        dxabs = 1e30; 
871        break;
872        }
873
874    if (trace) 
875      {
876        cout << scientific << setprecision(1)
877             << setw(3) << n_iter << " err = " 
878             << setw(7) << dxabs << "/" << setw(7) 
879             << eps << setprecision(5)
880             << "  x0 =" << setw(13) << x0 << endl;
881      }
882    } while ((dxabs >= eps) && (n_iter <= imax));
883
884 
885  //if momentum compact factor is negative, then swtich the fixed point(SFP->UFP) 
886  h_RF = Cell[Elem_GetPos(globval.cav, 1)].Elem.C->Ph;
887  if(globval.Alphac < 0) 
888     x0[ct_] =  -0.5*Cell[globval.Cell_nLoc].S/h_RF- x0[ct_]; //since the value of x0[ct_] is negative, so use "-0.5"
889  if(0)
890    cout<< " c*tau0 is: " << x0[ct_] <<endl;
891 
892  //get the one turn map for the close orbit
893  map.identity(); 
894  map += x0;
895  Cell_Pass(0, globval.Cell_nLoc, map, lastpos); 
896 
897 
898  status.codflag = dxabs < eps;
899
900  if (status.codflag) 
901    {
902      globval.CODvect = x0; 
903      getlinmat(6, map, globval.OneTurnMat);
904      Cell_Pass(0, globval.Cell_nLoc, x0, lastpos);
905    } 
906  else 
907    {
908       cout << "Cell_getCOD: failed to converge after " << n_iter << " iterations"
909            << ", dP=" << setw(13) << dP
910            << ", particle lost at element " << lastpos << endl;
911       cout << scientific << setprecision(5)
912            << " x0=" << setw(13) << x0 << endl;
913       cout << scientific << setprecision(5)
914            << "  x=" << setw(13) << map.cst() << endl;
915    }
916
917  return status.codflag;
918}
919
920/********************************************************
921bool GetCOD(long imax, double eps, double dP, long &lastpos)
922
923  Purpose: Get the closed orbit for the particle with energy spread dP
924 
925  Input:
926       imax number of iteration for cod search
927       dP   particle energy offset
928       eps  accuracy for cod search
929       
930  Output:
931 
932 
933  Return:
934    lastpos     last position when the particle is stale
935        cod     the bool status whether the COD is found or not     
936       
937   
938********************************************************/
939
940bool GetCOD(long imax, double eps, double dP, long &lastpos)
941{
942  bool  cod;
943
944  if (globval.MatMeth)
945    cod = Cell_GetCOD_M(imax, eps, dP, lastpos);
946  else
947    cod = Cell_getCOD(imax, eps, dP, lastpos);
948
949  return cod;
950}
951
952/************************************************
953void Cell_Init(void)
954
955  Purpose:
956    n....
957************************************************/
958void Cell_Init(void)
959{
960  long         i;
961  double       Stotal;
962  ElemFamType  *elemfamp;
963  elemtype     *elemp;
964
965  char  first_name[] = "begin          ";
966
967  if (debug)
968    printf("**  Cell_Init\n");
969
970  SI_init();  /* Initializes the constants for symplectic integrator */
971
972  memcpy(Cell[0].Elem.PName, first_name, sizeof(first_name));
973
974  for (i = 1; i <= globval.Elem_nFam; i++) 
975    {
976      elemfamp  = &ElemFam[i-1]; /* Get 1 of all elements stored in ElemFam
977                                  array */
978      elemp = &elemfamp->ElemF; // For switch structure: choice on element type
979      if (debug)
980        printf("Cell_Init, i:=%3ld: %*s\n", i, SymbolLength, elemp->PName);
981       
982      switch (elemp->Pkind) 
983      {
984        case drift:
985          Drift_Init(i);
986          break;
987     
988        case Mpole:
989          Mpole_Init(i);
990          break;
991     
992        case Wigl:
993          Wiggler_Init(i);
994          break;
995     
996        case FieldMap:
997          FieldMap_Init(i);
998          break;
999     
1000        case Insertion:
1001          Insertion_Init(i);
1002          break;
1003
1004        case Cavity:
1005          Cav_Init(i);
1006          break;
1007
1008        case marker:
1009          Marker_Init(i);
1010          break;
1011
1012        case Spreader:
1013          Spreader_Init(i);
1014          break;
1015
1016        case Recombiner:
1017          Recombiner_Init(i);
1018          break;
1019
1020        case Solenoid:
1021          Solenoid_Init(i);
1022          break;
1023
1024        default:
1025          printf("Cell_Init: undefined type\n");
1026          break;
1027      }
1028    }
1029 
1030 /* Computes s-location of each element in the structure */
1031  Stotal = 0.0;
1032  for (i = 0; i <= globval.Cell_nLoc; i++) 
1033    {
1034      Stotal += Cell[i].Elem.PL; 
1035      Cell[i].S = Stotal;
1036    }
1037}
Note: See TracBrowser for help on using the repository browser.