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

Last change on this file since 11 was 11, checked in by zhangj, 11 years ago
  • Property svn:executable set to *
File size: 18.1 KB
Line 
1/* Tracy-2
2
3   J. Bengtsson, CBP, LBL      1990 - 1994   Pascal version
4                 SLS, PSI      1995 - 1997
5   M. Boege      SLS, PSI      1998          C translation
6   L. Nadolski   SOLEIL        2002          Link to NAFF, Radia field maps
7   J. Bengtsson  NSLS-II, BNL  2004 -       
8
9*/
10/*
11 Current revision $Revision: 1.3 $
12 On branch $Name:  $
13 Latest change $Date: 2010/12/01 15:29:26 $ by $Author: nadolski $
14*/
15#include "complexeheader.h"
16
17double  pi = M_PI;
18
19/****************************************************************************/
20/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
21                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
22
23   Purpose:
24       Single particle tracking around the closed orbit for NTURN turns
25       The 6D phase trajectory is saved in a array, but the
26       tracked dp and ctau is not about the COD.
27
28   Input:
29       x, px, y, py 4 transverse coordinates
30       dp           energy offset
31       nmax         number of turns
32       pos          starting position for tracking
33       aperture     global physical aperture
34
35   Output:
36      lastn         last n (should be nmax if  not lost)
37      lastpos       last position in the ring
38      Tx            6xNTURN matrix of phase trajectory
39
40   Return:
41       none
42
43   Global variables:
44       NTURN number of turn for tracking
45       globval
46
47   Specific functions:
48       Cell_Pass
49
50   Comments:
51       useful for connection with NAFF
52       19/01/03 tracking around the closed orbit
53****************************************************************************/
54void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
55                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
56{
57  bool             lostF = false; /* Lost particle Flag */
58  ss_vect<double>  x1;  /* Tracking coordinates */
59  long             lastpos = globval.Cell_nLoc;
60  long             lastn = 0;
61  Vector2          aperture = {1.0, 1.0};
62
63  *status2 = true; /* stable */
64
65  if (globval.MatMeth) Cell_Concat(dp);
66
67  /* Get closed orbit */
68
69  getcod(dp, lastpos);
70
71 
72  if (!trace && status.codflag) 
73    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
74             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
75
76  /* Tracking coordinates around the closed orbit */
77  x1[0] =  x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
78  x1[2] =  y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
79  x1[4] = dp; x1[5] = ctau;
80
81  Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
82  Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
83  Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
84  lastn++;
85
86  do
87  { /* tracking through the ring */
88    if ((lastpos == globval.Cell_nLoc) &&
89        (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
90    {
91      if (globval.MatMeth)
92        Cell_fPass(x1, lastpos);
93      else
94        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
95      Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
96      Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
97      Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
98    }
99    else
100    {
101      printf("Trac_Simple: Particle lost \n");
102      fprintf(stdout, "%6ld plane: %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", 
103         lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
104      lostF = true;
105      *status2 = false;
106    }
107    lastn++;
108  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
109
110 
111  if (lastpos != globval.Cell_nLoc)
112  { /* Particle lost: Error message section */
113    *status2 = false;
114    printf("Trac_Simple: Particle lost \n");
115    fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
116             status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
117  }
118}
119
120
121/****************************************************************************/
122/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
123                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
124
125   Purpose:
126       Single particle tracking around the closed orbit for NTURN turns
127       The 6D phase trajectory is saved in a array, but the
128       tracked dp and ctau is not about the COD.
129
130   Input:
131       x, px, y, py 4 transverse coordinates
132       dp           energy offset
133       nmax         number of turns
134       pos          starting position for tracking
135       aperture     global physical aperture
136
137   Output:
138      lastn         last n (should be nmax if  not lost)
139      lastpos       last position in the ring
140      Tx            6xNTURN matrix of phase trajectory
141
142   Return:
143       none
144
145   Global variables:
146       NTURN number of turn for tracking
147       globval
148
149   Specific functions:
150       Cell_Pass
151
152   Comments:
153       useful for connection with NAFF
154       19/01/03 tracking around the closed orbit
155       
156       11/06/2013  Modified by Jianfeng Zhang @ LAL
157       The same feature as function 
158       Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
159                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
160       but add the feature to extract the position of the location of
161       the lost particle;
162       called by function
163       fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
164          double energy, bool diffusion, bool loss)       
165       in soleillib.cc.
166       
167****************************************************************************/
168void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
169                 double ctau, long nmax, double Tx[][NTURN], long& lastn, long& lastpos, ss_vect<double>& x1, bool *status2)
170{
171  bool             lostF = false; /* Lost particle Flag */
172  Vector2          aperture = {1.0, 1.0};
173
174   lastn = 0;
175   lastpos = globval.Cell_nLoc;
176  *status2 = true; /* stable */
177   x1[0]=0,x1[1]=0,x1[2]=0,x1[3]=0,x1[4]=0,x1[5]=0;
178 
179  if (globval.MatMeth) Cell_Concat(dp);
180
181  /* Get closed orbit */
182
183  getcod(dp, lastpos);
184
185 
186  if (!trace && status.codflag) 
187    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
188             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
189
190  /* Tracking coordinates around the closed orbit */
191  x1[0] =  x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
192  x1[2] =  y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
193  x1[4] = dp; x1[5] = ctau;
194
195  Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
196  Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
197  Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
198  lastn++;
199
200  do
201  { /* tracking through the ring */
202    if ((lastpos == globval.Cell_nLoc) &&
203        (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
204    {
205      if (globval.MatMeth)
206        Cell_fPass(x1, lastpos);
207      else
208        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
209      Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
210      Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
211      Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
212    }
213    else
214    {
215      printf("Trac_Simple: Particle lost \n");
216      fprintf(stdout, "%6ld plane: %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", 
217         lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
218      lostF = true;
219      *status2 = false;
220    }
221    lastn++;
222  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
223
224 
225  if (lastpos != globval.Cell_nLoc)
226  { /* Particle lost: Error message section */
227    *status2 = false;
228    printf("Trac_Simple: Particle lost \n");
229    fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
230             status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
231  }
232}
233
234
235/****************************************************************************/
236/* void Trac_Simple6DCOD(double x, double px, double y, double py, double dp, long nmax,
237                 double Tx[][NTURN], bool *status2)
238
239   Purpose:
240       Single particle tracking around the 6D closed orbit for NTURN turns
241       The 6D phase trajectory is saved in a array
242
243   Input:
244       x, px, y, py 4 transverses coordinates
245       dp           energy offset
246       nmax         number of turns
247       pos          starting position for tracking
248       aperture     global physical aperture
249
250   Output:
251      lastn         last n (should be nmax if  not lost)
252      lastpos       last position in the ring
253      Tx            6xNTURN matrix of phase trajectory
254
255   Return:
256       none
257
258   Global variables:
259       NTURN number of turn for tracking
260       globval
261
262   Specific functions:
263       Cell_Pass
264
265   Comments:
266       useful for connection with NAFF
267       19/01/03 tracking around the closed orbit
268****************************************************************************/
269void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
270                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
271{
272  bool             lostF = false; /* Lost particle Flag */
273  ss_vect<double>  x1;  /* Tracking coordinates */
274  long             lastpos = globval.Cell_nLoc;
275  long             lastn = 0;
276  Vector2          aperture = {1.0, 1.0};
277
278  *status2 = true; /* stable */
279
280  if (globval.MatMeth) Cell_Concat(dp);
281
282  /* Get closed orbit */
283
284  getcod(dp, lastpos);
285
286
287  if (trace && status.codflag)
288    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
289             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
290
291  /* Tracking coordinates around the closed orbit */
292  x1[0] =  x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
293  x1[2] =  y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
294  x1[4] = dp + globval.CODvect[4]; x1[5] = ctau + + globval.CODvect[5];
295
296  Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
297  Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
298  Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
299  lastn++;
300
301  do
302  { /* tracking through the ring */
303    if ((lastpos == globval.Cell_nLoc) &&
304        (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
305    {
306      if (globval.MatMeth)
307        Cell_fPass(x1, lastpos);
308      else
309        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
310      Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
311      Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
312      Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
313    }
314    else
315    {
316      printf("Trac_Simple: Particle lost \n");
317      fprintf(stdout, "%6ld plane: %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n",
318         lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
319      lostF = true;
320      *status2 = false;
321    }
322    lastn++;
323  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
324
325  if (lastpos != globval.Cell_nLoc)
326  { /* Particle lost: Error message section */
327    *status2 = false;
328    printf("Trac_Simple: Particle lost \n");
329    fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
330             status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
331  }
332}
333
334/****************************************************************************/
335/* void Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
336              double *fx, double *fz, int nb_freq[2])
337
338   Purpose:
339       Compute quasiperiodic approximation of a phase space trajectory
340       using NAFF Algorithm ((c) Laskar, IMCCE)
341
342   Input:
343       nterm number of frequencies to look for
344             if not multiple of 6, truncated to lower value
345       ndata size of the data to analyse
346       T     6D vector to analyse
347
348   Output:
349       fx frequencies found in the H-plane
350       fz frequencies found in the V-plane
351       nb_freq number of frequencies found out in each plane
352
353   Return:
354       none
355
356   Global variables:
357       g_NAFVariable  see modnaff.c
358       M_2_PI defined in math.h
359       trace ON or TRUE  for debugging
360
361   Specific functions:
362       naf_initnaf, naf_cleannaf
363       naf_mftnaf
364
365   Comments:
366       none
367
368****************************************************************************/
369/* Frequency Map Analysis */
370/* Analyse en Frequence */
371void Get_NAFF(int nterm, long ndata, double Tab[DIM][NTURN],
372              double *fx, double *fz, int nb_freq[2])
373{
374  /* Test whether ndata is divisible by 6 -- for NAFF -- */
375  /* Otherwise truncate ndata to lower value */
376  long r = 0; /* remainder of the euclidian division of ndata by 6 */
377  int i;
378
379  if ((r = ndata % 6) != 0) {
380    printf("Get_NAFF: Warning ndata = %ld, \n", ndata);
381    ndata -= r;
382    printf("New value for NAFF ndata = %ld \n", ndata);
383  }
384
385  g_NAFVariable.DTOUR      = M_2_PI;    /* size of one "cadran" */
386  g_NAFVariable.XH         = M_2_PI;    /* step */
387  g_NAFVariable.T0         = 0.0;       /* time t0 */
388  g_NAFVariable.NTERM      = nterm;     /* max term to find */
389  g_NAFVariable.KTABS      = ndata;     /* number of data: must be a multiple of 6 */
390  g_NAFVariable.m_pListFen = NULL;      /* no window */
391  g_NAFVariable.TFS        = NULL;      /* will contain frequency */
392  g_NAFVariable.ZAMP       = NULL;      /* will contain amplitude */
393  g_NAFVariable.ZTABS      = NULL;      /* will contain data to analyze */
394
395  /****************************************************/
396  /*               internal use in naf                */
397  g_NAFVariable.NERROR            = 0;
398  g_NAFVariable.ICPLX             = 1;
399  g_NAFVariable.IPRT              = -1;     /* 1 for diagnostics */
400  g_NAFVariable.NFPRT             = stdout; /* NULL   */
401  g_NAFVariable.NFS               = 0;
402  g_NAFVariable.IW                = 1;
403  g_NAFVariable.ISEC              = 1;
404  g_NAFVariable.EPSM              = 0;
405  g_NAFVariable.UNIANG            = 0;
406  g_NAFVariable.FREFON            = 0;
407  g_NAFVariable.ZALP              = NULL;
408  g_NAFVariable.m_iNbLineToIgnore = 1;      /* unused */
409  g_NAFVariable.m_dneps           = 1.e10;
410  g_NAFVariable.m_bFSTAB          = FALSE;  /* unused */
411  /*             end of interl use in naf             */
412  /****************************************************/
413
414  /* NAFF initialization */
415  naf_initnaf();
416
417  /**********************/
418  /* Analyse in H-plane */
419  /**********************/
420
421  /* fills up complexe vector for NAFF analysis */
422  for(i = 0; i < ndata; i++) {
423    g_NAFVariable.ZTABS[i].reel = Tab[0][i]; /* x  */
424    g_NAFVariable.ZTABS[i].imag = Tab[1][i]; /* xp */
425  }
426
427  /* Get out the mean value */
428  naf_smoy(g_NAFVariable.ZTABS);
429
430  naf_prtabs(g_NAFVariable.KTABS,g_NAFVariable.ZTABS, 20);
431//  trace=on;
432  naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
433
434  /* fill up H-frequency vector */
435  for (i = 1; i <= g_NAFVariable.NFS; i++)  {
436    fx[i-1] = g_NAFVariable.TFS[i];
437  }
438
439  nb_freq[0] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
440
441  if (trace)   /* print out results */
442  {
443    printf("(x,x') phase space: NFS=%d\n",g_NAFVariable.NFS);
444    for (i = 1; i <= g_NAFVariable.NFS; i++) {
445      printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
446              g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
447              i_compl_module(g_NAFVariable.ZAMP[i]),
448              i_compl_angle(g_NAFVariable.ZAMP[i]),
449              g_NAFVariable.TFS[i]);
450    }
451  }
452
453  /**********************/
454  /* Analyse in V-plane */
455  /**********************/
456
457  /* fill up complexe vector for NAFF analysis */
458  for (i = 0; i < ndata; i++) {
459    g_NAFVariable.ZTABS[i].reel = Tab[2][i];  /* z */
460    g_NAFVariable.ZTABS[i].imag = Tab[3][i];  /*zp */
461  }
462
463  naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
464
465  /* fills up V-frequency vector */
466  for (i = 1; i <= g_NAFVariable.NFS; i++) {
467    fz[i-1] =  g_NAFVariable.TFS[i];
468  }
469
470  nb_freq[1] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
471
472  if (trace)    /* print out results */
473  {
474    printf("(z,z') phase space: NFS=%d\n",g_NAFVariable.NFS);
475    for (i = 1; i <= g_NAFVariable.NFS; i++) {
476      printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
477              g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
478              i_compl_module(g_NAFVariable.ZAMP[i]),
479              i_compl_angle(g_NAFVariable.ZAMP[i]),
480              g_NAFVariable.TFS[i]);
481    }
482  }
483
484  /* free out memory used by NAFF */
485  naf_cleannaf();
486}
487
488/****************************************************************************/
489/* void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN],
490   long nbturn, long nshift)
491
492   Purpose:   used by fmap
493       Store in Tab0 values of Tab shifted by nshift turn
494
495   Input:
496       Tab    tracking tabular
497       nshift nb of turns to shift
498       nbturn nb of turns
499
500   Output:
501       Tab0  tracking tabular
502
503   Return:
504       none
505
506   Global variables:
507       none
508
509   Specific functions:
510       none
511
512   Comments:
513       
514
515****************************************************************************/
516void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN], long nbturn, long nshift)
517{
518  long k = 0L, k1 = 0L;
519
520  for (k = 0L; k < nbturn ; k++){
521    k1 = k + nshift;
522    Tab0[0][k] = Tab[0][k1];
523    Tab0[1][k] = Tab[1][k1];
524    Tab0[2][k] = Tab[2][k1];
525    Tab0[3][k] = Tab[3][k1];
526  }
527
528}
529
530/****************************************************************************/
531/* void Get_freq(double *fx, double *fz, double *nux, double *nuz)
532
533   Purpose:   used by fmap
534       Looks for tunes from frequency vectors output from NAFF
535
536   Input:
537       fx vector of horizontal frequencies
538       fz vector of verticcal frequencies
539
540   Output:
541       nux H tune
542       nuz V tune
543
544   Return:
545       none
546
547   Global variables:
548       none
549
550   Specific functions:
551       none
552
553   Comments:
554
555
556****************************************************************************/
557void Get_freq(double *fx, double *fz, double *nux, double *nuz)
558{
559  const double eps0 = 1e-4;
560  const double eps1 = 1e-6;
561 
562  // case of nux
563  if (fabs(fx[0]) < eps0){
564    *nux = fx[1];
565  }
566  else {
567    *nux = fx[0];
568  }
569
570  // case of nuz
571  if (fabs(fz[0]) < eps0) {
572     if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
573       if (fabs(fabs(fz[2]) - fabs(*nux)) < eps1) {
574          *nuz = fz[3];
575       }
576       else *nuz = fz[2];
577     }
578     else *nuz = fz[1];
579  }
580  else{
581    if (fabs(fabs(fz[0]) - fabs(*nux)) < eps1) {
582      if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
583         *nuz = fz[2];
584      }
585      else *nuz = fz[1];
586    }
587    else *nuz = fz[0];
588  }
589  *nuz = fabs(*nuz);  *nux = fabs(*nux);
590}
Note: See TracBrowser for help on using the repository browser.