From c6733bb8d2512227facf353df313e25211b09747 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Sat, 1 Apr 2023 21:32:54 +0200 Subject: [PATCH] test bend fringe --- atintegrators/ExactSectorBendPass.c | 4 +- atintegrators/exactbendfringe.c | 68 ++++++++++++++++++++++++++++- 2 files changed, 69 insertions(+), 3 deletions(-) diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 9fba0006c..e97fdf6c6 100644 --- 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 index d2bb6da44..20de3c14e 100644 --- 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; + +}