source: TRACY3/trunk/tracy/tracy/src/physlib.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: 99.0 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/* Current revision $Revision: 1.20 $
10 On branch $Name:  $
11 Latest change by $Author: zhang $
12*/
13
14/**************************/
15/* Routines for printing  */
16/**************************/
17
18/*******************************************************************************
19 *
20 *
21 *
22 *
23 ******************************************************************************/
24/**** same as asctime in C without the \n at the end****/
25char *asctime2(const struct tm *timeptr) {
26    // terminated with \0.
27    static char wday_name[7][4] = { "Sun", "Mon", "Tue", "Wed", "Thu", "Fri",
28            "Sat" };
29    // terminated with \0.
30    static char mon_name[12][4] = { "Jan", "Feb", "Mar", "Apr", "May", "Jun",
31            "Jul", "Aug", "Sep", "Oct", "Nov", "Dec" };
32    static char result[26];
33
34    sprintf(result, "%.3s %.3s%3d %.2d:%.2d:%.2d %d",
35            wday_name[timeptr->tm_wday], mon_name[timeptr->tm_mon],
36            timeptr->tm_mday, timeptr->tm_hour, timeptr->tm_min,
37            timeptr->tm_sec, 1900 + timeptr->tm_year);
38    return result;
39}
40/*******************************************************************************
41 *
42 *
43 *
44 *
45 ******************************************************************************/
46/** Get time and date **/
47struct tm* GetTime() {
48    struct tm *whattime;
49    /* Get time and date */
50    time_t aclock;
51    time(&aclock); /* Get time in seconds */
52    whattime = localtime(&aclock); /* Convert time to struct */
53    return whattime;
54}
55/****************************************************************************
56 * uint32_t stampstart()
57
58
59 Purpose: record time in millliseconds
60
61 Input:
62 none
63
64 Output:
65 non
66
67 Return:
68 time in milliseconds
69
70 Global variables:
71 none
72
73 specific functions:
74 none
75
76 Comments:
77 to be used with stampstop()
78
79 ****************************************************************************/
80
81uint32_t stampstart(void) {
82    struct timeval tv;
83    struct timezone tz;
84    struct tm *tm;
85    uint32_t start;
86    const bool timedebug = false;
87
88    // get the time
89    gettimeofday(&tv, &tz);
90    tm = localtime(&tv.tv_sec);
91
92    // print detailed time in milliseconds
93    if (timedebug)
94        printf("TIMESTAMP-START\t  %d:%02d:%02d:%d (~%d ms)\n", tm->tm_hour,
95                tm->tm_min, tm->tm_sec, tv.tv_usec, tm->tm_hour * 3600 * 1000
96                        + tm->tm_min * 60 * 1000 + tm->tm_sec * 1000
97                        + tv.tv_usec / 1000);
98
99    start = tm->tm_hour * 3600 * 1000 + tm->tm_min * 60 * 1000 + tm->tm_sec
100            * 1000 + tv.tv_usec / 1000;
101
102    return (start);
103
104}
105
106// compute time elapsed since start
107/****************************************************************************
108 * uint32_t stampstop(uint32_t start)
109
110 Purpose: compute time elapsed since start time
111
112 Input:
113 start starting time i millisecond
114
115 Output:
116 none
117
118 Return:
119 none
120
121 Global variables:
122 none
123
124 specific functions:
125 none
126
127 Comments:
128 to be used with stampstart
129
130 ****************************************************************************/
131uint32_t stampstop(uint32_t start) {
132
133    struct timeval tv;
134    struct timezone tz;
135    struct tm *tm;
136    uint32_t stop;
137    const bool timedebug = false;
138    bool prt = true;
139    // get the time
140    gettimeofday(&tv, &tz);
141    tm = localtime(&tv.tv_sec);
142
143    stop = tm->tm_hour * 3600 * 1000 + tm->tm_min * 60 * 1000 + tm->tm_sec
144            * 1000 + tv.tv_usec / 1000;
145
146    // print detailed time in milliseconds
147    if (timedebug) {
148        printf("TIMESTAMP-END\t  %d:%02d:%02d:%d (~%d ms) \n", tm->tm_hour,
149                tm->tm_min, tm->tm_sec, tv.tv_usec, tm->tm_hour * 3600 * 1000
150                        + tm->tm_min * 60 * 1000 + tm->tm_sec * 1000
151                        + tv.tv_usec / 1000);
152        printf("ELAPSED\t  %d ms\n", stop - start);
153    }
154
155    uint32_t delta, hour, minute, second, millisecond;
156    delta = stop - start;
157    hour = delta / 3600000;
158    minute = (delta - 3600000 * hour) / 60000;
159    second = (delta - 3600000 * hour - minute * 60000) / 1000;
160    millisecond = delta - 3600000 * hour - minute * 60000 - second * 1000;
161
162    if (prt)
163        printf("ELAPSED\t  %d h %d min %d s %d ms\n", hour, minute, second,
164                millisecond);
165
166    return (stop);
167
168}
169/****************************************************************************/
170/* void printglob(void)
171
172 Purpose:
173 Print global variables on screen
174 Print tunes and chromaticities
175 Print Oneturn matrix
176
177 Input:
178 none
179
180 Output:
181 output on the screen
182
183 Return:
184 none
185
186 Global variables:
187 globval
188
189 Specific functions:
190 none
191
192 Comments:
193 26/03/03 Oneturn matrix added
194 26/03/03 RF acceptance added
195 10/05/03 Momentum compaction factor added
196 16/05/03 Correction for a asymmetrical vaccum vessel
197 20/06/03 Add corrector, skew quad and bpm number
198 27/10/03 Add flag for radiation and chambre
199
200 Comments copied from Tracy 2.7(soleil),Written by L.Nadolski.
201 
202 03/06/2013 Add feature to print the summary on the sreen and an external file
203
204 ****************************************************************************/
205void printglob(void) {
206    printf("\n***************************************************************"
207        "***************\n");
208    printf("\n");
209    printf("  dPcommon     =  %9.3e  dPparticle   =  %9.3e"
210        "  Energy [GeV] = %.3f\n", globval.dPcommon, globval.dPparticle,
211            globval.Energy);
212    printf("  MaxAmplx [m] = %9.3e   MaxAmply [m] = %9.3e"
213        "   RFAccept [%%] = +/- %4.2f\n", Cell[0].maxampl[X_][1],
214            Cell[0].maxampl[Y_][1], globval.delta_RF * 1e2);
215    printf("  MatMeth      =  %s     ", globval.MatMeth ? "TRUE " : "FALSE");
216    printf(" Cavity_On    =  %s    ", globval.Cavity_on ? "TRUE " : "FALSE");
217    printf("  Radiation_On = %s     \n", globval.radiation ? "TRUE " : "FALSE");
218   
219    if(globval.bpm == 0)
220      printf("  bpm          =  0"); 
221    else
222      printf("  bpm          =  %3d", GetnKid(globval.bpm));
223    if(globval.qt == 0)
224      printf("                             qt           = 0           \n"); 
225    else
226      printf("                             qt           = %3d         \n", GetnKid(globval.qt));   
227    if(globval.hcorr == 0)
228      printf("  hcorr         =  0"); 
229    else
230      printf("  hcorr         = %3d", GetnKid(globval.hcorr));
231    if(globval.vcorr == 0)
232      printf("                             vcorr        =  0        \n"); 
233    else
234      printf("                             vcorr        = %3d      \n", GetnKid(globval.vcorr));
235                                           
236    printf(" Chambre_On   = %s     \n", globval.Aperture_on ? "TRUE " : "FALSE");
237 
238    printf("  alphac       =   %8.4e\n", globval.Alphac);
239    printf("  nux          =   %8.6f      nuy  =  %8.6f",
240            globval.TotalTune[X_], globval.TotalTune[Y_]);
241    if (globval.Cavity_on)
242        printf("  omega  = %13.9f\n", globval.Omega);
243    else {
244        printf("\n");
245        printf("  ksix         =    %8.6f      ksiy  =   %8.6f\n",
246                globval.Chrom[X_], globval.Chrom[Y_]);
247    }
248    printf("\n");
249    printf("  OneTurn matrix:\n");
250    printf("\n");
251    prtmat(ss_dim, globval.OneTurnMat);
252
253    printf("\nTwiss parameters at entrance:\n");
254    printf(
255            "   Betax [m] = % 9.3e  Alphax = % 9.3e Etax [m] = % 9.3e Etaxp = % 9.3e\n"
256                "   Betay [m] = % 9.3e  Alphay = % 9.3e Etay [m] = % 9.3e Etayp = % 9.3e\n\n",
257            Cell[1].Beta[X_], Cell[1].Alpha[X_], Cell[1].Eta[X_],
258            Cell[1].Etap[X_], Cell[1].Beta[Y_], Cell[1].Alpha[Y_],
259            Cell[1].Eta[Y_], Cell[1].Etap[Y_]);
260
261    fflush( stdout);
262}
263
264/****************************************************************************/
265/* void printglob2file(const char fic[])
266
267 Purpose:
268 Print global variables,
269       tunes and chromaticities, and
270       Oneturn matrix,
271       to an external file.
272
273 Input:
274 none
275
276 Output:
277 output on the screen
278
279 Return:
280 none
281
282 Global variables:
283 globval
284
285 Specific functions:
286 none
287
288 Comments:
289 03/06/2013  by Jianfeng Zhang @ LAL
290 The same feature as the file printglob(void), but
291 added feature to print the lattice summary to an external file
292
293 
294 ****************************************************************************/
295void printglob2file(const char fic[]) {
296  FILE *fout;
297  int i=0,j=0;
298 
299  fout = fopen(fic,"w");
300 
301    fprintf(fout, "\n***************************************************************"
302        "***************\n");
303    fprintf(fout,"\n");
304    fprintf(fout,"  dPcommon     =  %9.3e  dPparticle   =  %9.3e"
305        "  Energy [GeV] = %.3f\n", globval.dPcommon, globval.dPparticle,
306            globval.Energy);
307    fprintf(fout,"  MaxAmplx [m] = %9.3e   MaxAmply [m] = %9.3e"
308        "   RFAccept [%%] = +/- %4.2f\n", Cell[0].maxampl[X_][1],
309            Cell[0].maxampl[Y_][1], globval.delta_RF * 1e2);
310    fprintf(fout,"  MatMeth      =  %s     ", globval.MatMeth ? "TRUE " : "FALSE");
311    fprintf(fout," Cavity_On    =  %s    ", globval.Cavity_on ? "TRUE " : "FALSE");
312    fprintf(fout,"  Radiation_On = %s     \n", globval.radiation ? "TRUE " : "FALSE");
313   
314    if(globval.bpm == 0)
315      fprintf(fout,"  bpm          =  0"); 
316    else
317      fprintf(fout,"  bpm          =  %3d", GetnKid(globval.bpm));
318    if(globval.qt == 0)
319      fprintf(fout,"                             qt           = 0           \n"); 
320    else
321      fprintf(fout,"                             qt           = %3d         \n", GetnKid(globval.qt));   
322    if(globval.hcorr == 0)
323      fprintf(fout,"  hcorr         =  0"); 
324    else
325      fprintf(fout,"  hcorr         = %3d", GetnKid(globval.hcorr));
326    if(globval.vcorr == 0)
327      fprintf(fout,"                             vcorr        =  0        \n"); 
328    else
329      fprintf(fout,"                             vcorr        = %3d      \n", GetnKid(globval.vcorr));
330                                           
331    fprintf(fout," Chambre_On   = %s     \n", globval.Aperture_on ? "TRUE " : "FALSE");
332 
333    fprintf(fout,"  alphac       =   %8.4e\n", globval.Alphac);
334    fprintf(fout,"  nux          =   %8.6f      nuy  =  %8.6f",
335            globval.TotalTune[X_], globval.TotalTune[Y_]);
336    if (globval.Cavity_on)
337        fprintf(fout,"  omega  = %13.9f\n", globval.Omega);
338    else {
339        fprintf(fout,"\n");
340        fprintf(fout,"  ksix         =    %8.6f      ksiy  =   %8.6f\n",
341                globval.Chrom[X_], globval.Chrom[Y_]);
342    }
343    fprintf(fout,"\n");
344   
345    fprintf(fout,"  OneTurn matrix:\n");
346    for(i=0;i<ss_dim;i++){
347      fprintf(fout,"\n");
348      for(j=0;j<ss_dim;j++)
349        fprintf(fout,"%14.6e",globval.OneTurnMat[i][j]);
350    }
351
352    fprintf(fout,"\n\nTwiss parameters at entrance:\n");
353    fprintf(fout,
354            "   Betax [m] = % 9.3e  Alphax = % 9.3e Etax [m] = % 9.3e Etaxp = % 9.3e\n"
355                "   Betay [m] = % 9.3e  Alphay = % 9.3e Etay [m] = % 9.3e Etayp = % 9.3e\n\n",
356            Cell[1].Beta[X_], Cell[1].Alpha[X_], Cell[1].Eta[X_],
357            Cell[1].Etap[X_], Cell[1].Beta[Y_], Cell[1].Alpha[Y_],
358            Cell[1].Eta[Y_], Cell[1].Etap[Y_]);
359
360    fclose(fout);
361}
362
363/****************************************************************************/
364/* void printlatt(const char fic[])
365
366 Purpose:
367 Print twiss parameters into file linlat.out
368 name, position, alpha, beta, phase, dispersion and  its derivative
369
370 Input:
371 none
372
373 Output:
374 none
375
376 Return:
377 none
378
379 Global variables:
380 globval
381
382 Specific functions:
383 getelem
384
385 Comments:
386 28/04/03 outfile end with .out extension instead of .dat
387 Twiss parameters are computed at the end of elements
388
389 ****************************************************************************/
390//void printlatt(void)
391void printlatt(const char fic[]) {
392    long int i = 0;
393    FILE *outf;
394    struct tm *newtime;
395
396    /* Get time and date */
397    newtime = GetTime();
398   
399    if ((outf = fopen(fic, "w")) == NULL) {
400        fprintf(stdout, "printlatt: Error while opening file %s \n", fic);
401        exit(1);
402    }
403
404    fprintf(outf, "# TRACY III  -- %s -- %s \n", fic,
405            asctime2(newtime));
406    fprintf(
407            outf,
408            "#       name           s      alphax  betax   nux   etax   etapx  alphay  betay   nuy     etay        etapy\n");
409    fprintf(
410            outf,
411            "#                     [m]              [m]           [m]                   [m]             [m]\n");
412    fprintf(outf, "#                     exit        \n");
413
414    for (i = 1L; i <= globval.Cell_nLoc; i++) {
415        fprintf(
416                outf,
417                "%4ld:%.*s% 8.4f% 8.3f% 7.3f% 7.3f% 7.3f% 7.3f% 8.3f% 7.3f% 7.3f% 12.3e% 12.3e\n",
418                i, SymbolLength, Cell[i].Elem.PName, Cell[i].S,
419                Cell[i].Alpha[X_], Cell[i].Beta[X_], Cell[i].Nu[X_],
420                Cell[i].Eta[X_], Cell[i].Etap[X_], Cell[i].Alpha[Y_],
421                Cell[i].Beta[Y_], Cell[i].Nu[Y_], Cell[i].Eta[Y_],
422                Cell[i].Etap[Y_]);
423    }
424    fclose(outf);
425}
426/*******************************************************************************
427 *
428 *
429 *
430 *
431 ******************************************************************************/
432double int_curly_H1(long int n) {
433    /* Integration with Simpson's Rule */
434
435    double curly_H;
436    Vector2 alpha[3], beta[3], nu[3], eta[3], etap[3];
437
438    // only works for matrix style calculations
439    get_twiss3(n, alpha, beta, nu, eta, etap);
440
441    curly_H = (get_curly_H(alpha[0][X_], beta[0][X_], eta[0][X_], etap[0][X_])
442            + 4.0 * get_curly_H(alpha[1][X_], beta[1][X_], eta[1][X_],
443                    etap[1][X_]) + get_curly_H(alpha[2][X_], beta[2][X_],
444            eta[2][X_], etap[2][X_])) / 6.0;
445
446    return curly_H;
447}
448/*******************************************************************************
449 *
450 *
451 *
452 *
453 ******************************************************************************/
454void prt_sigma(void) {
455    long int i;
456    double code = 0.0;
457    FILE *outf;
458
459    outf = file_write("../out/sigma.out");
460
461    fprintf(outf, "#  name     s   sqrt(sx)   sqrt(sx')  sqrt(sy)  sqrt(sy')\n");
462    fprintf(outf, "#           [m]   [mm]       [mrad]     [mm]     [mrad]\n");
463    fprintf(outf, "#\n");
464
465    for (i = 0; i <= globval.Cell_nLoc; i++) {
466        switch (Cell[i].Elem.Pkind) {
467        case drift:
468            code = 0.0;
469            break;
470        case Mpole:
471            if (Cell[i].Elem.M->Pirho != 0)
472                code = 0.5;
473            else if (Cell[i].Elem.M->PBpar[Quad + HOMmax] != 0)
474                code = sgn(Cell[i].Elem.M->PBpar[Quad + HOMmax]);
475            else if (Cell[i].Elem.M->PBpar[Sext + HOMmax] != 0)
476                code = 1.5 * sgn(Cell[i].Elem.M->PBpar[Sext + HOMmax]);
477            else if (Cell[i].Fnum == globval.bpm)
478                code = 2.0;
479            else
480                code = 0.0;
481            break;
482        default:
483            code = 0.0;
484            break;
485        }
486        fprintf(outf, "%4ld %.*s %6.2f %4.1f %9.3e %9.3e %9.3e %9.3e\n", i,
487                SymbolLength, Cell[i].Elem.PName, Cell[i].S, code, 1e3 * sqrt(
488                        Cell[i].sigma[x_][x_]), 1e3 * sqrt(fabs(
489                        Cell[i].sigma[x_][px_])), 1e3 * sqrt(
490                        Cell[i].sigma[y_][y_]), 1e3 * sqrt(fabs(
491                        Cell[i].sigma[y_][py_])));
492    }
493
494    fclose(outf);
495}
496
497void recalc_S(void) {
498    long int k;
499    double S_tot;
500
501    S_tot = 0.0;
502    for (k = 0; k <= globval.Cell_nLoc; k++) {
503        S_tot += Cell[k].Elem.PL;
504        Cell[k].S = S_tot;
505    }
506}
507/**************************************************************************
508getcod(double dP, long &lastpos)
509
510   Purpose:
511       Find the closed orbit for the particle with energy offset dP
512       
513   Input:
514       imax number of iteration for cod search
515       dP   particle energy offset
516       eps  accuracy for cod search
517       
518   Output:
519   
520   Return:
521   
522   
523   Comments:
524     
525***************************************************************************/       
526bool getcod(double dP, long &lastpos) {
527    return GetCOD(globval.CODimax, globval.CODeps, dP, lastpos);
528}
529
530void get_twiss3(long int loc, Vector2 alpha[], Vector2 beta[], Vector2 nu[],
531        Vector2 eta[], Vector2 etap[]) {
532    /* Get Twiss functions at magnet entrance-, center-, and exit. */
533
534    long int j, k;
535    Vector2 dnu;
536    Matrix Ascr0, Ascr1;
537
538    // Lattice functions at the magnet entrance
539    for (k = 0; k <= 1; k++) {
540        alpha[0][k] = Cell[loc - 1].Alpha[k];
541        beta[0][k] = Cell[loc - 1].Beta[k];
542        nu[0][k] = Cell[loc - 1].Nu[k];
543        eta[0][k] = Cell[loc - 1].Eta[k];
544        etap[0][k] = Cell[loc - 1].Etap[k];
545    }
546
547    UnitMat(5, Ascr0);
548    for (k = 0; k <= 1; k++) {
549        Ascr0[2 * k][2 * k] = sqrt(Cell[loc - 1].Beta[k]);
550        Ascr0[2 * k][2 * k + 1] = 0.0;
551        Ascr0[2 * k + 1][2 * k] = -Cell[loc - 1].Alpha[k] / Ascr0[2 * k][2 * k];
552        Ascr0[2 * k + 1][2 * k + 1] = 1 / Ascr0[2 * k][2 * k];
553    }
554    Ascr0[0][4] = eta[0][X_];
555    Ascr0[1][4] = etap[0][X_];
556    Ascr0[2][4] = eta[1][Y_];
557    Ascr0[3][4] = etap[1][Y_];
558    CopyMat(5, Ascr0, Ascr1);
559    MulLMat(5, Cell[loc].Elem.M->AU55, Ascr1);
560    dnu[0]
561            = (atan(Ascr1[0][1] / Ascr1[0][0])
562                    - atan(Ascr0[0][1] / Ascr0[0][0])) / (2.0 * M_PI);
563    dnu[1]
564            = (atan(Ascr1[2][3] / Ascr1[2][2])
565                    - atan(Ascr0[2][3] / Ascr0[2][2])) / (2.0 * M_PI);
566
567    // Lattice functions at the magnet center
568    for (k = 0; k <= 1; k++) {
569        alpha[1][k] = -Ascr1[2 * k][2 * k] * Ascr1[2 * k + 1][2 * k] - Ascr1[2
570                * k][2 * k + 1] * Ascr1[2 * k + 1][2 * k + 1];
571        beta[1][k] = pow(Ascr1[2 * k][2 * k], 2.0) + pow(
572                Ascr1[2 * k][2 * k + 1], 2);
573        nu[1][k] = nu[0][k] + dnu[k];
574        j = 2 * k + 1;
575        eta[1][k] = Ascr1[j - 1][4];
576        etap[1][k] = Ascr1[j][4];
577    }
578
579    CopyMat(5, Ascr1, Ascr0);
580    MulLMat(5, Cell[loc].Elem.M->AD55, Ascr1);
581    dnu[0]
582            = (atan(Ascr1[0][1] / Ascr1[0][0])
583                    - atan(Ascr0[0][1] / Ascr0[0][0])) / (2.0 * M_PI);
584    dnu[1]
585            = (atan(Ascr1[2][3] / Ascr1[2][2])
586                    - atan(Ascr0[2][3] / Ascr0[2][2])) / (2.0 * M_PI);
587
588    // Lattice functions at the magnet exit
589    for (k = 0; k <= 1; k++) {
590        alpha[2][k] = -Ascr1[2 * k][2 * k] * Ascr1[2 * k + 1][2 * k] - Ascr1[2
591                * k][2 * k + 1] * Ascr1[2 * k + 1][2 * k + 1];
592        beta[2][k] = pow(Ascr1[2 * k][2 * k], 2.0) + pow(
593                Ascr1[2 * k][2 * k + 1], 2);
594        nu[2][k] = nu[1][k] + dnu[k];
595        j = 2 * k + 1;
596        eta[2][k] = Ascr1[j - 1][4];
597        etap[2][k] = Ascr1[j][4];
598    }
599}
600/*******************************************************************************
601 *
602 *
603 *
604 *
605 ******************************************************************************/
606void getabn(Vector2 &alpha, Vector2 &beta, Vector2 &nu) {
607    Vector2 gamma;
608    Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
609}
610/*******************************************************************************
611 *
612 *
613 *
614 *
615 ******************************************************************************/
616void TraceABN(long i0, long i1, const Vector2 &alpha, const Vector2 &beta,
617        const Vector2 &eta, const Vector2 &etap, const double dP) {
618    long i=0, j=0;
619    double sb=0.0;
620    ss_vect<tps> Ascr;
621
622    UnitMat(6, globval.Ascr);
623    for (i = 1; i <= 2; i++) {
624        sb = sqrt(beta[i - 1]);
625        j = i * 2 - 1;
626        globval.Ascr[j - 1][j - 1] = sb;
627        globval.Ascr[j - 1][j] = 0.0;
628        globval.Ascr[j][j - 1] = -(alpha[i - 1] / sb);
629        globval.Ascr[j][j] = 1 / sb;
630    }
631    globval.Ascr[0][4] = eta[0];
632    globval.Ascr[1][4] = etap[0];
633    globval.Ascr[2][4] = eta[1];
634    globval.Ascr[3][4] = etap[1];
635
636    for (i = 0; i < 6; i++)
637        globval.CODvect[i] = 0.0;
638    globval.CODvect[4] = dP;
639
640    if (globval.MatMeth)
641        Cell_Twiss_M(i0, i1, globval.Ascr, false, false, dP);
642    else {
643        for (i = 0; i <= 5; i++) {
644            Ascr[i] = tps(globval.CODvect[i]);
645            for (j = 0; j <= 5; j++)
646                Ascr[i] += globval.Ascr[i][j] * tps(0.0, j + 1);
647        }
648        Cell_Twiss(i0, i1, Ascr, false, false, dP);
649    }
650
651}
652
653/****************************************************************************/
654/* void FitTune(long qf, long qd, double nux, double nuy)
655
656 Purpose:
657 Fit tunes to the target values using quadrupoles "qf" and "qd"
658 Input:
659 qf : tuned quadrupole
660 qd : tuned quadrupole
661 nux: target horizontal tune
662 nuy: target vertical tune
663 Output:
664 none
665
666 Return:
667 none
668
669 Global variables:
670
671 specific functions:
672
673 Comments:
674
675 ****************************************************************************/
676void FitTune(long qf, long qd, double nux, double nuy) {
677    long i;
678    iVector2 nq = { 0, 0 };
679    Vector2 nu = { 0.0, 0.0 };
680    fitvect qfbuf, qdbuf;
681
682    /* Get elements for the first quadrupole family */
683    nq[X_] = GetnKid(qf);
684    for (i = 1; i <= nq[X_]; i++)
685        qfbuf[i - 1] = Elem_GetPos(qf, i);
686
687    /* Get elements for the second quadrupole family */
688    nq[Y_] = GetnKid(qd);
689    for (i = 1; i <= nq[Y_]; i++)
690        qdbuf[i - 1] = Elem_GetPos(qd, i);
691
692    nu[X_] = nux;
693    nu[Y_] = nuy;
694
695    /* fit tunes */
696    Ring_Fittune(nu, nueps, nq, qfbuf, qdbuf, nudkL, nuimax);
697}
698
699/****************************************************************************/
700/* void FitChrom(long sf, long sd, double ksix, double ksiy)
701
702 Purpose:
703 Fit chromaticities to the target values using sextupoles "sf" and "sd"
704 Input:
705 sf  : tuned sextupole
706 sd  : tuned sextupole
707 ksix: target horizontal chromaticity
708 ksiy: target vertical chromaticity
709 Output:
710 none
711
712 Return:
713 none
714
715 Global variables:
716
717 specific functions:
718
719 Comments:
720
721 ****************************************************************************/
722
723void FitChrom(long sf, long sd, double ksix, double ksiy) {
724    long i;
725    iVector2 ns = { 0, 0 };
726    fitvect sfbuf, sdbuf;
727    Vector2 ksi = { 0.0, 0.0 };
728
729    /* Get elements for the first sextupole family */
730    ns[X_] = GetnKid(sf);
731    for (i = 1; i <= ns[X_]; i++)
732        sfbuf[i - 1] = Elem_GetPos(sf, i);
733
734    /* Get elements for the second sextupole family */
735    ns[Y_] = GetnKid(sd);
736    for (i = 1; i <= ns[Y_]; i++)
737        sdbuf[i - 1] = Elem_GetPos(sd, i);
738
739    ksi[X_] = ksix;
740    ksi[Y_] = ksiy;
741
742    /* Fit chromaticities */
743    /*    Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, 1.0, 1);*/
744    Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, ksidkpL, ksiimax);
745}
746/*******************************************************************************
747 *
748 *
749 *
750 *
751 ******************************************************************************/
752void FitDisp(long q, long pos, double eta) {
753    long i=0L, nq=0L;
754    fitvect qbuf;
755
756    /* Get elements for the quadrupole family */
757    nq = GetnKid(q);
758    for (i = 1; i <= nq; i++)
759        qbuf[i - 1] = Elem_GetPos(q, i);
760
761    Ring_FitDisp(pos, eta, dispeps, nq, qbuf, dispdkL, dispimax);
762}
763/*******************************************************************************
764 *
765 *
766 *
767 *
768 ******************************************************************************/
769#define nfloq           4
770
771void getfloqs(Vector &x) {
772    // Transform to Floquet space
773    LinTrans(nfloq, globval.Ascrinv, x);
774}
775
776#undef nfloq
777
778/*******************************************************************************
779 *
780 *
781 *
782 *
783 ******************************************************************************/
784#define ntrack          4
785
786// 4D tracking in normal or Floquet space over nmax turns
787
788void track(const char *file_name, double ic1, double ic2, double ic3,
789        double ic4, double dp, long int nmax, long int &lastn,
790        long int &lastpos, int floqs, double f_rf) {
791    /* Single particle tracking around closed orbit:
792
793     Output                floqs
794
795     Phase Space               0     [x, px, y, py, delta, ct]
796     Floquet Space             1     [x^, px^, y^, py^, delta, ct]
797     Action-Angle Variables    2     [2Jx, phx, 2Jy, phiy, delta, ct]
798
799     */
800    long int i=0L;
801    double twoJx=0.0, twoJy=0.0, phix=0.0, phiy=0.0, scl_1 = 1.0, scl_2 = 1.0;
802    Vector x0, x1, x2, xf;
803    FILE *outf;
804
805    bool prt = false;
806
807    if ((floqs == 0)) {
808        scl_1 = 1e3;
809        scl_2 = 1e3;
810        x0[x_] = ic1;
811        x0[px_] = ic2;
812        x0[y_] = ic3;
813        x0[py_] = ic4;
814    } else if ((floqs == 1)) {
815        scl_1 = 1.0;
816        scl_2 = 1.0;
817        x0[x_] = ic1;
818        x0[px_] = ic2;
819        x0[y_] = ic3;
820        x0[py_] = ic4;
821        LinTrans(4, globval.Ascr, x0);
822    } else if (floqs == 2) {
823        scl_1 = 1e6;
824        scl_2 = 1.0;
825        x0[x_] = sqrt(ic1) * cos(ic2);
826        x0[px_] = -sqrt(ic1) * sin(ic2);
827        x0[y_] = sqrt(ic3) * cos(ic4);
828        x0[py_] = -sqrt(ic3) * sin(ic4);
829        LinTrans(4, globval.Ascr, x0);
830    }
831
832    outf = file_write(file_name);
833    fprintf(outf, "# Tracking with TRACY");
834    getcod(dp, lastpos);
835    if (floqs == 0)
836        fprintf(outf, "\n");
837    else if (floqs == 1) {
838        Ring_GetTwiss(false, dp);
839        fprintf(outf, "# (Floquet space)\n");
840    } else if (floqs == 2) {
841        Ring_GetTwiss(false, dp);
842        fprintf(outf, "# (Action-Angle variables)\n");
843    }
844    fprintf(outf, "#\n");
845    fprintf(outf, "#%3d%6ld% .1E% .1E% .1E% .1E% 7.5f% 7.5f\n", 1, nmax, 1e0,
846            1e0, 0e0, 0e0, globval.TotalTune[0], globval.TotalTune[1]);
847    if (floqs == 0) {
848        fprintf(outf,
849                "#    N       x            p_x            y            p_y");
850        fprintf(outf, "          delta          cdt\n");
851        fprintf(outf, "#           [mm]         [mrad]"
852            "         [mm]         [mrad]");
853    } else if (floqs == 1) {
854        fprintf(outf, "#    N     x^          px^          y^          py^");
855        fprintf(outf, "          delta          cdt");
856        fprintf(outf, "#                              "
857            "                            ");
858    } else if (floqs == 2) {
859        fprintf(outf,
860                "#    N     2Jx          phi_x          2Jy          phi_y");
861        fprintf(outf, "          delta          cdt\n");
862        fprintf(outf, "#                              "
863            "                            ");
864    }
865    if (f_rf == 0.0) {
866        fprintf(outf, "         [%%]           [mm]\n");
867        fprintf(outf, "#\n");
868        fprintf(outf, "%4d %23.16e %23.16e %23.16e %23.16e %23.16e %23.16e\n",
869                0, scl_1 * ic1, scl_2 * ic2, scl_1 * ic3, scl_2 * ic4,
870                1e2 * dp, 1e3 * globval.CODvect[ct_]);
871    } else {
872        fprintf(outf, "         [%%]           [deg]\n");
873        fprintf(outf, "#\n");
874        fprintf(outf, "%4d %23.16e %23.16e %23.16e %23.16e %23.16e %23.16e\n",
875                0, scl_1 * ic1, scl_2 * ic2, scl_1 * ic3, scl_2 * ic4,
876                1e2 * dp, 2.0 * f_rf * 180.0 * globval.CODvect[ct_] / c0);
877    }
878    x2[x_] = x0[x_] + globval.CODvect[x_];
879    x2[px_] = x0[px_] + globval.CODvect[px_];
880    x2[y_] = x0[y_] + globval.CODvect[y_];
881    x2[py_] = x0[py_] + globval.CODvect[py_];
882    if (globval.Cavity_on) {
883        x2[delta_] = dp + globval.CODvect[delta_];
884        x2[ct_] = globval.CODvect[ct_];
885    } else {
886        x2[delta_] = dp;
887        x2[ct_] = 0.0;
888    }
889
890    lastn = 0;
891
892    if (prt) {
893        printf("\n");
894        printf("track:\n");
895        printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn, 1e3
896                * x2[x_], 1e3 * x2[px_], 1e3 * x2[y_], 1e3 * x2[py_], 1e2
897                * x2[delta_], 1e3 * x2[ct_]);
898    }
899
900    if (globval.MatMeth)
901        Cell_Concat(dp);
902
903    do {
904        (lastn)++;
905        for (i = 0; i < nv_; i++)
906            x1[i] = x2[i];
907
908        if (globval.MatMeth) {
909            Cell_fPass(x2, lastpos);
910        } else {
911            Cell_Pass(0, globval.Cell_nLoc, x2, lastpos);
912        }
913
914        for (i = x_; i <= py_; i++)
915            xf[i] = x2[i] - globval.CODvect[i];
916
917        for (i = delta_; i <= ct_; i++)
918            if (globval.Cavity_on && (i != ct_))
919                xf[i] = x2[i] - globval.CODvect[i];
920            else
921                xf[i] = x2[i];
922
923        if (floqs == 1)
924            getfloqs(xf);
925        else if (floqs == 2) {
926            getfloqs(xf);
927            twoJx = pow(xf[x_], 2) + pow(xf[px_], 2);
928            twoJy = pow(xf[y_], 2) + pow(xf[py_], 2);
929            phix = atan2(xf[px_], xf[x_]);
930            phiy = atan2(xf[py_], xf[y_]);
931            xf[x_] = twoJx;
932            xf[px_] = phix;
933            xf[y_] = twoJy;
934            xf[py_] = phiy;
935        }
936        if (f_rf == 0.0)
937            fprintf(
938                    outf,
939                    "%4ld %23.16le %23.16le %23.16le %23.16le %23.16le %23.16le\n",
940                    lastn, scl_1 * xf[0], scl_2 * xf[1], scl_1 * xf[2], scl_2
941                            * xf[3], 1e2 * xf[4], 1e3 * xf[5]);
942        else
943            fprintf(
944                    outf,
945                    "%4ld %23.16le %23.16le %23.16le %23.16le %23.16le %23.16le\n",
946                    lastn, scl_1 * xf[0], scl_2 * xf[1], scl_1 * xf[2], scl_2
947                            * xf[3], 1e2 * xf[4], 2.0 * f_rf * 180.0 * xf[5]
948                            / c0);
949    } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc));
950
951    if (globval.MatMeth)
952        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
953
954    fclose(outf);
955}
956
957#undef ntrack
958
959/*******************************************************************************
960 *
961 *
962 *
963 *
964 ******************************************************************************/
965#define step            0.1
966#define px              0.0
967#define py              0.0
968void track_(double r, struct LOC_getdynap *LINK) {
969    long i=0L, lastn=0L, lastpos=0L;
970    Vector x;
971
972    x[0] = r * cos(LINK->phi);
973    x[1] = px;
974    x[2] = r * sin(LINK->phi);
975    x[3] = py;
976    x[4] = LINK->delta;
977    x[5] = 0.0;
978    /* transform to phase space */
979    if (LINK->floqs) {
980        LinTrans(5, globval.Ascr, x);
981    }
982    for (i = 0; i <= 3; i++)
983        x[i] += globval.CODvect[i];
984    lastn = 0;
985    do {
986        lastn++;
987        if (globval.MatMeth) {
988            Cell_fPass(x, lastpos);
989        } else {
990            Cell_Pass(0, globval.Cell_nLoc, x, lastpos);
991        }
992    } while (lastn != LINK->nturn && lastpos == globval.Cell_nLoc);
993    LINK->lost = (lastn != LINK->nturn);
994}
995#undef step
996#undef px
997#undef py
998
999/****************************************************************************/
1000/* void Trac(double x, double px, double y, double py, double dp, double ctau,
1001 long nmax, long pos, long &lastn, long &lastpos, FILE *outf1)
1002
1003 Purpose:
1004 Single particle tracking w/ respect to the chamber centrum
1005 Based on the version of tracy 2 in soleillib.c
1006 Input:
1007 x, px, y, py 4 transverses coordinates
1008 ctau         c*tau
1009 dp           energy offset
1010 nmax         number of turns
1011 pos          starting position for tracking
1012 aperture     global physical aperture
1013
1014 Output:
1015 lastn       last n (should be nmax if  not lost)
1016 lastpos     last position in the ring
1017 outf1       output file with 6D phase at different pos
1018
1019 Return:
1020 lastn:  last tracking turn that particle is not lost
1021 lastpos:   last element position that particle is not lost
1022
1023 Global variables:
1024 globval
1025
1026 specific functions:
1027 Cell_Pass
1028
1029 Comments:
1030 Absolute TRACKING with respect to the center of the vacuum vessel
1031 BUG: last printout is wrong because not at pos but at the end of
1032 the ring
1033 26/04/03 print output for phase space is for position pos now
1034 01/12/03 tracking from 0 to pos -1L instead of 0 to pos
1035 (wrong observation point)
1036
1037 23/07/10 change the call variable in Cell_Pass( ): pos-1L --> pos (L704, L717).
1038 since the Cell_Pass( ) is tracking from element i0 to i1(tracy 3), and
1039 the Cell_Pass( ) is tracking from element i0+1L to i1(tracy 2).
1040
1041 ****************************************************************************/
1042void Trac(double x, double px, double y, double py, double dp, double ctau,
1043        long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) {
1044 
1045    bool lostF = true; /* Lost particle Flag */
1046    Vector x1; /* tracking coordinates */
1047    Vector2 aperture;
1048
1049    /* Compute closed orbit: useful if insertion devices */
1050
1051    aperture[0] = 1e0;
1052    aperture[1] = 1e0;
1053
1054    x1[0] = x;
1055    x1[1] = px;
1056    x1[2] = y;
1057    x1[3] = py;
1058    x1[4] = dp;
1059    x1[5] = ctau;
1060
1061    lastn = 0L;
1062    (lastpos) = pos;
1063
1064    if (trace)
1065        fprintf(outf1, "\n");
1066    if (trace)
1067        fprintf(outf1,
1068                "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n",
1069                lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
1070
1071    //  Cell_Pass(pos -1L, globval.Cell_nLoc, x1, lastpos);
1072    Cell_Pass(pos, globval.Cell_nLoc, x1, lastpos);
1073
1074    if (lastpos == globval.Cell_nLoc)
1075        do {
1076            (lastn)++; /* 3) continue tracking for nmax turns*/
1077            Cell_Pass(0L, pos - 1L, x1, lastpos); /* 1) check whether particle is stable at element from 0 to pos-1L*/
1078
1079            if (trace)
1080                fprintf(
1081                        outf1,
1082                        "%6ld %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e %+10.5e \n",
1083                        lastn, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
1084
1085            if (lastpos == pos - 1L)
1086                Cell_Pass(pos, globval.Cell_nLoc, x1, lastpos); /* 2) check particle is stable at element from pos to end*/
1087            //   Cell_Pass(pos-1L,globval.Cell_nLoc, x1, lastpos);
1088
1089        } while (((lastn) < nmax) && ((lastpos) == globval.Cell_nLoc));
1090
1091    if (lastpos == globval.Cell_nLoc)
1092        Cell_Pass(0L, pos, x1, lastpos);
1093
1094    if (lastpos != pos) {
1095        printf("Trac: Particle lost \n");
1096        fprintf(stdout, "turn:%6ld plane: %1d"
1097            " %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn,
1098                status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
1099    }
1100}
1101
1102/****************************************************************************/
1103/*bool chk_if_lost(double x0, double y0, double delta,
1104 long int nturn, bool floqs)
1105
1106 Purpose:
1107 Binary search for dynamical aperture in Floquet space.
1108
1109 Input:
1110 none
1111
1112 Output:
1113 none
1114
1115 Return:
1116 none
1117
1118 Global variables:
1119 px_0, py_0
1120
1121 Specific functions:
1122 chk_if_lost
1123
1124 Comments:
1125 none
1126
1127 ****************************************************************************/
1128
1129#define nfloq     4
1130bool chk_if_lost(double x0, double y0, double delta, long int nturn, bool floqs) {
1131 
1132    long int i=0L, lastn=0L, lastpos=0L;
1133    Vector x;
1134
1135    bool prt = false;
1136
1137    x[x_] = x0;
1138    x[px_] = px_0;
1139    x[y_] = y0;
1140    x[py_] = py_0;
1141    x[delta_] = delta;
1142    x[ct_] = 0.0;
1143    if (floqs)
1144        // transform to phase space
1145        LinTrans(nfloq, globval.Ascr, x);
1146    for (i = 0; i <= 3; i++)
1147        x[i] += globval.CODvect[i];
1148
1149    lastn = 0;
1150    if (prt) {
1151        printf("\n");
1152        printf("chk_if_lost:\n");
1153        printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn,
1154                1e3 * x[x_], 1e3 * x[px_], 1e3 * x[y_], 1e3 * x[py_], 1e2
1155                        * x[delta_], 1e3 * x[ct_]);
1156    }
1157    do {
1158        lastn++;
1159        if (globval.MatMeth)
1160            Cell_fPass(x, lastpos);
1161        else
1162            Cell_Pass(0, globval.Cell_nLoc, x, lastpos);
1163        if (prt)
1164            printf("%4ld %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f\n", lastn, 1e3
1165                    * x[x_], 1e3 * x[px_], 1e3 * x[y_], 1e3 * x[py_], 1e2
1166                    * x[delta_], 1e3 * x[ct_]);
1167    } while ((lastn != nturn) && (lastpos == globval.Cell_nLoc));
1168    return (lastn != nturn);
1169}
1170#undef nfloq
1171
1172/****************************************************************************/
1173/* void getdynap(double *r, double phi, double delta, double eps,
1174 int nturn, bool floqs)
1175
1176 Purpose:
1177 Binary search for dynamical aperture in Floquet space.
1178
1179
1180 Input:
1181 none
1182
1183 Output:
1184 none
1185
1186 Return:
1187 none
1188
1189 Global variables:
1190 none
1191
1192 Specific functions:
1193 chk_if_lost
1194
1195 Comments:
1196 none
1197
1198 ****************************************************************************/
1199void getdynap(double &r, double phi, double delta, double eps, int nturn,
1200        bool floqs) {
1201    /* Determine dynamical aperture by binary search. */
1202    double rmin = 0.0, rmax = r;
1203
1204    const bool prt = false;
1205    const double r_reset = 1e-3, r0 = 10e-3;
1206
1207    if (prt)
1208        printf("\n");
1209
1210    if (globval.MatMeth)
1211        Cell_Concat(delta);
1212    while (!chk_if_lost(rmax * cos(phi), rmax * sin(phi), delta, nturn, floqs)) {
1213        if (rmax < r_reset)
1214            rmax = r0;
1215        rmax *= 2.0;
1216    }
1217    while (rmax - rmin >= eps) {
1218        r = rmin + (rmax - rmin) / 2.0;
1219        if (prt)
1220            printf("getdynap: %6.3f %6.3f %6.3f\n", 1e3 * rmin, 1e3 * rmax, 1e3
1221                    * r);
1222        if (!chk_if_lost(r * cos(phi), r * sin(phi), delta, nturn, floqs))
1223            rmin = r;
1224        else
1225            rmax = r;
1226        if (rmin > rmax) {
1227            printf("getdynap: rmin > rmax\n");
1228            exit_(0);
1229        }
1230
1231    }
1232    r = rmin;
1233}
1234
1235/****************************************************************************/
1236/* void getcsAscr(void)
1237
1238 Purpose:
1239 Get Courant-Snyder Ascr
1240
1241
1242 Input:
1243 none
1244
1245 Output:
1246 none
1247
1248 Return:
1249 none
1250
1251 Global variables:
1252 none
1253
1254 Specific functions:
1255 none
1256
1257 Comments:
1258 none
1259
1260 ****************************************************************************/
1261void getcsAscr(void) {
1262    long i=0L, j=0L;
1263    double phi=0.0;
1264    Matrix R;
1265
1266    UnitMat(6, R);
1267    for (i = 1; i <= 2; i++) {
1268        phi = -atan2(globval.Ascr[i * 2 - 2][i * 2 - 1],
1269                globval.Ascr[i * 2 - 2][i * 2 - 2]);
1270        R[i * 2 - 2][i * 2 - 2] = cos(phi);
1271        R[i * 2 - 1][i * 2 - 1] = R[i * 2 - 2][i * 2 - 2];
1272        R[i * 2 - 2][i * 2 - 1] = sin(phi);
1273        R[i * 2 - 1][i * 2 - 2] = -R[i * 2 - 2][i * 2 - 1];
1274    }
1275    MulRMat(6, globval.Ascr, R);
1276    for (i = 1; i <= 2; i++) {
1277        if (globval.Ascr[i * 2 - 2][i * 2 - 2] < 0.0) {
1278            for (j = 0; j <= 5; j++) {
1279                globval.Ascr[j][i * 2 - 2] = -globval.Ascr[j][i * 2 - 2];
1280                globval.Ascr[j][i * 2 - 1] = -globval.Ascr[j][i * 2 - 1];
1281            }
1282        }
1283    }
1284    if (!InvMat(6, globval.Ascrinv))
1285        printf("  *** Ascr is singular\n");
1286}
1287
1288/****************************************************************************
1289 void dynap(FILE *fp, double r, const double delta, const double eps,
1290        const int npoint, const int nturn, double x[], double y[],
1291        const bool floqs, const bool print)
1292
1293 Purpose:
1294 Determine the dynamical aperture by tracking using polar coordinates,
1295 and sampling in phase.
1296 Assumes mid-plane symmetry
1297 Write the dynamic aperture to file "fp"
1298 
1299 Input:
1300 r initial guess
1301 delta     off momentum energy
1302 eps       precision for binary search
1303 npoint    sample number for phase coordinate
1304 nturn     number of turn for computing da
1305 x[]       horizontal dynamics aperture
1306 y[]       vertical dynamics aperture
1307 floqs     true means Floquet space
1308 print     true means Print out to the screen
1309
1310 Output:
1311
1312
1313 Return:
1314 none
1315
1316 Global variables:
1317 none
1318
1319 Specific functions:
1320 getdynap
1321
1322 Comments:
1323 none
1324
1325 ****************************************************************************/
1326void dynap(FILE *fp, double r, const double delta, const double eps,
1327        const int npoint, const int nturn, double x[], double y[],
1328        const bool floqs, const bool print)
1329
1330{
1331    /* Determine the dynamical aperture by tracking.
1332     Assumes mid-plane symmetry.                    */
1333
1334    long int i=0L, lastpos=0L;
1335    double phi=0.0, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0;
1336
1337    getcod(delta, lastpos);
1338    if (floqs) {
1339        Ring_GetTwiss(false, delta);
1340        if (print) {
1341            printf("\n");
1342            printf("Dynamical Aperture (Floquet space):\n");
1343            printf("     x^         y^\n");
1344            printf("\n");
1345        }
1346        fprintf(fp, "# Dynamical Aperture (Floquet space):\n");
1347        fprintf(fp, "#      x^         y^\n");
1348        fprintf(fp, "#\n");
1349    } else {
1350        if (print) {
1351            printf("\n");
1352            printf("Dynamical Aperture:\n");
1353            printf("     x      y\n");
1354            printf("    [mm]   [mm]\n");
1355            printf("\n");
1356        }
1357        fprintf(fp, "# Dynamical Aperture:\n");
1358        fprintf(fp, "#    x      y\n");
1359        fprintf(fp, "#   [mm]   [mm]\n");
1360        fprintf(fp, "#\n");
1361    }
1362
1363    x_min = 0.0;
1364    x_max = 0.0;
1365    y_min = 0.0;
1366    y_max = 0.0;
1367    for (i = 0; i < npoint; i++) {
1368        phi = i * M_PI / (npoint - 1);
1369        if (i == 0)
1370            phi = 1e-3;
1371        else if (i == npoint - 1)
1372            phi -= 1e-3;
1373        getdynap(r, phi, delta, eps, nturn, floqs);
1374        x[i] = r * cos(phi);
1375        y[i] = r * sin(phi);
1376        x_min = min(x[i], x_min);
1377        x_max = max(x[i], x_max);
1378        y_min = min(y[i], y_min);
1379        y_max = max(y[i], y_max);
1380        if (!floqs) {
1381            if (print)
1382                printf("  %6.2f %6.2f\n", 1e3 * x[i], 1e3 * y[i]);
1383            fprintf(fp, "  %6.2f %6.2f\n", 1e3 * x[i], 1e3 * y[i]);
1384        } else {
1385            if (print)
1386                printf("  %10.3e %10.3e\n", x[i], y[i]);
1387            fprintf(fp, "  %10.3e %10.3e\n", x[i], y[i]);
1388        }
1389        fflush(fp);
1390    }
1391
1392    if (print) {
1393        printf("\n");
1394        printf("  x^: %6.2f - %5.2f y^: %6.2f - %5.2f mm\n", 1e3 * x_min, 1e3
1395                * x_max, 1e3 * y_min, 1e3 * y_max);
1396    }
1397}
1398
1399/****************************************************************************/
1400/* double get_aper(int n, double x[], double y[])
1401
1402 Purpose:
1403
1404
1405 Input:
1406 none
1407
1408 Output:
1409 none
1410
1411 Return:
1412 none
1413
1414 Global variables:
1415 none
1416
1417 Specific functions:
1418 none
1419
1420 Comments:
1421 none
1422
1423 ****************************************************************************/
1424double get_aper(int n, double x[], double y[]) {
1425    int i=0;
1426    double A=0.0;
1427
1428    A = 0.0;
1429    for (i = 2; i <= n; i++)
1430        A += x[i - 2] * y[i - 1] - x[i - 1] * y[i - 2];
1431    A += x[n - 1] * y[0] - x[0] * y[n - 1];
1432    // x2 from mid-plane symmetry
1433    A = fabs(A);
1434    //  printf("\n");
1435    //  printf("  Dyn. Aper.: %5.1f mm^2\n", 1e6*A);
1436    return (A);
1437}
1438
1439/**************************************************************************
1440 * void GetTrack(const char *file_name, long *n, double x[], double px[],
1441        double y[], double py[])
1442 *
1443 *
1444 *
1445 **************************************************************************/
1446void GetTrack(const char *file_name, long *n, double x[], double px[],
1447        double y[], double py[]) {
1448    int k=0;
1449    char line[200];
1450    FILE *inf;
1451
1452    inf = file_read(file_name);
1453
1454    do {
1455        fgets(line, 200, inf);
1456    } while (strstr(line, "#") != NULL);
1457
1458    // skip initial conditions
1459    fgets(line, 200, inf);
1460
1461    do {
1462        sscanf(line, "%d", &k);
1463        sscanf(line, "%*d %lf %lf %lf %lf", &x[k - 1], &px[k - 1], &y[k - 1],
1464                &py[k - 1]);
1465    } while (fgets(line, 200, inf) != NULL);
1466
1467    *n = k;
1468
1469    fclose(inf);
1470}
1471
1472/****************************************************************************/
1473/* void Getj(long n, double *x, double *px, double *y, double *py)
1474
1475 Purpose:
1476 Calculates the linear invariant
1477
1478 Input:
1479 none
1480
1481 Output:
1482 none
1483
1484 Return:
1485 none
1486
1487 Global variables:
1488 none
1489
1490 Specific functions:
1491 none
1492
1493 Comments:
1494 none
1495
1496 ****************************************************************************/
1497void Getj(long n, double *x, double *px, double *y, double *py) {
1498    long int i=0L;
1499
1500    for (i = 0; i < n; i++) {
1501        x[i] = (pow(x[i], 2) + pow(px[i], 2)) / 2.0;
1502        y[i] = (pow(y[i], 2) + pow(py[i], 2)) / 2.0;
1503    }
1504}
1505
1506/****************************************************************************/
1507/* double GetArg(double x, double px, double nu)
1508
1509 Purpose:
1510 get argument of x
1511
1512 Input:
1513 none
1514
1515 Output:
1516 none
1517
1518 Return:
1519 none
1520
1521 Global variables:
1522 none
1523
1524 Specific functions:
1525 none
1526
1527 Comments:
1528 17/07/03 use M_PI instead of pi
1529
1530 ****************************************************************************/
1531double GetArg(double x, double px, double nu) {
1532   
1533  double phi=0.0, val=0.0;
1534
1535    phi = GetAngle(x, px);
1536    if (phi < 0.0)
1537        phi += 2.0 * M_PI;
1538    val = phi + Fract(nu) * 2.0 * M_PI;
1539    if (val < 0.0)
1540        val += 2.0 * M_PI;
1541    return val;
1542}
1543
1544/****************************************************************************/
1545/* void GetPhi(long n, double *x, double *px, double *y, double *py)
1546
1547 Purpose:
1548 get linear phases
1549
1550 Input:
1551 none
1552
1553 Output:
1554 none
1555
1556 Return:
1557 none
1558
1559 Global variables:
1560 none
1561
1562 Specific functions:
1563 none
1564
1565 Comments:
1566 none
1567
1568 ****************************************************************************/
1569void GetPhi(long n, double *x, double *px, double *y, double *py) {
1570    /* Calculates the linear phase */
1571    long i=0L;
1572
1573    for (i = 1; i <= n; i++) {
1574        x[i - 1] = GetArg(x[i - 1], px[i - 1], i * globval.TotalTune[0]);
1575        y[i - 1] = GetArg(y[i - 1], py[i - 1], i * globval.TotalTune[1]);
1576    }
1577}
1578
1579/*********************************/
1580/* Routines for Fourier analysis */
1581/*********************************/
1582
1583void Sinfft(int n, double xr[]) {
1584    /* DFT with sine window */
1585    int i=0;
1586    double xi[n];
1587
1588    for (i = 0; i < n; i++) {
1589        xr[i] = sin((double) i / (double) n * M_PI) * xr[i];
1590        xi[i] = 0.0;
1591    }
1592    FFT(n, xr, xi);
1593    for (i = 0; i < n; i++)
1594        xr[i] = sqrt(xr[i] * xr[i] + xi[i] * xi[i]);
1595}
1596
1597void sin_FFT(int n, double xr[]) {
1598    /* DFT with sine window */
1599    int i;
1600    double *xi;
1601
1602    xi = dvector(1, 2 * n);
1603
1604    for (i = 1; i <= n; i++) {
1605        xi[2 * i - 1] = sin((double) i / n * M_PI) * xr[i - 1];
1606        xi[2 * i] = 0.0;
1607    }
1608    dfour1(xi, (unsigned long) n, 1);
1609    for (i = 1; i <= n; i++)
1610        xr[i - 1] = sqrt(pow(xi[2 * i - 1], 2) + pow(xi[2 * i], 2)) * 2.0 / n;
1611
1612    free_dvector(xi, 1, 2 * n);
1613}
1614/*******************************************************************************
1615 *
1616 *
1617 *
1618 *
1619 ******************************************************************************/
1620void sin_FFT(int n, double xr[], double xi[]) {
1621    /* DFT with sine window */
1622    int i=0;
1623    double *xri;
1624
1625    xri = dvector(1, 2 * n);
1626
1627    for (i = 1; i <= n; i++) {
1628        xri[2 * i - 1] = sin((double) i / n * M_PI) * xr[i - 1];
1629        xri[2 * i] = sin((double) i / n * M_PI) * xi[i - 1];
1630    }
1631    dfour1(xri, (unsigned long) n, 1);
1632    for (i = 1; i <= n; i++) {
1633        xr[i - 1] = sqrt(pow(xri[2 * i - 1], 2) + pow(xri[2 * i], 2)) * 2.0 / n;
1634        xi[i - 1] = atan2(xri[2 * i], xri[2 * i - 1]);
1635    }
1636
1637    free_dvector(xri, 1, 2 * n);
1638}
1639/*******************************************************************************
1640 *
1641 *
1642 *
1643 *
1644 ******************************************************************************/
1645void GetInd(int n, int k, int *ind1, int *ind3) {
1646    if (k == 1) {
1647        *ind1 = 2;
1648        *ind3 = 2;
1649    } else if (k == n / 2 + 1) {
1650        *ind1 = n / 2;
1651        *ind3 = n / 2;
1652    } else {
1653        *ind1 = k - 1;
1654        *ind3 = k + 1;
1655    }
1656}
1657/*******************************************************************************
1658 *
1659 *
1660 *
1661 *
1662 ******************************************************************************/
1663void GetInd1(int n, int k, int *ind1, int *ind3) {
1664    if (k == 1) {
1665        *ind1 = 2;
1666        *ind3 = 2;
1667    } else if (k == n) {
1668        *ind1 = n - 1;
1669        *ind3 = n - 1;
1670    } else {
1671        *ind1 = k - 1;
1672        *ind3 = k + 1;
1673    }
1674}
1675/*******************************************************************************
1676 *
1677 *
1678 *
1679 *
1680 ******************************************************************************/
1681void GetPeak(int n, double *x, int *k) {
1682    /* Locate peak in DFT spectrum */
1683    int ind1=0, ind2=0, ind3=0;
1684    double peak=0.0;
1685
1686    *k = 1;
1687    for (ind2 = 1; ind2 <= n / 2 + 1; ind2++) {
1688        GetInd(n, ind2, &ind1, &ind3);
1689        if (x[ind2 - 1] > peak && x[ind1 - 1] < x[ind2 - 1] && x[ind3 - 1]
1690                < x[ind2 - 1]) {
1691            peak = x[ind2 - 1];
1692            *k = ind2;
1693        }
1694    }
1695}
1696/*******************************************************************************
1697 *
1698 *
1699 *
1700 *
1701 ******************************************************************************/
1702void GetPeak1(int n, double *x, int *k) {
1703    /* Locate peak in DFT spectrum */
1704    int ind1=0, ind2=0, ind3=0;
1705    double peak=0.0;
1706
1707    *k = 1;
1708    for (ind2 = 1; ind2 <= n; ind2++) {
1709        GetInd1(n, ind2, &ind1, &ind3);
1710        if (x[ind2 - 1] > peak && x[ind1 - 1] < x[ind2 - 1] && x[ind3 - 1]
1711                < x[ind2 - 1]) {
1712            peak = x[ind2 - 1];
1713            *k = ind2;
1714        }
1715    }
1716}
1717/*******************************************************************************
1718 *
1719 *
1720 *
1721 *
1722 ******************************************************************************/
1723double Int2snu(int n, double *x, int k) {
1724    /* Get frequency by nonlinear interpolation with two samples
1725     for sine window. The interpolation is:
1726
1727     1              2 A(k)       1
1728     nu = - [ k - 1 + ------------- - - ] ,      k-1 <= N nu <= k
1729     N           A(k-1) + A(k)   2
1730     */
1731    int ind=0, ind1=0, ind3=0;
1732    double ampl1=0.0, ampl2=0.0;
1733
1734    GetInd(n, k, &ind1, &ind3);
1735    if (x[ind3 - 1] > x[ind1 - 1]) {
1736        ampl1 = x[k - 1];
1737        ampl2 = x[ind3 - 1];
1738        ind = k;
1739    } else {
1740        ampl1 = x[ind1 - 1];
1741        ampl2 = x[k - 1];
1742        /* Interpolate in right direction for 0 frequency */
1743        if (k != 1)
1744            ind = ind1;
1745        else
1746            ind = 0;
1747    }
1748    /* Avoid division by zero */
1749    if (ampl1 + ampl2 != 0.0)
1750        return ((ind - 1 + 2 * ampl2 / (ampl1 + ampl2) - 0.5) / n);
1751    else
1752        return 0.0;
1753}
1754
1755/****************************************************************************/
1756/* double Sinc(double omega)
1757
1758 Purpose:
1759
1760
1761 Input:
1762 none
1763
1764 Output:
1765 none
1766
1767 Return:
1768 none
1769
1770 Global variables:
1771 none
1772
1773 Specific functions:
1774 none
1775
1776 Comments:
1777 zero test to be changed numericallywise
1778
1779 ****************************************************************************/
1780double Sinc(double omega) {
1781    /*  Function to calculate:
1782
1783     sin( omega )
1784     ------------
1785     omega
1786     */
1787    if (omega != 0.0)
1788        return (sin(omega) / omega);
1789    else
1790        return 1.0;
1791}
1792
1793/****************************************************************************/
1794/* double intsampl(long n, double *x, double nu, long k)
1795
1796 Purpose:
1797
1798
1799 Input:
1800 none
1801
1802 Output:
1803 none
1804
1805 Return:
1806 none
1807
1808 Global variables:
1809 none
1810
1811 Specific functions:
1812 none
1813
1814 Comments:
1815 17/07/03 use M_PI instead of pi
1816
1817 ****************************************************************************/
1818double intsampl(int n, double *x, double nu, int k) {
1819    /* Get amplitude by nonlinear interpolation for sine window. The
1820     distribution is given by:
1821
1822     1    sin pi ( k + 1/2 )     sin pi ( k - 1/2 )
1823     F(k) =  - ( -------------------- + -------------------- )
1824     2      pi ( k + 1/2 )          pi ( k - 1/2 )
1825     */
1826    double corr=0.0;
1827
1828    corr = (Sinc(M_PI * (k - 1 + 0.5 - nu * n)) + Sinc(M_PI * (k - 1 - 0.5 - nu
1829            * n))) / 2;
1830    return (x[k - 1] / corr);
1831}
1832
1833/****************************************************************************/
1834/* double linint(long n, long k, double nu, double *x)
1835
1836 Purpose:
1837
1838
1839 Input:
1840 none
1841
1842 Output:
1843 none
1844
1845 Return:
1846 none
1847
1848 Global variables:
1849 none
1850
1851 Specific functions:
1852 none
1853
1854 Comments:
1855 17/07/03 use M_PI instead of pi
1856
1857 ****************************************************************************/
1858double linint(int n, int k, double nu, double *x) {
1859    /* Get phase by linear interpolation for rectangular window
1860     with -pi <= phi <= pi */
1861    int i=0;
1862    double phi=0.0;
1863    double xr[n], xi[n];
1864
1865    for (i = 0; i < n; i++) {
1866        xr[i] = x[i];
1867        xi[i] = 0.0;
1868    }
1869    FFT(n, xr, xi);
1870    phi = GetAngle(xr[k - 1], xi[k - 1]) - (n * nu - k + 1) * M_PI;
1871    if (phi > M_PI)
1872        phi -= 2.0 * M_PI;
1873    else if (phi < -M_PI)
1874        phi += 2.0 * M_PI;
1875    return phi;
1876}
1877
1878/****************************************************************************/
1879/* void FndRes(struct LOC_findres *LINK)
1880
1881 Purpose:
1882
1883
1884 Input:
1885 none
1886
1887 Output:
1888 none
1889
1890 Return:
1891 none
1892
1893 Global variables:
1894 none
1895
1896 Specific functions:
1897 none
1898
1899 Comments:
1900 none
1901
1902 ****************************************************************************/
1903void FndRes(struct LOC_findres *LINK) {
1904    int i=0, j=0, FORLIM=0, FORLIM1=0;
1905    double delta=0.0;
1906
1907    FORLIM = LINK->n;
1908    for (i = 0; i <= FORLIM; i++) {
1909        FORLIM1 = LINK->n;
1910        for (j = -LINK->n; j <= FORLIM1; j++) {
1911            delta = fabs(i * LINK->nux + j * LINK->nuy);
1912            delta -= (int) delta;
1913            if (delta > 0.5)
1914                delta = 1 - delta;
1915            delta = fabs(delta - LINK->f);
1916            delta -= (int) delta;
1917            if (delta > 0.5)
1918                delta = 1 - delta;
1919            if (delta < LINK->eps) {
1920                if (abs(i) + abs(j) < LINK->n && (i != 0 || j >= 0)) {
1921                    LINK->found = true;
1922                    *LINK->nx = i;
1923                    *LINK->ny = j;
1924                }
1925            }
1926        }
1927    }
1928}
1929
1930/****************************************************************************/
1931/* void FindRes(long n_, double nux_, double nuy_, double f_,
1932 long *nx_, long *ny_)
1933
1934 Purpose:
1935
1936
1937 Input:
1938 none
1939
1940 Output:
1941 none
1942
1943 Return:
1944 none
1945
1946 Global variables:
1947 none
1948
1949 Specific functions:
1950 none
1951
1952 Comments:
1953 none
1954
1955 ****************************************************************************/
1956void FindRes(int n_, double nux_, double nuy_, double f_, int *nx_, int *ny_) {
1957    /* Match f by a linear combination of nux and nuy */
1958    struct LOC_findres V;
1959
1960    V.n = n_;
1961    V.nux = nux_;
1962    V.nuy = nuy_;
1963    V.f = f_;
1964    V.nx = nx_;
1965    V.ny = ny_;
1966    V.found = false;
1967    V.eps = 0.5e-6;
1968    do {
1969        V.eps = 10 * V.eps;
1970        FndRes(&V);
1971    } while (!V.found);
1972}
1973/*******************************************************************************
1974 *
1975 *
1976 *
1977 *
1978 ******************************************************************************/
1979void GetPeaks(int n, double *x, int nf, double *nu, double *A) {
1980    int i=0, k=0, ind1=0, ind3=0;
1981
1982    for (i = 0; i < nf; i++) {
1983        GetPeak(n, x, &k);
1984        nu[i] = Int2snu(n, x, k);
1985        A[i] = intsampl(n, x, nu[i], k);
1986        /* Make peak flat to allow for new call */
1987        GetInd(n, k, &ind1, &ind3);
1988        if (x[ind1 - 1] > x[ind3 - 1])
1989            x[k - 1] = x[ind1 - 1];
1990        else
1991            x[k - 1] = x[ind3 - 1];
1992    }
1993}
1994/*******************************************************************************
1995 *
1996 *
1997 *
1998 *
1999 ******************************************************************************/
2000void GetPeaks1(int n, double *x, int nf, double *nu, double *A) {
2001    int i=0, k=0, ind1=0, ind3=0;
2002
2003    for (i = 0; i < nf; i++) {
2004        GetPeak1(n, x, &k);
2005        nu[i] = Int2snu(n, x, k);
2006        A[i] = intsampl(n, x, nu[i], k);
2007        /* Make peak flat to allow for new call */
2008        GetInd1(n, k, &ind1, &ind3);
2009        if (x[ind1 - 1] > x[ind3 - 1])
2010            x[k - 1] = x[ind1 - 1];
2011        else
2012            x[k - 1] = x[ind3 - 1];
2013    }
2014}
2015
2016/*******************************/
2017/* Routines for magnetic error */
2018/*******************************/
2019
2020void SetTol(int Fnum, double dxrms, double dyrms, double drrms) {
2021    int i=0;
2022    long k=0L;
2023
2024    for (i = 1; i <= GetnKid(Fnum); i++) {
2025        k = Elem_GetPos(Fnum, i);
2026        Cell[k].Elem.M->PdSrms[X_] = dxrms;
2027        Cell[k].Elem.M->PdSrnd[X_] = normranf();
2028        Cell[k].Elem.M->PdSrms[Y_] = dyrms;
2029        Cell[k].Elem.M->PdSrnd[Y_] = normranf();
2030        Cell[k].Elem.M->PdTrms = drrms;
2031        Cell[k].Elem.M->PdTrnd = normranf();
2032        Mpole_SetdS(Fnum, i);
2033        Mpole_SetdT(Fnum, i);
2034    }
2035}
2036/*******************************************************************************
2037 *
2038 *
2039 *
2040 *
2041 ******************************************************************************/
2042void Scale_Tol(int Fnum, double dxrms, double dyrms, double drrms) {
2043    int Knum=0;
2044    long int loc=0L;
2045
2046    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
2047        loc = Elem_GetPos(Fnum, Knum);
2048        Cell[loc].Elem.M->PdSrms[X_] = dxrms;
2049        Cell[loc].Elem.M->PdSrms[Y_] = dyrms;
2050        Cell[loc].Elem.M->PdTrms = drrms;
2051        Mpole_SetdS(Fnum, Knum);
2052        Mpole_SetdT(Fnum, Knum);
2053    }
2054}
2055
2056/****************************************************************************/
2057/* void SetaTol(int Fnum, int Knum, double dx, double dy, double dr)
2058
2059 Purpose:
2060 Set a known random multipole displacement error
2061
2062 Input:
2063 none
2064
2065 Output:
2066 none
2067
2068 Return:
2069 none
2070
2071 Global variables:
2072 none
2073
2074 Specific functions:
2075 none
2076
2077 Comments:
2078 none
2079
2080 ****************************************************************************/
2081void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) {
2082    long int loc=0L;
2083
2084    loc = Elem_GetPos(Fnum, Knum);
2085    Cell[loc].Elem.M->PdSrms[0] = dx;
2086    Cell[loc].Elem.M->PdSrnd[0] = 1e0;
2087    Cell[loc].Elem.M->PdSrms[1] = dy;
2088    Cell[loc].Elem.M->PdSrnd[1] = 1e0;
2089    Cell[loc].Elem.M->PdTrms = dr;
2090    Cell[loc].Elem.M->PdTrnd = 1e0;
2091    Mpole_SetdS(Fnum, Knum);
2092    Mpole_SetdT(Fnum, Knum);
2093}
2094/*******************************************************************************
2095 *
2096 *
2097 *
2098 *
2099 ******************************************************************************/
2100void ini_aper(const double Dxmin, const double Dxmax, const double Dymin,
2101        const double Dymax) {
2102    int k=0;
2103
2104    for (k = 0; k <= globval.Cell_nLoc; k++) {
2105        Cell[k].maxampl[X_][0] = Dxmin;
2106        Cell[k].maxampl[X_][1] = Dxmax;
2107        Cell[k].maxampl[Y_][0] = Dymin;
2108        Cell[k].maxampl[Y_][1] = Dymax;
2109    }
2110}
2111/*******************************************************************************
2112 *
2113 *
2114 *
2115 *
2116 ******************************************************************************/
2117void set_aper(const int Fnum, const double Dxmin, const double Dxmax,
2118        const double Dymin, const double Dymax) {
2119    int i=0;
2120    long int loc=0L;
2121
2122    for (i = 1; i <= GetnKid(Fnum); i++) {
2123        loc = Elem_GetPos(Fnum, i);
2124        Cell[loc].maxampl[X_][0] = Dxmin;
2125        Cell[loc].maxampl[X_][1] = Dxmax;
2126        Cell[loc].maxampl[Y_][0] = Dymin;
2127        Cell[loc].maxampl[Y_][1] = Dymax;
2128    }
2129}
2130/*************************************************
2131 * void LoadApertures(const char *ChamberFileName)
2132 *
2133 *
2134 *
2135 **************************************************/
2136void LoadApertures(const char *ChamberFileName) {
2137    char line[128], FamName[32];
2138    long Fnum=0L;
2139    double Xmin=0.0, Xmax=0.0, Ymin=0.0, Ymax=0.0;
2140    FILE *ChamberFile;
2141
2142    ChamberFile = file_read(ChamberFileName);
2143
2144    do
2145        fgets(line, 128, ChamberFile);
2146    while (strstr(line, "#") != NULL);
2147
2148    do {
2149        sscanf(line, "%s %lf %lf %lf %lf", FamName, &Xmin, &Xmax, &Ymin, &Ymax);
2150        Fnum = ElemIndex(FamName);
2151        if (Fnum > 0)
2152            set_aper(Fnum, Xmin, Xmax, Ymin, Ymax);
2153    } while (fgets(line, 128, ChamberFile) != NULL);
2154
2155    fclose(ChamberFile);
2156}
2157/**********************************************************************************
2158 * void LoadTolerances(const char *TolFileName)
2159 *
2160 *
2161 *
2162 **********************************************************************************/
2163// Load tolerances from the file
2164void LoadTolerances(const char *TolFileName) {
2165    char line[128], FamName[32];
2166    int Fnum=0;
2167    double dx=0.0, dy=0.0, dr=0.0;
2168    FILE *tolfile;
2169
2170    tolfile = file_read(TolFileName);
2171    if(tolfile == NULL){
2172      printf("LoadTolerances(): Error! Failure to read file: %s \n",TolFileName);
2173      exit_(1);
2174    }
2175
2176    do
2177        fgets(line, 128, tolfile);
2178    while (strstr(line, "#") != NULL);
2179
2180    do {
2181        if (strstr(line, "#") == NULL) {
2182            sscanf(line, "%s %lf %lf %lf", FamName, &dx, &dy, &dr);
2183            Fnum = ElemIndex(FamName);
2184            if (Fnum > 0) {
2185                SetTol(Fnum, dx, dy, dr);
2186            } else {
2187                printf("LoadTolerances: undefined element %s\n", FamName);
2188                exit_(1);
2189            }
2190        }
2191    } while (fgets(line, 128, tolfile) != NULL);
2192
2193    fclose(tolfile);
2194}
2195
2196/**************************************************************************************
2197 * void ScaleTolerances(const char *TolFileName, const double scl)
2198 *
2199 *
2200 *
2201 * ************************************************************************************/
2202// Load tolerances from the file
2203void ScaleTolerances(const char *TolFileName, const double scl) {
2204    char line[128], FamName[32];
2205    int Fnum=0;
2206    double dx=0.0, dy=0.0, dr=0.0;
2207    FILE *tolfile;
2208
2209    tolfile = file_read(TolFileName);
2210    if(tolfile == NULL){
2211      printf("LoadTolerances(): Error! Failure to read file: %s \n",TolFileName);
2212      exit_(1);
2213    }
2214   
2215    do
2216        fgets(line, 128, tolfile);
2217    while (strstr(line, "#") != NULL);
2218
2219    do {
2220        if (strstr(line, "#") == NULL) {
2221            sscanf(line, "%s %lf %lf %lf", FamName, &dx, &dy, &dr);
2222            Fnum = ElemIndex(FamName);
2223            if (Fnum > 0) {
2224                Scale_Tol(Fnum, scl * dx, scl * dy, scl * dr);
2225            } else {
2226                printf("ScaleTolerances: undefined element %s\n", FamName);
2227                exit_(1);
2228            }
2229        }
2230    } while (fgets(line, 128, tolfile) != NULL);
2231    fclose(tolfile);
2232}
2233/*******************************************************************************
2234 *
2235 *
2236 *
2237 *
2238 ******************************************************************************/
2239void SetKpar(int Fnum, int Knum, int Order, double k) {
2240
2241    Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order + HOMmax] = k;
2242    Mpole_SetPB(Fnum, Knum, Order);
2243}
2244/*******************************************************************************
2245 *
2246 *
2247 *
2248 *
2249 ******************************************************************************/
2250void SetL(int Fnum, int Knum, double L) {
2251
2252    Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L;
2253}
2254/*******************************************************************************
2255 *
2256 *
2257 *
2258 *
2259 ******************************************************************************/
2260void SetL(int Fnum, double L) {
2261    int i=0;
2262
2263    for (i = 1; i <= GetnKid(Fnum); i++)
2264        Cell[Elem_GetPos(Fnum, i)].Elem.PL = L;
2265}
2266/*******************************************************************************
2267 *
2268 *
2269 *
2270 *
2271 ******************************************************************************/
2272void SetdKpar(int Fnum, int Knum, int Order, double dk) {
2273
2274    Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order + HOMmax] += dk;
2275    Mpole_SetPB(Fnum, Knum, Order);
2276}
2277/*******************************************************************************
2278 *
2279 *
2280 *
2281 *
2282 ******************************************************************************/
2283void SetKLpar(int Fnum, int Knum, int Order, double kL) {
2284    long int loc=0L;
2285
2286    if (abs(Order) > HOMmax) {
2287        printf("SetKLPar: Error!....Multipole Order %d  > HOMmax %d\n", Order,
2288                HOMmax);
2289        exit_(1);
2290    }
2291
2292    loc = Elem_GetPos(Fnum, Knum);
2293    if (Cell[loc].Elem.PL != 0e0)
2294        Cell[loc].Elem.M->PBpar[Order + HOMmax] = kL / Cell[loc].Elem.PL;
2295    else
2296        Cell[loc].Elem.M->PBpar[Order + HOMmax] = kL;
2297    Mpole_SetPB(Fnum, Knum, Order);
2298}
2299/*******************************************************************************
2300 *
2301 *
2302 *
2303 *
2304 ******************************************************************************/
2305void SetdKLpar(int Fnum, int Knum, int Order, double dkL) {
2306    long int loc=0L;
2307
2308    loc = Elem_GetPos(Fnum, Knum);
2309    if (Cell[loc].Elem.PL != 0e0)
2310        Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkL / Cell[loc].Elem.PL;
2311    else
2312        Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkL;
2313    Mpole_SetPB(Fnum, Knum, Order);
2314}
2315
2316/*************************************************************
2317 void SetdKrpar(int Fnum, int Knum, int Order, double dkrel)
2318 
2319   Purpose:
2320     Increase or reduce the strength of element by a certain scale.
2321     
2322   Input:
2323     Fnum           family number
2324     Knum           kid number
2325     order          field strength order to be modified
2326     bnr            scale of the field strength with order "order" 
2327   
2328   Output:
2329   
2330   Comments:
2331   
2332     
2333**************************************************************/
2334void SetdKrpar(int Fnum, int Knum, int Order, double dkrel) {
2335    long int loc=0L;
2336
2337    loc = Elem_GetPos(Fnum, Knum);
2338    if (Order == Dip && Cell[loc].Elem.M->Pthick == thick)
2339        Cell[loc].Elem.M->PBpar[Dip + HOMmax] += dkrel
2340                * Cell[loc].Elem.M->Pirho;
2341    else
2342        Cell[loc].Elem.M->PBpar[Order + HOMmax] += dkrel
2343                * Cell[loc].Elem.M->PBpar[Order + HOMmax];
2344    Mpole_SetPB(Fnum, Knum, Order);
2345}
2346/*******************************************************************************
2347 *
2348 *
2349 *
2350 *
2351 ******************************************************************************/
2352void Setbn(int Fnum, int order, double bn) {
2353    int i=0;
2354
2355    for (i = 1; i <= GetnKid(Fnum); i++)
2356        SetKpar(Fnum, i, order, bn);
2357}
2358/*******************************************************************************
2359 *
2360 *
2361 *
2362 *
2363 ******************************************************************************/
2364void SetbnL(int Fnum, int order, double bnL) {
2365    int i=0;
2366
2367    for (i = 1; i <= GetnKid(Fnum); i++)
2368        SetKLpar(Fnum, i, order, bnL);
2369}
2370/*******************************************************************************
2371 *
2372 *
2373 *
2374 *
2375 ******************************************************************************/
2376void Setdbn(int Fnum, int order, double dbn) {
2377    int i=0;
2378
2379    for (i = 1; i <= GetnKid(Fnum); i++)
2380        SetdKpar(Fnum, i, order, dbn);
2381}
2382/*******************************************************************************
2383 *
2384 *
2385 *
2386 *
2387 ******************************************************************************/
2388void SetdbnL(int Fnum, int order, double dbnL) {
2389    int i=0;
2390
2391    for (i = 1; i <= GetnKid(Fnum); i++) {
2392        SetdKLpar(Fnum, i, order, dbnL);
2393    }
2394}
2395
2396/*************************************************************
2397 void Setbnr(int Fnum, int order, double bnr)
2398 
2399   Purpose:
2400     Increase or reduce the strength of element family by a certain
2401     scale.
2402     
2403   Input:
2404     Fnum           family number
2405     order          field strength order to be modified
2406     bnr            scale of the field strength with order "order" 
2407   
2408   Output:
2409   
2410   Comments:
2411     Jianfeng Zhang  07/04/2011
2412     Fix the bug in the type definition of the function parameter 'order'
2413     from 'long' to 'int'.
2414     
2415**************************************************************/
2416//void Setbnr(int Fnum, long order, double bnr) {
2417void Setbnr(int Fnum, int order, double bnr) {
2418    int i=0;
2419
2420    for (i = 1; i <= GetnKid(Fnum); i++)
2421        SetdKrpar(Fnum, i, order, bnr);
2422}
2423/*******************************************************************************
2424 *
2425 *
2426 *
2427 *
2428 ******************************************************************************/
2429void SetbnL_sys(int Fnum, int Order, double bnL_sys) {
2430    int Knum=0;
2431    long int loc=0L;
2432
2433    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
2434        loc = Elem_GetPos(Fnum, Knum);
2435        if (Cell[loc].Elem.PL != 0.0)
2436            Cell[loc].Elem.M->PBsys[Order + HOMmax] = bnL_sys
2437                    / Cell[loc].Elem.PL;
2438        else
2439            Cell[loc].Elem.M->PBsys[Order + HOMmax] = bnL_sys;
2440        Mpole_SetPB(Fnum, Knum, Order);
2441    }
2442}
2443/*******************************************************************************
2444 *
2445 *
2446 *
2447 *
2448 ******************************************************************************/
2449void set_dbn_rel(const int type, const int n, const double dbn_rel) {
2450    long int j=0L;
2451    double dbn=0.0;
2452
2453    printf("\n");
2454    printf("Setting Db_%d/b_%d = %6.1e for:\n", n, type, dbn_rel);
2455    printf("\n");
2456    for (j = 0; j <= globval.Cell_nLoc; j++)
2457        if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) {
2458            printf("%s\n", Cell[j].Elem.PName);
2459            dbn = dbn_rel * Cell[j].Elem.M->PBpar[type + HOMmax];
2460            Cell[j].Elem.M->PBrms[n + HOMmax] = dbn;
2461            Cell[j].Elem.M->PBrnd[n + HOMmax] = normranf();
2462            Mpole_SetPB(Cell[j].Fnum, Cell[j].Knum, n);
2463        }
2464}
2465
2466/********************************************************************
2467 double GetKpar(int Fnum, int Knum, int Order)
2468
2469 Purpose:
2470 Return the n-th order design field strength of the element
2471
2472 Input:
2473 Fnum     family index
2474 Knum     kid index
2475 Order     design field strength
2476
2477 Ouput:
2478 None
2479
2480 Return:
2481 n-th order design field strength
2482
2483 *********************************************************************/
2484double GetKpar(int Fnum, int Knum, int Order) {
2485    return (Cell[Elem_GetPos(Fnum, Knum)].Elem.M->PBpar[Order + HOMmax]);
2486}
2487
2488/********************************************************************
2489 double GetL(int Fnum, int Knum)
2490
2491 Purpose:
2492 Return the length of the element with "Fnum" and "Knum"
2493
2494 Input:
2495 Fnum     family index
2496 Knum     kid index
2497
2498
2499 Ouput:
2500 None
2501
2502 Return:
2503
2504
2505 *********************************************************************/
2506double GetL(int Fnum, int Knum) {
2507    return (Cell[Elem_GetPos(Fnum, Knum)].Elem.PL);
2508}
2509
2510/********************************************************************
2511 double GetKLpar(int Fnum, int Knum, int Order)
2512
2513 Purpose:
2514 Return the n-th order design integrated field strength of the element
2515
2516 Input:
2517 Fnum     family index
2518 Knum     kid index
2519 Order     design field strength
2520
2521 Ouput:
2522 None
2523
2524 Return:
2525 n-th order design integrated field strength
2526
2527 *********************************************************************/
2528
2529double GetKLpar(int Fnum, int Knum, int Order) {
2530    long int loc = 0L;
2531
2532    loc = Elem_GetPos(Fnum, Knum);
2533    if (Cell[loc].Elem.PL != 0e0)
2534        return (Cell[loc].Elem.M->PBpar[Order + HOMmax] * Cell[loc].Elem.PL);
2535    else
2536        return (Cell[loc].Elem.M->PBpar[Order + HOMmax]);
2537}
2538/*******************************************************************************
2539 *
2540 *
2541 *
2542 *
2543 ******************************************************************************/
2544void SetdKLrms(int Fnum, int Order, double dkLrms) {
2545    long int Knum=0L, loc=0L;
2546
2547    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
2548        loc = Elem_GetPos(Fnum, Knum);
2549        if (Cell[loc].Elem.PL != 0e0)
2550            Cell[loc].Elem.M->PBrms[Order + HOMmax] = dkLrms
2551                    / Cell[loc].Elem.PL;
2552        else
2553            Cell[loc].Elem.M->PBrms[Order + HOMmax] = dkLrms;
2554        Cell[loc].Elem.M->PBrnd[Order + HOMmax] = normranf();
2555        Mpole_SetPB(Fnum, Knum, Order);
2556    }
2557}
2558/*******************************************************************************
2559 *
2560 *
2561 *
2562 *
2563 ******************************************************************************/
2564void Setdkrrms(int Fnum, int Order, double dkrrms) {
2565    long int Knum=0L, loc=0L;
2566
2567    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
2568        loc = Elem_GetPos(Fnum, Knum);
2569        if (Order == Dip && Cell[loc].Elem.M->Pthick == thick)
2570            Cell[loc].Elem.M->PBrms[Dip + HOMmax] = dkrrms
2571                    * Cell[loc].Elem.M->Pirho;
2572        else
2573            Cell[loc].Elem.M->PBrms[Order + HOMmax] = dkrrms
2574                    * Cell[loc].Elem.M->PBpar[Order + HOMmax];
2575        Cell[loc].Elem.M->PBrnd[Order + HOMmax] = normranf();
2576        Mpole_SetPB(Fnum, Knum, Order);
2577    }
2578}
2579/*******************************************************************************
2580 *
2581 *
2582 *
2583 *
2584 ******************************************************************************/
2585void SetKL(int Fnum, int Order) {
2586    long int Knum=0L;
2587
2588    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++)
2589        Mpole_SetPB(Fnum, Knum, Order);
2590}
2591/*******************************************************************************
2592 *
2593 *
2594 *
2595 *
2596 ******************************************************************************/
2597void set_dx(const int type, const double sigma_x, const double sigma_y) {
2598    long int j=0L;
2599
2600    printf("\n");
2601    printf("Setting sigma_x,y = (%6.1e, %6.1e) for b_%d:\n", sigma_x, sigma_y,
2602            type);
2603    printf("\n");
2604    for (j = 0; j <= globval.Cell_nLoc; j++)
2605        if ((Cell[j].Elem.Pkind == Mpole) && (Cell[j].Elem.M->n_design == type)) {
2606            printf("%s\n", Cell[j].Elem.PName);
2607            Cell[j].Elem.M->PdSrms[X_] = sigma_x;
2608            Cell[j].Elem.M->PdSrms[Y_] = sigma_y;
2609            Cell[j].Elem.M->PdSrnd[X_] = normranf();
2610            Cell[j].Elem.M->PdSrnd[Y_] = normranf();
2611            Mpole_SetdS(Cell[j].Fnum, Cell[j].Knum);
2612        }
2613}
2614/*******************************************************************************
2615 *
2616 *
2617 *
2618 *
2619 ******************************************************************************/
2620void SetBpmdS(int Fnum, double dxrms, double dyrms) {
2621    long int Knum, loc=0L;
2622
2623    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
2624        loc = Elem_GetPos(Fnum, Knum);
2625        Cell[loc].dS[X_] = normranf() * dxrms;
2626        Cell[loc].dS[Y_] = normranf() * dyrms;
2627    }
2628}
2629
2630/******************************************************************************
2631void codstat(double *mean, double *sigma, double *xmax, long lastpos, bool all)
2632
2633  Purpose:
2634     Routines for closed orbit correction
2635     Return the mean orbit, rms orbit and maximum orbit, based on the orbits at
2636     all lattice elements or all bpm postion.
2637 
2638  Input:
2639     mean           mean value of the orbit, horizontal or vertical
2640     sigma          rms value of the orbit, horizontal or vertical
2641     xmax           maximum value of the orbit, horizontal or vertical
2642     lastpos        last element index in the lattice
2643     all            true, then do statistics on the orbit at all elements
2644                    false, ...............................at all bpm
2645****************************************************************************/
2646void codstat(double *mean, double *sigma, double *xmax, long lastpos, bool all) {
2647    long i=0L, j=0L, n=0L;
2648    Vector2 sum, sum2;
2649    double TEMP=0.0;
2650
2651    n = 0;
2652    for (j = 0; j <= 1; j++) {
2653        sum[j] = 0.0;
2654        sum2[j] = 0.0;
2655        xmax[j] = 0.0;
2656    }
2657    for (i = 0; i <= lastpos; i++) {
2658        if (all || Cell[i].Fnum == globval.bpm) {//get the sum and max orbit at all elements or all bpm
2659            n++;
2660            for (j = 1; j <= 2; j++) {
2661                sum[j - 1] += Cell[i].BeamPos[j * 2 - 2];
2662                TEMP = Cell[i].BeamPos[j * 2 - 2];
2663                sum2[j - 1] += TEMP * TEMP;
2664                xmax[j - 1]
2665                        = max(xmax[j - 1], fabs(Cell[i].BeamPos[j * 2 - 2]));
2666            }
2667        }
2668    }
2669    for (j = 0; j <= 1; j++) {
2670        if (n != 0)
2671            mean[j] = sum[j] / n; //mean value of the orbit
2672        else
2673            mean[j] = 0.0;
2674        if (n != 0 && n != 1) {
2675            TEMP = sum[j];
2676            sigma[j] = (n * sum2[j] - TEMP * TEMP) / (n * (n - 1.0));
2677        } else
2678            sigma[j] = 0.0;
2679        if (sigma[j] >= 0.0)
2680            sigma[j] = sqrt(sigma[j]);
2681        else
2682            sigma[j] = 0.0;
2683    }
2684}
2685
2686/****************************************************************************/
2687/* void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos,
2688 long bpmdis[mnp])
2689
2690 Purpose:
2691 Get statistics for  closed orbit
2692
2693 Input:
2694 none
2695
2696 Output:
2697 none
2698
2699 Return:
2700 none
2701
2702 Global variables:
2703 none
2704
2705 Specific functions:
2706 none
2707
2708 Comments:
2709 none
2710
2711 ****************************************************************************/
2712void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos,
2713        long bpmdis[mnp]) {
2714    long i=0L, j=0L, m=0L, n=0L;
2715    Vector2 sum, sum2;
2716    double TEMP=0.0;
2717
2718    m = n = 0;
2719    for (j = 0; j <= 1; j++) {
2720        sum[j] = 0.0;
2721        sum2[j] = 0.0;
2722        xmax[j] = 0.0;
2723    }
2724
2725    for (i = 0; i <= lastpos; i++) {
2726        if (Cell[i].Fnum == globval.bpm) {
2727            if (!bpmdis[m]) {
2728                for (j = 1; j <= 2; j++) {
2729                    sum[j - 1] += Cell[i].BeamPos[j * 2 - 2];
2730                    TEMP = Cell[i].BeamPos[j * 2 - 2];
2731                    sum2[j - 1] += TEMP * TEMP;
2732                    xmax[j - 1] = max(xmax[j - 1], fabs(Cell[i].BeamPos[j * 2
2733                            - 2]));
2734                }
2735                n++;
2736            }
2737            m++;
2738        }
2739    }
2740    for (j = 0; j <= 1; j++) {
2741        if (n != 0)
2742            mean[j] = sum[j] / n;
2743        else
2744            mean[j] = 0.0;
2745        if (n != 0 && n != 1) {
2746            TEMP = sum[j];
2747            sigma[j] = (n * sum2[j] - TEMP * TEMP) / (n * (n - 1.0));
2748        } else
2749            sigma[j] = 0.0;
2750        if (sigma[j] >= 0.0)
2751            sigma[j] = sqrt(sigma[j]);
2752        else
2753            sigma[j] = 0.0;
2754    }
2755}
2756
2757/****************************************************************************/
2758/* double digitize(double x, double maxkick, double maxsamp)
2759
2760 Purpose:
2761 Map x onto the integer interval (-maxsamp ... maxsamp) where maxsamp
2762 corresponds maxkick.
2763
2764
2765 Input:
2766 none
2767
2768 Output:
2769 none
2770
2771 Return:
2772 none
2773
2774 Global variables:
2775 none
2776
2777 Specific functions:
2778 none
2779
2780 Comments:
2781 none
2782
2783 ****************************************************************************/
2784double digitize(double x, double maxkick, double maxsamp) {
2785    if (maxkick > 0.)
2786        if (maxsamp > 1.)
2787            return Sgn(x) * maxkick / maxsamp * min(floor(fabs(x) / maxkick
2788                    * maxsamp), maxsamp - 1.);
2789        else {
2790            return Sgn(x) * min(fabs(x), maxkick);
2791        }
2792    else
2793        return x;
2794}
2795
2796/****************************************************************************/
2797/* double digitize2(long plane, long inum, double x, double maxkick, double maxsamp)
2798
2799 Purpose:
2800
2801
2802 Input:
2803 none
2804
2805 Output:
2806 none
2807
2808 Return:
2809 none
2810
2811 Global variables:
2812 none
2813
2814 Specific functions:
2815 none
2816
2817 Comments:
2818 none
2819
2820 ****************************************************************************/
2821svdarray xmemo[2];
2822
2823double digitize2(long plane, long inum, double x, double maxkick,
2824        double maxsamp) {
2825    double xint;
2826
2827    if (maxkick > 0.)
2828        if (maxsamp > 1.) {
2829            xint = min(floor(fabs(x) / maxkick * maxsamp), maxsamp - 1.);
2830
2831            if (fabs(xint - xmemo[inum][plane]) >= 1.) {
2832                xmemo[inum][plane] = xint;
2833            } else {
2834                xmemo[inum][plane] += 0.1;
2835                xint = min(xmemo[inum][plane], maxsamp - 1.);
2836            }
2837            return Sgn(x) * maxkick / maxsamp * xint;
2838        } else {
2839            return Sgn(x) * min(fabs(x), maxkick);
2840        }
2841    else
2842        return x;
2843}
2844
2845// MATH ROUTINE a mettre dans mathlib.c
2846
2847/****************************************************************************/
2848/*  void GetMean(n, x)
2849
2850 Purpose:
2851 Get out the mean value of vector x
2852
2853 Input:
2854 n vector size
2855 x vector to get out the mean value
2856
2857 Output:
2858 none
2859
2860 Return:
2861 none
2862
2863 Global variables:
2864 none
2865
2866 Specific functions:
2867 none
2868
2869 Comments:
2870 to be moved in mathlib
2871
2872 ****************************************************************************/
2873void GetMean(long n, double *x) {
2874    long i=0L;
2875    double mean = 0e0;
2876
2877    if (n < 1) {
2878        fprintf(stdout, "GetMean: error wrong vector size n=%ld\n", n);
2879        exit_(1);
2880    }
2881    for (i = 0; i < n; i++)
2882        mean += x[i];
2883    mean /= n;
2884    for (i = 0; i < n; i++)
2885        x[i] = x[i] - mean;
2886}
2887
2888/****************************************************************************/
2889/* double Fract(double x)
2890
2891 Purpose:
2892 Gets fractional part of x
2893
2894 Input:
2895 none
2896
2897 Output:
2898 none
2899
2900 Return:
2901 none
2902
2903 Global variables:
2904 none
2905
2906 Specific functions:
2907 none
2908
2909 Comments:
2910 none
2911
2912 ****************************************************************************/
2913double Fract(double x) {
2914    return (x - (long int) x);
2915}
2916
2917/****************************************************************************/
2918/* double Sgn (double x)
2919
2920 Purpose:
2921 Gets sign of x
2922
2923 Input:
2924 none
2925
2926 Output:
2927 0  if zero
2928 1  if positive
2929 -1 if negative
2930
2931 Return:
2932 none
2933
2934 Global variables:
2935 none
2936
2937 Specific functions:
2938 none
2939
2940 Comments:
2941 none
2942
2943 ****************************************************************************/
2944double Sgn(double x) {
2945    return (x == 0.0 ? 0.0 : (x > 0.0 ? 1.0 : -1.0));
2946}
2947
2948
2949
2950/*************************************************************************/
2951/*void PrintCh(void)
2952 
2953 Purpose:
2954 Output vacuum chamber limitation at each element to file "chambre.out"
2955
2956 Input:
2957 none
2958
2959 Output:
2960 none
2961
2962 Return:
2963 none
2964
2965 Global variables:
2966 none
2967
2968 Specific functions:
2969 none
2970
2971 Comments:
2972 none
2973
2974
2975 *************************************************************************/
2976void PrintCh(void) {
2977    long i = 0;
2978    struct tm *newtime;
2979    FILE *f;
2980
2981    const char *fic = "chambre.out";
2982
2983    newtime = GetTime();
2984
2985    f = file_write(fic);
2986    fprintf(f, "# TRACY III Synchrotron SOLEIL -- %s -- %s \n", fic, asctime2(
2987            newtime));
2988    fprintf(f,
2989            "#  i  name               s    -xch(mm) +xch(mm)  -ych(mm) +ych(mm)\n#\n");
2990
2991    for (i = 1; i <= globval.Cell_nLoc; i++)
2992        fprintf(f, "%4ld  %15s %6.2f  %7.3f  %7.3f  %7.3f %7.3f\n", i,
2993                Cell[i].Elem.PName, Cell[i].S, Cell[i].maxampl[X_][0] * 1E3,
2994                Cell[i].maxampl[X_][1] * 1E3, Cell[i].maxampl[Y_][0] * 1E3,
2995                Cell[i].maxampl[Y_][1] * 1E3);
2996
2997    fclose(f);
2998}
2999
3000
3001/****************************************************************************/
3002/* void GetChromTrac(long Nb, long Nbtour, double emax,
3003 double *xix, double *xiz)
3004
3005 Purpose:
3006 Computes chromaticities by tracking
3007
3008 Input:
3009 Nb      point number
3010 Nbtour  turn number
3011 emax    energy step
3012
3013 Output:
3014 xix horizontal chromaticity
3015 xiz vertical chromaticity
3016
3017 Return:
3018 none
3019
3020 Global variables:
3021 trace
3022
3023 Specific functions:
3024 Trac_Simple, Get_NAFF
3025
3026 Comments:
3027 27/04/03 chromaticities are now output arguments
3028 07/10/10 add test if unstable
3029
3030 ****************************************************************************/
3031#define nterm  2
3032#define ZERO 1E-8
3033
3034void GetChromTrac(long Nb, long Nbtour, double emax, double *xix, double *xiy) {
3035    bool status = true;
3036    int nb_freq[2] = { 0, 0 }; /* frequency number to look for */
3037    int i = 0;
3038    double Tab[6][NTURN], fx[nterm], fy[nterm], nux1, nux2, nuy1, nuy2;
3039
3040    double x  = 1e-6, xp  = 0.0, y  = 1e-6, yp  = 0.0;
3041    double x0 = 1e-6, xp0 = 0.0, y0 = 1e-6, yp0 = 0.0;
3042    //double xixExtra = 0.0, xizExtra= 0.0, xixhalf= 0.0, xizhalf= 0.0;
3043    //double nux3 = 0.0, nux4 = 0.0, nuz3 = 0.0, nuz4 = 0.0;
3044
3045    /* initializations */
3046    for (i = 0; i < nterm; i++) {
3047        fx[i] = 0.0;
3048        fy[i] = 0.0;
3049    }
3050    /* end init */
3051
3052    /* Tracking for delta = emax and computing tunes */
3053    x  = x0;
3054    xp = xp0;
3055    y  = y0;
3056    yp = yp0;
3057
3058    Trac_Simple4DCOD(x, xp, y, yp, emax, 0.0, Nbtour, Tab, &status);
3059    if (status){
3060    Get_NAFF(nterm, Nbtour, Tab, fx, fy, nb_freq);
3061    nux1 = (fabs(fx[0]) > ZERO ? fx[0] : fx[1]);
3062    nuy1 = fy[0];}
3063    else{
3064        nux1=999; nuy1=999;
3065    }
3066
3067    if (trace)
3068        fprintf(stdout,
3069                "\nGetChromTrac: Entering routine for chroma using tracking\n");
3070    if (trace)
3071        fprintf(stdout, "emax= % 10.6e nux1=% 10.6e nuy1= % 10.6e\n", emax,
3072                nux1, nuy1);
3073
3074    /* Tracking for delta = -emax and computing tunes */
3075    x  = x0;
3076    xp = xp0;
3077    y  = y0;
3078    yp = yp0;
3079
3080    Trac_Simple4DCOD(x, xp, y, yp, -emax, 0.0, Nbtour, Tab, &status);
3081    if (status){
3082        Get_NAFF(nterm, Nbtour, Tab, fx, fy, nb_freq);
3083    if (trace)
3084        fprintf(stdout, "nturn=%6ld x=% 10.5g xp=% 10.5g z=% 10.5g zp=% 10.5g"
3085            " delta=% 10.5g ctau=% 10.5g \n", Nbtour, Tab[0][Nbtour - 1],
3086                Tab[1][Nbtour - 1], Tab[2][Nbtour - 1], Tab[3][Nbtour - 1],
3087                Tab[4][Nbtour - 1], Tab[5][Nbtour - 1]);
3088
3089    nux2 = (fabs(fx[0]) > ZERO ? fx[0] : fx[1]);
3090    nuy2 = fy[0];
3091
3092    if (trace)
3093        fprintf(stdout, "emax= % 10.6e nux2= % 10.6e nuy2= % 10.6e\n", -emax,
3094                nux2, nuy2);
3095
3096    /* Computing chromaticities */
3097    *xix = (nux2 - nux1) * 0.5 / emax;
3098    *xiy = (nuy2 - nuy1) * 0.5 / emax;
3099
3100    if (trace)
3101        fprintf(stdout,
3102                "GetChromTrac: Exiting routine for chroma using tracking\n\n");
3103    }
3104    else{ // unstable
3105        *xix = -99999;
3106        *xiy = -99999;
3107    }
3108
3109    /*
3110     // Compute for half a step to diagnose precision
3111     Trac_Simple(x, xp, z, zp, emax*0.5, 0.0, Nbtour, Tab, &status);
3112     Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
3113     nux3 = (fabs (fx[0]) > 1e-8 ? fx[0] : fx[1]); nuz3 = fz[0];
3114
3115     Trac_Simple(x, xp, z, zp, -emax*0.5, 0.0, Nbtour, Tab, &status);
3116     Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
3117     nux4 = (fabs(fx[0]) > 1e-8 ? fx[0] : fx[1]); nuz4 = fz[0];
3118
3119     xixhalf = (nux4-nux3)/emax; xizhalf = (nuz4-nuz3)/emax;
3120
3121     // Richardson extrapolation
3122     xixExtra = (4.0*xixhalf-*xix)/3.0;
3123     xizExtra = (4.0*xizhalf-*xiz)/3.0;
3124
3125     fprintf(stdout, "chroma evaluated at +/- %6.2g, xix = % f xiz = % f\n",
3126     emax, *xix, *xiz);
3127     fprintf(stdout, "chroma evaluated at +/- %6.2g, xix = % f xiz = % f\n",
3128     emax/2, xixhalf, xizhalf);
3129     fprintf(stdout, "chroma evaluated from Richardson Extrapolation, xix = % f xiz = % f\n",
3130     xixExtra, xizExtra);
3131     */
3132}
3133#undef nterm
3134#undef ZERO
3135
3136/****************************************************************************/
3137/* void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz)
3138
3139 Purpose:
3140 Computes tunes by tracking
3141
3142 Input:
3143 Nb      point number
3144 Nbtour  turn number
3145 emax    energy step
3146
3147 Output:
3148 nux horizontal tune (0.0 if unstable)
3149 nuz vertical tune (0.0 if unstable)
3150
3151 Return:
3152 none
3153
3154 Global variables:
3155 trace
3156
3157 Specific functions:
3158 Trac_Simple, Get_NAFF
3159
3160 Comments:
3161   Add test on stability
3162
3163 ****************************************************************************/
3164#define nterm  2
3165#define ZERO 1E-8
3166void GetTuneTrac(long Nbtour, double emax, double *nux, double *nuz) {
3167    double Tab[6][NTURN], fx[nterm], fz[nterm];
3168    int nb_freq[2];
3169    bool status;
3170
3171    double x = 1e-6, xp = 0.0, z = 1e-6, zp = 0.0;
3172
3173    Trac_Simple4DCOD(x, xp, z, zp, emax, 0.0, Nbtour, Tab, &status);
3174    if (status){
3175       Get_NAFF(nterm, Nbtour, Tab, fx, fz, nb_freq);
3176        *nux = (fabs(fx[0]) > ZERO ? fx[0] : fx[1]);
3177        *nuz = fz[0];}
3178    else{ // particle unstable
3179        *nux = 0.0;
3180        *nuz = 0.0;}
3181}
3182#undef nterm
3183#undef ZERO
3184
3185/****************************************************************************/
3186/* void TransTwiss(double *alpha, double *beta, double *eta, double *etap, double *codvect)
3187
3188 Purpose: high level application
3189 Calculate Twiss functions for a transport line
3190
3191 Input:
3192 alpha   alpha functions at the line entrance
3193 beta    beta functions at the line entrance
3194 eta     dispersion functions at the line entrance
3195 etap    dispersion derivatives functions at the line entrance
3196 codvect closed orbit functions at the line entrance
3197
3198 Output:
3199 none
3200
3201 Return:
3202 none
3203
3204 Global variables:
3205
3206
3207 Specific functions:
3208 TransTrace
3209
3210 Comments:
3211 redundant with ttwiss
3212
3213 ****************************************************************************/
3214void TransTwiss(Vector2 &alpha, Vector2 &beta, Vector2 &eta, Vector2 &etap,
3215        Vector &codvect) {
3216    TransTrace(0, globval.Cell_nLoc, alpha, beta, eta, etap, codvect);
3217}
3218
3219/****************************************************************************/
3220/* void ttwiss(double *alpha, double *beta, double *eta, double *etap, double dP)
3221
3222 Purpose:
3223 Calculate Twiss functions for transport line
3224
3225 Input:
3226 none
3227
3228 Output:
3229 none
3230
3231 Return:
3232 none
3233
3234 Global variables:
3235 none
3236
3237 Specific functions:
3238 none
3239
3240 Comments:
3241 redundant with TransTwiss
3242
3243 ****************************************************************************/
3244void ttwiss(const Vector2 &alpha, const Vector2 &beta, const Vector2 &eta,
3245        const Vector2 &etap, const double dP) {
3246    TraceABN(0, globval.Cell_nLoc, alpha, beta, eta, etap, dP);
3247}
3248
3249/****************************************************************************/
3250/* void findcodS(double dP)
3251
3252 Purpose:
3253 Search for the closed orbit using a numerical method
3254 Algo: Newton_Raphson method
3255 Quadratic convergence
3256 May need a guess starting point
3257 Simple precision algorithm
3258
3259 Input:
3260 dP energy offset
3261
3262 Output:
3263 none
3264
3265 Return:
3266 none
3267
3268 Global variables:
3269 none
3270
3271 specific functions:
3272 Newton_Raphson
3273
3274 Comments:
3275 Method introduced because of bad convergence of da for ID using RADIA maps
3276
3277 ****************************************************************************/
3278void findcodS(double dP) {
3279    double *vcod;
3280    Vector x0;
3281    const int ntrial = 40; // maximum number of trials for closed orbit
3282    const double tolx = 1e-8; // numerical precision
3283    int k=0;
3284    int dim=0; // 4D or 6D tracking
3285    long lastpos=0L;
3286
3287    vcod = dvector(1, 6);
3288
3289    // starting point
3290    for (k = 1; k <= 6; k++)
3291        vcod[k] = 0.0;
3292
3293    vcod[5] = dP; // energy offset
3294
3295    if (globval.Cavity_on) {
3296        dim = 6; /* 6D tracking*/
3297        fprintf(stdout, "Error looking for cod in 6D\n");
3298        exit_(1);
3299    } else {
3300        dim = 4; /* 4D tracking */
3301        vcod[1] = Cell[0].Eta[0] * dP;
3302        vcod[2] = Cell[0].Etap[0] * dP;
3303        vcod[3] = Cell[0].Eta[1] * dP;
3304        vcod[4] = Cell[0].Etap[1] * dP;
3305    }
3306
3307    Newton_RaphsonS(ntrial, vcod, dim, tolx);
3308
3309    if (status.codflag == false)
3310        fprintf(stdout, "Error No COD found\n");
3311    if (trace) {
3312        for (k = 1; k <= 6; k++)
3313            x0[k - 1] = vcod[k];
3314        fprintf(stdout, "Before cod % .5e % .5e % .5e % .5e % .5e % .5e \n",
3315                x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]);
3316        Cell_Pass(0, globval.Cell_nLoc, x0, lastpos);
3317        fprintf(stdout, "After  cod % .5e % .5e % .5e % .5e % .5e % .5e \n",
3318                x0[0], x0[1], x0[2], x0[3], x0[4], x0[5]);
3319        Cell_Pass(0, globval.Cell_nLoc, x0, lastpos);
3320    }
3321    free_dvector(vcod, 1, 6);
3322}
3323
3324/****************************************************************************/
3325/* void findcod(double dP)
3326
3327 Purpose:
3328 Search for the closed orbit using a numerical method
3329 Algo: Newton_Raphson method
3330 Quadratic convergence
3331 May need a guess starting point
3332 Simple precision algorithm
3333 4D
3334 Starting point: linear closed orbit
3335
3336 6D
3337 Starting point: zero
3338 if radiation on : x[5] is the synchroneous phase (equilibrium RF phase)
3339 off: x[5] is zero
3340
3341 Input:
3342 dP energy offset
3343
3344 Output:
3345 none
3346
3347 Return:
3348 vcod:  6-D closed orbit
3349
3350 Global variables:
3351 none
3352
3353 specific functions:
3354 Newton_Raphson
3355
3356 Comments:
3357 Method introduced because of bad convergence of da for ID
3358 using RADIA maps
3359
3360 ****************************************************************************/
3361void findcod(double dP) {
3362    Vector vcod;
3363    const int ntrial = 40; // maximum number of trials for closed orbit
3364    const double tolx = 1e-10; // numerical precision
3365    int k=0, dim = 0;
3366    long lastpos=0L;
3367
3368    // initializations
3369    for (k = 0; k <= 5; k++)
3370        vcod[k] = 0.0;
3371
3372    if (globval.Cavity_on) {
3373        fprintf(stdout, "warning looking for cod in 6D\n");
3374        dim = 6;
3375    } else { // starting point linear closed orbit
3376        dim = 4;
3377        vcod[0] = Cell[0].Eta[0] * dP;
3378        vcod[1] = Cell[0].Etap[0] * dP;
3379        vcod[2] = Cell[0].Eta[1] * dP;
3380        vcod[3] = Cell[0].Etap[1] * dP;
3381        vcod[4] = dP; // energy offset
3382    }
3383
3384    Newton_Raphson(dim, vcod, ntrial, tolx);
3385
3386    if (status.codflag == false)
3387        fprintf(stdout, "Error No COD found\n");
3388
3389    CopyVec(6, vcod, globval.CODvect); // save closed orbit at the ring entrance
3390
3391    if (trace) {
3392        fprintf(stdout, "Before cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n",
3393                vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]);
3394        Cell_Pass(0, globval.Cell_nLoc, vcod, lastpos);
3395        fprintf(stdout, "After  cod2 % .5e % .5e % .5e % .5e % .5e % .5e \n",
3396                vcod[0], vcod[1], vcod[2], vcod[3], vcod[4], vcod[5]);
3397    }
3398}
3399/****************************************************************************/
3400/* void computeFandJS(double *x, int n, double **fjac, double *fvect)
3401
3402 Purpose:
3403 Simple precision algo
3404 Tracks x over one turn. And computes the Jacobian matrix of the
3405 transformation by numerical differentiation.
3406 using forward difference formula : faster but less accurate
3407 using symmetric difference formula
3408
3409 Input:
3410 x vector for evaluation
3411 n dimension 4 or 6
3412
3413 Output:
3414 fvect transport of x over one turn
3415 fjac  Associated jacobian matrix
3416
3417 Return:
3418 none
3419
3420 Global variables:
3421 none
3422
3423 specific functions:
3424 none
3425
3426 Comments:
3427 none
3428
3429 ****************************************************************************/
3430
3431void computeFandJS(double *x, int n, double **fjac, double *fvect) {
3432    int i=0, k=0;
3433    long lastpos = 0L;
3434    Vector x0, fx, fx1, fx2;
3435
3436    const double deps = 1e-8; //stepsize for numerical differentiation
3437
3438    for (i = 1; i <= 6; i++)
3439        x0[i - 1] = x[i];
3440
3441    Cell_Pass(0, globval.Cell_nLoc, x0, lastpos);
3442
3443    for (i = 1; i <= n; i++) {
3444        fvect[i] = x0[i - 1];
3445        fx[i - 1] = x0[i - 1];
3446    }
3447
3448    // compute Jacobian matrix by numerical differentiation
3449    for (k = 0; k < n; k++) {
3450        for (i = 1; i <= 6; i++)
3451            x0[i - 1] = x[i];
3452        x0[k] += deps; // differential step in coordinate k
3453
3454        Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring
3455        for (i = 1; i <= 6; i++)
3456            fx1[i - 1] = x0[i - 1];
3457
3458        for (i = 1; i <= 6; i++)
3459            x0[i - 1] = x[i];
3460        x0[5] = 0.0;
3461        x0[k] -= deps; // differential step in coordinate k
3462
3463        Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring
3464        for (i = 1; i <= 6; i++)
3465            fx2[i - 1] = x0[i - 1];
3466
3467        for (i = 1; i <= n; i++) // symmetric difference formula
3468            fjac[i][k + 1] = 0.5 * (fx1[i - 1] - fx2[i - 1]) / deps;
3469        //~ for (i = 1; i <= n; i++) // forward difference formula
3470        //~ fjac[i][k + 1] = (float) ((x0[i - 1] - fx[i - 1]) / deps);
3471    }
3472}
3473
3474/****************************************************************************/
3475/* void computeFand(int n, float *x, float **fjac, float *fvect)
3476
3477 Purpose:
3478 Tracks x over one turn. And computes the Jacobian matrix of the
3479 transformation by numerical differentiation.
3480 using symmetric difference formula
3481 double precision algorithm
3482
3483 Input:
3484 x vector for evaluation
3485
3486 Output:
3487 fvect transport of x over one turn
3488 fjac  Associated jacobian matrix
3489
3490 Return:
3491 none
3492
3493 Global variables:
3494 none
3495
3496 specific functions:
3497 none
3498
3499 Comments:
3500 none
3501
3502 ****************************************************************************/
3503void computeFandJ(int n, Vector &x, Matrix &fjac, Vector &fvect) {
3504    int i=0, k=0;
3505    long lastpos = 0;
3506    Vector x0, fx1, fx2;
3507
3508    const double deps = 1e-8; //stepsize for numerical differentiation
3509
3510    CopyVec(6, x, x0);
3511
3512    Cell_Pass(0, globval.Cell_nLoc, x0, lastpos);
3513    CopyVec(n, x0, fvect);
3514
3515    // compute Jacobian matrix by numerical differentiation
3516    for (k = 0; k < n; k++) {
3517        CopyVec(6L, x, x0);
3518        x0[k] += deps; // differential step in coordinate k
3519
3520        Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring
3521        CopyVec(6L, x0, fx1);
3522
3523        CopyVec(6L, x, x0);
3524        x0[k] -= deps; // differential step in coordinate k
3525
3526        Cell_Pass(0, globval.Cell_nLoc, x0, lastpos); // tracking along the ring
3527        CopyVec(6L, x0, fx2);
3528
3529        for (i = 0; i < n; i++) // symmetric difference formula
3530            fjac[i][k] = 0.5 * (fx1[i] - fx2[i]) / deps;
3531    }
3532}
3533
3534/****************************************************************************/
3535/* void Newton_RaphsonS(int ntrial,double x[],int n,double tolx, double tolf)
3536 
3537 Purpose:
3538 Newton_Rapson algorithm from Numerical Recipes
3539 single precision algorithm
3540 Robustess: quadratic convergence
3541 Hint: for n-dimensional problem, the algo can be stuck on local minimum
3542 In this case, it should be enough to provide a resonable starting
3543 point.
3544
3545 Method:
3546 look for closed orbit solution of f(x) = x
3547 This problems is equivalent to finding the zero of g(x)= f(x) - x
3548 g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h)
3549 Then at first order we solve h:
3550 h = - inverse(Jacobian(f) -Id) * (f(x)-x)
3551 the new guess is then xnew = x + h
3552 By iteration, this converges quadratically.
3553
3554 The algo is stopped whenever  |x -xnew| < tolx
3555
3556 f(x) is computes by tracking over one turn
3557 Jacobian(f) is computed numerically by numerical differentiation
3558 These two operations are provided by the function computeFandJ
3559
3560 Input:
3561 ntrial number of iterations for closed zero search
3562 n number of dimension 4 or 6
3563 x intial guess for the closed orbit
3564 tolx tolerance over the solution x
3565 tolf tolerance over the evalution f(x)
3566
3567 Output:
3568 x closed orbit
3569
3570 Return:
3571 none
3572
3573 Global variables:
3574 status
3575
3576 specific functions:
3577 computeFandJS
3578 ludcmp,lubksb
3579
3580 Comments:
3581 none
3582
3583 ****************************************************************************/
3584
3585void Newton_RaphsonS(int ntrial, double x[], int n, double tolx) {
3586    int k=0, i=0, *indx;
3587    double errx=0.0, d=0.0, *bet, *fvect, **alpha;
3588
3589    errx = 0.0;
3590    // NR arrays start from 1 and not 0 !!!
3591    indx = ivector(1, n);
3592    bet = dvector(1, n);
3593    fvect = dvector(1, n);
3594    alpha = dmatrix(1, n, 1, n);
3595
3596    for (k = 1; k <= ntrial; k++) { // loop over number of iterations
3597        // supply function values at x in fvect and Jacobian matrix in fjac
3598        computeFandJS(x, n, alpha, fvect);
3599
3600        // Jacobian -Id
3601        for (i = 1; i <= n; i++)
3602            alpha[i][i] -= 1.0;
3603        for (i = 1; i <= n; i++)
3604            bet[i] = x[i] - fvect[i]; // right side of linear equation
3605        // solve linear equations using LU decomposition using NR routines
3606        dludcmp(alpha, n, indx, &d);
3607        dlubksb(alpha, n, indx, bet);
3608        errx = 0.0; // check root convergence
3609        for (i = 1; i <= n; i++) { // update solution
3610            errx += fabs(bet[i]);
3611            x[i] += bet[i];
3612        }
3613
3614        if (trace)
3615            fprintf(
3616                    stdout,
3617                    "%02d: cod % .5e % .5e % .5e % .5e % .5e % .5e  errx =% .5e\n",
3618                    k, x[1], x[2], x[3], x[4], x[5], x[6], errx);
3619        if (errx <= tolx) {
3620            status.codflag = true;
3621            break;
3622        }
3623    }
3624    // check whenver closed orbit found out
3625    if ((k >= ntrial) && (errx >= tolx * 100))
3626        status.codflag = false;
3627
3628    free_dmatrix(alpha, 1, n, 1, n);
3629    free_dvector(bet, 1, n);
3630    free_dvector(fvect, 1, n);
3631    free_ivector(indx, 1, n);
3632}
3633
3634/****************************************************************************/
3635/* int Newton_Raphson(int n, double x[], int ntrial, double tolx)
3636 
3637 Purpose:
3638 Newton_Rapson algorithm from Numerical Recipes
3639 double precision algorithm
3640 Robustess: quadratic convergence
3641 Hint: for n-dimensional problem, the algo can be stuck on local minimum
3642 In this case, it should be enough to provide a resonable starting
3643 point.
3644
3645 Method:
3646 look for closed orbit solution of f(x) = x
3647 This problems is equivalent to finding the zero of g(x)= f(x) - x
3648 g(x+h) ~= f(x) - x + (Jacobian(f) -Id) h + O(h*h)
3649 Then at first order we solve h:
3650 h = - inverse(Jacobian(f) -Id) * (f(x)-x)
3651 the new guess is then xnew = x + h
3652 By iteration, this converges quadratically.
3653
3654 The algo is stopped whenever  |x -xnew| < tolx
3655
3656 f(x) is computes by tracking over one turn
3657 Jacobian(f) is computed numerically by numerical differentiation
3658 These two operations are provided by the function computeFandJ
3659
3660 Input:
3661 ntrial number of iterations for closed zero search
3662 x intial guess for the closed orbit
3663 tolx tolerance over the solution x
3664 tolf tolerance over the evalution f(x)
3665
3666 Output:
3667 x closed orbit
3668
3669 Return:
3670 none
3671
3672 Global variables:
3673 status
3674
3675 specific functions:
3676 computeFandJ
3677 InvMat, LinTrans
3678
3679 Comments:
3680 none
3681
3682 ****************************************************************************/
3683int Newton_Raphson(int n, Vector &x, int ntrial, double tolx) {
3684    int k=0, i=0;
3685    double errx=0.0;
3686    Vector bet, fvect;
3687    Matrix alpha;
3688
3689    errx = 0.0;
3690
3691    for (k = 1; k <= ntrial; k++) { // loop over number of iterations
3692        // supply function values at x in fvect and Jacobian matrix in fjac
3693        computeFandJ(n, x, alpha, fvect);
3694
3695        // Jacobian - Id
3696        for (i = 0; i < n; i++)
3697            alpha[i][i] -= 1.0;
3698        for (i = 0; i < n; i++)
3699            bet[i] = x[i] - fvect[i]; // right side of linear equation
3700        // inverse matrix using gauss jordan method from Tracy (from NR)
3701        if (!InvMat((long) n, alpha))
3702            fprintf(stdout, "Matrix non inversible ...\n");
3703        LinTrans((long) n, alpha, bet); // bet = alpha*bet
3704        errx = 0.0; // check root convergence
3705        for (i = 0; i < n; i++) { // update solution
3706            errx += fabs(bet[i]);
3707            x[i] += bet[i];
3708        }
3709
3710        if (trace)
3711            fprintf(
3712                    stdout,
3713                    "%02d: cod2 % .5e % .5e % .5e % .5e % .5e % .5e  errx =% .5e\n",
3714                    k, x[0], x[1], x[2], x[3], x[4], x[5], errx);
3715        if (errx <= tolx) {
3716            status.codflag = true;
3717            return 1;
3718        }
3719    }
3720    // check whever closed orbit found out
3721    if ((k >= ntrial) && (errx >= tolx)) {
3722        status.codflag = false;
3723        return 1;
3724    }
3725    return 0;
3726}
3727/*******************************************************************************
3728 *
3729 *
3730 *
3731 *
3732 ******************************************************************************/
3733void rm_mean(long int n, double x[]) {
3734    long int i=0L;
3735    double mean=0.0;
3736
3737    for (i = 0; i < n; i++)
3738        mean += x[i];
3739    mean /= n;
3740    for (i = 0; i < n; i++)
3741        x[i] -= mean;
3742}
Note: See TracBrowser for help on using the repository browser.