test bend fringe

atcollab · Apr 1, 2023 · c6733bb · c6733bb
1 parent 8083146
commit c6733bb
Show file tree

Hide file tree

Showing 2 changed files with 69 additions and 3 deletions.
diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c
@@ -89,7 +89,7 @@ static void ExactSectorBend(
 
       /*  integrator  */
       Yrot(r6, entrance_angle);
-      // bend_fringe(r6, irho, gK);
+      bend_fringe(r6, irho, gK);
       if (do_fringe)
         multipole_fringe(r6, le, A, B, max_order, 1.0, 1);
       bend_edge(r6, irho, -entrance_angle);
@@ -112,7 +112,7 @@ static void ExactSectorBend(
       bend_edge(r6, irho, -exit_angle);
       if (do_fringe)
         multipole_fringe(r6, le, A, B, max_order, -1.0, 1);
-      // bend_fringe(r6, -irho, gK);
+      bend_fringe(r6, -irho, gK);
       Yrot(r6, exit_angle);
 
       /* Check physical apertures at the exit of the magnet */

diff --git a/atintegrators/exactbendfringe.c b/atintegrators/exactbendfringe.c
@@ -30,7 +30,7 @@ static void Yrot(double *r6, double phi)
     }
 }
 
-void bend_fringe(double *r6, double irho, double gK)
+void bend_fri(double *r6, double irho, double gK)
 {
     /* Forest 13.13, bend fringe in the hard-edge limit */
 
@@ -116,3 +116,69 @@ static void bend_edge(double *r6, double rhoinv, double theta)
         r6[ct_] += dct;
     }
 }
+
+#define Power pow
+#define ArcTan atan
+#define pow2(u) ((u)*(u))
+
+double get_pz(double * x)
+{
+    return sqrt(pow2(1 + x[delta_]) - pow2(x[px_]) - pow2(x[py_]));
+}
+
+void bend_fringe(double * x, double irho, double gK)
+{
+
+    double dpx, dpy, dd, b0, px, py, pz, g, K, d, phi, xp, yp, yf, xf, lf, pyf;
+
+    b0 = irho;
+
+    /* gK always multiplied together so put everything in g and set K to one */
+
+    K = 1.0;
+    g = gK;
+
+    pz = get_pz(x);
+    px = x[px_];
+    py = x[py_];
+    d  = x[delta_];
+    xp = px / pz;
+    yp = py / pz;
+
+    phi = -b0 * tan( b0 * g * K * (1 +    pow2(xp)*(2 + pow2(yp)))*pz - atan(xp / (1 + pow2(yp))));
+
+    /* these are the partial derivatives of phi with respect to px, py and delta
+       total horror from Mathematica. This could benefit from some mini-TPSA */
+
+    dpx = -((b0*(Power(px,2)*Power(pz,4)*(Power(py,2) - Power(pz,2)) - Power(pz,6)*(Power(py,2) + Power(pz,2)) +
+                 b0*g*K*px*(Power(pz,2)*Power(Power(py,2) + Power(pz,2),2)*(2*Power(py,2) + 3*Power(pz,2)) + Power(px,4)*(3*Power(py,2)*Power(pz,2) + 2*Power(pz,4)) +
+                            Power(px,2)*(3*Power(py,6) + 8*Power(py,4)*Power(pz,2) + 9*Power(py,2)*Power(pz,4) + 5*Power(pz,6))))*
+             Power(Sec((b0*g*K*(Power(pz,4) + Power(px,2)*(Power(py,2) + 2*Power(pz,2))))/Power(pz,3) - ArcTan((px*pz)/(Power(py,2) + Power(pz,2)))),2))/
+            (Power(pz,5)*(Power(py,4) + Power(px,2)*Power(pz,2) + 2*Power(py,2)*Power(pz,2) + Power(pz,4))));
+
+
+    dpy = -((b0*py*(px*Power(pz,4)*(Power(py,2) + Power(pz,2)) + b0*g*K*(-(Power(pz,4)*Power(Power(py,2) + Power(pz,2),2)) +
+                                                                         Power(px,4)*(3*Power(py,2)*Power(pz,2) + 4*Power(pz,4)) +
+                                                                         Power(px,2)*(3*Power(py,6) + 10*Power(py,4)*Power(pz,2) + 11*Power(py,2)*Power(pz,4) + 3*Power(pz,6))))*
+             Power(Sec((b0*g*K*(Power(pz,4) + Power(px,2)*(Power(py,2) + 2*Power(pz,2))))/Power(pz,3) - ArcTan((px*pz)/(Power(py,2) + Power(pz,2)))),2))/
+            (Power(pz,5)*(Power(py,4) + Power(px,2)*Power(pz,2) + 2*Power(py,2)*Power(pz,2) + Power(pz,4))));
+
+    dd = (b0*(1 + d)*(px*Power(pz,4)*(Power(py,2) - Power(pz,2)) + b0*g*K*
+                      (-(Power(pz,4)*Power(Power(py,2) + Power(pz,2),2)) + Power(px,4)*(3*Power(py,2)*Power(pz,2) + 2*Power(pz,4)) +
+                       Power(px,2)*(3*Power(py,6) + 8*Power(py,4)*Power(pz,2) + 7*Power(py,2)*Power(pz,4) + Power(pz,6))))*
+          Power(Sec((b0*g*K*(Power(pz,4) + Power(px,2)*(Power(py,2) + 2*Power(pz,2))))/Power(pz,3) - ArcTan((px*pz)/(Power(py,2) + Power(pz,2)))),2))/
+        (Power(pz,5)*(Power(py,4) + Power(px,2)*Power(pz,2) + 2*Power(py,2)*Power(pz,2) + Power(pz,4)));
+
+    /* solve quadratic equation in yf (Forest fringe_part_I.pdf) */
+
+    yf = (2 * x[y_]) / (1 + sqrt(1 - 2 * dpy * x[y_]));
+    xf = x[x_] + 0.5 * dpx * pow2(yf);
+    lf = x[ct_] - 0.5 * dd * pow2(yf);
+    pyf = py - phi * yf;
+
+    x[y_]  = yf;
+    x[x_]  = xf;
+    x[py_] = pyf;
+    x[ct_] = lf;
+
+}