Quad curvature correction in bending magnets

atintegrators/BndMPoleSymplectic4Pass.c

@@ -2,7 +2,8 @@
 #include "atelem.c"
 #include "atlalib.c"
 #include "atphyslib.c"
-#include "driftkick.c"		/* fastdrift and bndthinkick */
+#include "drift_fast.h"
+#include "kick_h_k0h_k1h_kn.h"
 #include "quadfringe.c"		/* QuadFringePassP, QuadFringePassN */
 
 struct elem
@@ -55,12 +56,12 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl
     double K2 = SL*KICK2;
     bool useLinFrEleEntrance = (fringeIntM0 != NULL && fringeIntP0 != NULL  && FringeQuadEntrance==2);
     bool useLinFrEleExit = (fringeIntM0 != NULL && fringeIntP0 != NULL  && FringeQuadExit==2);
-    double B0 = B[0];
-    double A0 = A[0];
+    double B0 = 0.0;
+    double A0 = 0.0;
 
-    if (KickAngle) {   /* Convert corrector component to polynomial coefficients */
-        B[0] -= sin(KickAngle[0])/le;
-        A[0] += sin(KickAngle[1])/le;
+    if (KickAngle) { /* Convert corrector component to polynomial coefficients */
+        B0 = -sin(KickAngle[0]) / le;
+        A0 = sin(KickAngle[1]) / le;
     }
 
     #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \
@@ -95,13 +96,13 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl
             }
             /* integrator */
             for (m=0; m < num_int_steps; m++) { /* Loop over slices */
-                fastdrift(r6, NormL1);
-                bndthinkick(r6, A, B, K1, irho, max_order);
-                fastdrift(r6, NormL2);
-                bndthinkick(r6, A, B, K2, irho, max_order);
-                fastdrift(r6, NormL2);
-                bndthinkick(r6, A, B, K1, irho, max_order);
-                fastdrift(r6, NormL1);
+                drift(r6, NormL1);
+                kick(r6, A0, B0, A, B, max_order, K1, irho);
+                drift(r6, NormL2);
+                kick(r6, A0, B0, A, B, max_order, K2, irho);
+                drift(r6, NormL2);
+                kick(r6, A0, B0, A, B, max_order, K1, irho);
+                drift(r6, NormL1);
             }
             /* quadrupole gradient fringe */
             if (FringeQuadExit && B[1]!=0) {
@@ -122,10 +123,6 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl
             if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling);
         }
     }
-    if (KickAngle) {  /* Remove corrector component in polynomial coefficients */
-        B[0] = B0;
-        A[0] = A0;
-    }
 }
 
 #if defined(MATLAB_MEX_FILE) || defined(PYAT)

atintegrators/BndMPoleSymplectic4RadPass.c

@@ -1,8 +1,8 @@
 #include "atelem.c"
 #include "atlalib.c"
-#include "diff_bend_fringe.c"
-#include "diff_bnd_kick.c"
 #include "diff_drift.c"
+#include "diff_h_k0h_k1h_kn.c"
+#include "diff_bend_fringe.c"
 #include "quadfringe.c"		/* QuadFringePassP, QuadFringePassN */
 
 struct elem
@@ -57,15 +57,16 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do
     double K2 = SL*KICK2;
     bool useLinFrEleEntrance = (fringeIntM0 != NULL && fringeIntP0 != NULL  && FringeQuadEntrance==2);
     bool useLinFrEleExit = (fringeIntM0 != NULL && fringeIntP0 != NULL  && FringeQuadExit==2);
-    double B0 = B[0];
-    double A0 = A[0];
     double rad_const = RAD_CONST*pow(gamma, 3);
     double diff_const = DIF_CONST*pow(gamma, 5);
+    double B0 = 0.0;
+    double A0 = 0.0;
 
-    if (KickAngle) {   /* Convert corrector component to polynomial coefficients */
-        B[0] -= sin(KickAngle[0])/le;
-        A[0] += sin(KickAngle[1])/le;
+    if (KickAngle) { /* Convert corrector component to polynomial coefficients */
+        B0 = -sin(KickAngle[0]) / le;
+        A0 = sin(KickAngle[1]) / le;
     }
+
     #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \
     shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,bdiff,\
     irho,gap,A,B,L1,L2,K1,K2,max_order,num_int_steps,rad_const, diff_const,scaling,\
@@ -95,11 +96,11 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do
             /* integrator */
             for (m=0; m < num_int_steps; m++) { /* Loop over slices */
                 diff_drift(r6, L1, bdiff);
-                diff_bnd_kick(r6, A, B, max_order, K1, irho, rad_const, diff_const, bdiff);
+                diff_kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, bdiff);
                 diff_drift(r6, L2, bdiff);
-                diff_bnd_kick(r6, A, B, max_order, K2, irho, rad_const, diff_const, bdiff);
+                diff_kick(r6, A0, B0, A, B, max_order, K2, irho, rad_const, diff_const, bdiff);
                 diff_drift(r6, L2, bdiff);
-                diff_bnd_kick(r6, A, B, max_order, K1, irho, rad_const, diff_const, bdiff);
+                diff_kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, bdiff);
                 diff_drift(r6, L1, bdiff);
             }
             /* quadrupole gradient fringe */
@@ -121,10 +122,6 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do
             if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling);
         }
     }
-    if (KickAngle) {  /* Remove corrector component in polynomial coefficients */
-        B[0] = B0;
-        A[0] = A0;
-    }
 }
 
 #if defined(MATLAB_MEX_FILE) || defined(PYAT)

atintegrators/ExactMultipoleRadPass.c

@@ -139,11 +139,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem,
     }
 
     /* Check energy */
-    Energy = atEnergy(Param->energy, Energy);
-    if (Energy == 0) {
-        atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
-        check_error();
-    }
+    Energy = atEnergy(Param->energy, Energy); check_error();
 
     Elem = (struct elem *)atMalloc(sizeof(struct elem));
     Elem->Length = Length;

atintegrators/ExactRectBendPass.c

@@ -22,7 +22,8 @@ struct elem
     int FringeBendExit;
     int FringeQuadEntrance;
     int FringeQuadExit;
-    double gK;
+    double gK_entrance;
+    double gK_exit;
     double x0ref;
     double refdz;
     double *R1;
@@ -40,34 +41,43 @@ static void ExactRectangularBend(double *r, double le, double bending_angle,
         double entrance_angle, double exit_angle,
         int FringeBendEntrance, int FringeBendExit,
         int FringeQuadEntrance, int FringeQuadExit,
-        double gK, double x0ref, double refdz,
+        double gK_entrance, double gK_exit,
+        double x0ref, double refdz,
         double *T1, double *T2,
         double *R1, double *R2,
         double *RApertures, double *EApertures,
         double *KickAngle, double scaling, int num_particles)
 {
     double irho = bending_angle / le;
     double phi2 = 0.5 * bending_angle;
-    double LR = phi2 < 1.e-10 ? le : le *sin(phi2) / phi2;
+    double phi_entrance = phi2-entrance_angle;
+    double phi_exit = phi2-exit_angle;
+    double LR = fabs(phi2) < 1.e-10 ? le : le *sin(phi2) / phi2;
     double SL = LR/num_int_steps;
     double L1 = SL*DRIFT1;
     double L2 = SL*DRIFT2;
     double K1 = SL*KICK1;
     double K2 = SL*KICK2;
     double B0 = B[0];
     double A0 = A[0];
+    double k1_entrance_angle = 0.0;
+    double k1_exit_angle = 0.0;
 
-    if (KickAngle) {   /* Convert corrector component to polynomial coefficients */
+    if (KickAngle) { /* Convert corrector component to polynomial coefficients */
         B[0] -= sin(KickAngle[0])/le;
         A[0] += sin(KickAngle[1])/le;
     }
+    if (max_order >= 1) {
+        k1_entrance_angle = B[1] * entrance_angle;
+        k1_exit_angle = B[1] * exit_angle;
+    }
 
     #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \
     shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,\
-    irho,gK,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\
-    entrance_angle,exit_angle,x0ref,refdz,\
+    irho,gK_entrance,gK_exit,A0,B0,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\
+    k1_entrance_angle,k1_exit_angle,entrance_angle,exit_angle,x0ref,refdz,\
     FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\
-    LR,le,phi2)
+    LR,le,phi_entrance,phi_exit)
     for (int c = 0; c<num_particles; c++) { /* Loop over particles */
         double *r6 = r + 6*c;
         if (!atIsNaN(r6[0])) {
@@ -85,14 +95,17 @@ static void ExactRectangularBend(double *r, double le, double bending_angle,
             if (RApertures) checkiflostRectangularAp(r6,RApertures);
             if (EApertures) checkiflostEllipticalAp(r6,EApertures);
 
-            /* edge focus */
+            /* Entry face */
+            r6[0] += x0ref;
             if (FringeBendEntrance)
-                bend_fringe(r6, irho, gK);
+                bend_fringe(r6, irho, gK_entrance);
             if (FringeQuadEntrance)
                 multipole_fringe(r6, le, A, B, max_order, 1.0, 1);
-            bend_edge(r6, irho, phi2-entrance_angle);
+            if (phi_entrance != 0.0) {
+                if (k1_entrance_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_entrance_angle);
+                bend_edge(r6, irho, phi_entrance);
+            }
 
-            r6[0] += x0ref;
             if (num_int_steps == 0) {
                 exact_straight_bend(r6, irho, LR);
             }
@@ -107,17 +120,20 @@ static void ExactRectangularBend(double *r, double le, double bending_angle,
                     exact_straight_bend(r6, irho, L1);
                 }
             }
-            r6[0] -= x0ref;
 
             /* Convert absolute path length to path lengthening */
             r6[5] -= (le+refdz);
 
-            /* edge focus */
-            bend_edge(r6, irho, phi2-exit_angle);
+            /* Exit face */
+            if (phi_exit != 0.0) {
+                bend_edge(r6, irho, phi_exit);
+                if (k1_exit_angle != 0.0 && FringeQuadExit) quad_wedge(r6, -k1_exit_angle);
+            }
             if (FringeQuadExit)
                 multipole_fringe(r6, le, A, B, max_order, -1.0, 1);
             if (FringeBendExit)
-                bend_fringe(r6, -irho, gK);
+                bend_fringe(r6, -irho, gK_exit);
+            r6[0] -= x0ref;
 
             /* Check physical apertures at the exit of the magnet */
             if (RApertures) checkiflostRectangularAp(r6, RApertures);
@@ -158,7 +174,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
         int FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error();
         int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error();
         int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error();
-        double gK=atGetOptionalDouble(ElemData,"gK", 0.0); check_error();
+        double FullGap=atGetOptionalDouble(ElemData,"FullGap",0.0); check_error();
+        double FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1",0.0); check_error();
+        double FringeInt2=atGetOptionalDouble(ElemData,"FringeInt2",0.0); check_error();
         double x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error();
         double refdz=atGetOptionalDouble(ElemData,"RefDZ", 0.0); check_error();
         double *R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error();
@@ -192,7 +210,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
         Elem->FringeBendExit=FringeBendExit;
         Elem->FringeQuadEntrance=FringeQuadEntrance;
         Elem->FringeQuadExit=FringeQuadExit;
-        Elem->gK=gK;
+        Elem->gK_entrance=FullGap*FringeInt1;
+        Elem->gK_exit=FullGap*FringeInt2;
         Elem->x0ref=x0ref;
         Elem->refdz=refdz;
         Elem->R1=R1;
@@ -207,7 +226,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
             Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle,
             Elem->FringeBendEntrance,Elem->FringeBendExit,
             Elem->FringeQuadEntrance, Elem->FringeQuadExit,
-            Elem->gK,Elem->x0ref,Elem->refdz,
+            Elem->gK_entrance, Elem->gK_exit,
+            Elem->x0ref,Elem->refdz,
             Elem->T1, Elem->T2, Elem->R1, Elem->R2,
             Elem->RApertures, Elem->EApertures,
             Elem->KickAngle, Elem->Scaling, num_particles);
@@ -241,7 +261,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         int FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error();
         int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error();
         int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error();
-        double gK=atGetOptionalDouble(ElemData,"gK", 0.0); check_error();
+        double FullGap=atGetOptionalDouble(ElemData,"FullGap",0.0); check_error();
+        double FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1",0.0); check_error();
+        double FringeInt2=atGetOptionalDouble(ElemData,"FringeInt2",0.0); check_error();
         double x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error();
         double refdz=atGetOptionalDouble(ElemData,"RefDZ", 0.0); check_error();
         double *R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error();
@@ -263,11 +285,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         /* ALLOCATE memory for the output array of the same size as the input  */
         plhs[0] = mxDuplicateArray(prhs[1]);
         r_in = mxGetDoubles(plhs[0]);
+
         ExactRectangularBend(r_in, Length, BendingAngle, PolynomA, PolynomB,
             MaxOrder, NumIntSteps, EntranceAngle, ExitAngle,
             FringeBendEntrance, FringeBendExit,
             FringeQuadEntrance, FringeQuadExit,
-            gK, x0ref, refdz,
+            FullGap*FringeInt1, FullGap*FringeInt2,
+            x0ref, refdz,
             T1, T2, R1, R2, RApertures, EApertures,
             KickAngle, Scaling, num_particles);
     } else if (nrhs == 0) {
@@ -284,12 +308,14 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
         mxSetCell(plhs[0], i0++, mxCreateString("ExitAngle"));
         if (nlhs>1) {    /* list of optional fields */
             int i1 = 0;
-            plhs[1] = mxCreateCellMatrix(15, 1);
+            plhs[1] = mxCreateCellMatrix(17, 1);
             mxSetCell(plhs[1], i1++, mxCreateString("FringeBendEntrance"));
             mxSetCell(plhs[1], i1++, mxCreateString("FringeBendExit"));
             mxSetCell(plhs[1], i1++, mxCreateString("FringeQuadEntrance"));
             mxSetCell(plhs[1], i1++, mxCreateString("FringeQuadExit"));
-            mxSetCell(plhs[1], i1++, mxCreateString("gK"));
+            mxSetCell(plhs[1], i1++, mxCreateString("FullGap"));
+            mxSetCell(plhs[1], i1++, mxCreateString("FringeInt1"));
+            mxSetCell(plhs[1], i1++, mxCreateString("FringeInt2"));
             mxSetCell(plhs[1], i1++, mxCreateString("X0ref"));
             mxSetCell(plhs[1], i1++, mxCreateString("RefDZ"));
             mxSetCell(plhs[1], i1++, mxCreateString("T1"));