Changeset 23 in TRACY3


Ignore:
Timestamp:
Dec 6, 2013, 5:12:43 PM (10 years ago)
Author:
zhangj
Message:

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

Location:
branches/tracy3-3.10.1b/tracy
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • branches/tracy3-3.10.1b/tracy/tools/soltracy.cc

    r11 r23  
    1 /************************************
     1/*
    22 Current revision $Revision$
    33 On branch $Name$
    44 Latest change $Date$ by $Author$
    5 *************************************/
    6 #define ORDER 1   
    7 //#define ORDER 4   
     5*/
     6#define ORDER 1         
    87
    98int no_tps = ORDER; // arbitrary TPSA order is defined locally
    109
     10
     11
    1112extern bool freq_map;
    1213
    13 // #if HAVE_CONFIG_H
    14 // #include <config.h>
    15 // #endif
    16 
    17 // #if MPI_EXEC
    18 //   //library of parallel computation,must be before "stdio.h"
    19 //   #include <mpi.h>
    20 // #endif
     14#if HAVE_CONFIG_H
     15#include <config.h>
     16#endif
     17
     18/*
     19#if MPI_EXEC
     20  //library of parallel computation,must be before "stdio.h"
     21  #include <mpi.h>
     22#endif
     23*/
    2124
    2225  #include "tracy_lib.h"
    … …  
    4144  globval.H_exact = false;
    4245
    43  
    44   //output files
    45   char fic_twiss[S_SIZE + 4] = "";   //twiss file
    46   char fic_summary[S_SIZE + 4]="";   //summary file
    47  
    48  
    4946  /* parameters to read the user input script .prm */
    5047  long i=0L; //initialize the for loop to read command string
    … …  
    6461  UserCommand UserCommandFlag[NCOMMAND];
    6562 
    66 // #if MPI_EXEC
    67 //   //Initialize parallel computation
    68 //     MPI_Init(&argc, &argv);
    69 // #endif
     63  /*
     64#if MPI_EXEC
     65  //Initialize parallel computation
     66    MPI_Init(&argc, &argv);
     67#endif
     68  */
     69
    7070  /************************************************************************
    7171   read in files and flags
    … …  
    7676    read_script(argv[1], true,CommandNo, UserCommandFlag);
    7777  } else {
    78     fprintf(stdout, "Not enough parameters\n Syntax is program parameter file\n");
     78    fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
    7979    exit_(1);
    8080  }
    … …  
    105105  else
    106106    globval.ge = ElemIndex(ge_name);
    107  
    108  
     107 
    109108//  globval.g = ElemIndex("g");  /* get family index of  girder*/
    110  strcpy(fic_twiss,lat_file);
    111   strcpy(fic_summary,lat_file);
    112   strcat(fic_twiss,".twi");
    113   strcat(fic_summary,".sum");
    114  
    115     //generate the twiss and summary files with proper appendix 
    116   /* print the summary of the element in lattice to the screen and an external file*/ 
     109   
     110  /* print the summary of the element in lattice */ 
    117111  printglob();   
    118   printglob2file(fic_summary);
    119  //print the twiss file
    120    /* print out lattice functions, with all the optical information for the lattice with design values */
    121   printlatt(fic_twiss); 
    122  
    123112  /************************************************************************
    124113    print files, very important file for debug
    … …  
    263252   }
    264253   
    265    else if(strcmp(CommandStr,"PrintTrackElemFlag") == 0) {
    266       cout << "\n";
    267       cout << "print the tracked coordinates to file: "
    268            << UserCommandFlag[i]._PrintTrackElem_track_file << "\n";
    269        
    270       PrintTrackElem(UserCommandFlag[i]._PrintTrackElem_track_file,
    271                  UserCommandFlag[i]._PrintTrackElem_x, UserCommandFlag[i]._PrintTrackElem_px,
    272                  UserCommandFlag[i]._PrintTrackElem_y, UserCommandFlag[i]._PrintTrackElem_py,
    273                  UserCommandFlag[i]._PrintTrackElem_delta,UserCommandFlag[i]._PrintTrackElem_ctau,
    274                  UserCommandFlag[i]._PrintTrackElem_nelem1,UserCommandFlag[i]._PrintTrackElem_nelem2); 
    275    }
    276254    //print the girder
    277255  // else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
    … …  
    468446  else if(strcmp(CommandStr,"FmapFlag") == 0) {
    469447    printf("\n begin Fmap calculation for on momentum particles: \n");
    470    
    471   //  #if MPI_EXEC
     448 
     449  //   /* 
     450 //   #if MPI_EXEC
    472451
    473452 //    //initialization for parallel computing
    … …  
    490469 //        UserCommandFlag[i]._FmapFlag_delta,
    491470 //        UserCommandFlag[i]._FmapFlag_diffusion,
    492  //        UserCommandFlag[i]._FmapFlag_printloss,
    493471 //        numprocs,myid);
    494472
    … …  
    543521           UserCommandFlag[i]._FmapFlag_ymax,
    544522           UserCommandFlag[i]._FmapFlag_delta,
    545            UserCommandFlag[i]._FmapFlag_diffusion,
    546            UserCommandFlag[i]._FmapFlag_printloss);
     523           UserCommandFlag[i]._FmapFlag_diffusion);
    547524      //#endif
     525   
    548526  }
    549527
    … …  
    551529  else if(strcmp(CommandStr,"FmapdpFlag") == 0) {
    552530    printf("\n begin Fmap calculation for off momentum particles: \n");
    553    
     531   
     532  //   /*
    554533 //   #if MPI_EXEC
    555534
    … …  
    572551 //          UserCommandFlag[i]._FmapdpFlag_z,
    573552 //          UserCommandFlag[i]._FmapdpFlag_diffusion,
    574  //            UserCommandFlag[i]._FmapdpFlag_printloss,
    575553 //          numprocs,myid);
    576554
    … …  
    625603              UserCommandFlag[i]._FmapdpFlag_emax,
    626604              UserCommandFlag[i]._FmapdpFlag_z,
    627               UserCommandFlag[i]._FmapdpFlag_diffusion,
    628                UserCommandFlag[i]._FmapdpFlag_printloss);
    629        //  #endif
     605              UserCommandFlag[i]._FmapdpFlag_diffusion);
     606       // #endif
     607   
    630608  }
    631609
    … …  
    657635    printf("\n Calculate momentum acceptance: \n");
    658636
    659 //  #if MPI_EXEC
     637   
     638 // #if MPI_EXEC
    660639   
    661640//     /* calculate momentum acceptance*/
    … …  
    789768
    790769// #else
     770
    791771     MomentumAcceptance(UserCommandFlag[i]._MomentumAccFlag_momacc_file,
    792772                        UserCommandFlag[i]._MomentumAccFlag_istart,
    … …  
    800780                        UserCommandFlag[i]._MomentumAccFlag_nturn,
    801781                        UserCommandFlag[i]._MomentumAccFlag_zmax);
    802      //  #endif   
     782     /*
     783  #endif   
     784     */
    803785
    804786    /* restore the initial values*/
    … …  
    863845          UserCommandFlag[i]._Phase_ctau,
    864846          UserCommandFlag[i]._Phase_nturn);
    865     printf("6D phase space tracking: \n the simulation time for phase space in tracy 3 is \n");
     847    printf("the simulation time for phase space in tracy 3 is \n");
    866848    stop = stampstop(start);
    867849
    … …  
    10211003
    10221004 
    1023  // #if MPI_EXEC
    1024  //  MPI_Finalize();   //finish parallel computation
    1025  // #endif
    1026  
     1005  /*
     1006 #if MPI_EXEC
     1007  MPI_Finalize();   //finish parallel computation
     1008 #endif
     1009  */
     1010
    10271011  return 0;
    10281012}//end of main()

  • branches/tracy3-3.10.1b/tracy/tracy/inc/naffutils.h

    r11 r23  
    2727                 double ctau, long nmax, double Tx[][NTURN], bool *status);
    2828
    29 void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
    30                  double ctau, long nmax, double Tx[][NTURN], long& lastn, long& lastpos,ss_vect<double>& x1,bool *status2);
    31 
    3229void Trac_Simple6DCOD(double x, double px, double y, double py, double dp,
    3330         double ctau, long nmax, double Tx[][NTURN], bool *status);

  • branches/tracy3-3.10.1b/tracy/tracy/inc/physlib.h

    r11 r23  
    5555
    5656void printglob(void);
    57 
    58 void printglob2file(const char fic[]);
    5957
    6058void printlatt(const char fic[]);

  • branches/tracy3-3.10.1b/tracy/tracy/inc/read_script.h

    r11 r23  
    3737   long  _PrintTrack_nmax;
    3838   
    39    //start track particle coordinates; PrintTrackFlag
    40    char _PrintTrackElem_track_file[max_str];
    41    double _PrintTrackElem_x, _PrintTrackElem_px, _PrintTrackElem_y,
    42           _PrintTrackElem_py, _PrintTrackElem_delta, _PrintTrackElem_ctau;
    43    long  _PrintTrackElem_nelem1, _PrintTrackElem_nelem2;
    44    
    45    
    4639   //twiss file
    4740   char _PrintTwiss_twiss_file[max_str];
    … …  
    6558   long _FmapFlag_nxpoint, _FmapFlag_nypoint, _FmapFlag_nturn;
    6659   double _FmapFlag_xmax, _FmapFlag_ymax, _FmapFlag_delta;
    67    bool _FmapFlag_diffusion;
    68    bool _FmapFlag_printloss;
    69    
     60   bool _FmapFlag_diffusion;
    7061 
    7162 //extern bool FmapdpFlag;
    … …  
    7364   long _FmapdpFlag_nxpoint, _FmapdpFlag_nepoint, _FmapdpFlag_nturn;
    7465   double _FmapdpFlag_xmax, _FmapdpFlag_emax, _FmapdpFlag_z;
    75    bool _FmapdpFlag_diffusion;
    76    bool _FmapdpFlag_printloss;
     66   bool _FmapdpFlag_diffusion; 
    7767
    7868   
    … …  
    120110   _FmapFlag_nxpoint=31L, _FmapFlag_nypoint=21L, _FmapFlag_nturn=516L;
    121111   _FmapFlag_xmax=0.025, _FmapFlag_ymax=0.005, _FmapFlag_delta=0.0;
    122    _FmapFlag_diffusion = true, _FmapFlag_printloss = false;
     112   _FmapFlag_diffusion = true;
    123113
    124114/*fmap for off momentum particle*/
    … …  
    126116 _FmapdpFlag_nxpoint=31L, _FmapdpFlag_nepoint=21L, _FmapdpFlag_nturn=516L;
    127117 _FmapdpFlag_xmax=0.025, _FmapdpFlag_emax=0.005, _FmapdpFlag_z=0.0;
    128  _FmapdpFlag_diffusion = true, _FmapdpFlag_printloss = false;
    129  
     118 _FmapdpFlag_diffusion = true;                 
     119
    130120/* tune shift with amplitude*/
    131121 strcpy(_AmplitudeTuneShift_nudx_file,"nudx.out");
    … …  
    172162
    173163/***** file names************/
    174 extern char    lat_file[max_str];
    175 
    176 //files with multipole errors; for soleil lattice
     164 //files with multipole errors; for soleil lattice
    177165 extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
    178166 extern char multipole_file[max_str];

  • branches/tracy3-3.10.1b/tracy/tracy/inc/soleillib.h

    r11 r23  
    5858//              double z, bool diffusion, bool matlab);
    5959void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    60           double energy, bool diffusion);       
    61 void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    62           double energy, bool diffusion, bool printloss);       
     60          double energy, bool diffusion);         
    6361void fmap_p(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    64           double energy, bool diffusion, bool printloss, int numprocs, int myid);         
     62          double energy, bool diffusion, int numprocs, int myid);         
    6563void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
    66               double z, bool diffusion, bool printloss);
     64              double z, bool diffusion);
    6765void fmapdp_p(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
    68             double z, bool diffusion, bool printloss, int numprocs, int myid);
     66            double z, bool diffusion, int numprocs, int myid);
    6967void Nu_Naff(void);
    7068
    … …  
    130128/* fit tunes for soleil lattice, in which each quadrupole is cut into two parts*/                       
    131129void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy);                   
    132 //print the coordinates at lattice elements tracked for n turns
     130//print the coordinates at lattice elements
    133131void PrintTrack(const char *TrackFile, double x, double px,double y,double py,
    134132                double delta, double ctau, long int nmax);
    135 //print the tracked coordinates after a lattice element
    136 void PrintTrackElem(const char *TrackFile, double x, double px,double y,double py,
    137                 double delta, double ctau, long int nelem1,long int nelem2);
     133               
    138134
    139135

  • branches/tracy3-3.10.1b/tracy/tracy/inc/t2elem.h

    r11 r23  
    3838
    3939template<typename T>
    40 void Drift(double L,double h_bend, ss_vect<T> &x);
    41 
    42 template<typename T>
    4340void Drift(double L, ss_vect<T> &x);
    4441
    … …  
    4744
    4845template<typename T>
    49 void bend_fringe(double hb, double gap, ss_vect<T> &x);
    50 
    51 template<typename T>
    52 static void bendCurvature(double irho, double H, ss_vect<T> &x);
    53 
    54 template<typename T>
    55 static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x, bool Entrance);
     46static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x);
    5647
    5748template<typename T>
    … …  
    6152template<typename T>
    6253void Drift_Pass(CellType &Cell, ss_vect<T> &x);
    63 
    64 template<typename T>
    65 void dipole_kick(double L, double h_ref, double h_bend, ss_vect<T> &x);
    66 template<typename T>
    67 void multipole_kick(int Order, double MB[], double L, double h_bend, ss_vect<T> &x);
    6854
    6955template<typename T>

  • branches/tracy3-3.10.1b/tracy/tracy/inc/tpsa_lin.h

    r3 r23  
    7575void dacos(const tps_buf &x, tps_buf &z);
    7676
    77 void datan(const tps_buf &x, tps_buf &z);
    78 
    79 void dactan(const tps_buf &x, tps_buf &z);
    80 
    8177void dasinh(const tps_buf &x, tps_buf &z);
    8278
    8379void dacosh(const tps_buf &x, tps_buf &z);
    8480
    85 void daarcsin(const tps_buf &x, tps_buf &z);
    86 
    87 void daarccos(const tps_buf &x, tps_buf &z);
     81void datan(const tps_buf &x, tps_buf &z);
    8882
    8983void daarctan(const tps_buf &x, tps_buf &z);
    90 
    91 void daarcctan(const tps_buf &x, tps_buf &z);
    9284
    9385void dafun_(const char *fun, const tps_buf &x, tps_buf &z);

  • branches/tracy3-3.10.1b/tracy/tracy/inc/tracy.h

    r11 r23  
    110110  bool tuneflag, chromflag, codflag, mapflag, passflag, overflag, chambre;
    111111  int lossplane; /* lost in: horizontal    1
    112                              vertical      2
    113                              longitudinal  3
    114                               not lost in any of the above plane: 0*/
     112                                     vertical      2
     113                                         longitudinal  3 */
    115114} statusrec;
    116115

  • branches/tracy3-3.10.1b/tracy/tracy/inc/tracy_global.h

    r11 r23  
    7979  int         Pmethod;   // Integration Method
    8080  int         PN;        // Number of integration steps
    81   long        dipEdge_effect1; // dipole edge effect at the entrance
    82   long        dipEdge_effect2; // dipole edge effect at the exit
    8381  long        quadFF1;         /* Entrance quadrupole Fringe field flag */
    8482  long        quadFF2;         /* Exit quadrupole Fringe field flag */
    … …  
    10098  double      PdTrnd;    // (normal)random scale factor to rotation error PdTrms
    10199  // Multipole strengths
    102   mpolArray   PBpar;     // design field gradient; bn, and an
    103   mpolArray   PBsys;     // systematic multipole errors gradient, bn and an
    104   mpolArray   PBrms;     // rms multipole field errors gradient, bn and an
    105   mpolArray   PBrnd;     // random scale factor of rms error PBrms gradient, bn and an
    106   mpolArray   PB;        // total field strength(design,sys,rms) gradient, bn and an
     100  mpolArray   PBpar;     // design field strength
     101  mpolArray   PBsys;     // systematic multipole errors
     102  mpolArray   PBrms;     // rms multipole field errors
     103  mpolArray   PBrnd;     // random scale factor of rms error PBrms
     104  mpolArray   PB;        // total field strength(design,sys,rms)
    107105  int         Porder;    // The highest order in PB
    108106  int         n_design;  // multipole order (design, All = 0, Dip = 1, Quad = 2, Sext = 3, Oct = 4, Dec = 5, Dodec = 6)
    … …  
    111109  double PTx1;           // horizontal entrance angle [deg]
    112110  double PTx2;           // horizontal exit angle [deg]
    113   double PH1;            // bending curvature of the entrance pole face of dipole, see P116 SAC-75.
    114   double PH2;            // bending curvature of the exit pole face of dipole, see P116 SAC-75.
    115111  double Pgap;           // total magnet gap [m]
    116112  double Pirho;          // angle([radian], but in lattice definition, angle is with unit degree)/length of the dipole, 1/rho [1/m]

  • branches/tracy3-3.10.1b/tracy/tracy/src/naffutils.cc

    r11 r23  
    1818
    1919/****************************************************************************/
    20 /* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp,
    21                  double ctau, long nmax, double Tx[][NTURN], bool *status2)
     20/* void Trac_Simple4DCOD(double x, double px, double y, double py, double dp, long nmax,
     21                 double Tx[][NTURN], bool *status)
    2222
    2323   Purpose:
    2424       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.
     25       The 6D phase trajectory is saved in a array
    2726
    2827   Input:
    … …  
    7069
    7170 
    72   if (!trace && status.codflag)
     71  if (trace && status.codflag)
    7372    printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
    7473             dp*1e2, globval.CODvect[0]*1e3, globval.CODvect[2]*1e3);
    … …  
    108107  } while ((lastn < nmax) && (lastpos == globval.Cell_nLoc) && (lostF == false));
    109108
    110  
    111109  if (lastpos != globval.Cell_nLoc)
    112110  { /* Particle lost: Error message section */
    … …  
    118116}
    119117
    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 ****************************************************************************/
    168 void 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 
    235118/****************************************************************************/
    236119/* void Trac_Simple6DCOD(double x, double px, double y, double py, double dp, long nmax,
    237                  double Tx[][NTURN], bool *status2)
     120                 double Tx[][NTURN], bool *status)
    238121
    239122   Purpose:

  • branches/tracy3-3.10.1b/tracy/tracy/src/physlib.cc

    r11 r23  
    1616/**************************/
    1717
    18 /*******************************************************************************
    19  *
    20  *
    21  *
    22  *
    23  ******************************************************************************/
     18
    2419/**** same as asctime in C without the \n at the end****/
    2520char *asctime2(const struct tm *timeptr) {
    … …  
    3833    return result;
    3934}
    40 /*******************************************************************************
    41  *
    42  *
    43  *
    44  *
    45  ******************************************************************************/
     35
    4636/** Get time and date **/
    4737struct tm* GetTime() {
    … …  
    199189
    200190 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
    203191
    204192 ****************************************************************************/
    … …  
    263251
    264252/****************************************************************************/
    265 /* void printglob2file(void)
    266 
    267  Purpose:
    268  Print global variables on screen
    269  Print tunes and chromaticities
    270  Print Oneturn matrix
    271 
    272  Input:
    273  none
    274 
    275  Output:
    276  output on the screen
    277 
    278  Return:
    279  none
    280 
    281  Global variables:
    282  globval
    283 
    284  Specific functions:
    285  none
    286 
    287  Comments:
    288  03/06/2013  by Jianfeng Zhang @ LAL
    289  The same feature as the file printglob(void), but
    290  added feature to print the lattice summary to an external file
    291 
    292  
    293  ****************************************************************************/
    294 void printglob2file(const char fic[]) {
    295   FILE *fout;
    296   int i=0,j=0;
    297  
    298   fout = fopen(fic,"w");
    299  
    300     fprintf(fout, "\n***************************************************************"
    301         "***************\n");
    302     fprintf(fout,"\n");
    303     fprintf(fout,"  dPcommon     =  %9.3e  dPparticle   =  %9.3e"
    304         "  Energy [GeV] = %.3f\n", globval.dPcommon, globval.dPparticle,
    305             globval.Energy);
    306     fprintf(fout,"  MaxAmplx [m] = %9.3e   MaxAmply [m] = %9.3e"
    307         "   RFAccept [%%] = +/- %4.2f\n", Cell[0].maxampl[X_][1],
    308             Cell[0].maxampl[Y_][1], globval.delta_RF * 1e2);
    309     fprintf(fout,"  MatMeth      =  %s     ", globval.MatMeth ? "TRUE " : "FALSE");
    310     fprintf(fout," Cavity_On    =  %s    ", globval.Cavity_on ? "TRUE " : "FALSE");
    311     fprintf(fout,"  Radiation_On = %s     \n", globval.radiation ? "TRUE " : "FALSE");
    312    
    313     if(globval.bpm == 0)
    314       fprintf(fout,"  bpm          =  0"); 
    315     else
    316       fprintf(fout,"  bpm          =  %3d", GetnKid(globval.bpm));
    317     if(globval.qt == 0)
    318       fprintf(fout,"                             qt           = 0           \n");
    319     else
    320       fprintf(fout,"                             qt           = %3d         \n", GetnKid(globval.qt));   
    321     if(globval.hcorr == 0)
    322       fprintf(fout,"  hcorr         =  0"); 
    323     else
    324       fprintf(fout,"  hcorr         = %3d", GetnKid(globval.hcorr));
    325     if(globval.vcorr == 0)
    326       fprintf(fout,"                             vcorr        =  0        \n"); 
    327     else
    328       fprintf(fout,"                             vcorr        = %3d      \n", GetnKid(globval.vcorr));
    329                                            
    330     fprintf(fout," Chambre_On   = %s     \n", globval.Aperture_on ? "TRUE " : "FALSE");
    331  
    332     fprintf(fout,"  alphac       =   %8.4e\n", globval.Alphac);
    333     fprintf(fout,"  nux          =   %8.6f      nuy  =  %8.6f",
    334             globval.TotalTune[X_], globval.TotalTune[Y_]);
    335     if (globval.Cavity_on)
    336         fprintf(fout,"  omega  = %13.9f\n", globval.Omega);
    337     else {
    338         fprintf(fout,"\n");
    339         fprintf(fout,"  ksix         =    %8.6f      ksiy  =   %8.6f\n",
    340                 globval.Chrom[X_], globval.Chrom[Y_]);
    341     }
    342     fprintf(fout,"\n");
    343    
    344     fprintf(fout,"  OneTurn matrix:\n");
    345     for(i=0;i<ss_dim;i++){
    346       fprintf(fout,"\n");
    347       for(j=0;j<ss_dim;j++)
    348         fprintf(fout,"%14.6e",globval.OneTurnMat[i][j]);
    349     }
    350 
    351     fprintf(fout,"\n\nTwiss parameters at entrance:\n");
    352     fprintf(fout,
    353             "   Betax [m] = % 9.3e  Alphax = % 9.3e Etax [m] = % 9.3e Etaxp = % 9.3e\n"
    354                 "   Betay [m] = % 9.3e  Alphay = % 9.3e Etay [m] = % 9.3e Etayp = % 9.3e\n\n",
    355             Cell[1].Beta[X_], Cell[1].Alpha[X_], Cell[1].Eta[X_],
    356             Cell[1].Etap[X_], Cell[1].Beta[Y_], Cell[1].Alpha[Y_],
    357             Cell[1].Eta[Y_], Cell[1].Etap[Y_]);
    358 
    359     fclose(fout);
    360 }
    361 
    362 /****************************************************************************/
    363 /* void printlatt(const char fic[])
     253/* void printlatt(void)
    364254
    365255 Purpose:
    … …  
    423313    fclose(outf);
    424314}
    425 /*******************************************************************************
    426  *
    427  *
    428  *
    429  *
    430  ******************************************************************************/
     315
    431316double int_curly_H1(long int n) {
    432317    /* Integration with Simpson's Rule */
    … …  
    445330    return curly_H;
    446331}
    447 /*******************************************************************************
    448  *
    449  *
    450  *
    451  *
    452  ******************************************************************************/
     332
    453333void prt_sigma(void) {
    454334    long int i;
    … …  
    597477    }
    598478}
    599 /*******************************************************************************
    600  *
    601  *
    602  *
    603  *
    604  ******************************************************************************/
     479
    605480void getabn(Vector2 &alpha, Vector2 &beta, Vector2 &nu) {
    606481    Vector2 gamma;
    607482    Cell_GetABGN(globval.OneTurnMat, alpha, beta, gamma, nu);
    608483}
    609 /*******************************************************************************
    610  *
    611  *
    612  *
    613  *
    614  ******************************************************************************/
     484
    615485void TraceABN(long i0, long i1, const Vector2 &alpha, const Vector2 &beta,
    616486        const Vector2 &eta, const Vector2 &etap, const double dP) {
    … …  
    743613    Ring_Fitchrom(ksi, ksieps, ns, sfbuf, sdbuf, ksidkpL, ksiimax);
    744614}
    745 /*******************************************************************************
    746  *
    747  *
    748  *
    749  *
    750  ******************************************************************************/
     615
    751616void FitDisp(long q, long pos, double eta) {
    752     long i=0L, nq=0L;
     617    long i, nq;
    753618    fitvect qbuf;
    754619
    … …  
    760625    Ring_FitDisp(pos, eta, dispeps, nq, qbuf, dispdkL, dispimax);
    761626}
    762 /*******************************************************************************
    763  *
    764  *
    765  *
    766  *
    767  ******************************************************************************/
     627
    768628#define nfloq           4
    769629
    … …  
    775635#undef nfloq
    776636
    777 /*******************************************************************************
    778  *
    779  *
    780  *
    781  *
    782  ******************************************************************************/
    783637#define ntrack          4
    784638
    … …  
    797651
    798652     */
    799     long int i=0L;
    800     double twoJx=0.0, twoJy=0.0, phix=0.0, phiy=0.0, scl_1 = 1.0, scl_2 = 1.0;
     653    long int i;
     654    double twoJx, twoJy, phix, phiy, scl_1 = 1.0, scl_2 = 1.0;
    801655    Vector x0, x1, x2, xf;
    802656    FILE *outf;
    … …  
    956810#undef ntrack
    957811
    958 /*******************************************************************************
    959  *
    960  *
    961  *
    962  *
    963  ******************************************************************************/
    964812#define step            0.1
    965813#define px              0.0
    966814#define py              0.0
    967815void track_(double r, struct LOC_getdynap *LINK) {
    968     long i=0L, lastn=0L, lastpos=0L;
     816    long i, lastn, lastpos;
    969817    Vector x;
    970818
    … …  
    1041889void Trac(double x, double px, double y, double py, double dp, double ctau,
    1042890        long nmax, long pos, long &lastn, long &lastpos, FILE *outf1) {
    1043  
    1044     bool lostF = true; /* Lost particle Flag */
     891    bool lostF; /* Lost particle Flag */
    1045892    Vector x1; /* tracking coordinates */
    1046893    Vector2 aperture;
    … …  
    1060907    lastn = 0L;
    1061908    (lastpos) = pos;
     909    lostF = true;
    1062910
    1063911    if (trace)
    … …  
    1128976#define nfloq     4
    1129977bool chk_if_lost(double x0, double y0, double delta, long int nturn, bool floqs) {
    1130  
    1131     long int i=0L, lastn=0L, lastpos=0L;
     978    long int i, lastn, lastpos;
    1132979    Vector x;
    1133980
    … …  
    12591106 ****************************************************************************/
    12601107void getcsAscr(void) {
    1261     long i=0L, j=0L;
    1262     double phi=0.0;
     1108    long i, j;
     1109    double phi;
    12631110    Matrix R;
    12641111
    … …  
    13311178     Assumes mid-plane symmetry.                    */
    13321179
    1333     long int i=0L, lastpos=0L;
    1334     double phi=0.0, x_min=0.0, x_max=0.0, y_min=0.0, y_max=0.0;
     1180    long int i, lastpos;
     1181    double phi, x_min, x_max, y_min, y_max;
    13351182
    13361183    getcod(delta, lastpos);
    … …  
    14221269 ****************************************************************************/
    14231270double get_aper(int n, double x[], double y[]) {
    1424     int i=0;
    1425     double A=0.0;
     1271    int i;
     1272    double A;
    14261273
    14271274    A = 0.0;
    … …  
    14361283}
    14371284
    1438 /**************************************************************************
    1439  * void GetTrack(const char *file_name, long *n, double x[], double px[],
    1440         double y[], double py[])
    1441  *
    1442  *
    1443  *
    1444  **************************************************************************/
    14451285void GetTrack(const char *file_name, long *n, double x[], double px[],
    14461286        double y[], double py[]) {
    1447     int k=0;
     1287    int k;
    14481288    char line[200];
    14491289    FILE *inf;
    … …  
    14951335 ****************************************************************************/
    14961336void Getj(long n, double *x, double *px, double *y, double *py) {
    1497     long int i=0L;
     1337    long int i;
    14981338
    14991339    for (i = 0; i < n; i++) {
    … …  
    15291369 ****************************************************************************/
    15301370double GetArg(double x, double px, double nu) {
    1531    
    1532   double phi=0.0, val=0.0;
     1371    double phi, val;
    15331372
    15341373    phi = GetAngle(x, px);
    … …  
    15681407void GetPhi(long n, double *x, double *px, double *y, double *py) {
    15691408    /* Calculates the linear phase */
    1570     long i=0L;
     1409    long i;
    15711410
    15721411    for (i = 1; i <= n; i++) {
    … …  
    15821421void Sinfft(int n, double xr[]) {
    15831422    /* DFT with sine window */
    1584     int i=0;
     1423    int i;
    15851424    double xi[n];
    15861425
    … …  
    16111450    free_dvector(xi, 1, 2 * n);
    16121451}
    1613 /*******************************************************************************
    1614  *
    1615  *
    1616  *
    1617  *
    1618  ******************************************************************************/
     1452
    16191453void sin_FFT(int n, double xr[], double xi[]) {
    16201454    /* DFT with sine window */
    1621     int i=0;
     1455    int i;
    16221456    double *xri;
    16231457
    … …  
    16361470    free_dvector(xri, 1, 2 * n);
    16371471}
    1638 /*******************************************************************************
    1639  *
    1640  *
    1641  *
    1642  *
    1643  ******************************************************************************/
     1472
    16441473void GetInd(int n, int k, int *ind1, int *ind3) {
    16451474    if (k == 1) {
    … …  
    16541483    }
    16551484}
    1656 /*******************************************************************************
    1657  *
    1658  *
    1659  *
    1660  *
    1661  ******************************************************************************/
     1485
    16621486void GetInd1(int n, int k, int *ind1, int *ind3) {
    16631487    if (k == 1) {
    … …  
    16721496    }
    16731497}
    1674 /*******************************************************************************
    1675  *
    1676  *
    1677  *
    1678  *
    1679  ******************************************************************************/
     1498
    16801499void GetPeak(int n, double *x, int *k) {
    16811500    /* Locate peak in DFT spectrum */
    1682     int ind1=0, ind2=0, ind3=0;
    1683     double peak=0.0;
    1684 
     1501    int ind1, ind2, ind3;
     1502    double peak;
     1503
     1504    peak = 0.0;
    16851505    *k = 1;
    16861506    for (ind2 = 1; ind2 <= n / 2 + 1; ind2++) {
    … …  
    16931513    }
    16941514}
    1695 /*******************************************************************************
    1696  *
    1697  *
    1698  *
    1699  *
    1700  ******************************************************************************/
     1515
    17011516void GetPeak1(int n, double *x, int *k) {
    17021517    /* Locate peak in DFT spectrum */
    1703     int ind1=0, ind2=0, ind3=0;
    1704     double peak=0.0;
    1705 
     1518    int ind1, ind2, ind3;
     1519    double peak;
     1520
     1521    peak = 0.0;
    17061522    *k = 1;
    17071523    for (ind2 = 1; ind2 <= n; ind2++) {
    … …  
    17141530    }
    17151531}
    1716 /*******************************************************************************
    1717  *
    1718  *
    1719  *
    1720  *
    1721  ******************************************************************************/
     1532
    17221533double Int2snu(int n, double *x, int k) {
    17231534    /* Get frequency by nonlinear interpolation with two samples
    … …  
    17281539     N           A(k-1) + A(k)   2
    17291540     */
    1730     int ind=0, ind1=0, ind3=0;
    1731     double ampl1=0.0, ampl2=0.0;
     1541    int ind, ind1, ind3;
     1542    double ampl1, ampl2;
    17321543
    17331544    GetInd(n, k, &ind1, &ind3);
    … …  
    18231634     2      pi ( k + 1/2 )          pi ( k - 1/2 )
    18241635     */
    1825     double corr=0.0;
     1636    double corr;
    18261637
    18271638    corr = (Sinc(M_PI * (k - 1 + 0.5 - nu * n)) + Sinc(M_PI * (k - 1 - 0.5 - nu
    … …  
    18581669    /* Get phase by linear interpolation for rectangular window
    18591670     with -pi <= phi <= pi */
    1860     int i=0;
    1861     double phi=0.0;
     1671    int i;
     1672    double phi;
    18621673    double xr[n], xi[n];
    18631674
    … …  
    19011712 ****************************************************************************/
    19021713void FndRes(struct LOC_findres *LINK) {
    1903     int i=0, j=0, FORLIM=0, FORLIM1=0;
    1904     double delta=0.0;
     1714    int i, j, FORLIM, FORLIM1;
     1715    double delta;
    19051716
    19061717    FORLIM = LINK->n;
    … …  
    19701781    } while (!V.found);
    19711782}
    1972 /*******************************************************************************
    1973  *
    1974  *
    1975  *
    1976  *
    1977  ******************************************************************************/
     1783
    19781784void GetPeaks(int n, double *x, int nf, double *nu, double *A) {
    1979     int i=0, k=0, ind1=0, ind3=0;
     1785    int i, k, ind1, ind3;
    19801786
    19811787    for (i = 0; i < nf; i++) {
    … …  
    19911797    }
    19921798}
    1993 /*******************************************************************************
    1994  *
    1995  *
    1996  *
    1997  *
    1998  ******************************************************************************/
     1799
    19991800void GetPeaks1(int n, double *x, int nf, double *nu, double *A) {
    2000     int i=0, k=0, ind1=0, ind3=0;
     1801    int i, k, ind1, ind3;
    20011802
    20021803    for (i = 0; i < nf; i++) {
    … …  
    20181819
    20191820void SetTol(int Fnum, double dxrms, double dyrms, double drrms) {
    2020     int i=0;
    2021     long k=0L;
     1821    int i;
     1822    long k;
    20221823
    20231824    for (i = 1; i <= GetnKid(Fnum); i++) {
    … …  
    20331834    }
    20341835}
    2035 /*******************************************************************************
    2036  *
    2037  *
    2038  *
    2039  *
    2040  ******************************************************************************/
     1836
    20411837void Scale_Tol(int Fnum, double dxrms, double dyrms, double drrms) {
    2042     int Knum=0;
    2043     long int loc=0L;
     1838    int Knum;
     1839    long int loc;
    20441840
    20451841    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
    … …  
    20791875 ****************************************************************************/
    20801876void SetaTol(int Fnum, int Knum, double dx, double dy, double dr) {
    2081     long int loc=0L;
     1877    long int loc;
    20821878
    20831879    loc = Elem_GetPos(Fnum, Knum);
    … …  
    20911887    Mpole_SetdT(Fnum, Knum);
    20921888}
    2093 /*******************************************************************************
    2094  *
    2095  *
    2096  *
    2097  *
    2098  ******************************************************************************/
     1889
    20991890void ini_aper(const double Dxmin, const double Dxmax, const double Dymin,
    21001891        const double Dymax) {
    2101     int k=0;
     1892    int k;
    21021893
    21031894    for (k = 0; k <= globval.Cell_nLoc; k++) {
    … …  
    21081899    }
    21091900}
    2110 /*******************************************************************************
    2111  *
    2112  *
    2113  *
    2114  *
    2115  ******************************************************************************/
     1901
    21161902void set_aper(const int Fnum, const double Dxmin, const double Dxmax,
    21171903        const double Dymin, const double Dymax) {
    2118     int i=0;
    2119     long int loc=0L;
     1904    int i;
     1905    long int loc;
    21201906
    21211907    for (i = 1; i <= GetnKid(Fnum); i++) {
    … …  
    21271913    }
    21281914}
    2129 /*************************************************
    2130  * void LoadApertures(const char *ChamberFileName)
    2131  *
    2132  *
    2133  *
    2134  **************************************************/
     1915
    21351916void LoadApertures(const char *ChamberFileName) {
    21361917    char line[128], FamName[32];
    2137     long Fnum=0L;
    2138     double Xmin=0.0, Xmax=0.0, Ymin=0.0, Ymax=0.0;
     1918    long Fnum;
     1919    double Xmin, Xmax, Ymin, Ymax;
    21391920    FILE *ChamberFile;
    21401921
    … …  
    21541935    fclose(ChamberFile);
    21551936}
    2156 /**********************************************************************************
    2157  * void LoadTolerances(const char *TolFileName)
    2158  *
    2159  *
    2160  *
    2161  **********************************************************************************/
     1937
    21621938// Load tolerances from the file
    21631939void LoadTolerances(const char *TolFileName) {
    21641940    char line[128], FamName[32];
    2165     int Fnum=0;
    2166     double dx=0.0, dy=0.0, dr=0.0;
     1941    int Fnum;
     1942    double dx, dy, dr;
    21671943    FILE *tolfile;
    21681944
    21691945    tolfile = file_read(TolFileName);
    2170     if(tolfile == NULL){
    2171       printf("LoadTolerances(): Error! Failure to read file: %s \n",TolFileName);
    2172       exit_(1);
    2173     }
    21741946
    21751947    do
    … …  
    21931965}
    21941966
    2195 /**************************************************************************************
    2196  * void ScaleTolerances(const char *TolFileName, const double scl)
    2197  *
    2198  *
    2199  *
    2200  * ************************************************************************************/
    22011967// Load tolerances from the file
    22021968void ScaleTolerances(const char *TolFileName, const double scl) {
    22031969    char line[128], FamName[32];
    2204     int Fnum=0;
    2205     double dx=0.0, dy=0.0, dr=0.0;
     1970    int Fnum;
     1971    double dx, dy, dr;
    22061972    FILE *tolfile;
    22071973
    22081974    tolfile = file_read(TolFileName);
    2209     if(tolfile == NULL){
    2210       printf("LoadTolerances(): Error! Failure to read file: %s \n",TolFileName);
    2211       exit_(1);
    2212     }
    2213    
     1975
    22141976    do
    22151977        fgets(line, 128, tolfile);
    … …  
    22301992    fclose(tolfile);
    22311993}
    2232 /*******************************************************************************
    2233  *
    2234  *
    2235  *
    2236  *
    2237  ******************************************************************************/
     1994
    22381995void SetKpar(int Fnum, int Knum, int Order, double k) {
    22391996
    … …  
    22411998    Mpole_SetPB(Fnum, Knum, Order);
    22421999}
    2243 /*******************************************************************************
    2244  *
    2245  *
    2246  *
    2247  *
    2248  ******************************************************************************/
     2000
    22492001void SetL(int Fnum, int Knum, double L) {
    22502002
    22512003    Cell[Elem_GetPos(Fnum, Knum)].Elem.PL = L;
    22522004}
    2253 /*******************************************************************************
    2254  *
    2255  *
    2256  *
    2257  *
    2258  ******************************************************************************/
     2005
    22592006void SetL(int Fnum, double L) {
    2260     int i=0;
     2007    int i;
    22612008
    22622009    for (i = 1; i <= GetnKid(Fnum); i++)
    22632010        Cell[Elem_GetPos(Fnum, i)].Elem.PL = L;
    22642011}
    2265 /*******************************************************************************
    2266  *
    2267  *
    2268  *
    2269  *
    2270  ******************************************************************************/
     2012
    22712013void SetdKpar(int Fnum, int Knum, int Order, double dk) {
    22722014
    … …  
    22742016    Mpole_SetPB(Fnum, Knum, Order);
    22752017}
    2276 /*******************************************************************************
    2277  *
    2278  *
    2279  *
    2280  *
    2281  ******************************************************************************/
     2018
    22822019void SetKLpar(int Fnum, int Knum, int Order, double kL) {
    2283     long int loc=0L;
     2020    long int loc;
    22842021
    22852022    if (abs(Order) > HOMmax) {
    … …  
    22962033    Mpole_SetPB(Fnum, Knum, Order);
    22972034}
    2298 /*******************************************************************************
    2299  *
    2300  *
    2301  *
    2302  *
    2303  ******************************************************************************/
     2035
    23042036void SetdKLpar(int Fnum, int Knum, int Order, double dkL) {
    2305     long int loc=0L;
     2037    long int loc;
    23062038
    23072039    loc = Elem_GetPos(Fnum, Knum);
    … …  
    23322064**************************************************************/
    23332065void SetdKrpar(int Fnum, int Knum, int Order, double dkrel) {
    2334     long int loc=0L;
     2066    long int loc;
    23352067
    23362068    loc = Elem_GetPos(Fnum, Knum);
    … …  
    23432075    Mpole_SetPB(Fnum, Knum, Order);
    23442076}
    2345 /*******************************************************************************
    2346  *
    2347  *
    2348  *
    2349  *
    2350  ******************************************************************************/
     2077
    23512078void Setbn(int Fnum, int order, double bn) {
    2352     int i=0;
     2079    int i;
    23532080
    23542081    for (i = 1; i <= GetnKid(Fnum); i++)
    23552082        SetKpar(Fnum, i, order, bn);
    23562083}
    2357 /*******************************************************************************
    2358  *
    2359  *
    2360  *
    2361  *
    2362  ******************************************************************************/
     2084
    23632085void SetbnL(int Fnum, int order, double bnL) {
    2364     int i=0;
     2086    int i;
    23652087
    23662088    for (i = 1; i <= GetnKid(Fnum); i++)
    23672089        SetKLpar(Fnum, i, order, bnL);
    23682090}
    2369 /*******************************************************************************
    2370  *
    2371  *
    2372  *
    2373  *
    2374  ******************************************************************************/
     2091
    23752092void Setdbn(int Fnum, int order, double dbn) {
    2376     int i=0;
     2093    int i;
    23772094
    23782095    for (i = 1; i <= GetnKid(Fnum); i++)
    23792096        SetdKpar(Fnum, i, order, dbn);
    23802097}
    2381 /*******************************************************************************
    2382  *
    2383  *
    2384  *
    2385  *
    2386  ******************************************************************************/
     2098
    23872099void SetdbnL(int Fnum, int order, double dbnL) {
    2388     int i=0;
     2100    int i;
    23892101
    23902102    for (i = 1; i <= GetnKid(Fnum); i++) {
    … …  
    24152127//void Setbnr(int Fnum, long order, double bnr) {
    24162128void Setbnr(int Fnum, int order, double bnr) {
    2417     int i=0;
     2129    int i;
    24182130
    24192131    for (i = 1; i <= GetnKid(Fnum); i++)
    24202132        SetdKrpar(Fnum, i, order, bnr);
    24212133}
    2422 /*******************************************************************************
    2423  *
    2424  *
    2425  *
    2426  *
    2427  ******************************************************************************/
     2134
    24282135void SetbnL_sys(int Fnum, int Order, double bnL_sys) {
    2429     int Knum=0;
    2430     long int loc=0L;
     2136    int Knum;
     2137    long int loc;
    24312138
    24322139    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
    … …  
    24402147    }
    24412148}
    2442 /*******************************************************************************
    2443  *
    2444  *
    2445  *
    2446  *
    2447  ******************************************************************************/
     2149
    24482150void set_dbn_rel(const int type, const int n, const double dbn_rel) {
    2449     long int j=0L;
    2450     double dbn=0.0;
     2151    long int j;
     2152    double dbn;
    24512153
    24522154    printf("\n");
    … …  
    25352237        return (Cell[loc].Elem.M->PBpar[Order + HOMmax]);
    25362238}
    2537 /*******************************************************************************
    2538  *
    2539  *
    2540  *
    2541  *
    2542  ******************************************************************************/
     2239
    25432240void SetdKLrms(int Fnum, int Order, double dkLrms) {
    2544     long int Knum=0L, loc=0L;
     2241    long int Knum, loc;
    25452242
    25462243    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
    … …  
    25552252    }
    25562253}
    2557 /*******************************************************************************
    2558  *
    2559  *
    2560  *
    2561  *
    2562  ******************************************************************************/
     2254
    25632255void Setdkrrms(int Fnum, int Order, double dkrrms) {
    2564     long int Knum=0L, loc=0L;
     2256    long int Knum, loc;
    25652257
    25662258    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
    … …  
    25762268    }
    25772269}
    2578 /*******************************************************************************
    2579  *
    2580  *
    2581  *
    2582  *
    2583  ******************************************************************************/
     2270
    25842271void SetKL(int Fnum, int Order) {
    2585     long int Knum=0L;
     2272    long int Knum;
    25862273
    25872274    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++)
    25882275        Mpole_SetPB(Fnum, Knum, Order);
    25892276}
    2590 /*******************************************************************************
    2591  *
    2592  *
    2593  *
    2594  *
    2595  ******************************************************************************/
     2277
    25962278void set_dx(const int type, const double sigma_x, const double sigma_y) {
    2597     long int j=0L;
     2279    long int j;
    25982280
    25992281    printf("\n");
    … …  
    26112293        }
    26122294}
    2613 /*******************************************************************************
    2614  *
    2615  *
    2616  *
    2617  *
    2618  ******************************************************************************/
     2295
    26192296void SetBpmdS(int Fnum, double dxrms, double dyrms) {
    2620     long int Knum, loc=0L;
     2297    long int Knum, loc;
    26212298
    26222299    for (Knum = 1; Knum <= GetnKid(Fnum); Knum++) {
    … …  
    26442321****************************************************************************/
    26452322void codstat(double *mean, double *sigma, double *xmax, long lastpos, bool all) {
    2646     long i=0L, j=0L, n=0L;
     2323    long i, j, n;
    26472324    Vector2 sum, sum2;
    2648     double TEMP=0.0;
     2325    double TEMP;
    26492326
    26502327    n = 0;
    … …  
    27112388void CodStatBpm(double *mean, double *sigma, double *xmax, long lastpos,
    27122389        long bpmdis[mnp]) {
    2713     long i=0L, j=0L, m=0L, n=0L;
     2390    long i, j, m, n;
    27142391    Vector2 sum, sum2;
    2715     double TEMP=0.0;
     2392    double TEMP;
    27162393
    27172394    m = n = 0;
    … …  
    28712548 ****************************************************************************/
    28722549void GetMean(long n, double *x) {
    2873     long i=0L;
     2550    long i;
    28742551    double mean = 0e0;
    28752552
    … …  
    32802957    const int ntrial = 40; // maximum number of trials for closed orbit
    32812958    const double tolx = 1e-8; // numerical precision
    3282     int k=0;
    3283     int dim=0; // 4D or 6D tracking
    3284     long lastpos=0L;
     2959    int k;
     2960    int dim; // 4D or 6D tracking
     2961    long lastpos;
    32852962
    32862963    vcod = dvector(1, 6);
    … …  
    33623039    const int ntrial = 40; // maximum number of trials for closed orbit
    33633040    const double tolx = 1e-10; // numerical precision
    3364     int k=0, dim = 0;
    3365     long lastpos=0L;
     3041    int k, dim = 0;
     3042    long lastpos;
    33663043
    33673044    // initializations
    … …  
    34293106
    34303107void computeFandJS(double *x, int n, double **fjac, double *fvect) {
    3431     int i=0, k=0;
     3108    int i, k;
    34323109    long lastpos = 0L;
    34333110    Vector x0, fx, fx1, fx2;
    … …  
    35013178 ****************************************************************************/
    35023179void computeFandJ(int n, Vector &x, Matrix &fjac, Vector &fvect) {
    3503     int i=0, k=0;
     3180    int i, k;
    35043181    long lastpos = 0;
    35053182    Vector x0, fx1, fx2;
    … …  
    35833260
    35843261void Newton_RaphsonS(int ntrial, double x[], int n, double tolx) {
    3585     int k=0, i=0, *indx;
    3586     double errx=0.0, d=0.0, *bet, *fvect, **alpha;
     3262    int k, i, *indx;
     3263    double errx, d, *bet, *fvect, **alpha;
    35873264
    35883265    errx = 0.0;
    … …  
    36813358 ****************************************************************************/
    36823359int Newton_Raphson(int n, Vector &x, int ntrial, double tolx) {
    3683     int k=0, i=0;
    3684     double errx=0.0;
     3360    int k, i;
     3361    double errx;
    36853362    Vector bet, fvect;
    36863363    Matrix alpha;
    … …  
    37243401    return 0;
    37253402}
    3726 /*******************************************************************************
    3727  *
    3728  *
    3729  *
    3730  *
    3731  ******************************************************************************/
     3403
    37323404void rm_mean(long int n, double x[]) {
    3733     long int i=0L;
    3734     double mean=0.0;
    3735 
     3405    long int i;
     3406    double mean;
     3407
     3408    mean = 0.0;
    37363409    for (i = 0; i < n; i++)
    37373410        mean += x[i];

  • branches/tracy3-3.10.1b/tracy/tracy/src/prtmfile.cc

    r11 r23  
    153153
    154154 Purpose:
    155  Print out flatfile;
    156  print the lattice parameters in an external file.
     155 Print out flatfile
    157156
    158157 Input:
    … …  
    182181            Cell[i].dS[Y_], Cell[i].Elem.M->PdTpar, Cell[i].Elem.M->PdTsys
    183182                + Cell[i].Elem.M->PdTrms * Cell[i].Elem.M->PdTrnd);
    184         fprintf(mfile, " %23.16e %23.16e %23.16e %23.16e %23.16e  %23.16e %23.16e\n",
     183        fprintf(mfile, " %23.16e %23.16e %23.16e %23.16e %23.16e\n",
    185184            Cell[i].Elem.PL, Cell[i].Elem.M->Pirho, Cell[i].Elem.M->PTx1,
    186             Cell[i].Elem.M->PTx2, Cell[i].Elem.M->PH1,Cell[i].Elem.M->PH2,Cell[i].Elem.M->Pgap);
     185            Cell[i].Elem.M->PTx2, Cell[i].Elem.M->Pgap);
    187186        prtHOM(mfile, Cell[i].Elem.M->n_design, Cell[i].Elem.M->PB,
    188187            Cell[i].Elem.M->Porder);

  • branches/tracy3-3.10.1b/tracy/tracy/src/rdmfile.cc

    r11 r23  
    217217                if (prt)
    218218                    printf("%s\n", line);
    219                 sscanf(line, "%lf %lf %lf %lf %lf  %lf %lf", &Cell[i].Elem.PL,
     219                sscanf(line, "%lf %lf %lf %lf %lf", &Cell[i].Elem.PL,
    220220                        &Cell[i].Elem.M->Pirho, &Cell[i].Elem.M->PTx1,
    221                         &Cell[i].Elem.M->PTx2, &Cell[i].Elem.M->PH1, &Cell[i].Elem.M->PH2, &Cell[i].Elem.M->Pgap);
     221                        &Cell[i].Elem.M->PTx2, &Cell[i].Elem.M->Pgap);
    222222                if (Cell[i].Elem.M->Pirho != 0.0)
    223223                    Cell[i].Elem.M->Porder = 1;

  • branches/tracy3-3.10.1b/tracy/tracy/src/read_script.cc

    r11 r23  
    2525   Comments:
    2626       Written by Jianfeng Zhang 03/2011 soleil
    27        
    28        (1) 15/10/2013 by Jianfeng Zhang @ LAL
    29            Fix the bugs of fgets():
    30               to read lines with arbritracy length,
    31               to read the empty space at the end of the line.
    32              
    33            Remove the rontine to move the address of "line" if the
    34            contents of "line" is empty space; otherwise free() can't
    35            find the orignial address of "line".
    36            
    3727
    3828****************************************************************************/
    … …  
    4030/* files */
    4131
    42 char    lat_file[max_str]="voidlattice";
    4332// girder file
    4433char girder_file[max_str];
    … …  
    7867 
    7968 
    80   char    str[max_str]="voidstring", dummy[max_str]="voidstring",dummy2[max_str]="voidstring", nextpara[max_str]="voidpara";
     69  char    str[max_str]="voidstring", dummy[max_str]="voidstring", nextpara[max_str]="voidpara";
    8170  char    in[max_str];  //temporary line with preceding white space
    82   char    *line = NULL; //line to store the command without preceding white space
    83   size_t   len = 0;
    84   ssize_t  read;
     71  char    *line; //line to store the command without preceding white space
    8572  char    name[max_str]="voidname";
    86  // char    lat_file[max_str]="voidlattice";
     73  char    lat_file[max_str]="voidlattice";
    8774  char    EndName[]="void";
    88 
     75 
    8976  FILE    *inf;
    9077  const bool  prt = false; // for debugging printout each line of input file
    … …  
    10087  // Manipulation of the parameter file
    10188  strcpy(dummy, param_file_name);
    102 
    103   sprintf(full_param_file_name,"%s",dummy);
    104   if (!prt) printf("\n reading in script file: %s\n",full_param_file_name);
     89  pch = strstr (dummy,".prm"); // search for extension .prm
     90  /* remove additional .prm if exist */
     91  if (pch != NULL) strncpy(pch,"\0",1);
     92
     93  sprintf(full_param_file_name,"%s.prm",dummy);
     94  if (prt) printf("\n reading in script file: %s\n",full_param_file_name);
    10595
    10696  //
    10797  inf = file_read(full_param_file_name);
    108    if(inf == NULL){
    109      printf("read_script(): Error! Read File %s Failure!\n",full_param_file_name);
    110      exit_(1); 
    111   }
     98
    11299 
    113100  // read parameter file, line by line
    114    while ((read = getline(&line, &len, inf)) != -1) {
    115      
    116      LineNum++;
    117      if(prt){
    118        printf("Line # %ld \n",LineNum);
    119        printf("Retrieved line of length %zu : \n",read);
    120        cout << line << endl;
    121      }
    122                
    123      
     101 //  while (fgets(line, max_str, inf) != NULL) {
     102   while (line = fgets(in, max_str, inf)) {
     103   
     104     /* kill preceding whitespace generated by "table" key
     105        or "space" key, but leave \n
     106        so we're guaranteed to have something*/
     107     while(*line == ' ' || *line == '\t') {
     108       line++;
     109     }       
     110   
     111    if(prt)
     112      printf("cfg: %s",line);
     113     
     114      LineNum++;
    124115           
    125116    /* read the end of line symbol '\n','\r' or '\r\n' at different operation system*/
    126     if (strstr(line, "#") == NULL && line[0] != '\n' &&
    127         line[0] != '\r' && strcmp(line, "\r\n") != 0 ){
     117    if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0 &&
     118        strcmp(line,"\r") != 0 && strcmp(line,"\r\n") != 0) {
    128119      // get initial command token
    129120      sscanf(line, "%s", name);
    … …  
    260251          strcpy(UserCommandFlag[CommNo].CommandStr,name);
    261252      } 
    262       //print the cooridinates using tracking at each element
    263       else if (strcmp("PrintTrackElemFlag", name) == 0){
    264           sscanf(line, "%*s %s",nextpara);
    265          
    266           if(strcmp(nextpara,"voidpara")!=0)
    267             sscanf(line, "%*s %s %lf %lf %lf %lf %lf %lf %ld %ld",
    268                         UserCommandFlag[CommNo]._PrintTrackElem_track_file,
    269                         &(UserCommandFlag[CommNo]._PrintTrackElem_x),
    270                         &(UserCommandFlag[CommNo]._PrintTrackElem_px),
    271                         &(UserCommandFlag[CommNo]._PrintTrackElem_y),
    272                         &(UserCommandFlag[CommNo]._PrintTrackElem_py),
    273                         &(UserCommandFlag[CommNo]._PrintTrackElem_delta),
    274                         &(UserCommandFlag[CommNo]._PrintTrackElem_ctau),
    275                         &(UserCommandFlag[CommNo]._PrintTrackElem_nelem1),
    276                         &(UserCommandFlag[CommNo]._PrintTrackElem_nelem2));
    277           strcpy(UserCommandFlag[CommNo].CommandStr,name);
    278       }
    279253      //print close orbit(COD) flag
    280254      else if (strcmp("PrintGirderFlag", name) == 0){
    … …  
    293267      else if (strcmp("FmapFlag", name) == 0){       
    294268          strcpy(dummy, "");
    295           strcpy(dummy2,"");
     269         
    296270          sscanf(line, "%*s %s",nextpara);
    297271         
    298           //if no definition of loss flag in the lattice file, then "dummy2" is empty string.
    299272          if(strcmp(nextpara,"voidpara")!=0)
    300             sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s",
     273            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s",
    301274                          UserCommandFlag[CommNo]._FmapFlag_fmap_file,
    302275                          &(UserCommandFlag[CommNo]._FmapFlag_nxpoint),
    … …  
    306279                          &(UserCommandFlag[CommNo]._FmapFlag_ymax),
    307280                          &(UserCommandFlag[CommNo]._FmapFlag_delta),
    308                           dummy,dummy2);
     281                          dummy);
    309282         
    310283        if(strcmp(dummy, "true") == 0)
    … …  
    313286          UserCommandFlag[CommNo]._FmapFlag_diffusion = false;
    314287       
    315         if(strcmp(dummy2, "true") == 0)
    316           UserCommandFlag[CommNo]._FmapFlag_printloss = true;
    317         else if(strcmp(dummy2, "false") == 0)
    318           UserCommandFlag[CommNo]._FmapFlag_printloss = false;
    319         else
    320           UserCommandFlag[CommNo]._FmapFlag_printloss = false;
    321        
    322        
    323         //cout << "debug:      " << line << endl;
    324         //cout << "dummy2 =  " << dummy2 << "     lossflag =  " << UserCommandFlag[CommNo]._FmapFlag_printloss << endl;
    325        
    326288       // FmapFlag = true;
    327289         strcpy(UserCommandFlag[CommNo].CommandStr,name);
    … …  
    330292      else if (strcmp("FmapdpFlag", name) == 0){         
    331293          strcpy(dummy, "");
    332           strcpy(dummy2,"");
    333294          sscanf(line, "%*s %s",nextpara);
    334           
     295         
    335296          if(strcmp(nextpara,"voidpara")!=0)
    336             sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s %s",
     297            sscanf(line, "%*s %s %ld %ld %ld %lf %lf %lf %s",
    337298                        UserCommandFlag[CommNo]._FmapdpFlag_fmapdp_file,
    338299                        &(UserCommandFlag[CommNo]._FmapdpFlag_nxpoint),
    … …  
    342303                        &(UserCommandFlag[CommNo]._FmapdpFlag_emax),
    343304                        &(UserCommandFlag[CommNo]._FmapdpFlag_z),
    344                         dummy,dummy2);
     305                        dummy);
    345306         
    346307        if(strcmp(dummy, "true") == 0)
    … …  
    349310          UserCommandFlag[CommNo]._FmapdpFlag_diffusion = false;
    350311
    351         if(strcmp(dummy2, "true") == 0)
    352           UserCommandFlag[CommNo]._FmapdpFlag_printloss = true;
    353         else if(strcmp(dummy2, "false") == 0)
    354           UserCommandFlag[CommNo]._FmapdpFlag_printloss = false;
    355         else
    356           UserCommandFlag[CommNo]._FmapdpFlag_printloss = false;
    357        
    358        
    359         cout << "debug:      " << line << endl;
    360         cout << "dummy2 =  " << dummy2 << "     lossflag =  " << UserCommandFlag[CommNo]._FmapdpFlag_printloss << endl;
    361        
    362312          //  FmapdpFlag = true;
    363313         strcpy(UserCommandFlag[CommNo].CommandStr,name);
    … …  
    585535        exit_(1);
    586536      }
    587   }
     537    }
    588538    /* continue read in the line */
    589539    else
    590      continue;
    591   }//end of while loop
    592  
    593   free(line);
    594  
     540      continue;
     541  }
    595542  fclose(inf);
    596543} 

  • branches/tracy3-3.10.1b/tracy/tracy/src/soleillib.cc

    r11 r23  
    2525       none
    2626
    27        
    2827   Global variables:
    2928       trace
    … …  
    3837void Get_Disp_dp(void)
    3938{
    40   long i=0L;
     39  long i;
    4140//  long lastpos = 0;
    4241  const char nomfic[] = "dispersion.out";
    … …  
    9695{
    9796  Vector        x1;     /* tracking coordinates */
    98   long          i = 0L, k = 0L, imax = 50L;
     97  long          i = 0L, k = 0L, imax = 50;
    9998  FILE *        outf;
    10099  double        dP = 0.0, dP20 = 0.0, dpmax = 0.06;
    101100  Vector2       amp = {0.0, 0.0}, H = {0.0, 0.0};
    102101  const char    nomfic[] = "amp_ind.out";
    103   long          lastpos = 0L;
     102  long          lastpos = 0;
    104103  CellType      Celldebut, Cell;
    105104  Vector        codvector[Cell_nLocMax];
    … …  
    201200{
    202201  CellType Cell;
    203   long i=0L;
     202  long i;
    204203Vector2 H;
    205204
    … …  
    252251{
    253252  CellType Cell;
    254   long i=0L;
     253  long i;
    255254  long lastpos = 1L;
    256255  Vector2 H;
    … …  
    462461       See also LoadApers in nsrl-ii.cc
    463462       J.Zhang 07/10 soleil
    464        
    465       (1) 15/10/2013 Jianfeng Zhang @ LAL
    466          Replace fgets() by getline() to fix the bug
    467          to read arbritrary line length.
    468463****************************************************************************/
    469464void ReadCh(const char *AperFile)
    470465{
    471  char    Name1[max_str], Name2[max_str];
    472  char    *line = NULL;
    473  size_t  len = 0;
    474  ssize_t read;
     466 char    in[max_str], Name1[max_str], Name2[max_str];
     467 char    *line;
    475468  int     Fnum1=0, Fnum2=0, Kidnum1=0, Kidnum2=0, k1=0, k2=0;
    476469  int     i=0, j=0,LineNum=0;
    … …  
    480473
    481474  fp = file_read(AperFile);
    482   if(fp == NULL){
    483     printf("ReadCh(): Error! Failure to read file %s\n", AperFile);
    484     exit_(1);
    485   }
    486  
    487   printf("\n Loading and setting vacuum apertures to lattice elements...\n");
    488 
    489   while ((read=getline(&line, &len, fp)) != -1) {
    490    
    491       /* count the line number that has been read*/
    492     LineNum++;
    493     if(!prt){
    494       printf("Line # %ld : \n",LineNum);
    495       printf("Retrieved line of length %zu :\n",read);
    496       printf("%s\n",line);     
    497     }
    498    
     475
     476  printf("\n");
     477  printf("Loading and setting vacuum apertures to lattice elements...\n");
     478
     479  while (line=fgets(in, max_str, fp)) {
     480  /* kill preceding whitespace generated by "table" key
     481        or "space" key, but leave \n
     482        so we're guaranteed to have something*/
     483     while(*line == ' ' || *line == '\t') {
     484       line++;
     485     }       
     486    /* count the line number that has been read*/
     487    LineNum++; 
    499488    /* NOT read comment line or blank line with the end of line symbol '\n','\r' or '\r\n'*/
    500     if (strstr(line, "#") == NULL && line[0] != '\n' &&
    501         line[0] != '\r' &&strcmp(line,"\r\n") != 0)
     489    if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0 &&
     490        strcmp(line,"\r") != 0 &&strcmp(line,"\r\n") != 0)
    502491    /* read the aperture setting */
    503492    {
    … …  
    570559 //     printf("%s", line);
    571560  }
    572   free(line);
    573561  fclose(fp);
    574562  // turn on the global flag for CheckAmpl()
    … …  
    617605  bool lostF = true; /* Lost particle Flag */
    618606  Vector x1;            /* tracking coordinates */
    619   long i=0L;
     607  long i;
    620608  Vector2  aperture;
    621609  aperture[0] = 1e0; aperture[1] = 1e0;
    … …  
    811799#undef nterm
    812800
    813 /*******************************************************************************
    814  *
    815  *
    816  *
    817  *
    818  ******************************************************************************/
     801
    819802double get_D(const double df_x, const double df_y)
    820803{
    821   double  D=0.0;
     804  double  D;
    822805
    823806  const double D_min = -2.0, D_max = -10.0;
    … …  
    851834   Input:
    852835       FmapFile file to save calculated frequency map analysis
    853        Nbx                horizontal step number
    854        Nby                vertical step number
    855        xmax               horizontal maximum amplitude
    856        zmax              vertical maximum amplitude
    857        Nbtour            number of turn for tracking
    858        energy            particle energy offset
    859        diffusion        flag to calculate tune diffusion/ tune
    860                         difference between the first Nbtour and last Nbtour
     836       Nbx    horizontal step number
     837       Nby    vertical step number
     838       xmax   horizontal maximum amplitude
     839       zmax   vertical maximum amplitude
     840       Nbtour number of turn for tracking
     841       energy particle energy offset
    861842       matlab  set file print format for matlab post-process; specific for nsrl-ii
    862843       
    863844   Output:
    864        status2 true if stable
     845       status true if stable
    865846              false otherwise
    866847
    … …  
    895876 int nb_freq[2] = {0, 0};
    896877 long nturn = Nbtour;
    897  bool status2 = true;
     878 bool status = true;
    898879 struct tm *newtime;
    899880
    … …  
    946927 //    z  = z0 + sgn(j)*sqrt((double)abs(j))*zstep;
    947928     // tracking around closed orbit
    948      Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
    949      if (status2) { // if trajectory is stable
     929     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
     930     if (status) { // if trajectory is stable
    950931       // gets frequency vectors
    951932       Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    … …  
    993974#undef NTERM2
    994975
    995 /****************************************************************************/
    996 /* void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    997           double energy, bool diffusion, bool loss)
    998 
    999    Purpose:
    1000        
    1001          
    1002        Compute a frequency map of Nbx x Nbz points
    1003        For each set of initial conditions the particle is tracked over
    1004        Nbtour for an energy offset dp
    1005 
    1006        Frequency map is based on fixed beam energy, trace x versus z,
    1007        or, tracking transverse dynamic aperture for fixed momentum
    1008        (usually, on-momentum) particle.
    1009        
    1010        The stepsize follows a square root law
    1011 
    1012        Results in fmap.out
    1013 
    1014    Input:
    1015        FmapFile file to save calculated frequency map analysis
    1016        Nbx                horizontal step number
    1017        Nby                vertical step number
    1018        xmax               horizontal maximum amplitude
    1019        zmax              vertical maximum amplitude
    1020        Nbtour            number of turn for tracking
    1021        energy            particle energy offset
    1022        diffusion        flag to calculate tune diffusion/ tune
    1023                         difference between the first Nbtour and last Nbtour
    1024        matlab  set file print format for matlab post-process; specific for nsrl-ii
    1025        
    1026    Output:
    1027        status2 true if stable
    1028               false otherwise
    1029 
    1030    Return:
    1031        none
    1032 
    1033    Global variables:
    1034        none
    1035 
    1036    Specific functions:
    1037        Trac_Simple, Get_NAFF
    1038 
    1039    Comments:
    1040        15/10/03 run for the diffusion: nasty patch for retrieving the closed orbit
    1041        16/02/03 patch removed
    1042        19/07/11 add interface of file defined by user which is used to save calculated
    1043                 frequency map analysis
    1044                
    1045        11/06/2013  Modified by Jianfeng Zhang @ LAL
    1046        The same as famp(onst char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    1047           double energy, bool diffusion); but with a flag to print/not print the final information of
    1048           the particle; if the final turn is not the "Nbtour", then the particle is lost.
    1049 ****************************************************************************/
    1050 #define NTERM2  10
    1051 void fmap(const char *FmapFile, long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    1052           double energy, bool diffusion, bool printloss)         
    1053 {
    1054  FILE * outf;   //file with the tunes at the grid point
    1055  FILE *outf_loss; //file with the information of the lost particle
    1056  long i = 0L, j = 0L;
    1057  double Tab[DIM][NTURN], Tab0[DIM][NTURN];
    1058  double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
    1059  double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
    1060  double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
    1061  const double ctau = 0.0;
    1062  double xstep = 0.0, zstep = 0.0;
    1063  double nux1 = 0.0, nuz1 = 0.0, nux2 = 0.0, nuz2 = 0.0;
    1064  int nb_freq[2] = {0, 0};
    1065  long nturn = Nbtour;
    1066  bool status2 = true;
    1067  
    1068  struct tm *newtime;
    1069  
    1070 //variables of the returned the tracked particle
    1071  long lastn = 0;
    1072  long lastpos = 1;
    1073  ss_vect<double> x1;
    1074  
    1075  
    1076  //if(loss){
    1077  char FmapLossFile[S_SIZE+5]=" ";
    1078  strcpy(FmapLossFile,FmapFile);
    1079  strcat(FmapLossFile,".loss");
    1080  //}
    1081  
    1082  /* Get time and date */
    1083  time_t aclock;
    1084  time(&aclock);                 /* Get time in seconds */
    1085  newtime = localtime(&aclock);  /* Convert time to struct */
    1086 
    1087  if (trace) printf("Entering fmap ... results in %s\n\n",FmapFile);
    1088 
    1089  /* Opening file */
    1090  if ((outf = fopen(FmapFile, "w")) == NULL) {
    1091    fprintf(stdout, "fmap: error while opening file %s\n", FmapFile);
    1092    exit_(1);
    1093  }
    1094 
    1095 
    1096  fprintf(outf,"# TRACY III -- %s -- %s \n", FmapFile, asctime2(newtime));
    1097  fprintf(outf,"# nu = f(x) \n");
    1098 // fprintf(outf,"#    x[mm]          z[mm]           fx             fz"
    1099 //       "            dfx            dfz      D=log_10(sqrt(df_x^2+df_y^2))\n");
    1100 //
    1101          fprintf(outf,"#    x[m]          z[m]           fx             fz"
    1102          "            dfx            dfz\n");
    1103 
    1104 
    1105 //file with lost particle information
    1106 if(printloss){
    1107 
    1108   if ((outf_loss = fopen(FmapLossFile, "w")) == NULL) {
    1109       fprintf(stdout, "fmap: error while opening file %s\n", FmapLossFile);
    1110       exit_(1);
    1111    }
    1112 
    1113  fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
    1114  fprintf(outf_loss,"# information of the lost particle: "
    1115                    "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
    1116 
    1117         fprintf(outf_loss,"#   x[m]    z[m]   Nturn   lostp    S[m]       x[m]"
    1118          "       xp[rad]    z[m]      zp[rad]     delta      ctau\n");
    1119  }
    1120          
    1121  if ((Nbx < 1) || (Nbz < 1))
    1122    fprintf(stdout,"fmap: Error Nbx=%ld Nbz=%ld\n",Nbx,Nbz);
    1123  
    1124  // steps in both planes
    1125  xstep = xmax/sqrt((double)Nbx);
    1126  zstep = zmax/sqrt((double)Nbz);
    1127 
    1128  // double number of turn if diffusion to compute
    1129  if (diffusion) nturn = 2*Nbtour;
    1130 
    1131  // px and pz zeroed
    1132  xp = xp0;
    1133  zp = zp0;
    1134 
    1135 // Tracking part + NAFF
    1136  for (i = 0; i <= Nbx; i++) {
    1137  x  = x0 + sqrt((double)i)*xstep;
    1138 // for (i = -Nbx; i <= Nbx; i++) {
    1139  //  x  = x0 + sgn(i)*sqrt((double)abs(i))*xstep;
    1140 //   if (!matlab) fprintf(outf,"\n");
    1141    fprintf(stdout,"\n");
    1142    for (j = 0; j<= Nbz; j++) {
    1143    z  = z0 + sqrt((double)j)*zstep;
    1144  //  for (j = -Nbz; j<= Nbz; j++) {
    1145  //    z  = z0 + sgn(j)*sqrt((double)abs(j))*zstep;
    1146  
    1147    //print out the lost information
    1148    if(printloss)
    1149      Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
    1150      else
    1151      // tracking around closed orbit
    1152      Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
    1153 
    1154    if (status2) { // if trajectory is stable
    1155        // gets frequency vectors
    1156        Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    1157        Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
    1158        if (diffusion) { // diffusion
    1159          // shift data for second round NAFF
    1160          Get_Tabshift(Tab,Tab0,Nbtour,Nbtour);
    1161          // gets frequency vectors
    1162          Get_NAFF(NTERM2, Nbtour, Tab0, fx2, fz2, nb_freq);
    1163          Get_freq(fx2,fz2,&nux2,&nuz2); // gets nux and nuz
    1164        }
    1165      } // unstable trajectory
    1166      else { //zeroing output
    1167       nux1 = 0.0; nuz1 = 0.0;
    1168       nux2 = 0.0; nuz2 = 0.0;
    1169      }
    1170      
    1171      // printout value
    1172      if (!diffusion){
    1173 //       fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n",
    1174 //             1e3*x, 1e3*z, nux1, nuz1);
    1175 //       fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
    1176 //             1e3*x, 1e3*z, nux1, nuz1);
    1177                fprintf(outf,"%10.6e %10.6e %10.6e %10.6e\n",
    1178                x, z, nux1, nuz1);
    1179        fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n",
    1180                x, z, nux1, nuz1);
    1181      }
    1182      else {
    1183        dfx = nux1 - nux2; dfz = nuz1 - nuz2;
    1184        fprintf(outf,"%10.6e %10.6e %10.6e %10.6e %10.6e %10.6e\n",
    1185                x, z, nux1, nuz1, dfx, dfz);
    1186        fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n",
    1187                x, z, nux1, nuz1, dfx, dfz);
    1188      }
    1189    
    1190      //print out the information of the lost particle
    1191      if(printloss)
    1192         fprintf(outf_loss, "%6.3e %6.3e %8ld %3d %+12.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
    1193                             x, z, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
    1194                             x1[2], x1[3], x1[4], x1[5]);
    1195 
    1196    }
    1197  }
    1198 
    1199  fclose(outf);
    1200  if(printloss)
    1201    fclose(outf_loss);
    1202 }
    1203 #undef NTERM2
    1204 
    1205976
    1206977/****************************************************************************/
    … …  
    12531024#define NTERM2  10
    12541025void fmap_p(const char *FmapFile_p, long Nbx, long Nbz, long Nbtour, double xmax,
    1255             double zmax, double energy, bool diffusion, bool printloss, int numprocs, int myid)   
     1026            double zmax, double energy, bool diffusion, int numprocs, int myid)   
    12561027{
    12571028 FILE * outf;
    1258  FILE *outf_loss; //file with the information of the lost particle
    12591029 long i = 0L, j = 0L;
    12601030 double Tab[DIM][NTURN], Tab0[DIM][NTURN];
    1261  double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
     1031 double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
    12621032 double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
    12631033 double x0 = 1e-6, xp0 = 0.0, z0 = 1e-6, zp0 = 0.0;
    … …  
    12671037 int nb_freq[2] = {0, 0};
    12681038 long nturn = Nbtour;
    1269  bool status2 = true;
     1039 bool status = true;
    12701040 struct tm *newtime;
    12711041
    … …  
    12741044 strcat(FmapFile,FmapFile_p);
    12751045 printf("%s\n",FmapFile);
    1276 
    1277 //variables of the returned the tracked particle
    1278  long lastn = 0;
    1279  long lastpos = 1;
    1280  ss_vect<double> x1;
    1281 
    1282  char FmapLossFile[S_SIZE+5]=" ";
    1283  strcpy(FmapLossFile,FmapFile);
    1284  strcat(FmapLossFile,".loss");
    1285 
    12861046
    12871047 /* Get time and date */
    … …  
    12991059 }
    13001060
    1301 //file with lost particle information
    1302  if(printloss){
    1303    if ((outf_loss = fopen(FmapLossFile, "w")) == NULL) {
    1304      fprintf(stdout, "fmap: error while opening file %s\n", FmapLossFile);
    1305      exit_(1);
    1306    }
    1307  }
    1308 
    13091061 if(myid==0)
    13101062 {
    … …  
    13121064  fprintf(outf,"# nu = f(x) \n");
    13131065  fprintf(outf,"#    x[m]          z[m]           fx             fz            dfx            dfz\n");
    1314  
    1315   if(printloss){
    1316     fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapLossFile, asctime2(newtime));
    1317     fprintf(outf_loss,"# information of the lost particle: "
    1318             "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
    1319    
    1320     fprintf(outf_loss,"#   x[m]    z[m]   Nturn   lostp    S[m]       x[m]"
    1321             "       xp[rad]    z[m]      zp[rad]     delta      ctau\n");
    1322   }
    1323 }
     1066 }
    13241067 
    13251068 if ((Nbx < 1) || (Nbz < 1))
    … …  
    13701113 
    13711114     // tracking around closed orbit
    1372      if(printloss)
    1373        Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
    1374      else
    1375        Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status2);
    1376 
    1377 
    1378      if (status2) // if trajectory is stable
     1115     Trac_Simple4DCOD(x,xp,z,zp,energy,ctau,nturn,Tab,&status);
     1116     if (status) // if trajectory is stable
    13791117     {
    13801118       // gets frequency vectors
    … …  
    14101148     fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n", x, z, nux1, nuz1, dfx, dfz);
    14111149     }
    1412 
    1413  //print out the information of the lost particle
    1414      if(printloss)
    1415         fprintf(outf_loss, "%6.3e %6.3e %8ld %3d %+12.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
    1416                             x, z, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1],
    1417                             x1[2], x1[3], x1[4], x1[5]);
    14181150   }
    14191151 }
    14201152
    14211153 fclose(outf);
    1422 
    1423  if(printloss)
    1424    fclose(outf_loss);
    1425 
    14261154}
    14271155#undef NTERM2
    … …  
    14531181       emax             maximum energy
    14541182       z                vertical amplitude
    1455        diffusion        flag to calculate tune diffusion/ tune
    1456                         difference between the first Nbtour and last Nbtour
     1183       diffusion        flag to calculate tune diffusion
    14571184       matlab  set file print format for matlab post-process; specific for nsrl-ii
    14581185   Output:
    1459        status2 true if stable
     1186       status true if stable
    14601187              false otherwise
    14611188
    … …  
    14741201       is negative for negative enrgy offset since this is true for the cod
    14751202       19/07/11  add features to save calculated fmapdp in the user defined file
    1476 
    1477        18/06/2013   by Jianfeng Zhang @ LAL
    1478 
    1479       Add bool flag to print out the last information      of the tracked particle
    14801203****************************************************************************/
    14811204#define NTERM2  10
    14821205void fmapdp(const char *FmapdpFile, long Nbx, long Nbe, long Nbtour, double xmax, double emax,
    1483               double z, bool diffusion, bool printloss)
     1206              double z, bool diffusion)
    14841207{
    14851208 FILE * outf;
    1486 FILE *outf_loss; //file with the information of the lost particle
    14871209 long i = 0L, j = 0L;
    14881210 double Tab[DIM][NTURN], Tab0[DIM][NTURN];
    1489  double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
     1211 double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
    14901212 double x = 0.0, xp = 0.0, zp = 0.0, dp = 0.0, ctau = 0.0;
    14911213 double x0 = 1e-6, xp0 = 0.0, zp0 = 0.0;
    … …  
    14951217 int nb_freq[2] = {0, 0};
    14961218 long nturn = Nbtour;
    1497  bool status2=true;
     1219 bool status=true;
    14981220 struct tm *newtime;
    1499 
    1500  //variables of the returned the tracked particle
    1501  long lastn = 0;
    1502  long lastpos = 1;
    1503  ss_vect<double> x1;
    1504  
    1505  char FmapdpLossFile[S_SIZE+5]=" ";
    1506  strcpy(FmapdpLossFile,FmapdpFile);
    1507  strcat(FmapdpLossFile,".loss");
    15081221
    15091222 /* Get time and date */
    … …  
    15141227 if (diffusion && globval.Cavity_on == false) nturn = 2*Nbtour;
    15151228
    1516  if (trace) printf("Entering fmapdp ... results in %s\n\n",FmapdpFile);
     1229 if (trace) printf("Entering fmap ... results in %s\n\n",FmapdpFile);
    15171230
    15181231 /* Opening file */
    … …  
    15271240fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx           dfz\n");
    15281241 
    1529 
    1530 //file with lost particle information
    1531 if(printloss){
    1532 
    1533   if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL) {
    1534       fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpLossFile);
    1535       exit_(1);
    1536    }
    1537 
    1538  fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
    1539  fprintf(outf_loss,"# information of the lost particle: "
    1540                    "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
    1541 
    1542         fprintf(outf_loss,"#  dp      x[m]    Nturn   lostp     S[m]       x[m]"
    1543          "      xp[rad]     z[m]     zp[rad]     delta      ctau\n");
    1544  }
    1545 
    15461242 if ((Nbx <= 1) || (Nbe <= 1))
    15471243   fprintf(stdout,"fmapdp: Error Nbx=%ld Nbe=%ld\n",Nbx,Nbe);
    … …  
    15661262        diffusion = false;
    15671263     }   
    1568      else{
     1264     else
    15691265      // x  = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
    15701266         x  = x0 + sqrt((double)j)*xstep;
    1571      }
    1572 
    1573 if(printloss)
    1574      Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
    1575      else
    1576      Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status2);
    1577 
    1578      if (status2) {
     1267     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
     1268     if (status) {
    15791269       Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    15801270       Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
    … …  
    16101300        dp, x, nux1, nuz2, dfx, dfz);
    16111301     }
    1612    
    1613      if(printloss)
    1614        fprintf(outf_loss," %-+8.3f %-+8.3f %-8ld %2d %+12.3f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
    1615                dp, x, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
    1616 
    16171302   }
    16181303 }
    16191304
    16201305 fclose(outf);
    1621 
    1622 if(printloss)
    1623    fclose(outf_loss);
    16241306}
    16251307#undef NTERM2
    … …  
    16561338
    16571339   Output:
    1658        status2 true if stable
     1340       status true if stable
    16591341              false otherwise
    16601342
    … …  
    16711353       14/11/2011  add features to parallel calculate fmapdp.
    16721354                   Merged with the version written by Mao-sen Qiu at Taiwan light source.
    1673                    
    1674        18/06/2013   by Jianfeng Zhang @ LAL
    1675         Add bool flag to print out the last information of the tracked particle           
    16761355****************************************************************************/
    16771356#define NTERM2  10
    16781357void fmapdp_p(const char *FmapdpFile_p, long Nbx, long Nbe, long Nbtour, double xmax,
    1679               double emax, double z, bool diffusion, bool printloss, int numprocs, int myid)
     1358              double emax, double z, bool diffusion, int numprocs, int myid)
    16801359{
    16811360 FILE * outf;
    1682 FILE *outf_loss; //file with the information of the lost particle
    16831361 long i = 0L, j = 0L;
    16841362 double Tab[DIM][NTURN], Tab0[DIM][NTURN];
    1685  double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx=0.0, dfz=0.0;
     1363 double fx[NTERM2], fz[NTERM2], fx2[NTERM2], fz2[NTERM2], dfx, dfz;
    16861364 double x = 0.0, xp = 0.0, zp = 0.0, dp = 0.0, ctau = 0.0;
    16871365 double x0 = 1e-6, xp0 = 0.0, zp0 = 0.0;
    … …  
    16911369 int nb_freq[2] = {0, 0};
    16921370 long nturn = Nbtour;
    1693  bool status2=true;
     1371 bool status=true;
    16941372 struct tm *newtime;
    16951373
    … …  
    16991377 printf("%s\n",FmapdpFile);
    17001378
    1701 //variables of the returned the tracked particle
    1702  long lastn = 0;
    1703  long lastpos = 1;
    1704  ss_vect<double> x1;
    1705 
    1706  char FmapdpLossFile[S_SIZE+5]=" ";
    1707  strcpy(FmapdpLossFile,FmapdpFile);
    1708  strcat(FmapdpLossFile,".loss");
    1709 
    17101379 /* Get time and date */
    17111380 time_t aclock;
    … …  
    17211390   fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpFile);
    17221391   exit_(1);
    1723  }
    1724 
    1725 if(printloss){
    1726   if ((outf_loss = fopen(FmapdpLossFile, "w")) == NULL) {
    1727       fprintf(stdout, "fmapdp: error while opening file %s\n", FmapdpLossFile);
    1728       exit_(1);
    1729    }
    17301392 }
    17311393
    … …  
    17361398     // fprintf(outf,"#    dp[%%]         x[mm]          fx            fz           dfx           dfz\n");
    17371399     fprintf(outf,"#    dp[m]         x[m]           fx            fz           dfx           dfz\n");
    1738  
    1739 
    1740      if(printloss){
    1741        fprintf(outf_loss,"# TRACY III -- %s -- %s \n", FmapdpLossFile, asctime2(newtime));
    1742        fprintf(outf_loss,"# information of the lost particle: "
    1743                "  lostp = 0 (no particle lost /1(horizontal) /2(vertical) /3(longitudinal) plane) \n");
    1744 
    1745        fprintf(outf_loss,"#  dp      x[m]    Nturn   lostp     S[m]       x[m]"
    1746                "      xp[rad]     z[m]     zp[rad]     delta      ctau\n");
    1747      }
    1748   }
     1400   }
    17491401 
    17501402if ((Nbx <= 1) || (Nbe <= 1))
    … …  
    17961448        diffusion = false;
    17971449     }   
    1798      else{
     1450     else
    17991451      // x  = x0 + sgn(j)*sqrt((double)abs(j))*xstep;
    18001452         x  = x0 + sqrt((double)j)*xstep;
    1801      }
    1802 
    1803      if(printloss)
    1804        Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,lastn, lastpos, x1, &status2);
    1805      else
    1806        Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status2);
    1807 
    1808      if (status2) {
     1453     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,nturn,Tab,&status);
     1454     if (status) {
    18091455       Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    18101456       Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
    … …  
    18401486        dp, x, nux1, nuz2, dfx, dfz);
    18411487     }
    1842 
    1843      if(printloss)
    1844        fprintf(outf_loss," %-+8.3f %-+8.3f %-8ld %2d %+12.3f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f %+10.6f \n",
    1845                dp, x, lastn,status.lossplane, Cell[lastpos].S, x1[0],x1[1], x1[2], x1[3], x1[4], x1[5]);
    1846 
    18471488   }
    18481489 }
    18491490
    18501491 fclose(outf);
    1851 
    1852  if(printloss)
    1853    fclose(outf_loss);
    18541492}
    18551493#undef NTERM2
    … …  
    18971535  double nux1 = 0.0, nuz1 = 0.0;
    18981536  int nb_freq[2] = {0, 0};
    1899   bool status2 = true;
     1537  bool status = true;
    19001538  struct tm *newtime;
    19011539
    … …  
    19291567    ctau = ctau0;
    19301568   
    1931     Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,Nbtour,Tab,&status2); // tracking around closed orbit
    1932     if (status2) {
     1569    Trac_Simple4DCOD(x,xp,z,zp,dp,ctau,Nbtour,Tab,&status); // tracking around closed orbit
     1570    if (status) {
    19331571       Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq); // get frequency vectors
    19341572       Get_freq(fx,fz,&nux1,&nuz1);  // gets nux and nuz
    … …  
    19361574    else {
    19371575       nux1 = 0.0; nuz1 = 0.0;
    1938        status2 = true;
     1576       status = true;
    19391577    }
    19401578   
    … …  
    19911629  FILE *outf;
    19921630  const char *fic = phasefile;
    1993   int i=0;
    1994   bool status2=false;
     1631  int i;
     1632  bool status;
    19951633  struct tm *newtime;
    19961634
    … …  
    20251663  }
    20261664 
    2027   Trac_Simple6DCOD(x,xp,y,yp,energy,ctau,Nbtour,Tab,&status2);
    2028 
     1665  Trac_Simple6DCOD(x,xp,y,yp,energy,ctau,Nbtour,Tab,&status);
    20291666  for (i = 0; i < Nbtour; i++) {
    20301667    fprintf(outf,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
    20311668            Tab[0][i],Tab[1][i],Tab[2][i],Tab[3][i],Tab[4][i],Tab[5][i]);
    20321669  }
    2033   //  cout << "status is: " << status2 << endl;
    20341670  fclose(outf);
    20351671}
    … …  
    20691705  FILE *outf;
    20701706  const char  *fic="phasepoly.out";
    2071   long        lastpos = 0L,lastn = 0L;
    2072   int         i=0,j=0;
    2073   double      x=0.0, z=0.0, px=0.0, pz=0.0, delta=0.0, ctau=0.0;
     1707  long        lastpos = 0,lastn = 0;
     1708  int         i,j;
     1709  double      x, z, px, pz, delta, ctau;
    20741710  double      ex = 1368E-9, el = 1.78E-4;
    20751711  double      betax = 9.0, /*betaz = 8.2, */betal = 45.5;
    … …  
    21631799  double start = 0.0, step = 0.0;
    21641800  double x = 0.0, px = 0.0, z = 0.0, pz = 0.0, delta = 0.0, ctau = 0.0;
    2165   bool status2 = true;
     1801  bool status = true;
    21661802  struct tm *newtime;
    21671803
    … …  
    22221858   fprintf(stdout,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
    22231859            x,px,z,pz,delta,ctau);
    2224     Trac_Simple4DCOD(x,px,z,pz,delta,ctau,Nbtour,Tab,&status2);
     1860    Trac_Simple4DCOD(x,px,z,pz,delta,ctau,Nbtour,Tab,&status);
    22251861   for (i = 0; i < Nbtour; i++) {
    22261862      fprintf(outf,"% .5e % .5e % .5e % .5e % .5e % .5e\n",
    … …  
    22651901  FILE *outf;
    22661902  const char fic[] = "check_ampl.out";
    2267   int i=0;
     1903  int i;
    22681904
    22691905  if ((outf = fopen(fic, "w")) == NULL)
    … …  
    23201956  FILE *outf;
    23211957  const char fic[] = "enveloppe.out";
    2322   int i=0,j=0 ;
     1958  int i,j ;
    23231959  CellType Cell;
    23241960
    … …  
    35623198  FILE *fi;
    35633199  const char fic_skew[] = "QT-solamor_2_3.dat";
    3564   int i=0;
     3200  int i;
    35653201  double theta[500]; /* array for skew quad tilt*/
    3566   double b2=0.0, mKL=0.0;
     3202  double b2, mKL;
    35673203  CellType Cell;
    3568   long mOrder=0L;
     3204  long mOrder;
    35693205
    35703206  int nquad = 0;  /* Number of skew quadrupoles */
    … …  
    37673403  struct tm     *newtime;  // for time
    37683404  Vector        codvector[Cell_nLocMax];
    3769   bool          cavityflag=false, radiationflag=false;
     3405  bool          cavityflag, radiationflag;
    37703406  bool          trace=true;
    37713407 
    … …  
    41653801  struct tm     *newtime;  // for time
    41663802  Vector        codvector[Cell_nLocMax];
    4167   bool          cavityflag=false, radiationflag=false;
     3803  bool          cavityflag, radiationflag;
    41683804  bool          trace=true;
    41693805 
    … …  
    46504286 const char xfic[] = "xspectrum.out";
    46514287 const char zfic[] = "zspectrum.out";
    4652  long i=0L, j=0L, k=0L;
     4288 long i, j, k;
    46534289 #define nterm2  20
    46544290 double Tab[6][NTURN], fx[nterm2], fz[nterm2], fx2[nterm2], fz2[nterm2];
    … …  
    46584294 int nb_freq[2] = {0, 0};
    46594295 long nturn = Nbtour;
    4660  bool status2=true;
     4296 bool status=true;
    46614297 struct tm *newtime;
    46624298
    … …  
    46994335   for (j = 0; j<= Nbz; j++) {
    47004336     z  = z0 + sqrt((double)j)*zstep;
    4701      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
    4702      if (status2) {
     4337     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
     4338     if (status) {
    47034339      Get_NAFF(nterm2, Nbtour, Tab, fx, fz, nb_freq);
    47044340     }
    … …  
    47794415            long nmax, long pos, long &lastn, long &lastpos, FILE *outf1)
    47804416{
    4781   bool lostF=false; /* Lost particle Flag */
     4417  bool lostF; /* Lost particle Flag */
    47824418  Vector x1;     /* tracking coordinates */
    47834419  Vector2  aperture;
    … …  
    47924428  getelem(pos-1,&Cell);
    47934429
    4794   if (trace) printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
     4430  if (!trace) printf("dp= % .5e %% xcod= % .5e mm zcod= % .5e mm \n",
    47954431             dp*1e2, Cell.BeamPos[0]*1e3, Cell.BeamPos[2]*1e3);
    47964432
    … …  
    48634499  CellType Cell;
    48644500  int qlist[320];
    4865   int nquad=0, i=0;
     4501  int nquad=0, i;
    48664502  double A = 0.0;
    48674503
    … …  
    48974533/****************************************************************************/
    48984534/* void fmapfull(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    4899    double energy, bool *status2)
     4535   double energy, bool *status)
    49004536
    49014537   Purpose:
    … …  
    49174553
    49184554   Output:
    4919        status2 true if stable
     4555       status true if stable
    49204556              false otherwise
    49214557
    … …  
    49394575 FILE * outf;
    49404576 const char fic[] = "fmapfull.out";
    4941  int i=0, j=0, k=0;
     4577 int i, j, k;
    49424578 double Tab[DIM][NTURN], Tab0[DIM][NTURN];
    49434579 double fx[NTERM], fz[NTERM], fx2[NTERM], fz2[NTERM];
    … …  
    49484584 double nux1[NTERM], nuz1[NTERM],nux2[NTERM], nuz2[NTERM];
    49494585 long nturn = Nbtour;
    4950  bool status2=true;
     4586 bool status=true;
    49514587 struct tm *newtime;
    49524588 char name[14];
    … …  
    50084644   for (j = 0; j<= Nbz; j++) {
    50094645     z  = z0 + sqrt((double)j)*zstep;
    5010      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
    5011 
    5012      if (status2) {
     4646     Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
     4647
     4648     if (status) {
    50134649       Get_NAFF(NTERM, Nbtour, Tab, fx, fz, nb_freq);
    50144650
    … …  
    50914727/****************************************************************************/
    50924728/* void Dyna(long Nbx, long Nbz, long Nbtour, double xmax, double zmax,
    5093    double energy, bool *status2)
     4729   double energy, bool *status)
    50944730
    50954731   Purpose:
    … …  
    51114747
    51124748   Output:
    5113        status2 true if stable
     4749       status true if stable
    51144750              false otherwise
    51154751
    … …  
    51334769  FILE * outf;
    51344770  const char fic[] = "dyna.out";
    5135   long i=0, j=0;
     4771  long i, j;
    51364772  double Tab[6][NTURN], fx[NTERM2], fz[NTERM2];
    51374773  double x = 0.0, xp = 0.0, z = 0.0, zp = 0.0;
    … …  
    51404776  int nb_freq[2] = {0, 0};
    51414777  long nturn = Nbtour;
    5142   bool status2=true;
     4778  bool status=true;
    51434779  struct tm *newtime;
    51444780
    … …  
    51744810    for (j = 0; j<= Nbz; j++) {
    51754811      z  = z0 + sqrt((double)j)*zstep;
    5176       Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
    5177       if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
     4812      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
     4813      if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    51784814      else {
    51794815       fx[0] = 0.0; fz[0] = 0.0;
    51804816      }
    5181       fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
    5182       fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
     4817      fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
     4818      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
    51834819      if (diffusion) {
    5184         if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    5185         fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
    5186         fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
     4820        if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
     4821        fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
     4822        fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
    51874823      }
    51884824    }
    … …  
    51964832    for (j = 0; j<= Nbz; j++) {
    51974833      z  = z0 + sqrt((double)j)*zstep;
    5198       Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status2);
    5199       if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
     4834      Trac_Simple4DCOD(x,xp,z,zp,energy,0.0,nturn,Tab,&status);
     4835      if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    52004836      else {
    52014837       fx[0] = 0.0; fz[0] =0.0;
    52024838      }
    5203       fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
    5204       fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status2);
     4839      fprintf(outf,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
     4840      fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e %d\n", x, z, fx[0], fz[0], status);
    52054841      if (diffusion) {
    5206         if (status2) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
     4842        if (status) Get_NAFF(NTERM2, Nbtour, Tab, fx, fz, nb_freq);
    52074843        fprintf(outf,"%14.6e %14.6e %14.6e %14.6e\n", x, z, fx[0], fz[0]);
    52084844        fprintf(stdout,"%14.6e %14.6e %14.6e %14.6e\n", x, z, fx[0], fz[0]);
    … …  
    52474883  FILE *outf;
    52484884  const char fic[] = "phase2.out";
    5249   long lastpos = 0L,lastn = 0L;
     4885  long lastpos = 0,lastn = 0;
    52504886  struct tm *newtime;
    52514887
    … …  
    53735009  int i,j;
    53745010  bool chroma=true, trace=false;
    5375   bool radiationflag=false, cavityflag=false;
     5011  bool radiationflag, cavityflag;
    53765012  double dP      = 0.0;
    53775013  double diffcmu = 0.0,   /* cos(mu1) - cos(mu2)*/
    … …  
    57955431void PhaseLongitudinalHamiltonien(void)
    57965432{
    5797   long i=0L,j=0L;
     5433  long i,j;
    57985434  const double t = T;        // To get a one turn map
    5799   double phi=0.0, delta=0.0, H0=0.0;
     5435  double phi, delta, H0;
    58005436  long imax = 1000L,         // turn number
    58015437       jmax = 25L;          // starting condition number
    … …  
    59935629}
    59945630
    5995 /*******************************************************************************
    5996  *
    5997  *
    5998  *
    5999  *
    6000  ******************************************************************************/
     5631
    60015632double EnergySmall(double *X, double irho)
    60025633{
    6003  double A=0.0, B=0.0;
     5634 double A, B;
    60045635 double h = irho;
    60055636
    … …  
    60085639 return (A+B);
    60095640}
    6010 /*******************************************************************************
    6011  *
    6012  *
    6013  *
    6014  *
    6015  ******************************************************************************/
     5641
    60165642double EnergyDrift(double *X)
    60175643{
    6018  double A=0.0;
     5644 double A;
    60195645
    60205646 A = (X[1]*X[1]+X[3]*X[3])/2.0/(1.0+X[4]);
    … …  
    60555681  FILE *outf;
    60565682  const char fic[] = "enveloppe2.out";
    6057   int i=0,j=0 ;
     5683  int i,j ;
    60585684  CellType Cell;
    60595685  /* Array for Enveloppes */
    … …  
    61855811void set_RFVoltage(const int Fnum, const double V_RF){
    61865812
    6187   int k=0, n = 0;
     5813  int k, n = 0;
    61885814 
    61895815 
    … …  
    62325858       01/2011  Fix the bug for reading the end of line symbol "\n" , "\r",'\r\n'
    62335859                at different operation system
    6234        04/2011  Change the set of 'seed' for rms error in file, now it's mandatory.
    6235        10/2013 Jianfeng Zhang @ LAL
    6236                Fix the bug of fgets() to read lines of arbritrary length.
     5860       04/2011  Change the set of 'seed' for rms error in file, now it's mandatory.     
    62375861*****************************************************************************************************/
    62385862void ReadFieldErr(const char *FieldErrorFile)
    62395863{ 
    6240   bool  rms=false, set_rnd = false;
    6241   char     name[max_str],keywrd[max_str], *prm;
    6242   char    *line = NULL;
    6243   size_t   len = 0;
    6244   ssize_t   read;
     5864  bool  rms, set_rnd = false;
     5865  char    in[max_str], name[max_str],keywrd[max_str], *prm;
     5866  char    *line;
    62455867  int     n = 0;    /* field error order*/
    62465868  int     LineNum = 0;
    6247   int     seed_val = 0; // random seed number for the rms error
     5869  int     seed_val; // random seed number for the rms error
    62485870  double  Bn = 0.0, An = 0.0, r0 = 0.0; /* field error components and radius when the field error is measured */
    62495871  /* conversion number from A to T.m for soleil*/
    … …  
    62545876
    62555877  inf = file_read(FieldErrorFile);
    6256   if(inf == NULL){
    6257     printf("ReadFieldErr(): Error! Failure to read file %s !\n",FieldErrorFile);
    6258     exit_(1);
    6259   }
    6260  
     5878
    62615879  printf("\n");
    62625880  /* read lines*/
    6263   while ((read=getline(&line, &len, inf)) != -1) { 
     5881  while (line=fgets(in, max_str, inf)) {
     5882 
     5883  /* kill preceding whitespace generated by "table" key
     5884        or "space" key, but leave \n
     5885        so we're guaranteed to have something*/
     5886     while(*line == ' ' || *line == '\t') {
     5887       line++;
     5888     }   
    62645889   
    62655890    /* count line number for debug*/
    62665891    LineNum++;
    6267     if(prt){
    6268       printf("Line # %ld \n",LineNum);
    6269       printf("Retrieved line of length %zu : \n",read);
    6270       printf("%s\n",line);
    6271     }
    62725892   
    62735893    /* check the line is whether comment line or null line*/
    … …  
    62875907        /*read and assign the key words and measure radius*/
    62885908          sscanf(line, " %*s %s %lf",keywrd, &r0);
    6289           if (!prt) printf("\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,r0);
     5909          if (prt) printf("\nsetting <%s> multipole error to: %-5s r0 = %7.1le\n",keywrd,name,r0);
    62905910         
    62915911          rms = (strcmp("rms", keywrd) == 0)? true : false;
    … …  
    63115931            sscanf(prm, "%lf", &An);
    63125932         
    6313             if (!prt)
     5933            if (prt)
    63145934              printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
    63155935         
    … …  
    63355955     // printf("%s", line);
    63365956  }
    6337   free(line);
     5957
    63385958  fclose(inf);
    63395959}
    … …  
    63936013      AddFieldValues_fam(Fnum,keywrd, r0, n, Bn, An);
    63946014    else
    6395       printf("AddFieldErrors(): Error! undefined lattice element %s !\n", name);
     6015      printf("SetFieldErrors: undefined element %s\n", name);
    63966016  }
    63976017}
    … …  
    66746294         which also functions as these correctors. 
    66756295                 
    6676        10/2010  Written by Jianfeng Zhang @ SOLEIL
     6296       10/2010  Written by Jianfeng Zhang
    66776297**********************************************************************/
    66786298void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const char *keywrd, const double r0,
    … …  
    66826302  double  bnL = 0.0, anL = 0.0;
    66836303  double  brho = 0.0, conv_strength = 0.0;
    6684   double  corr = 0.0;   /* skew quadrupole horizontal or vertical corrector error, read from a file*/
     6304  double  corr;   /* skew quadrupole horizontal or vertical corrector error, read from a file*/
    66856305  int    corrlistThick[120];   /* index of associated sextupole*/
    66866306 
    … …  
    67986418void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy)
    67996419{
    6800   long      i=0L;
     6420  long      i;
    68016421  iVector2  nq1 = {0,0},nq2 = {0,0}, nq={0,0};
    68026422  Vector2   nu = {0.0, 0.0};
    … …  
    68476467      delta            initial delta relative to closed orbit
    68486468      ctau             initial c*tau relative to closed orbit
    6849       nmax             maximum number of tracking turns
     6469      nmax             maximum number of turns
    68506470     
    68516471     
    … …  
    68596479{
    68606480   
    6861     long int i=0L, pos = 1L;
    6862     long int lastn = 0L, lastpos = 0L;
     6481    long int i,pos = 1;
     6482    long int lastn = 0, lastpos = 0;
    68636483    Vector x0, x1, x2, xf,codvector[Cell_nLocMax];
    68646484    FILE *outf;
    … …  
    69056525    do {
    69066526        pos = 1;//track from first element
    6907         (lastn)++; //number of turns
     6527        (lastn)++;
    69086528        for (i = 0; i < nv_; i++)  //nv_ = 6
    69096529            x1[i] = x2[i];
    69106530
    6911         //tracking for one turn
    69126531        while( pos <= globval.Cell_nLoc){ 
    69136532         
    … …  
    69336552                     pos, lastpos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_],
    69346553                     1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);
    6935           }
    6936            fprintf(outf,"i      name     S    x(mm)      px(mrad)     y(mm)     py(mrad)    delta(%)    z(mm)    Nturn\n");
    6937           fprintf(outf,"%6d  %s   %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
     6554          }     
     6555          fprintf(outf,"%6d  %s %8.4e %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e %4ld \n",
    69386556                     pos,Cell[pos].Elem.PName,Cell[pos].S, 1e3*xf[x_], 1e3*xf[px_],
    69396557                     1e3*xf[y_], 1e3*xf[py_],1e2*xf[delta_], 1e3*xf[ct_],lastn);           
    … …  
    69426560        }//finish of tracking and printing to file 
    69436561       
    6944     } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc)); //track for nmax turns
    6945 
    6946 //     if (globval.MatMeth)
    6947 //         Cell_Pass(0, globval.Cell_nLoc, x1, lastpos);
    6948 
    6949     fclose(outf);   
    6950 }
    6951 /**********************************************************************
    6952 void PrintTrackElem(const char *TrackFile, double x, double px, double y,double py,
    6953            double delta, double ctau, long int nelem1, long int nelem2)
    6954            
    6955   Purpose:         
    6956     Print the coordinates tracked around a lattice element by tracking around COD
    6957 
    6958   Input:
    6959       TrackFile        file to be print
    6960       x                initial x relative to closed orbit
    6961       px               initial px relative to closed orbit
    6962       y                initial y relative to closed orbit
    6963       py               initial py relative to closed orbit
    6964       delta            initial delta relative to closed orbit
    6965       ctau             initial c*tau relative to closed orbit
    6966       nelem1           start lattice index of the tracked element
    6967       nelem2           end lattice index of the tracked element
    6968      
    6969      
    6970   Output:
    6971  
    6972   Comments:
    6973     Written by Jianfeng Zhang  @ LAL 04/2013   
    6974 **********************************************************************/
    6975 void PrintTrackElem(const char *TrackFile, double x, double px, double y,double py,
    6976            double delta, double ctau, long int nelem1, long int nelem2)
    6977 {
    6978    
    6979     long int i=0L, pos = 1L;
    6980     Vector x0, x1;
    6981     FILE *outf;
    6982     struct    tm *newtime;
    6983    
    6984     bool prt = true;
    6985            
    6986     outf = file_write(TrackFile);
    6987     /* Get time and date */
    6988     newtime = GetTime();
    6989 
    6990     fprintf(outf, "# Element tracking using TRACY III-- %s -- %s\n",TrackFile, asctime2(newtime));
    6991     fprintf(outf, "#\n"); 
    6992    // fprintf(outf, "# work tunes: %7.5f %7.5f\n",globval.TotalTune[0], globval.TotalTune[1]);               
    6993    // fprintf(outf, "#    i  ElemName            S           x            p_x           y          p_y");
    6994    // fprintf(outf, "         delta         cdt     NElem \n");
    6995    // fprintf(outf, "#                          [m]        [mm]         [mrad]        [mm]        [mrad]");     
    6996    // fprintf(outf, "        [%%]          [mm]\n");
    6997 
    6998       fprintf(outf,"i      ElemName     S[m]    x(m)      px(rad)     y(m)     py(rad)    delta(-)    z(m)   \n");
    6999     //initial coordinates
    7000     x0[x_] = x;
    7001     x0[px_] = px;
    7002     x0[y_] = y;
    7003     x0[py_] = py;
    7004     x0[delta_] = delta;
    7005     x0[ct_] = ctau;         
    7006    
    7007      for (i = 0; i < nv_; i++)  //nv_ = 6
    7008             x1[i] = x0[i];
    7009 
    7010  //print the coordinates at the element   
    7011         pos = nelem1+1;//first element "1" of Tracy is a default "debut", it's diffirent from the real lattice index
    7012      
    7013         //tracking from the start element to the end element
    7014         while( pos <= nelem2+1){ 
    7015          
    7016           Elem_Pass(pos,x1);           
    7017        
    7018            
    7019           if (prt) {
    7020             printf("%4ld %s %6.4f %16.8f %16.8f %16.8f %16.8f %16.8f %16.8f \n",
    7021                      pos-1,Cell[pos].Elem.PName,Cell[pos].S, x1[x_], x1[px_],
    7022                      x1[y_], x1[py_],x1[delta_], x1[ct_]);
    7023           }
    7024          
    7025           fprintf(outf,"%6d  %s   %8.4e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e \n",
    7026                      pos-1,Cell[pos].Elem.PName,Cell[pos].S, x1[x_], x1[px_],
    7027                      x1[y_], x1[py_],x1[delta_], x1[ct_]);         
    7028        
    7029         pos++;               
    7030         }//finish of tracking and printing to file 
    7031        
    7032  
     6562    } while ((lastn != nmax) && (lastpos == globval.Cell_nLoc));
    70336563
    70346564//     if (globval.MatMeth)

  • branches/tracy3-3.10.1b/tracy/tracy/src/t2elem.cc

    r19 r23  
    1 /*Tracy-3
     1/* Tracy-3
    22
    33 J. Bengtsson, CBP, LBL      1990 - 1994   Pascal version
    … …  
    99 Element propagators.                                                      */
    1010
    11 
    1211double c_1, d_1, c_2, d_2, cl_rad, q_fluct;
    1312double I2, I4, I5, dcurly_H, dI4;
    … …  
    2019
    2120// Dynamic model
    22 
    23 
    2421
    2522/****************************************************************************/
    … …  
    8481
    8582 Purpose:
    86  Get the longitudinal momentum ps; for the exact/approximation Hamiltonian
     83 Get the longitudinal momentum ps
    8784 
    8885    for approximation Hamitonian:
    8986    ps = 1+delta   
    9087   
    91    
    9288    For exact hamitonian:
    9389   
    94     (1) ps = sqrt((1+delta)^2 - px_^2 - py_^2)
    95    
    96     (2) using the check of TEAPOT to check ps < 0
     90    ps
    9791 **********************************************************/
    9892template<typename T>
    … …  
    10195    T p_s, p_s2;
    10296
    103    if (!globval.H_exact)
     97    if (!globval.H_exact)
    10498    p_s = 1.0+x[delta_];
    10599    else {
    … …  
    109103        else {
    110104            // avoid compile warning
    111             //p_s = 0.0;
    112              p_s = 1.0e-20;
    113             cout << "        " << 1+x[delta_] << "   " << x[px_] << "   " << x[py_] << endl;
     105            p_s = 0.0;
    114106            printf("get_p_s: *** Speed of light exceeded!\n");
    115            // exit_(1);
    116            
     107            exit_(1);
    117108        }
    118109    }
    … …  
    124115
    125116 Purpose:
    126  Exact Drift pass in the curvilinear coordinate for the dipole
    127 
    128  See Eqn. (10.26) on page 307 in E. Forest's beam dynamics, a new attitude
    129  
    130  H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta                                 
     117 Drift pass
     118
     119 If H_exact = false, using approximation Hamiltonian:
     120
     121                                   px^2+py^2
     122              H(x,y,l,px,py,delta) = -------------
     123                                       2(1+delta)     
     124       
     125
     126 Otherwise, using exact Hamiltonian
     127
    131128
    132129 Input:
    … …  
    144141
    145142 Specific functions:
    146  
    147  
    148  
    149  Comments:
    150  
    151  Test version....rewritten by Jianfeng zhang @ LAL, 31/07/2013
    152 
    153  ****************************************************************************/
    154 
    155 
    156 template<typename T>
    157 void Drift(double L,double h_bend, ss_vect<T> &x) {
    158     T u;
    159     ss_vect<T> x0;
    160 
    161     x0=x;
    162    
    163     if (globval.H_exact) {
    164       if(get_p_s(x) == 0)
    165         return;
    166 
    167       x[x_] = x0[x_]/cos(L*h_bend)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
    168       x[px_]= x0[px_]*cos(L*h_bend) + get_p_s(x)*sin(L*h_bend);
    169       x[y_] = x0[y_]+x0[px_]*x0[x_]*tan(L*h_bend)/get_p_s(x)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
    170       x[py_]= x0[py_];
    171       x[ct_]= x0[ct_]+((1+x0[delta_])*x0[x_]*tan(L*h_bend))/get_p_s(x)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
    172  
    173    
    174     }
    175 }
    176 
    177 /****************************************************************************/
    178 /* Drift(double L, ss_vect<T> &x)
    179 
    180  Purpose:
    181  Drift pass in the Cartesian coordinate
    182 
    183  If H_exact = false, using approximation Hamiltonian (only valid for big ring):
    184 
    185                                        px^2+py^2
    186               H(x,y,cT,px,py,delta) = -------------
    187                                        2(1+delta)     
    188        
    189 
    190  Otherwise, using the exact Hamiltonian (valid for small ring)
    191  (J. Bengtsson: Tracy-2 user's mannual):
    192  
    193  H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta                                 
    194 
    195  Input:
    196  L:  drift length
    197  &x:  pointer to initial vector: x
    198 
    199  Output:
    200  none
    201 
    202  Return:
    203  none
    204 
    205  Global variables:
    206 
    207 
    208  Specific functions:
    209143
    210144 ****************************************************************************/
    … …  
    215149    T u;
    216150
    217     if (globval.H_exact) {
    218       if(get_p_s(x) == 0)
    219         return;
    220      
    221    u = L/get_p_s(x);
     151    if (!globval.H_exact) {
     152    u = L/(1.0+x[delta_]);
     153    x[ct_] += u*(sqr(x[px_])+sqr(x[py_]))/(2.0*(1.0+x[delta_]));
     154    } else {
     155    u = L/get_p_s(x);
    222156    x[ct_] += u*(1.0+x[delta_]) - L;
    223     } else {
    224      u = L/(1.0+x[delta_]);
    225     x[ct_] += u*(sqr(x[px_])+sqr(x[py_]))/(2.0*(1.0+x[delta_]));
    226     }
    227    
     157    }
    228158    x[x_] += x[px_] * u;
    229159    x[y_] += x[py_] * u;
    … …  
    671601    Correction for magnet gap (longitudinal fringe field)
    672602
    673     Input:
    674      irho:         h = 1/rho [1/m]
    675      phi:         dipole edge (entrance or exit) angle
    676      gap:        full gap between poles
     603     irho h = 1/rho [1/m]
     604     phi  dipole edge angle
     605     gap  full gap between poles
    677606
    678607                             2
    … …  
    683612     K1 is usually 1/2
    684613     K2 is zero here                                                 
    685      For the type of Tanh-like fringe field
    686      
    687 
    688      18/06/2012  Jianfeng Zhang@LAL
    689      
    690      Fix the bug to calculate psi for the sector magnets which has phi=0
     614
     615
    691616
    692617*********************************************************************/
    … …  
    697622    const double k1 = 0.5, k2 = 0.0;
    698623
    699 //     if (phi == 0.0)  //NOT hard edge, but sector magnets!!!
    700 //         psi = 0.0; 
    701 //     else
     624    if (phi == 0.0)  //hard edge
     625        psi = 0.0; 
     626    else
    702627        psi = k1 * gap * irho * (1.0 + sqr(sin(dtor(phi)))) / cos(dtor(phi))
    703628                * (1.0 - k2 * gap * irho * tan(dtor(phi)));
    704629
    705630    return psi;
    706 }
    707 
    708 /*****************************************************************
    709  * Exact map of the sector dipole
    710  *
    711  * Forest's beam dynamics P.360 Eqn. (12.18)
    712  * Also see Laurent's thesis: P219-220.
    713  *
    714  * Written by Jianfeng Zhang @ LAL, 01/08/2013
    715  *
    716  * Test version........................................
    717  *
    718  * ****************************************************************/
    719 template<typename T>
    720 void dipole_kick(double L, double h_ref, double h_bend, ss_vect<T> &x){
    721 
    722   ss_vect<T> x0;
    723   T ps, u, dpxf;
    724    
    725   x0 = x;
    726   ps = get_p_s(x0);
    727  
    728   u = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
    729  
    730   x[py_] = x0[py_];
    731  
    732   x[delta_] =x0[delta_];
    733  
    734   x[px_]= x0[px_]*cos(L*h_ref) + (ps -h_bend*(1/h_ref+x0[x_]))*sin(L*h_ref);
    735  
    736   dpxf = h_ref*sqrt((1+x[delta_])*(1+x[delta_]) - x[px_]*x[px_] - x[py_]*x[py_])-h_bend*(1+h_ref*x[x_]);
    737  
    738  
    739   x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
    740  
    741   x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x[delta_])-sqr(x[px_])-sqr(x[py_]))-dpxf - h_bend);
    742  
    743   x[ct_] += (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
    744  
    745 }
    746 
    747 /*****************************************************************
    748  * pure multipole kick
    749  *
    750  * Forest's beam dynamics P.360 Eqn. (12.18)
    751  * Also see Laurent's thesis: P219-220.
    752  *
    753  * Written by Jianfeng Zhang @ LAL, 01/08/2013
    754  *
    755  * Test version........................................
    756  *
    757  * ****************************************************************/
    758 template<typename T>
    759 void multipole_kick(int Order, double MB[], double L, double h_bend, ss_vect<T> &x){
    760 
    761   int j=0;
    762   ss_vect<T> x0;
    763  
    764   T ByoBrho = 0, BxoBrho = 0, ByoBrho1=0 ;
    765  
    766   if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) {
    767         x0 = x;
    768         /* compute the total magnetic field Bx and By with Horner's rule */
    769         ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e)
    770         BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e)
    771         for (j = Order - 1; j >= 1; j--) {
    772             ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax];
    773             BxoBrho = x0[y_] * ByoBrho + x0[x_] * BxoBrho + MB[HOMmax - j];
    774             ByoBrho = ByoBrho1;
    775         }
    776   }
    777  
    778   x[px_] -= L * (ByoBrho + h_bend*1);
    779         //vertical kick due to the magnets
    780         x[py_] += L * BxoBrho;
    781631}
    782632
    … …  
    786636
    787637  Purpose:
    788         Refer to Tracy 2 manual.
    789         Calculate multipole kick.
    790   ,
    791   (1)     The exact Hamiltonian is
     638        Calculate multipole kick. The Hamiltonian is
    792639
    793640            H = A + B where A and B are the kick part defined by
    794                  
    795                  (exact form of the H of the drift)
    796                  A(x,y,-l,px,py,dP) = -sqrt( (1 + delta)^2 - px^2 - py^2 )  + delta
    797                                      
    798                  (exact form of the H of the kick from the magnet)                                       
    799                  B(x,y,-l,px,py,dP) = -h*x*sqrt( (1 + delta)^2 - px^2 - py^2 ) + h x + h^2*x^2/2 + int(Re(By+iBx)/Brho)   
    800                  
    801                  so this is the exact Hamitonian for small ring.
    802                  
    803 
    804         The kick is given by
    805                                                                                                 e                     
    806          Dp_x = L* (h_ref*sqrt( (1 + delta)^2 - (px^2 + py^2) ) - h_bend - x*h_bend*h_ref -   ---- B_y )
    807                                                                                                p_0                 
    808          Dp_y = L* e/p_0 * B_x
    809          
    810          x    = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * p_x
    811          
    812          y    = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * p_y
    813          
    814          CT   = L* h*x_ref / sqrt( (1 + delta)^2 - px^2 - py^2) * (1+delta)
    815                                                                                            
     641                                         2    2
     642                                       px + py
     643                 A(x,y,-l,px,py,dP) = ---------
     644                                       2(1+dP)
     645                                                   2 2
     646                 B(x,y,-l,px,py,dP) = -h*x*dP + 0.5 h x + int(Re(By+iBx)/Brho)
     647
     648                 so this is the appproximation for large ring
     649                 the chromatic term is missing hx*A
     650
     651
     652                The kick is given by
     653
     654              e L       L delta    L x              e L
     655     Dp_x = - --- B_y + ------- - ----- ,    Dp_y = --- B_x
     656              p_0         rho     rho^2             p_0
     657
    816658    where
    817659        e      1
    … …  
    823665                           /
    824666                           ====
    825  
    826   (2)     The approximate Hamiltonian is
    827 
    828             H = A + B where A and B are defined by
    829                                          2    2
    830                                        px + py
    831                  A(x,y,-l,px,py,dP) = ---------                   (approximate form of the H of the drift)
    832                                       2(1+delta)
    833                                                         2 2
    834                  B(x,y,-l,px,py,dP) = -h*x*delta + 0.5 h x + int(Re(By+iBx)/Brho)   (approximate form of the H of the kick from the magnet)
    835 
    836                  so this is the appproximation for large ring
    837                  the chromatic term is missing hx*A
    838 
    839 
    840                 The kick is given by
    841 
    842               e L       L delta    L x              e L
    843      Dp_x =  --- B_y + ------- - ----- ,    Dp_y = --- B_x
    844               p_0         rho     rho^2             p_0
    845 
    846      //second order of the Hamiltonian motion
    847      x    =  h_ref * x_0 / (1+delta) * p_x
    848      
    849      y    =  h_ref * x_0 / (1+delta) * p_y
    850      
    851     CT    =  h_ref * x_0 / (1+delta) * (p_x^2 + p_y^2)/2/(1+delta)
    852      
    853     where
    854         e      1
    855       --- = -----
    856       p_0   B rho   
    857                            ====
    858                            \
    859       (B_y + iB_x) = B rho  >   (ia_n  + b_n ) (x + iy)^n-1
    860                            /
    861                            ====
    862 
    863                     = B rho [(b_n*x - a_n*y) + i*(bn*y + an*x) + b_(n-1) + i*a_(n-1)] * (x+iy)^n-2
    864                     =......
    865  
     667
    866668 Input:
    867669 Order  maximum non zero multipole component
    868670 MB     array of an and bn, magnetic field components
    869671 L      multipole length
    870  h_bend   B_dipole/Brho is the dipole field normalized by the reference momentum; one of the example is the combined dipole
    871  h_ref   1/rho [m^-1] is the curvature of the reference orbit in the dipoles which is used in curvilinear coordinate,
     672 h_bend   1/rho in curvilinear coordinate,
     673 0    in cartisian cooridinate.
     674 h_ref   1/rho [m^-1]
    872675 x      initial coordinates vector
    873676
    … …  
    885688
    886689 Comments:
    887  06/11/2012     Jianfeng Nadolski @ LAL
    888                 Fix the bug in the kick map of the exact Hamiltonian.
    889                 Add the second order of the approximate Hamiltonian.
     690 none
    890691
    891692 **************************************************************************/
    … …  
    895696    int j;
    896697    T BxoBrho, ByoBrho, ByoBrho1, B[3];
    897     T sqrtpx, dpx, ps2new, psnew;
    898     T n1, n2;
    899698    ss_vect<T> x0, cod;
    900     T u=0.0;
    901    
     699
    902700    if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) {
    903701        x0 = x;
    904         /* compute the total magnetic field Bx and By with Horner's rule */
     702        /* compute field with Horner's rule */
    905703        ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e)
    906704        BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e)
     705
    907706        for (j = Order - 1; j >= 1; j--) {
    908707            ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax];
    … …  
    911710        }
    912711
    913        
    914712        if (globval.radiation || globval.emittance) {
    915713            B[X_] = BxoBrho;
    … …  
    919717        }
    920718
    921         //curvilinear cocodinates, from the dipole components
    922719        if (h_ref != 0.0) {
    923          
    924           //exact Hamiltonian
    925           if (globval.H_exact) {
    926            
    927             if(get_p_s(x0)==0)
    928               return;
    929            
    930             //dipole kick map from the exact Hamiltonian; not symplectic
    931 //          u = L * h_ref * x0[x_] /get_p_s(x0);
    932 //             x[x_] += u * x0[px_];
    933 //             x[y_] += u * x0[py_];
    934 //          x[px_]-= L * (-h_ref * get_p_s(x0) + h_bend + h_ref * h_bend
    935 //                   * x0[x_] + ByoBrho);
    936 //             x[ct_] += u * (1.0+x0[delta_]);
    937              
    938               //kick map due to the dipole vector potential
    939              // x[px_] -= L*(h_bend + h_ref*h_bend*x0[x_]+ByoBrho);
    940  
    941 //          sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
    942 // //   
    943 //      x[px_] = x0[px_]*cos(L*h_ref) + (get_p_s(x0)-h_bend*(h_ref+x0[x_]))*sin(L*h_ref)-L*ByoBrho;
    944 //     
    945 //      dpx = -x0[px_]*h_ref*sin(L*h_ref)+(get_p_s(x0)-h_bend*(h_ref+x0[x_]))*h_ref*cos(L*h_ref);
    946 //     
    947 //      x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[py_]))-dpx - h_bend);
    948 //     
    949 //      x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(sqrtpx*x0[px_])-asin(sqrtpx*x[px_]));
    950 //     
    951 //      x[ct_]+= (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(sqrtpx*x0[px_])-asin(sqrtpx*x[px_]));
    952        
    953         //map of a wedge
    954 //          x[px_] = x0[px_]*cos(L*h_bend) + (get_p_s(x0) - h_ref*x0[x_])*sin(L*h_bend) -L*ByoBrho;
    955 //       
    956 //       sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
    957 //       n1 = x0[px_]*sqrtpx;
    958 //       n2 = x[px_]*sqrtpx;
    959 //       
    960 //          x[x_] = x0[x_]*cos(L*h_bend) + (x0[x_]*x0[px_]*sin(2*L*h_bend)+sqr(sin(L*h_bend))*(2*x0[x_]*get_p_s(x0)-h_ref*sqr(x0[x_])))
    961 //                  /(sqrt(sqr(1+x0[delta_])-sqr(x[px_])-sqr(x0[y_]))+get_p_s(x0)*cos(L*h_bend)-x0[px_]*sin(L*h_bend));
    962 //
    963 //          x[y_] += x0[py_]*L*h_bend*h_ref + x0[py_]*h_ref*(asin(n1) - asin(n2));
    964 //          x[ct_]+= (1+x0[delta_])*L*h_bend*h_ref + (1+x0[delta_])*h_ref*(asin(n1) - asin(n2));
    965 // 
    966           }//approximate Hamiltonian
    967           else {//first order map
    968720            x[px_] -= L * (ByoBrho + (h_bend - h_ref) / 2.0 + h_ref * h_bend
    969721                    * x0[x_] - h_ref * x0[delta_]);
    970722            x[ct_] += L * h_ref * x0[x_];
    971  
    972             //  second order, cubic term of H, H3; has problem when do tracking...works for THOMX!!!!!
    973             // not symplectic????
    974          if(0){
    975             u = L * h_ref * x0[x_] /(1.0+x0[delta_]);
    976             x[x_] += u * x0[px_];
    977             x[y_] += u * x0[py_];
    978          //   x[px_]+= 0.5*u/x0[x_]*(x0[px_]*x0[px_]+x0[py_]*x0[py_]);
    979             x[ct_] += u*(sqr(x0[px_])+sqr(x0[py_]))/(2.0*(1.0+x0[delta_]));
    980            
    981            // cout << "test......2.....  u = " << u<< endl;
    982           }     
    983       }
    984          
    985         }//from other straight magnets
    986         else
    987             x[px_] -= L * (ByoBrho + h_bend*1); //h_bend should be trigger on for combined dipoles
    988 
    989         //vertical kick due to the magnets
     723        } else
     724            x[px_] -= L * (ByoBrho + h_bend);
     725
    990726        x[py_] += L * BxoBrho;
    991      
    992         /*
    993         if(h_ref == 0 && h_bend !=0){
    994         sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
    995         ps2new = sqr(1+x0[delta_]) - sqr(x[px_]) -sqr(x0[py_]);
    996         psnew = sqrt(ps2new);
    997        
    998         x[x_] += 1/h_bend*(psnew - get_p_s(x0));
    999         x[y_] += x0[py_]/h_bend*(asin(x0[px_]*sqrtpx)-asin(x[px_]*sqrtpx));
    1000         x[ct_]+= (1+x0[delta_])/h_bend*(asin(x0[px_]*sqrtpx)-asin(x[px_]*sqrtpx));
    1001         }*/
    1002        
    1003        
    1004727    }
    1005728}
    … …  
    1012735 Compute edge focusing for a dipole
    1013736 There is no radiation coming from the edge
    1014  The standard formula used is : (SLAC-75, K. Brown's first order of fringe field but with the correction of 1/(1+delta) in R43)
     737 The standard formula used is :
    1015738
    1016739       px = px0 +  irho tan(phi) *x0
    … …  
    1029752 gap  = gap of the dipole for longitudinal fringing field (see psi)
    1030753 x    = input particle coordinates
    1031 Entrance:   bool flag, true, then calculate edge focus and fringe field at the entrance of the dipole
    1032                        false,...........................................at the exit fo the dipole
    1033  
     754
    1034755 Output:
    1035756 x    output particle coordinates
    … …  
    1053774
    1054775 April 2011, No energy dependence for H-plane according to SSC-141 note (E. Forest)
    1055  Agreement better with MAD8 (LNLS), SOLEIL small effect (.1). Modified by Laurent Nadolski @ soleil
    1056  Comments: Jianfeng Zhang @ LAL, 05/06/2013
    1057  This model only works for Soleil lattice, but not works for ThomX and SuperB DR lattice.
    1058  
    1059  18/06/2012  Jianfeng Zhang @ LAL
    1060  Add the contribution of x' to the x[py_] for small ring like Thom-X, using Alexandre Loulergue's correction.
    1061  This models works for ThomX, Soleil, and SuperB DR lattices; and
    1062  close to the Forest's model and K. Brown's models.
    1063  
     776 Agreement better with MAD8 (LNLS), SOLEIL small effect (.1)
    1064777 ****************************************************************************/
    1065778template<typename T>
    1066 static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x, bool Entrance) {
    1067  
    1068 
     779static void EdgeFocus(double irho, double phi, double gap, ss_vect<T> &x) {
    1069780    if (true) {
    1070      
    1071        // cout << "Dipole edge effect WITH Alex's  correction: " << endl;
    1072      
    1073       //with the contribution from the entrance angle of the particle at the edge
    1074 // Original model, written by L. Nadolski, this is a reasonal model!!!!
    1075 //    warning: => diverging Taylor map (see SSC-141)
    1076 //         x[px_] += irho * tan(dtor(phi)) * x[x_];
    1077 //         x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]
    1078 //                 / (1.0 + x[delta_]);
    1079 
    1080 
    1081 // K.Brown 2rd order transfer matrix of the dipole fringe field for sector dipoles (phi=0)
    1082 // compared the fringe field model used in Tracy
    1083 // results:  K. Brown's model is wrong for small ring like Thom-X!!!!
    1084 // The original model of tracy is correct for tracking!!!!!!!!!!
    1085   /*
    1086   double psi;
    1087     psi=gap*0.5*irho*(1+sin(dtor(phi))*sin(dtor(phi)))/cos(dtor(phi));
    1088  
    1089     if(Entrance){
    1090    x[x_] += irho/2*x[y_]*x[y_];
    1091    x[px_] += 0.5*x[x_]*x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5*x[y_]*x[y_];
    1092    x[py_] -= -irho*tan(-dtor(psi))*x[y_]+1*x[x_]*x[y_]+irho*x[px_]*x[y_]+ (irho*psi/cos(-dtor(psi))/cos(-dtor(psi))) * x[y_]
    1093                  / (1.0 + x[delta_]);
    1094     }
    1095     else
    1096     {
    1097        x[x_] -= irho/2*x[y_]*x[y_];
    1098    x[px_] += 0.5*x[x_]*x[x_] + irho * tan(dtor(phi)) * x[x_]-0.5*x[y_]*x[y_];
    1099    x[py_] -= -irho*tan(-dtor(psi))*x[y_]+1*x[x_]*x[y_]-irho*x[px_]*x[y_]+ (irho*psi/cos(-dtor(psi))/cos(-dtor(psi))) * x[y_]
    1100                  / (1.0 + x[delta_]);
    1101     }
    1102     */
    1103 
    1104     //add the contribution of the entrance momentum of the particle, from Alex Loulergue.
    1105     if(Entrance){
    1106        x[px_] += irho * tan(dtor(phi) +x[py_]/(1+x[delta_]*1)) * x[x_]/ (1.0 + x[delta_]*0);
    1107        x[py_] -= irho * tan(dtor(phi) + x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)/(1+x[delta_]*0)) * x[y_]
    1108                 / (1.0 + x[delta_]*0);
    1109     }else{
    1110        x[px_] += irho * tan(dtor(phi) - x[py_]/(1+x[delta_]*1)) * x[x_]/ (1.0 + x[delta_]*0);
    1111        x[py_] -= irho * tan(dtor(phi) - x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)/(1+x[delta_]*0)) * x[y_]
    1112                 / (1.0 + x[delta_]*0);
    1113     }
    1114      
    1115     } else {//original one
    1116 
    1117    //cout << "Dipole edge effect WITHOUT Alex's  correction: " << endl;
     781        // warning: => diverging Taylor map (see SSC-141)
    1118782        x[px_] += irho * tan(dtor(phi)) * x[x_];
    1119         x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]/ (1.0 + x[delta_]*1);
    1120     }
    1121 }
    1122 
    1123 
    1124 /******************************************************************
    1125  template<typename T>
    1126 static void BendCurvature(double irho, double H, ss_vect<T> &x)
    1127 
    1128 The leanding order term T211, T233, and T413 on page 117 & page 118
    1129 of the dipole edge effect in SLAC-75, using the symplectic
    1130 form in
    1131  E. Forest: "Fringe field in MAD, Part II: Bend curvature in MAD-X for the
    1132              mudule PTC".
    1133              P.10 Eq. (5) and (42).
    1134  
    1135  For an rectangular dipole, the bending curvature term is zero, since 1/R = 0.
    1136  
    1137  irho:   curvature of the central orbit inside the dipole
    1138  H:      curvature of the entrance/exit pole face of the dipole
    1139  
    1140  Comments:
    1141  Written by Jianfeng Zhang @ LAL, 01/10/2012.
    1142 ******************************************************************/
    1143 template<typename T>
    1144 static void BendCurvature(double irho,  double H, ss_vect<T> &x) {
    1145  
    1146   if(trace)
    1147   cout << "Forest's bend curvature...."<<endl;
    1148  
    1149   T pm2,u, coeff1, coeff2,coeff3;
    1150  
    1151    ss_vect<T> x0;
    1152    
    1153    x0 = x;
    1154  
    1155    pm2 = sqr(1+x0[delta_]) - sqr(x0[px_]);
    1156    u = irho*H/2/pm2;
    1157    
    1158    coeff1 = u * 2*x0[px_]*(1+x0[delta_])/pm2; //
    1159    coeff2 = u * (1+x0[delta_]); //
    1160    coeff3 = u * (pm2 - 2*(1+x0[delta_]))/pm2;
    1161    
    1162  
    1163    x[x_]  += x0[x_]/(1-coeff1*x0[y_]*x0[y_]);
    1164    x[px_] -= coeff2*x0[y_]*x0[y_] - irho*H/2*sqr(x0[x_]);
    1165    x[py_] -= 2*coeff2*x[x_]*x0[y_];
    1166    x[ct_] -= coeff3*x[x_]*x0[y_]*x0[y_];
    1167 }   
    1168    
    1169    
    1170      
    1171  
     783        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]
     784                / (1.0 + x[delta_]);
     785    } else {
     786        x[px_] += irho * tan(dtor(phi)) * x[x_];
     787        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_];
     788    }
     789}
    1172790
    1173791/****************************************************************
    … …  
    1176794
    1177795  Purpose:
    1178      Rotate from the beam coorinate to the dipole face.
    1179      This is a pure geometric operation, no physical meaning.
    1180      
    1181      See P.307 Eq. (10.26) in
    1182      E. Forest "Beam dynamics: A new attitude and framework".
    1183      
     796     the effect of edge focusing of the dipole
    1184797  Input:
    1185798     phi:  entrance or exit angle of the dipole edge
    1186799     x:    initial cooridinate of the particle
    1187      
    1188800  Output:
    1189801 
    … …  
    1207819}
    1208820
    1209 
    1210821/***************************************************************
    1211822template<typename T>
    … …  
    1213824
    1214825  Purpose:
    1215      the effect of longitudinal fringe field using exact Hamitonian.
    1216      This is a hard effect model, the fringe field and edge focus
    1217      happens at the entrance/exit point of the dipole pole face.
    1218 
    1219       include only the first order of irho.
    1220      Input:
     826     the effect of longitudinal fringe field using exact Hamitonian
     827  Input:
    1221828 
    1222829  Output:
    … …  
    1242849        x[ct_] = x1[ct_] - coeff * x1[px_] * sqr(x[y_]) * (1.0 + x1[delta_])
    1243850                / ps3;
    1244     } else {//give the error message and set the unstable x value during the tracking
    1245              //this x value will be discard as an unstable value during the tracking
     851    } else {
    1246852        printf("bend_fringe: *** Speed of light exceeded!\n");
    1247         x[y_] = 10;
    1248         x[x_] = 10;
    1249         x[py_] =10;
    1250         x[ct_] = 10;
    1251         //exit_(1);
    1252     }
    1253 }
    1254 
    1255 /***************************************************************
    1256 template<typename T>
    1257 void bend_fringe(double hb, double gap, ss_vect<T> &x)
    1258 
    1259   Purpose:
    1260      the effect of longitudinal fringe field using exact Hamitonian.
    1261      This is a hard effect model, the fringe field and edge focus
    1262      happens at the entrance/exit point of the dipole pole face.
    1263 
    1264       include the second order of irho.
    1265      
    1266       See E. Forest and et al.
    1267       "Fringe field in MAD Part I: Second order Fringe in MAD-X for the module PTC",
    1268       P. 8-9, Eq. (35) and (34).
    1269  
    1270    Input:
    1271  
    1272   Output:
    1273  
    1274   Return:
    1275  
    1276 ****************************************************************/
    1277 template<typename T>
    1278 void bend_fringe(double hb, double gap, ss_vect<T> &x) {
    1279  
    1280   bool track=false;
    1281  
    1282   if(track) cout << "bend_fringe(): Forest's model"<<endl;
    1283  
    1284     T coeff1,coeff2,coeff3,coeff4;
    1285     T ps, ps2, ps3, ps5;
    1286    
    1287     ss_vect<T> x1;
    1288  
    1289     gap = gap*0.5;
    1290    
    1291     x1 = x;
    1292     ps = get_p_s(x);
    1293    
    1294     if(ps==0)
    1295       return;
    1296    
    1297     ps2 = sqr(ps);
    1298     ps3 = ps * ps2;
    1299     ps5 = ps2*ps3;
    1300    
    1301     coeff1 = 1.0 - 4.0*sqr(hb)*gap*x1[py_]*x1[y_]/ps3;
    1302     coeff2 = sqr(hb)*gap*2*x1[px_]*(2*sqr(x1[px_])-sqr(1.0+x1[delta_]))/ps5 + 
    1303             hb*(ps/(sqr(1+x1[delta_])-sqr(x1[px_])) + 2*sqr(x1[px_])*ps/sqr(sqr(1+x1[delta_])-sqr(x1[px_])));
    1304              
    1305     coeff3 = hb*(x1[px_]/ps)/(1+sqr(x1[py_])/ps2) - sqr(hb)*gap*((ps2+sqr(x1[px_]))/ps3+sqr(x1[px_])/ps2*(ps2+x1[py_])/ps3);
    1306     coeff4 = -sqr(hb)*gap*4*sqr(x1[px_])*(1+x1[delta_])/ps5;
    1307        
    1308     if (coeff1 >= 0.0) {
    1309         x[y_] = 2.0 * x1[y_] / (1.0 + sqrt(coeff1));
    1310     } else {
    1311       //give the error message and set the unstable x value during the tracking
    1312              //this x value will be discard as an unstable value during the tracking
    1313         printf("bend_fringe: *** Speed of light exceeded!\n");
    1314         x[y_] = 10;
    1315         x[x_] = 10;
    1316         x[py_] =10;
    1317         x[ct_] = 10;
    1318      // exit_(1);
    1319     }
    1320 
    1321   x[x_] += 0.5*coeff2*sqr(x[y_]);
    1322   x[py_]-= coeff3*x[y_];
    1323   x[ct_] -= 0.5*coeff4*sqr(x[y_]);
     853        exit_(1);
     854    }
    1324855}
    1325856
    … …  
    1450981    elemtype *elemp;
    1451982    MpoleType *M;
    1452    
    1453     ss_vect<T> x1, x2;
    1454    
    1455983
    1456984    elemp = &Cell.Elem;
    … …  
    1460988    GtoL(x, Cell.dS, Cell.dT, M->Pc0, M->Pc1, M->Ps1);
    1461989
    1462     //entrance of the magnet
    1463990    if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) { /* fringe fields */
    1464991
    … …  
    1468995        }
    1469996
    1470         //fringe field and edge focusing of dipoles
    1471        if (M->Pirho != 0.0 && M->dipEdge_effect1 == 1){
    1472          if (!globval.H_exact) { //big ring */ modified linear term T43 on page 117 & 118 in SLAC-75.
    1473                        EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x,true);
    1474            } else {//small and big rings; Forest's model
    1475           if(M->PH1!=0){
    1476             // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0.
    1477            BendCurvature(M->Pirho, M->PH1, x);
    1478           }
    1479            p_rot(M->PTx1,  x);  //rotate from cartesian cooridnate to curlinear curved beam coordinate;
    1480              // since the map of a sector dipole is used in Tracy 3
    1481            bend_fringe(M->Pirho, M->Pgap, x); //fringe field
    1482          }
    1483        }
     997        if (!globval.H_exact) {
     998            if (M->Pirho != 0.0)
     999                EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x);
     1000        } else {
     1001            p_rot(M->PTx1, x);
     1002            bend_fringe(M->Pirho, x);
     1003        }
    14841004    }
    14851005
    14861006    if (M->Pthick == thick) {
    1487       if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates
    1488       // if (M->n_design == Dip || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates
     1007        if (!globval.H_exact || ((M->PTx1 == 0.0) && (M->PTx2 == 0.0))) {// polar coordinates,curvilinear coordinates
    14891008            h_ref = M->Pirho;
    14901009            dL = elemp->PL / M->PN;
    … …  
    14941013                dL = elemp->PL / M->PN;
    14951014            else
    1496                 dL = 1.0 / M->Pirho * (sin(dtor(M->PTx1)) + sin(dtor(M->PTx2)))  //straight length of the dipole
     1015                dL = 1.0 / M->Pirho * (sin(dtor(M->PTx1)) + sin(dtor(M->PTx2)))
    14971016                        / M->PN;
    14981017        }
    … …  
    15221041
    15231042    case Meth_Second:
    1524      
    1525       //  cout << "Mpole_Pass: Meth_Second not supported" << endl;
    1526        // exit_(0);
    1527        // break;
    1528 
    1529     //specially for the test of the dipole map,see Forest's beam dynamics, p361. eqn.(12.19)
    1530  
    1531     dL1 = 0.5*dL;
    1532     dL2 = dL;
    1533            
    1534             for (seg = 1; seg <= M->PN; seg++) {
    1535               dipole_kick(dL1,M->Pirho,h_ref,x);
    1536               multipole_kick(M->Porder, M->PB, dkL1, M->Pirho, x);
    1537               dipole_kick(dL1,M->Pirho,h_ref,x);
    1538             }
    1539 
    1540            
    1541        
    1542     case Meth_Fourth:  //forth order symplectic integrator
     1043        cout << "Mpole_Pass: Meth_Second not supported" << endl;
     1044        exit_(0);
     1045        break;
     1046
     1047    case Meth_Fourth:
    15431048        if (M->Pthick == thick) {
    15441049            dL1 = c_1 * dL;
    … …  
    15561061                }
    15571062
    1558                 Drift(dL1, x);
     1063                Drift(dL1, x);
    15591064                thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
    15601065                Drift(dL2, x);
    1561                 thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
     1066                thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
    15621067
    15631068                if (globval.emittance && (!globval.Cavity_on) && (M->Pirho
    … …  
    15661071                    dI4 += 4.0 * is_tps<tps>::get_dI4(x);
    15671072                }
    1568                
    1569                
     1073
    15701074                Drift(dL2, x);
    1571                 thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
    1572                 Drift(dL1, x);
     1075                thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
     1076                Drift(dL1, x);
    15731077
    15741078                if (globval.emittance && (!globval.Cavity_on) && (M->Pirho
    … …  
    15941098    }
    15951099
    1596     //exit of the magnets
    15971100    if ((M->Pmethod == Meth_Second) || (M->Pmethod == Meth_Fourth)) {
    1598        
    1599       /* dipole fringe fields */
    1600              if (M->Pirho != 0.0 && M->dipEdge_effect2 == 1){
    1601             if (!globval.H_exact) { //big ring, linear model correction
    1602                                EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x,false);
    1603                }else{
    1604              bend_fringe(-M->Pirho, M->Pgap,x); //Forest's fringe field of the dipole
    1605               p_rot(M->PTx2, x); //rotate back to the cartesian cooridinate
    1606              if(M->PH2!=0){
    1607             // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0.
    1608              BendCurvature(M->Pirho, M->PH2, x);
    1609              }
    1610          }
    1611      }
    1612         //quadrupole fringe field
     1101        /* fringe fields */
     1102        if (!globval.H_exact) {
     1103            if (M->Pirho != 0.0)
     1104                EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x);
     1105        } else {
     1106            bend_fringe(-M->Pirho, x);
     1107            p_rot(M->PTx2, x);
     1108        }
    16131109        if (globval.quad_fringe && (M->PB[Quad + HOMmax] != 0.0) && (M->quadFF2
    16141110                == 1))
    … …  
    29162412    mergeto4by5(M, ah, av);
    29172413}
    2918   /* ************************************************************** 
    2919    * Compute transfert matrix for a quadrupole magnet
     2414
     2415void Mpole_Setmatrix(int Fnum1, int Knum1, double K) {
     2416    /*   Compute transfert matrix for a quadrupole magnet
    29202417     the transfert matrix A is plitted into two part
    29212418     A = AD55xAU55
    … …  
    29292426     corresponding to a half magnet w/ an entrance
    29302427     angle and no exit angle.
    2931      The linear fringe field is taken into account               
    2932      **************************************************************/
    2933 void Mpole_Setmatrix(int Fnum1, int Knum1, double K) {
    2934  
     2428     The linear frindge field is taken into account                    */
    29352429
    29362430    CellType *cellp;
    … …  
    34262920        printf("Elem_GetPos: there are no kids in family %d (%s)\n", Fnum1,
    34272921                ElemFam[Fnum1 - 1].ElemF.PName);
    3428         exit_(1);
     2922        exit_(0);
    34292923    }
    34302924
    … …  
    34572951    d_2 = 1e0 - 2e0 * d_1;
    34582952
    3459    
    3460    
    3461    
    3462    
    34632953    // classical radiation
    34642954    //  C_gamma = 8.846056192e-05;
    … …  
    38533343    M->PTx1 = 0.0; /* Entrance angle */
    38543344    M->PTx2 = 0.0; /* Exit angle */
    3855     M->PH1 = 0.0; /* Entrance pole face curvature*/
    3856     M->PH2 = 0.0; /* Exit pole face curvature */
    38573345    M->Pgap = 0.0; /* Gap for fringe field ??? */
    38583346

  • branches/tracy3-3.10.1b/tracy/tracy/src/t2lat.cc

    r11 r23  
    4444  leq, pwrsym, lss, geo, gtr, period_, charcon, stringcon, ident, geq, lsym,
    4545  bobrhosym, bobrhovsym, bobrhohsym, kxvsym, kxhsym, phisym, ksym,
    46   tsym, t1sym, t2sym,h1sym,h2sym,
     46  tsym, t1sym, t2sym,
    4747  gapsym, thksym, invsym, thnsym,
    4848  endsym, tsksym, bemsym, corsym, prnsym, tblsym, possym, prmsym,
    … …  
    5151  nbdsym, frgsym, latsym, mpsym, dbnsym, kssym, homsym, lmdsym, dtsym, xytsym,
    5252  vrfsym, harnumsym, frqsym, gstsym, typsym, rollsym, idsym,
    53   fnamesym1, fnamesym2, scalingsym1, scalingsym2, fmsym, harmsym, sprsym, recsym, solsym,edge_effect1sym,edge_effect2sym,
     53  fnamesym1, fnamesym2, scalingsym1, scalingsym2, fmsym, harmsym, sprsym, recsym, solsym,
    5454  ff1sym, ff2sym, ffscalingsym, tiltsym
    5555} Lat_symbol;   /*\*/
    5656// idsym fnamesym1 fnamesym2 scalingsym added for insertion
    5757// ring sym added
    58   //edge_effect1sym, edge_effect2sym added, to turn on/off the edge effect of dipoles @ 05/10/2012 by Jianfeng Zhang
    5958// ff1sym ff2sym for quadrupole entrance and exit fringe field
    6059// ffscalingsym scaling factor for entrance and exit fringe field.  /*J.Zhang
    … …  
    658657      if (!strncmp(id, "t              ", sizeof(alfa_)))
    659658        *sym = tsym;
    660       else if (!strncmp(id, "h1            ", sizeof(alfa_)))
    661         *sym = h1sym;
    662       else if (!strncmp(id, "h2            ", sizeof(alfa_)))
    663         *sym = h2sym;
    664659      else if (!strncmp(id, "gap            ", sizeof(alfa_)))
    665660        *sym = gapsym;
    … …  
    10591054        else if (!strncmp(fname, "t2             ", sizeof(partsName)))
    10601055          x = WITH1->PTx2;
    1061         else if (!strncmp(fname, "h1            ", sizeof(partsName)))
    1062           x = WITH1->PH1;
    1063         else if (!strncmp(fname, "h2            ", sizeof(partsName)))
    1064           x = WITH1->PH2;
    10651056        else if (!strncmp(fname, "gap            ", sizeof(partsName)))
    10661057          x = WITH1->Pgap;
    … …  
    20262017  struct LOC_Lat_DealElement V;
    20272018  bool           Result = false;
    2028   double         t=0.0, t1=0.0, t2=0.0, h1=0,h2=0,gap=0.0, QL = 0.0, QK=0.0;
    2029   double         QKick=0.0;  /* kick angle of the corrector [rad]*/
    2030   int            Kplane=0; /* kick plane of the corrector, 1 for H plane, -1 for V plane*/
    2031   double         QKV=0.0, QKH=0.0, QKxV=0.0, QKxH=0.0, QPhi=0.0, QKS=0.0;
    2032   double         dt=0.0, Frf=0.0, Vrf=0.0;
    2033   long            harnum=0L;
    2034   int            k1=1; //number of cut slice of the element
    2035   int            k2=4; //integration method; Meth_Linear=0, Meth_First=1,Meth_Second=2,Meth_Fourth=4
    2036   int            edge_effect1=0,edge_effect2=0; //flag to turn on/off dipole fringe field
    2037   int            FF1=0, FF2=0;
    2038   double         FFscaling=0.0;
     2019  double         t, t1, t2, gap, QL = 0.0, QK;
     2020  double         QKick;  /* kick angle of the corrector [rad]*/
     2021  int            Kplane; /* kick plane of the corrector, 1 for H plane, -1 for V plane*/
     2022  double         QKV, QKH, QKxV, QKxH, QPhi, QKS;
     2023  double         dt, Frf, Vrf;
     2024  long k1, k2, harnum, FF1, FF2;
     2025  double FFscaling;
    20392026  Lat_symbol     sym1;
    20402027  symset         mysys, SET;
    … …  
    20892076    test__(P_expset(SET, 1 << ((long)semicolon)), "<;> expected", &V);
    20902077    GetSym__(&V);
    2091     QL = *V.rnum;
    20922078    globval.Elem_nFam++;
    20932079    if (globval.Elem_nFam <= Elem_nFamMax) {
    … …  
    20952081      WITH1 = &WITH->ElemF;
    20962082      memcpy(WITH1->PName, ElementName, sizeof(partsName));
    2097    //   WITH1->PL = *V.rnum;
    2098        WITH1->PL = QL;
     2083      WITH1->PL = *V.rnum;
    20992084      WITH1->Pkind = PartsKind(drift);
    21002085      Drift_Alloc(&WITH->ElemF);
    … …  
    21162101            T2     = <exit angle>, ( [degree] )
    21172102            gap    = <total magnet gap>, ( [m] )
    2118    edge_effect1    = 1/0, turn on/ff the entrance edge effect,
    2119                      and fringe field (if gap != 0)
    2120    edge_effect2    = 1/0, turn on/ff the exit edge effect,
    2121                      and fringe field (if gap != 0)
    21222103            K      = <K-value>, ( [m^-2] )
    21232104                     ( K > 0 : focusing in horizontal plane )
    … …  
    21352116    Example
    21362117
    2137       B: bending, L=0.70, T=10.0, T1:=5.0, T2:=5.0, H1=0, H2=0, gap=0.03, edge_effect1=1, edge_effect2=1, K=-1.0, N=8, Method=2;
     2118      B: bending, L=0.70, T=10.0, T1:=5.0, T2:=5.0, K=-1.0, N=8, Method=2;
    21382119
    21392120    *************************************************************************/
    … …  
    21432124    QL = 0.0;   /* L */
    21442125    QK = 0.0;   /* K, quadrupole components*/
    2145     k1 = 1;   /* N */
     2126    k1 = 0;   /* N */
    21462127    t  = 0.0;   /* T */
    21472128    t1 = 0.0;   /* T1 */
    21482129    t2 = 0.0;   /* T2 */
    2149     h1 = 0.0;
    2150     h2 = 0.0;
    21512130    gap = 0.0;   /* gap */
    2152     edge_effect1 = 0; /* edge_effect effect at the entrance */
    2153     edge_effect2 = 0; /* edge effect at the exit */
    2154     k2 = Meth_Fourth;   /* method */
     2131    k2 = Meth_Linear;   /* method */
    21552132    dt = 0.0;
    21562133    ClearHOMandDBN(&V);
    … …  
    21632140    P_addset(mysys, (long)t1sym);
    21642141    P_addset(mysys, (long)t2sym);
    2165     P_addset(mysys, (long)h1sym);
    2166     P_addset(mysys, (long)h2sym);
    21672142    P_addset(mysys, (long)gapsym);
    2168     P_addset(mysys, (long)edge_effect1sym);
    2169     P_addset(mysys, (long)edge_effect2sym);
    21702143    P_addset(mysys, (long)homsym);
    21712144    P_addset(mysys, (long)dbnsym);
    … …  
    22052178        break;
    22062179
    2207         case h1sym:
    2208         h1 = EVAL_(&V);
    2209         break;
    2210 
    2211       case h2sym:
    2212         h2 = EVAL_(&V);
    2213         break;
    2214        
    22152180      case gapsym:
    22162181        gap = EVAL_(&V);
    2217         break;
    2218 
    2219       case edge_effect1sym:
    2220         edge_effect1 = (long)EVAL_(&V);
    2221         break;
    2222 
    2223       case edge_effect2sym:
    2224         edge_effect2 = (long)EVAL_(&V);
    22252182        break;
    22262183
    … …  
    22662223      else
    22672224        WITH2->Pirho = t * M_PI / 180.0;
    2268       WITH2->PTx1 = t1; WITH2->PTx2 = t2; //entrance and exit angle of the dipoles (relative to rectangular magnets)
    2269       WITH2->PH1 = h1; WITH2->PH2 = h2; //entrance and exit angle of the dipoles
    2270       WITH2->Pgap = gap;
    2271       WITH2->dipEdge_effect1 = edge_effect1;//edge_effect effect at the entrance
    2272       WITH2->dipEdge_effect2 = edge_effect2;//edge effect at the exit
     2225      WITH2->PTx1 = t1; WITH2->PTx2 = t2; WITH2->Pgap = gap;
    22732226      WITH2->n_design = Dip;
    22742227      AssignHOM(globval.Elem_nFam, &V);
    … …  
    23022255
    23032256      Q1: Quadrupole, L=0.5, K=2.213455, N=4, Method=4;
    2304       Q1 : Quadrupole, L=0.5, K=2.213455, N=4, FF1 =1, FF2=0, FFscaling=1, Method=4;
     2257      Q1 : Quadrupole, L=0.5, K=2.213455, N=4, FF1 =1, FF2, FFscaling=1, Method=4;
    23052258
    23062259    **************************************************************/
    … …  
    23102263    QL = 0.0;   /* L */
    23112264    QK = 0.0;   /* K */
    2312     k1 = 1;   /* N */
     2265    k1 = 0;   /* N */
    23132266    k2 = Meth_Fourth;   /* method */
    23142267    dt = 0.0;
    … …  
    24382391    QL = 0.0;   /* L */
    24392392    QK = 0.0;   /* K */
    2440     k1 = 1;   /* N */
     2393    k1 = 0;   /* N */
    24412394    k2 = Meth_Fourth;   /* method */
    24422395    FF1 = 0;     /* Entrance Fringe field */
    … …  
    26532606    QKick = 0.0; /* kick angle of the corrector [rad]*/
    26542607    Kplane = 0;   /* 1 is horizontal corrector, -1 is vertical corrector */
    2655     k1 = 1;     /* N */
     2608    k1 = 0;     /* N */
    26562609    k2 = Meth_Linear;   /* method */
    26572610    dt = 0.0;
    … …  
    29312884    QL = 0.0;   /* L */
    29322885    QK = 0.0;   /* K */
    2933     k1 = 1;   /* N */
     2886    k1 = 0;   /* N */
    29342887    t = 0.0;   /* T */
    29352888    t1 = 0.0;   /* T1 */
    29362889    t2 = 0.0;   /* T2 */
    2937     h1 = 0.0;
    2938     h2 = 0.0;
    29392890    gap = 0.0;   /* gap */
    2940     k2 = Meth_Fourth;   /* method */
     2891    k2 = Meth_Linear;   /* method */
    29412892    dt = 0.0;
    29422893    ClearHOMandDBN(&V);
    … …  
    29472898    P_addset(mysys, (long)t1sym);
    29482899    P_addset(mysys, (long)t2sym);
    2949     P_addset(mysys, (long)h1sym);
    2950     P_addset(mysys, (long)h2sym);
    29512900    P_addset(mysys, (long)gapsym);
    29522901    P_addset(mysys, (long)rollsym);
    … …  
    29832932        break;
    29842933
    2985       case h1sym:
    2986         h1 = EVAL_(&V);
    2987         break;
    2988 
    2989       case h2sym:
    2990         h2 = EVAL_(&V);
    2991         break;
    2992        
    29932934      case gapsym:
    29942935        gap = EVAL_(&V);
    … …  
    30382979      }
    30392980      WITH2->PN = k1; WITH2->Pmethod = k2;
    3040       WITH2->PTx1 = t1; WITH2->PTx2 = t2;
    3041       WITH2->PH1 = h1; WITH2->PH2 = h2;
    3042       WITH2->Pgap = gap;
     2981      WITH2->PTx1 = t1; WITH2->PTx2 = t2; WITH2->Pgap = gap;
    30432982      WITH2->PdTpar = dt;
    30442983      AssignHOM(globval.Elem_nFam, &V);
    … …  
    30833022    getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
    30843023    GetSym__(&V);
    3085     QL = 0e0;
    3086     QK = 0e0;
    3087     QKV = 0e0;
    3088     QKH = 0e0;
    3089     QKxV = 0e0;
    3090     QKxH = 0e0;
    3091     QPhi = 0e0;
    3092     QKS = 0e0;
    3093     k1 = 1;
    3094     k2 = Meth_Fourth;
    3095     dt = 0e0;
     3024    QL = 0e0; QK = 0e0; QKV = 0e0; QKH = 0e0; QKxV = 0e0; QKxH = 0e0;
     3025    QPhi = 0e0; QKS = 0e0; k1 = 0; k2 = Meth_Linear; dt = 0e0;
    30963026    ClearHOMandDBN(&V);
    30973027    P_addset(P_expset(mysys, 0), (long)lsym);
    … …  
    32193149    getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
    32203150    GetSym__(&V);
    3221     QL = 0.0;
    3222     k1 = 1;
    3223     strcpy(str1, "");
    3224     strcpy(str2, "");
     3151    QL = 0.0; k1 = 0; strcpy(str1, ""); strcpy(str2, "");
    32253152    P_addset(P_expset(mysys, 0), (long)lsym);
    32263153    P_addset(mysys, (long)nsym);
    … …  
    35723499    QL = 0.0; /* L */
    35733500    QK = 0.0; /* K */
    3574     k1 = 1;   /* N */
     3501    k1 = 0;   /* N */
    35753502    P_addset(P_expset(mysys, 0), (long)lsym);
    35763503    P_addset(mysys, (long)bobrhosym);
    … …  
    37313658  Reg("drift          ", drfsym, &V);
    37323659  Reg("dt             ", dtsym, &V);
    3733   Reg("edge_effect1   ", edge_effect1sym, &V);     /* Jianfeng Zhang*/
    3734   Reg("edge_effect2   ", edge_effect2sym, &V);     /* Jianfeng Zhang */
    37353660  Reg("end            ", endsym, &V);
    37363661  Reg("ff1            ", ff1sym, &V);     /* Laurent */
    … …  
    37463671  Reg("gap            ", gapsym, &V);
    37473672  Reg("ghost          ", gstsym, &V);
    3748   Reg("h1             ", h1sym, &V);
    3749   Reg("h2             ", h2sym, &V);
    37503673  Reg("harm           ", harmsym, &V);
    37513674  Reg("harnum         ", harnumsym, &V);

  • branches/tracy3-3.10.1b/tracy/tracy/src/tpsa_lin.cc

    r3 r23  
    414414}
    415415
    416 /****************************
    417  * tan(x)
    418  * *************************/
     416
    419417void datan(const tps_buf &x, tps_buf &z)
    420418{
    … …  
    424422}
    425423
    426 /*************************
    427  *  cotan(x)
    428  * ***********************/
    429 void dactan(const tps_buf &x, tps_buf &z)
    430 {
    431   tps_buf  c, s;
    432 
    433   dacos(x, c); dasin(x, s); dadiv_(c, s, z);
    434 }
    435 
    436 /*************************************
    437  *
    438  * Purpose:
    439  *     Get the derivatives of the arcsin(x)
    440  *   
    441  *    d(arcsin(u))/dx = (1/sqrt(1-u^2))*du/dx;
    442  *
    443  *  Input:
    444  *         x: u
    445  *  Output:
    446  *        z: the differential of x.
    447  *
    448  *       
    449  * Comments:
    450  *           13/11/2012
    451  *          Written by Jianfeng Nadolski @ LAL
    452  *                     
    453  * ***********************************/
    454 void daarcsin(const tps_buf &x, tps_buf &z)
    455 {
    456   int     i;
    457   double  a;  //differential of x; a = sqrt(1/1-u^2)
    458 
    459   a = x[0]; z[0] = asin(a); a = 1.0/sqrt(1.0-sqr(a));
    460   for (i = 1; i <= ss_dim; i++)
    461     z[i] = a*x[i];
    462 }
    463 
    464 /*************************************
    465  *
    466  * Purpose:
    467  *     Get the derivatives of the arccos(x)
    468  *   
    469  *    d(arccos(u))/dx = -(1/sqrt(1-u^2))*du/dx;
    470  *
    471  *  Input:
    472  *         x: u
    473  *  Output:
    474  *        z: the differential of x.
    475  *
    476  *       
    477  * Comments:
    478  *           13/11/2012
    479  *          Written by Jianfeng Nadolski @ LAL
    480  *                     
    481  * ***********************************/
    482 void daarccos(const tps_buf &x, tps_buf &z)
    483 {
    484 int     i;
    485   double  a;  //differential of x; a = -sqrt(1/1-u^2)
    486 
    487   a = x[0]; z[0] = asin(a); a = -1.0/sqrt(1.0-sqr(a));
    488   for (i = 1; i <= ss_dim; i++)
    489     z[i] = a*x[i];
    490  
    491 }
    492 /*************************************
    493  *
    494  * Purpose:
    495  *     Get the derivatives of the arctan(x)
    496  *   
    497  *    d(arctan(u))/dx = (1/1+u^2)*du/dx;
    498  *
    499  *  Input:
    500  *         x: u
    501  *  Output:
    502  *        z: the differential of x.
    503  *
    504  *         
    505  * ***********************************/
     424
    506425void daarctan(const tps_buf &x, tps_buf &z)
    507426{
    508427  int     i;
    509   double  a;  //differential of x; a = (1/1+u^2)
     428  double  a;
    510429
    511430  a = x[0]; z[0] = atan(a); a = 1.0/(1.0+sqr(a));
    … …  
    514433}
    515434
    516 /*************************************
    517  *
    518  * Purpose:
    519  *     Get the derivatives of the arcctan(x)
    520  *   
    521  *    d(arcctan(u))/dx = -(1/1+u^2)*du/dx;
    522  *
    523  *  Input:
    524  *         x: u
    525  *  Output:
    526  *        z: the differential of x.
    527  *
    528  *         
    529  * ***********************************/
    530 void daarcctan(const tps_buf &x, tps_buf &z)
    531 {
    532   int     i;
    533   double  a;  //differential of x; a = -(1/1+u^2)
    534 
    535   a = x[0]; z[0] = atan(a); a = -1.0/(1.0+sqr(a));
    536   for (i = 1; i <= ss_dim; i++)
    537     z[i] = a*x[i];
    538 }
    539435
    540436void dafun_(const char *fun, const tps_buf &x, tps_buf &z)
    541437{
    542438  tps_buf  u;
    543  
     439
    544440  if (!strncmp(fun, "INV ", dafunlen))
    545441    dainv_(x, u);
    … …  
    554450  else if (!strncmp(fun, "COS ", dafunlen))
    555451    dacos(x, u);
    556   else if (!strncmp(fun, "TAN ", dafunlen))
    557     datan(x, u);
    558   else if (!strncmp(fun, "CTAN ", dafunlen))
    559     dactan(x, u);
    560452  else if (!strncmp(fun, "SINH ", dafunlen))
    561453    dasinh(x, u);
    562454  else if (!strncmp(fun, "COSH ", dafunlen))
    563455    dacosh(x, u);
    564   else if (!strncmp(fun, "ASIN", dafunlen))
    565     daarcsin(x, u);
    566   else if (!strncmp(fun, "ACOS", dafunlen))
    567     daarccos(x, u);
     456  else if (!strncmp(fun, "TAN ", dafunlen))
     457    datan(x, u);
    568458  else if (!strncmp(fun, "ATAN", dafunlen))
    569459    daarctan(x, u);
    570   else if (!strncmp(fun, "ACTAN", dafunlen))
    571     daarcctan(x, u);
    572460  else {
    573461    printf("dafun: illegal function %s\n", fun);

  • branches/tracy3-3.10.1b/tracy/tracy/src/tracy.cc

    r11 r23  
    9797template void radiate(ss_vect<tps> &, const double, const double, const tps[]);
    9898
    99 template void Drift(double, double, ss_vect<double> &);
    100 
    101 template void Drift(double, double, ss_vect<tps> &);
    102 
    10399template void Drift(double, ss_vect<double> &);
    104100
    … …  
    109105template void bend_fringe(double, ss_vect<tps> &);
    110106
    111 template void bend_fringe(double, double, ss_vect<double> &);
    112 
    113 template void bend_fringe(double, double, ss_vect<tps> &);
    114 
    115 template static void BendCurvature(double, double, ss_vect<double> &);
    116 
    117 template static void BendCurvature(double, double, ss_vect<tps> &);
    118 
    119 template static void EdgeFocus(double, double, double, ss_vect<double> &, bool);
    120 
    121 template static void EdgeFocus(double, double, double, ss_vect<tps> &,bool);
     107template static void EdgeFocus(double, double, double, ss_vect<double> &);
     108
     109template static void EdgeFocus(double, double, double, ss_vect<tps> &);
    122110
    123111template void quad_fringe(double, ss_vect<double> &);
    … …  
    128116
    129117template void Drift_Pass(CellType &, ss_vect<tps> &);
    130 
    131 
    132 template void dipole_kick(double, double, double, ss_vect<double> &);
    133 template void dipole_kick(double, double, double, ss_vect<tps> &);
    134 template void multipole_kick(int, double[], double, double, ss_vect<double> &);
    135 template void multipole_kick(int, double[], double, double, ss_vect<tps> &);
    136118
    137119template void thin_kick(int, double[], double, double, double,
Note: See TracChangeset for help on using the changeset viewer.