Changeset 11 in TRACY3 for trunk/tracy/tracy/src/t2elem.cc

Timestamp:

Oct 21, 2013, 10:40:43 AM (11 years ago)

Author:

zhangj

Message:

 

File:

• trunk/tracy/tracy/src/t2elem.cc (modified) (23 diffs)

trunk/tracy/tracy/src/t2elem.cc

@@ -113 +113 @@
             cout << "        " << 1+x[delta_] << "   " << x[px_] << "   " << x[py_] << endl;
             printf("get_p_s: *** Speed of light exceeded!\n");
-          //  exit_(1);
+           // exit_(1);
             
         }
@@ -124 +124 @@
 
  Purpose:
- Drift pass
+ Exact Drift pass in the curvilinear coordinate for the dipole
+
+ See Eqn. (10.26) on page 307 in E. Forest's beam dynamics, a new attitude
+ 
+ H(x,y,cT,px,py,delta) =(1 + h_ref*x) * sqrt((1+delta)^ - px^2 - py^2) + delta                                  
+
+ Input:
+ L:  drift length
+ &x:  pointer to initial vector: x
+
+ Output:
+ none
+
+ Return:
+ none
+
+ Global variables:
+
+
+ Specific functions:
+ 
+ 
+ 
+ Comments:
+ 
+ Test version....rewritten by Jianfeng zhang @ LAL, 31/07/2013
+
+ ****************************************************************************/
+
+
+template<typename T>
+void Drift(double L,double h_bend, ss_vect<T> &x) {
+    T u;
+    ss_vect<T> x0;
+
+    x0=x;
+    
+    if (globval.H_exact) {
+      if(get_p_s(x) == 0)
+        return;
+
+      x[x_] = x0[x_]/cos(L*h_bend)/(1-x0[px_]/get_p_s(x)*tan(L*h_bend));
+      x[px_]= x0[px_]*cos(L*h_bend) + get_p_s(x)*sin(L*h_bend);
+      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));
+      x[py_]= x0[py_];
+      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));
+ 
+    
+    }
+}
+
+/****************************************************************************/
+/* Drift(double L, ss_vect<T> &x)
+
+ Purpose:
+ Drift pass in the Cartesian coordinate
 
  If H_exact = false, using approximation Hamiltonian (only valid for big ring):
@@ -616 +671 @@
     Correction for magnet gap (longitudinal fringe field)
 
-     irho h = 1/rho [1/m]
-     phi  dipole edge angle
-     gap  full gap between poles
+    Input: 
+     irho:         h = 1/rho [1/m]
+     phi:         dipole edge (entrance or exit) angle
+     gap:        full gap between poles
 
                              2
@@ -650 +706 @@
 }
 
+/*****************************************************************
+ * Exact map of the sector dipole
+ * 
+ * Forest's beam dynamics P.360 Eqn. (12.18)
+ * Also see Laurent's thesis: P219-220.
+ * 
+ * Written by Jianfeng Zhang @ LAL, 01/08/2013
+ * 
+ * Test version........................................
+ * 
+ * ****************************************************************/
+template<typename T>
+void dipole_kick(double L, double h_ref, double h_bend, ss_vect<T> &x){
+
+  ss_vect<T> x0;
+  T ps, u, dpxf;
+    
+  x0 = x;
+  ps = get_p_s(x0);
+  
+  u = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
+  
+  x[py_] = x0[py_];
+  
+  x[delta_] =x0[delta_];
+  
+  x[px_]= x0[px_]*cos(L*h_ref) + (ps -h_bend*(1/h_ref+x0[x_]))*sin(L*h_ref);
+  
+  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_]);
+  
+  
+  x[y_] += h_ref/h_bend*x0[py_]*L + x0[py_]/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
+  
+  x[x_] = 1/(h_bend*h_ref)*(h_ref*sqrt(sqr(1+x[delta_])-sqr(x[px_])-sqr(x[py_]))-dpxf - h_bend);
+  
+  x[ct_] += (1+x0[delta_])*h_ref/h_bend*L + (1+x0[delta_])/h_bend*(asin(x0[px_]/u)-asin(x[px_]/u));
+ 
+}
+
+/*****************************************************************
+ * pure multipole kick
+ * 
+ * Forest's beam dynamics P.360 Eqn. (12.18)
+ * Also see Laurent's thesis: P219-220.
+ * 
+ * Written by Jianfeng Zhang @ LAL, 01/08/2013
+ * 
+ * Test version........................................
+ * 
+ * ****************************************************************/
+template<typename T>
+void multipole_kick(int Order, double MB[], double L, double h_bend, ss_vect<T> &x){
+
+  int j=0;
+  ss_vect<T> x0;
+  
+  T ByoBrho = 0, BxoBrho = 0, ByoBrho1=0 ;
+  
+  if ((h_bend != 0.0) || ((1 <= Order) && (Order <= HOMmax))) {
+        x0 = x;
+        /* compute the total magnetic field Bx and By with Horner's rule */
+        ByoBrho = MB[HOMmax + Order]; // normalized By, By/(p0/e)
+        BxoBrho = MB[HOMmax - Order]; // normalized Bx, Bx/(p0/e)
+        for (j = Order - 1; j >= 1; j--) {
+            ByoBrho1 = x0[x_] * ByoBrho - x0[y_] * BxoBrho + MB[j + HOMmax];
+            BxoBrho = x0[y_] * ByoBrho + x0[x_] * BxoBrho + MB[HOMmax - j];
+            ByoBrho = ByoBrho1; 
+        }
+  }
+  
+  x[px_] -= L * (ByoBrho + h_bend*1); 
+        //vertical kick due to the magnets
+        x[py_] += L * BxoBrho;
+}
+
 /****************************************************************************/
 /* template<typename T>
@@ -657 +788 @@
         Refer to Tracy 2 manual.
         Calculate multipole kick. 
-  
+  , 
   (1)     The exact Hamiltonian is
 
@@ -737 +868 @@
  MB     array of an and bn, magnetic field components
  L      multipole length
- h_bend   B_dipole/Brho is the dipole field normalized by the reference momentum 
+ h_bend   B_dipole/Brho is the dipole field normalized by the reference momentum; one of the example is the combined dipole
  h_ref   1/rho [m^-1] is the curvature of the reference orbit in the dipoles which is used in curvilinear coordinate,
  x      initial coordinates vector
@@ -797 +928 @@
               return;
             
-            u = L * h_ref * x0[x_] /get_p_s(x0);
-            x[x_] += u * x0[px_];
-            x[y_] += u * x0[py_];
-            x[px_]-= L * (-h_ref * get_p_s(x0) + h_bend + h_ref * h_bend
-                     * x0[x_] + ByoBrho);
-            x[ct_] += u * (1.0+x0[delta_]);
+            //dipole kick map from the exact Hamiltonian; not symplectic
+//          u = L * h_ref * x0[x_] /get_p_s(x0);
+//             x[x_] += u * x0[px_];
+//             x[y_] += u * x0[py_];
+//          x[px_]-= L * (-h_ref * get_p_s(x0) + h_bend + h_ref * h_bend
+//                   * x0[x_] + ByoBrho);
+//             x[ct_] += u * (1.0+x0[delta_]);
+             
+              //kick map due to the dipole vector potential
+             // x[px_] -= L*(h_bend + h_ref*h_bend*x0[x_]+ByoBrho);
  
 //          sqrtpx = sqrt(sqr(1+x0[delta_])-sqr(x0[py_]));
@@ -835 +970 @@
             x[ct_] += L * h_ref * x0[x_];
   
-            //  second order of the h_ref; has problem when do tracking...
-           u = L * h_ref * x0[x_] /(1.0+x[delta_]);
+            //  second order, cubic term of H, H3; has problem when do tracking...works for THOMX!!!!!
+            // not symplectic????
+         if(0){ 
+            u = L * h_ref * x0[x_] /(1.0+x0[delta_]);
             x[x_] += u * x0[px_];
             x[y_] += u * x0[py_];
             x[ct_] += u*(sqr(x0[px_])+sqr(x0[py_]))/(2.0*(1.0+x0[delta_]));
             
-            }
+           // cout << "test......2.....  u = " << u<< endl;
+          }     
+      }
           
         }//from other straight magnets 
         else 
-            x[px_] -= L * (ByoBrho + h_bend);
+            x[px_] -= L * (ByoBrho + h_bend*1); //h_bend should be trigger on for combined dipoles
 
         //vertical kick due to the magnets
@@ -872 +1011 @@
  Compute edge focusing for a dipole
  There is no radiation coming from the edge
- The standard formula used is : (Where these formula comes from?????)
+ 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)
 
        px = px0 +  irho tan(phi) *x0
@@ -913 +1052 @@
 
  April 2011, No energy dependence for H-plane according to SSC-141 note (E. Forest)
- Agreement better with MAD8 (LNLS), SOLEIL small effect (.1)
+ Agreement better with MAD8 (LNLS), SOLEIL small effect (.1). Modified by Laurent Nadolski @ soleil
+ Comments: Jianfeng Zhang @ LAL, 05/06/2013
+ This model only works for Soleil lattice, but not works for ThomX and SuperB DR lattice.
  
- Add the contribution of x' to the x[py_] for small ring like Thom-X.
  18/06/2012  Jianfeng Zhang @ LAL
- 
+ Add the contribution of x' to the x[py_] for small ring like Thom-X, using Alexandre Loulergue's correction.
+ This models works for ThomX, Soleil, and SuperB DR lattices; and 
+ close to the Forest's model and K. Brown's models.
  
  ****************************************************************************/
@@ -926 +1068 @@
     if (true) {
       
-        //cout << "With Alex's  correction: " << endl;
+       // cout << "Dipole edge effect WITH Alex's  correction: " << endl;
       
       //with the contribution from the entrance angle of the particle at the edge
@@ -939 +1081 @@
 // compared the fringe field model used in Tracy
 // results:  K. Brown's model is wrong for small ring like Thom-X!!!!
-// The original model of tracy is correct !!!!!!!!!!
+// The original model of tracy is correct for tracking!!!!!!!!!!
   /*
   double psi;
@@ -959 +1101 @@
     */
 
-    //add the contribution of the entrance momentum of the particle, from Alex Loulergue. This model is NOT
-    //correct, because it's not symplectic!!!!!!!!
+    //add the contribution of the entrance momentum of the particle, from Alex Loulergue. 
     if(Entrance){
-       x[px_] += irho * tan(dtor(phi)) * x[x_]/ (1.0 + x[delta_]*0);
-       x[py_] -= irho * tan(dtor(phi) + x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)) * x[y_]
-                / (1.0 + x[delta_]);
+       x[px_] += irho * tan(dtor(phi) +x[py_]/(1+x[delta_]*1)) * x[x_]/ (1.0 + x[delta_]*0);
+       x[py_] -= irho * tan(dtor(phi) + x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)/(1+x[delta_]*0)) * x[y_]
+                / (1.0 + x[delta_]*0);
     }else{
-       x[px_] += irho * tan(dtor(phi)) * x[x_]/ (1.0 + x[delta_]*0);
-       x[py_] -= irho * tan(dtor(phi) - x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)) * x[y_]
-                / (1.0 + x[delta_]);
+       x[px_] += irho * tan(dtor(phi) - x[py_]/(1+x[delta_]*1)) * x[x_]/ (1.0 + x[delta_]*0);
+       x[py_] -= irho * tan(dtor(phi) - x[px_]/(1+x[delta_]*1) - get_psi(irho, phi, gap)/(1+x[delta_]*0)) * x[y_]
+                / (1.0 + x[delta_]*0);
     }
       
     } else {//original one
 
-  //cout << "Without Alex's  correction: " << endl;
+   //cout << "Dipole edge effect WITHOUT Alex's  correction: " << endl;
         x[px_] += irho * tan(dtor(phi)) * x[x_];
-        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_];
+        x[py_] -= irho * tan(dtor(phi) - get_psi(irho, phi, gap)) * x[y_]/ (1.0 + x[delta_]*1);
     }
 }
@@ -981 +1122 @@
 
 /******************************************************************
- 
+ template<typename T>
+static void BendCurvature(double irho, double H, ss_vect<T> &x)
+
+The leanding order term T211, T233, and T413 on page 117 & page 118
+of the dipole edge effect in SLAC-75, using the symplectic
+form in 
  E. Forest: "Fringe field in MAD, Part II: Bend curvature in MAD-X for the 
              mudule PTC". 
              P.10 Eq. (5) and (42).
  
+ For an rectangular dipole, the bending curvature term is zero, since 1/R = 0.
+ 
+ irho:   curvature of the central orbit inside the dipole
+ H:      curvature of the entrance/exit pole face of the dipole
+ 
+ Comments:
+ Written by Jianfeng Zhang @ LAL, 01/10/2012.
 ******************************************************************/
 template<typename T>
-static void EdgeFocus(double irho, ss_vect<T> &x) {
+static void BendCurvature(double irho,  double H, ss_vect<T> &x) {
+  
+  if(trace)
+  cout << "Forest's bend curvature...."<<endl;
   
   T pm2,u, coeff1, coeff2,coeff3;
   
-   ss_vect<T> x1;
+   ss_vect<T> x0;
    
-   x1 = x;
+   x0 = x;
   
-   pm2 = sqr(1+x1[delta_]) - sqr(x1[px_]); 
-   u = irho*irho/2/pm2;
+   pm2 = sqr(1+x0[delta_]) - sqr(x0[px_]); 
+   u = irho*H/2/pm2;
    
-   coeff1 = u * 2*x1[px_]*(1+x1[delta_])/pm2; 
-   coeff2 = u * (1+x1[delta_]);
-   coeff3 = u *(-1) * (sqr(1+x1[delta_]) + sqr(x1[px_]))/pm2; 
+   coeff1 = u * 2*x0[px_]*(1+x0[delta_])/pm2; //
+   coeff2 = u * (1+x0[delta_]); //
+   coeff3 = u * (pm2 - 2*(1+x0[delta_]))/pm2; 
    
  
-   x[x_] = x1[x_]/(1-coeff1*x1[y_]*x1[y_]);
-   x[px_] -= coeff2*x1[y_]*x1[y_] - irho*irho/2*sqr(x1[x_]);
-   x[py_] -= 2*coeff2*x[x_]*x1[y_];
-   x[ct_] -= coeff3*x[x_]*x1[y_]*x1[y_];
+   x[x_]  += x0[x_]/(1-coeff1*x0[y_]*x0[y_]);
+   x[px_] -= coeff2*x0[y_]*x0[y_] - irho*H/2*sqr(x0[x_]);
+   x[py_] -= 2*coeff2*x[x_]*x0[y_];
+   x[ct_] -= coeff3*x[x_]*x0[y_]*x0[y_];
 }   
    
@@ -1313 +1469 @@
         //fringe field and edge focusing of dipoles
        if (M->Pirho != 0.0 && M->dipEdge_effect1 == 1){
-         if (!globval.H_exact) { //big ring */
-           //            EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x,true);
-           //} else {//small and big rings
-          if(M->PTx1!=0){
-            EdgeFocus(M->Pirho, x);
-           //p_rot(M->PTx1,  x); 
-             //rotate from cartesian cooridnate to curlinear curved beam coordinate; 
+         if (!globval.H_exact) { //big ring */ modified linear term T43 on page 117 & 118 in SLAC-75.
+                       EdgeFocus(M->Pirho, M->PTx1, M->Pgap, x,true);
+           } else {//small and big rings; Forest's model
+          if(M->PH1!=0){
+            // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0.
+           BendCurvature(M->Pirho, M->PH1, x); 
+          }
+           p_rot(M->PTx1,  x);  //rotate from cartesian cooridnate to curlinear curved beam coordinate; 
              // since the map of a sector dipole is used in Tracy 3 
-          }
-           p_rot(M->PTx1,  x); 
            bend_fringe(M->Pirho, M->Pgap, x); //fringe field
          }
@@ -1366 +1521 @@
 
     case Meth_Second:
-        cout << "Mpole_Pass: Meth_Second not supported" << endl;
-        exit_(0);
-        break;
-
+      
+      //  cout << "Mpole_Pass: Meth_Second not supported" << endl;
+       // exit_(0);
+       // break;
+
+    //specially for the test of the dipole map,see Forest's beam dynamics, p361. eqn.(12.19)
+ 
+    dL1 = 0.5*dL;
+    dL2 = dL;
+           
+            for (seg = 1; seg <= M->PN; seg++) {
+              dipole_kick(dL1,M->Pirho,h_ref,x);
+              multipole_kick(M->Porder, M->PB, dkL1, M->Pirho, x);
+              dipole_kick(dL1,M->Pirho,h_ref,x);
+            }
+
+            
+       
     case Meth_Fourth:  //forth order symplectic integrator
         if (M->Pthick == thick) {
@@ -1386 +1555 @@
                 }
 
-                Drift(dL1, x);
+                 Drift(dL1, x);
                 thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
                 Drift(dL2, x);
-                thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
+                thin_kick(M->Porder, M->PB, dkL2, M->Pirho, h_ref, x);
 
                 if (globval.emittance && (!globval.Cavity_on) && (M->Pirho
@@ -1396 +1565 @@
                     dI4 += 4.0 * is_tps<tps>::get_dI4(x);
                 }
-
+                
+                
                 Drift(dL2, x);
-                thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
-                Drift(dL1, x);
+                thin_kick(M->Porder, M->PB, dkL1, M->Pirho, h_ref, x);
+                Drift(dL1, x);
 
                 if (globval.emittance && (!globval.Cavity_on) && (M->Pirho
@@ -1428 +1598 @@
       /* dipole fringe fields */
              if (M->Pirho != 0.0 && M->dipEdge_effect2 == 1){
-            if (!globval.H_exact) { //big ring
-              //                 EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x,false);
-              // }else{
-              
-              
-             bend_fringe(-M->Pirho, M->Pgap,x); //fringe field of the dipole
-              p_rot(M->PTx2, x); 
-             if(M->PTx2!=0){
-              //p_rot(M->PTx2, x); //rotate back to the cartesian cooridinate
-            EdgeFocus(M->Pirho, x);
+            if (!globval.H_exact) { //big ring, linear model correction
+                               EdgeFocus(M->Pirho, M->PTx2, M->Pgap, x,false);
+               }else{
+             bend_fringe(-M->Pirho, M->Pgap,x); //Forest's fringe field of the dipole
+              p_rot(M->PTx2, x); //rotate back to the cartesian cooridinate
+             if(M->PH2!=0){
+            // curvature of the magnet pole face; for a sector/wedge/rectangular dipole, this term is 0.
+             BendCurvature(M->Pirho, M->PH2, x);
              }
-        
-
          }
      }
@@ -2762 +2928 @@
      corresponding to a half magnet w/ an entrance
      angle and no exit angle.
-     The linear frindge field is taken into account                
+     The linear fringe field is taken into account                
      **************************************************************/
 void Mpole_Setmatrix(int Fnum1, int Knum1, double K) {
@@ -3259 +3425 @@
         printf("Elem_GetPos: there are no kids in family %d (%s)\n", Fnum1,
                 ElemFam[Fnum1 - 1].ElemF.PName);
-        exit_(0);
+        exit_(1);
     }
 
@@ -3290 +3456 @@
     d_2 = 1e0 - 2e0 * d_1;
 
+    
+    
+    
+    
     // classical radiation
     //  C_gamma = 8.846056192e-05;
@@ -3682 +3852 @@
     M->PTx1 = 0.0; /* Entrance angle */
     M->PTx2 = 0.0; /* Exit angle */
+    M->PH1 = 0.0; /* Entrance pole face curvature*/
+    M->PH2 = 0.0; /* Exit pole face curvature */
     M->Pgap = 0.0; /* Gap for fringe field ??? */