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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.