source: TRACY3/branches/tracy3-3.10.1b/tracy/tracy/src/naffutils.cc @ 23

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

Clean version of Tracy: SoleilVersion at the end of 2011.Use this clean version to find the correct dipole fringe field to have the correct FMAP and FMAPDP of ThomX. Modified files: tpsa_lin.cc, soleillib.cc, prtmfile.cc, rdmfile.cc, read_script.cc, physlib.cc, tracy.cc, t2lat.cc, t2elem.cc, naffutils.cc in /tracy/tracy/src folder; naffutils.h, tracy_global.h, physlib.h, tracy.h, read_script.h, solielilib.h, t2elem.h in /tracy/tracy.inc folder; soltracy.cc in tracy/tools folder
  • Property svn:executable set to *
File size: 14.2 KB
Line 
1/* Tracy-2
2
3   J. Bengtsson, CBP, LBL      1990 - 1994   Pascal version
4                 SLS, PSI      1995 - 1997
5   M. Boege      SLS, PSI      1998          C translation
6   L. Nadolski   SOLEIL        2002          Link to NAFF, Radia field maps
7   J. Bengtsson  NSLS-II, BNL  2004 -       
8
9*/
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, long nmax,
21                 double Tx[][NTURN], bool *status)
22
23   Purpose:
24       Single particle tracking around the closed orbit for NTURN turns
25       The 6D phase trajectory is saved in a array
26
27   Input:
28       x, px, y, py 4 transverse coordinates
29       dp           energy offset
30       nmax         number of turns
31       pos          starting position for tracking
32       aperture     global physical aperture
33
34   Output:
35      lastn         last n (should be nmax if  not lost)
36      lastpos       last position in the ring
37      Tx            6xNTURN matrix of phase trajectory
38
39   Return:
40       none
41
42   Global variables:
43       NTURN number of turn for tracking
44       globval
45
46   Specific functions:
47       Cell_Pass
48
49   Comments:
50       useful for connection with NAFF
51       19/01/03 tracking around the closed orbit
52****************************************************************************/
53void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, 
54                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
55{
56  bool             lostF = false; /* Lost particle Flag */
57  ss_vect<double>  x1;  /* Tracking coordinates */
58  long             lastpos = globval.Cell_nLoc;
59  long             lastn = 0;
60  Vector2          aperture = {1.0, 1.0};
61
62  *status2 = true; /* stable */
63
64  if (globval.MatMeth) Cell_Concat(dp);
65
66  /* Get closed orbit */
67
68  getcod(dp, lastpos);
69
70 
71  if (trace && status.codflag) 
72    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
73             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
74
75  /* Tracking coordinates around the closed orbit */
76  x1[0] =  x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
77  x1[2] =  y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
78  x1[4] = dp; x1[5] = ctau;
79
80  Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
81  Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
82  Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
83  lastn++;
84
85  do
86  { /* tracking through the ring */
87    if ((lastpos == globval.Cell_nLoc) &&
88        (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
89    {
90      if (globval.MatMeth)
91        Cell_fPass(x1, lastpos);
92      else
93        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
94      Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
95      Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
96      Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
97    }
98    else
99    {
100      printf("Trac_Simple: Particle lost \n");
101      fprintf(stdout, "%6ld plane: %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", 
102         lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
103      lostF = true;
104      *status2 = false;
105    }
106    lastn++;
107  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
108
109  if (lastpos != globval.Cell_nLoc)
110  { /* Particle lost: Error message section */
111    *status2 = false;
112    printf("Trac_Simple: Particle lost \n");
113    fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
114             status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
115  }
116}
117
118/****************************************************************************/
119/* void Trac_Simple6DCOD(double x, double px, double y, double py, double dp, long nmax,
120                 double Tx[][NTURN], bool *status)
121
122   Purpose:
123       Single particle tracking around the 6D closed orbit for NTURN turns
124       The 6D phase trajectory is saved in a array
125
126   Input:
127       x, px, y, py 4 transverses coordinates
128       dp           energy offset
129       nmax         number of turns
130       pos          starting position for tracking
131       aperture     global physical aperture
132
133   Output:
134      lastn         last n (should be nmax if  not lost)
135      lastpos       last position in the ring
136      Tx            6xNTURN matrix of phase trajectory
137
138   Return:
139       none
140
141   Global variables:
142       NTURN number of turn for tracking
143       globval
144
145   Specific functions:
146       Cell_Pass
147
148   Comments:
149       useful for connection with NAFF
150       19/01/03 tracking around the closed orbit
151****************************************************************************/
152void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
153                 double ctau, long nmax, double Tx[][NTURN], bool *status2)
154{
155  bool             lostF = false; /* Lost particle Flag */
156  ss_vect<double>  x1;  /* Tracking coordinates */
157  long             lastpos = globval.Cell_nLoc;
158  long             lastn = 0;
159  Vector2          aperture = {1.0, 1.0};
160
161  *status2 = true; /* stable */
162
163  if (globval.MatMeth) Cell_Concat(dp);
164
165  /* Get closed orbit */
166
167  getcod(dp, lastpos);
168
169
170  if (trace && status.codflag)
171    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
172             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
173
174  /* Tracking coordinates around the closed orbit */
175  x1[0] =  x + globval.CODvect[0]; x1[1] = px + globval.CODvect[1];
176  x1[2] =  y + globval.CODvect[2]; x1[3] = py + globval.CODvect[3];
177  x1[4] = dp + globval.CODvect[4]; x1[5] = ctau + + globval.CODvect[5];
178
179  Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
180  Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
181  Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
182  lastn++;
183
184  do
185  { /* tracking through the ring */
186    if ((lastpos == globval.Cell_nLoc) &&
187        (fabs(x1[0]) < aperture[0]) && (fabs(x1[2]) < aperture[1]) && status.codflag)
188    {
189      if (globval.MatMeth)
190        Cell_fPass(x1, lastpos);
191      else
192        Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
193      Tx[0][lastn] = x1[0]; Tx[1][lastn] = x1[1];
194      Tx[2][lastn] = x1[2]; Tx[3][lastn] = x1[3];
195      Tx[4][lastn] = x1[4]; Tx[5][lastn] = x1[5];
196    }
197    else
198    {
199      printf("Trac_Simple: Particle lost \n");
200      fprintf(stdout, "%6ld plane: %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n",
201         lastn, status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
202      lostF = true;
203      *status2 = false;
204    }
205    lastn++;
206  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
207
208  if (lastpos != globval.Cell_nLoc)
209  { /* Particle lost: Error message section */
210    *status2 = false;
211    printf("Trac_Simple: Particle lost \n");
212    fprintf(stdout, "turn=%5ld plane= %1d %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g %+10.5g \n", lastn-1,
213             status.lossplane, x1[0], x1[1], x1[2], x1[3], x1[4], x1[5]);
214  }
215}
216
217/****************************************************************************/
218/* void Get_NAFF(int nterm, long ndata, double T[DIM][NTURN],
219              double *fx, double *fz, int nb_freq[2])
220
221   Purpose:
222       Compute quasiperiodic approximation of a phase space trajectory
223       using NAFF Algorithm ((c) Laskar, IMCCE)
224
225   Input:
226       nterm number of frequencies to look for
227             if not multiple of 6, truncated to lower value
228       ndata size of the data to analyse
229       T     6D vector to analyse
230
231   Output:
232       fx frequencies found in the H-plane
233       fz frequencies found in the V-plane
234       nb_freq number of frequencies found out in each plane
235
236   Return:
237       none
238
239   Global variables:
240       g_NAFVariable  see modnaff.c
241       M_2_PI defined in math.h
242       trace ON or TRUE  for debugging
243
244   Specific functions:
245       naf_initnaf, naf_cleannaf
246       naf_mftnaf
247
248   Comments:
249       none
250
251****************************************************************************/
252/* Frequency Map Analysis */
253/* Analyse en Frequence */
254void Get_NAFF(int nterm, long ndata, double Tab[DIM][NTURN],
255              double *fx, double *fz, int nb_freq[2])
256{
257  /* Test whether ndata is divisible by 6 -- for NAFF -- */
258  /* Otherwise truncate ndata to lower value */
259  long r = 0; /* remainder of the euclidian division of ndata by 6 */
260  int i;
261
262  if ((r = ndata % 6) != 0) {
263    printf("Get_NAFF: Warning ndata = %ld, \n", ndata);
264    ndata -= r;
265    printf("New value for NAFF ndata = %ld \n", ndata);
266  }
267
268  g_NAFVariable.DTOUR      = M_2_PI;    /* size of one "cadran" */
269  g_NAFVariable.XH         = M_2_PI;    /* step */
270  g_NAFVariable.T0         = 0.0;       /* time t0 */
271  g_NAFVariable.NTERM      = nterm;     /* max term to find */
272  g_NAFVariable.KTABS      = ndata;     /* number of data: must be a multiple of 6 */
273  g_NAFVariable.m_pListFen = NULL;      /* no window */
274  g_NAFVariable.TFS        = NULL;      /* will contain frequency */
275  g_NAFVariable.ZAMP       = NULL;      /* will contain amplitude */
276  g_NAFVariable.ZTABS      = NULL;      /* will contain data to analyze */
277
278  /****************************************************/
279  /*               internal use in naf                */
280  g_NAFVariable.NERROR            = 0;
281  g_NAFVariable.ICPLX             = 1;
282  g_NAFVariable.IPRT              = -1;     /* 1 for diagnostics */
283  g_NAFVariable.NFPRT             = stdout; /* NULL   */
284  g_NAFVariable.NFS               = 0;
285  g_NAFVariable.IW                = 1;
286  g_NAFVariable.ISEC              = 1;
287  g_NAFVariable.EPSM              = 0;
288  g_NAFVariable.UNIANG            = 0;
289  g_NAFVariable.FREFON            = 0;
290  g_NAFVariable.ZALP              = NULL;
291  g_NAFVariable.m_iNbLineToIgnore = 1;      /* unused */
292  g_NAFVariable.m_dneps           = 1.e10;
293  g_NAFVariable.m_bFSTAB          = FALSE;  /* unused */
294  /*             end of interl use in naf             */
295  /****************************************************/
296
297  /* NAFF initialization */
298  naf_initnaf();
299
300  /**********************/
301  /* Analyse in H-plane */
302  /**********************/
303
304  /* fills up complexe vector for NAFF analysis */
305  for(i = 0; i < ndata; i++) {
306    g_NAFVariable.ZTABS[i].reel = Tab[0][i]; /* x  */
307    g_NAFVariable.ZTABS[i].imag = Tab[1][i]; /* xp */
308  }
309
310  /* Get out the mean value */
311  naf_smoy(g_NAFVariable.ZTABS);
312
313  naf_prtabs(g_NAFVariable.KTABS,g_NAFVariable.ZTABS, 20);
314//  trace=on;
315  naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
316
317  /* fill up H-frequency vector */
318  for (i = 1; i <= g_NAFVariable.NFS; i++)  {
319    fx[i-1] = g_NAFVariable.TFS[i];
320  }
321
322  nb_freq[0] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
323
324  if (trace)   /* print out results */
325  {
326    printf("(x,x') phase space: NFS=%d\n",g_NAFVariable.NFS);
327    for (i = 1; i <= g_NAFVariable.NFS; i++) {
328      printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
329              g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
330              i_compl_module(g_NAFVariable.ZAMP[i]),
331              i_compl_angle(g_NAFVariable.ZAMP[i]),
332              g_NAFVariable.TFS[i]);
333    }
334  }
335
336  /**********************/
337  /* Analyse in V-plane */
338  /**********************/
339
340  /* fill up complexe vector for NAFF analysis */
341  for (i = 0; i < ndata; i++) {
342    g_NAFVariable.ZTABS[i].reel = Tab[2][i];  /* z */
343    g_NAFVariable.ZTABS[i].imag = Tab[3][i];  /*zp */
344  }
345
346  naf_mftnaf(nterm,fabs(g_NAFVariable.FREFON)/g_NAFVariable.m_dneps);
347
348  /* fills up V-frequency vector */
349  for (i = 1; i <= g_NAFVariable.NFS; i++) {
350    fz[i-1] =  g_NAFVariable.TFS[i];
351  }
352
353  nb_freq[1] = g_NAFVariable.NFS; /* nb of frequencies found out by NAFF */
354
355  if (trace)    /* print out results */
356  {
357    printf("(z,z') phase space: NFS=%d\n",g_NAFVariable.NFS);
358    for (i = 1; i <= g_NAFVariable.NFS; i++) {
359      printf("AMPL=%15.8E+i*%15.8E abs(AMPL)=%15.8E arg(AMPL)=%15.8E FREQ=%15.8E\n",
360              g_NAFVariable.ZAMP[i].reel,g_NAFVariable.ZAMP[i].imag,
361              i_compl_module(g_NAFVariable.ZAMP[i]),
362              i_compl_angle(g_NAFVariable.ZAMP[i]),
363              g_NAFVariable.TFS[i]);
364    }
365  }
366
367  /* free out memory used by NAFF */
368  naf_cleannaf();
369}
370
371/****************************************************************************/
372/* void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN],
373   long nbturn, long nshift)
374
375   Purpose:   used by fmap
376       Store in Tab0 values of Tab shifted by nshift turn
377
378   Input:
379       Tab    tracking tabular
380       nshift nb of turns to shift
381       nbturn nb of turns
382
383   Output:
384       Tab0  tracking tabular
385
386   Return:
387       none
388
389   Global variables:
390       none
391
392   Specific functions:
393       none
394
395   Comments:
396       
397
398****************************************************************************/
399void Get_Tabshift(double Tab[DIM][NTURN], double Tab0[DIM][NTURN], long nbturn, long nshift)
400{
401  long k = 0L, k1 = 0L;
402
403  for (k = 0L; k < nbturn ; k++){
404    k1 = k + nshift;
405    Tab0[0][k] = Tab[0][k1];
406    Tab0[1][k] = Tab[1][k1];
407    Tab0[2][k] = Tab[2][k1];
408    Tab0[3][k] = Tab[3][k1];
409  }
410
411}
412
413/****************************************************************************/
414/* void Get_freq(double *fx, double *fz, double *nux, double *nuz)
415
416   Purpose:   used by fmap
417       Looks for tunes from frequency vectors output from NAFF
418
419   Input:
420       fx vector of horizontal frequencies
421       fz vector of verticcal frequencies
422
423   Output:
424       nux H tune
425       nuz V tune
426
427   Return:
428       none
429
430   Global variables:
431       none
432
433   Specific functions:
434       none
435
436   Comments:
437
438
439****************************************************************************/
440void Get_freq(double *fx, double *fz, double *nux, double *nuz)
441{
442  const double eps0 = 1e-4;
443  const double eps1 = 1e-6;
444 
445  // case of nux
446  if (fabs(fx[0]) < eps0){
447    *nux = fx[1];
448  }
449  else {
450    *nux = fx[0];
451  }
452
453  // case of nuz
454  if (fabs(fz[0]) < eps0) {
455     if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
456       if (fabs(fabs(fz[2]) - fabs(*nux)) < eps1) {
457          *nuz = fz[3];
458       }
459       else *nuz = fz[2];
460     }
461     else *nuz = fz[1];
462  }
463  else{
464    if (fabs(fabs(fz[0]) - fabs(*nux)) < eps1) {
465      if (fabs(fabs(fz[1]) - fabs(*nux)) < eps1) {
466         *nuz = fz[2];
467      }
468      else *nuz = fz[1];
469    }
470    else *nuz = fz[0];
471  }
472  *nuz = fabs(*nuz);  *nux = fabs(*nux);
473}
Note: See TracBrowser for help on using the repository browser.