From 89a33ca2b73d15e97b6ca739ef1f6fa368a86f97 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 10:04:28 +0100 Subject: [PATCH 01/34] quad curvature correction --- atintegrators/BndMPoleSymplectic4Pass.c | 30 ++++----- atintegrators/BndMPoleSymplectic4RadPass.c | 24 +++---- atintegrators/ExactSectorBendPass.c | 21 +++--- atintegrators/ExactSectorBendRadPass.c | 21 +++--- atintegrators/a_fastdrift.c | 13 ++++ atintegrators/a_h_k0h_k1h_kn.c | 42 ++++++++++++ atintegrators/b_drift_kick_drift_expanded.c | 2 + .../diff_drift_kick_drift_expanded.c | 2 + .../{diff_bnd_kick.c => diff_h_k0h_k1h_kn.c} | 67 +++++-------------- 9 files changed, 118 insertions(+), 104 deletions(-) create mode 100644 atintegrators/a_fastdrift.c create mode 100644 atintegrators/a_h_k0h_k1h_kn.c create mode 100644 atintegrators/b_drift_kick_drift_expanded.c create mode 100644 atintegrators/diff_drift_kick_drift_expanded.c rename atintegrators/{diff_bnd_kick.c => diff_h_k0h_k1h_kn.c} (68%) diff --git a/atintegrators/BndMPoleSymplectic4Pass.c b/atintegrators/BndMPoleSymplectic4Pass.c index 1aef9d8e2..bb696e6d2 100644 --- a/atintegrators/BndMPoleSymplectic4Pass.c +++ b/atintegrators/BndMPoleSymplectic4Pass.c @@ -2,7 +2,7 @@ #include "atelem.c" #include "atlalib.c" #include "atphyslib.c" -#include "driftkick.c" /* fastdrift and bndthinkick */ +#include "b_drift_kick_drift_expanded.c" /* kick */ #include "quadfringe.c" /* QuadFringePassP, QuadFringePassN */ struct elem @@ -55,12 +55,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 +95,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, K1, irho, max_order); + drift(r6, NormL2); + kick(r6, A0, B0, A, B, K2, irho, max_order); + drift(r6, NormL2); + kick(r6, A0, B0, A, B, K1, irho, max_order); + drift(r6, NormL1); } /* quadrupole gradient fringe */ if (FringeQuadExit && B[1]!=0) { @@ -122,10 +122,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) diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index 897906ff9..4d361dc9e 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -1,8 +1,7 @@ #include "atelem.c" #include "atlalib.c" #include "diff_bend_fringe.c" -#include "diff_bnd_kick.c" -#include "diff_drift.c" +#include "diff_drift_kick_drift_expanded.c" #include "quadfringe.c" /* QuadFringePassP, QuadFringePassN */ struct elem @@ -57,15 +56,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 +95,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 +121,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) diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 44b910834..3dfd642ae 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -1,7 +1,7 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "driftkick.c" /* strthinkick.c */ +#include "a_strthinkick.c" /* kick */ #include "exactbend.c" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -50,12 +50,12 @@ static void ExactSectorBend(double *r, double le, double bending_angle, double L2 = SL*DRIFT2; double K1 = SL*KICK1; double K2 = SL*KICK2; - 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) \ @@ -90,11 +90,11 @@ static void ExactSectorBend(double *r, double le, double bending_angle, else { for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_bend(r6, irho, L1); - strthinkick(r6, A, B, K1, max_order); + kick(r6, A0, B0, A, B, K1, max_order); exact_bend(r6, irho, L2); - strthinkick(r6, A, B, K2, max_order); + kick(r6, A0, B0, A, B, K2, max_order); exact_bend(r6, irho, L2); - strthinkick(r6, A, B, K1, max_order); + kick(r6, A0, B0, A, B, K1, max_order); exact_bend(r6, irho, L1); } } @@ -122,9 +122,6 @@ static void ExactSectorBend(double *r, double le, double bending_angle, if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling); } } - /* Remove corrector component in polynomial coefficients */ - B[0] = B0; - A[0] = A0; } #if defined(MATLAB_MEX_FILE) || defined(PYAT) diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index d4ed0b06b..9d5088bce 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -1,7 +1,7 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "exactkickrad.c" +#include "exactbndkickrad.c" #include "exactbend.c" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -50,14 +50,14 @@ static void ExactSectorBendRad(double *r, double le, double irho, double L2 = SL*DRIFT2; double K1 = SL*KICK1; double K2 = SL*KICK2; - 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) \ @@ -89,11 +89,11 @@ static void ExactSectorBendRad(double *r, double le, double irho, for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_bend(r6, irho, L1); - ex_bndthinkickrad(r6, A, B, max_order, K1, irho, rad_const, diff_const, NULL); + kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, NULL); exact_bend(r6, irho, L2); - ex_bndthinkickrad(r6, A, B, max_order, K2, irho, rad_const, diff_const, NULL); + kick(r6, A0, B0, A, B, max_order, K2, irho, rad_const, diff_const, NULL); exact_bend(r6, irho, L2); - ex_bndthinkickrad(r6, A, B, max_order, K1, irho, rad_const, diff_const, NULL); + kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, NULL); exact_bend(r6, irho, L1); } @@ -120,9 +120,6 @@ static void ExactSectorBendRad(double *r, double le, double irho, if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling); } } - /* Remove corrector component in polynomial coefficients */ - B[0] = B0; - A[0] = A0; } #if defined(MATLAB_MEX_FILE) || defined(PYAT) diff --git a/atintegrators/a_fastdrift.c b/atintegrators/a_fastdrift.c new file mode 100644 index 000000000..6ebad8752 --- /dev/null +++ b/atintegrators/a_fastdrift.c @@ -0,0 +1,13 @@ + +static void drift(double* r, double NormL) + +/* NormL=(Physical Length)/(1+delta) is computed externally to speed up calculations + in the loop if momentum deviation (delta) does not change + such as in 4-th order symplectic integrator w/o radiation + */ + +{ + r[0] += NormL*r[1]; + r[2] += NormL*r[3]; + r[5] += NormL*(r[1]*r[1]+r[3]*r[3])/(2*(1+r[4])); +} diff --git a/atintegrators/a_h_k0h_k1h_kn.c b/atintegrators/a_h_k0h_k1h_kn.c new file mode 100644 index 000000000..c463d2fe8 --- /dev/null +++ b/atintegrators/a_h_k0h_k1h_kn.c @@ -0,0 +1,42 @@ +static void kick(double* r, double A0, double B0, double* A, double* B, double L, double irho, int max_order) +/***************************************************************************** +Calculate multipole kick in a curved elemrnt (bending magnet) +The reference coordinate system has the curvature given by the inverse +(design) radius irho. +IMPORTANT !!! +The magnetic field Bo that provides this curvature MUST NOT be included in the dipole term +PolynomB[1](MATLAB notation)(C: B[0] in this function) of the By field expansion + +The kick is given by + + e L L delta L x +theta = - --- B + ------- - ----- , + x p y rho 2 + 0 rho + + e L +theta = --- B + y p x + 0 + +*************************************************************************/ +{ + double ReSum = B[max_order]; + double ImSum = A[max_order]; + double ReSumTemp; + double x = r[0]; + double y = r[2]; + + /* recursively calculate the local transverse magnetic field */ + for (int i=max_order-1; i>=0; i--) { + ReSumTemp = ReSum*x - ImSum*y + B[i]; + ImSum = ImSum*x + ReSum*y + A[i]; + ReSum = ReSumTemp; + } + ReSum += B0; + ImSum += A0; + + r[1] -= L*(ReSum - (r[4]-x*irho)*irho + irho*B[1]*x*y); + r[3] += L*(ImSum - irho*B[1]*(x*x-y*y/2.0)); + r[5] += L*irho*x; /* pathlength */ +} diff --git a/atintegrators/b_drift_kick_drift_expanded.c b/atintegrators/b_drift_kick_drift_expanded.c new file mode 100644 index 000000000..e63bdf22e --- /dev/null +++ b/atintegrators/b_drift_kick_drift_expanded.c @@ -0,0 +1,2 @@ +#include "a_fastdrift.c" +#include "a_h_k0h_k1h_kn.c" \ No newline at end of file diff --git a/atintegrators/diff_drift_kick_drift_expanded.c b/atintegrators/diff_drift_kick_drift_expanded.c new file mode 100644 index 000000000..71b64cbea --- /dev/null +++ b/atintegrators/diff_drift_kick_drift_expanded.c @@ -0,0 +1,2 @@ +#include "diff_drift.c" +#include "diff_h_k0h_k1h_kn.c" diff --git a/atintegrators/diff_bnd_kick.c b/atintegrators/diff_h_k0h_k1h_kn.c similarity index 68% rename from atintegrators/diff_bnd_kick.c rename to atintegrators/diff_h_k0h_k1h_kn.c index 69ac38393..702468aa3 100644 --- a/atintegrators/diff_bnd_kick.c +++ b/atintegrators/diff_h_k0h_k1h_kn.c @@ -98,58 +98,29 @@ static double B2perp(double bx, double by, double irho, return (SQR(by*(1+x*irho)) + SQR(bx*(1+x*irho)) + SQR(bx*ypr - by*xpr))*v_norm2 ; } -static void diff_bnd_kick(double *r6, double *A, double *B, int max_order, - double L, double irho, double rad_const, double diff_const, double *bdiff) { - /* clang-format off */ -/***************************************************************************** -The design magnetic field Byo that provides this curvature By0 = irho * E0 /(c*e) -MUST NOT be included in the dipole term PolynomB(1)(MATLAB notation)(B[0] C notation) -of the By field expansion -HOWEVER!!! to calculate the effect of classical radiation the full field must be -used in the square of the |v x B|. -When calling B2perp(Bx, By, ...), use the By = ReSum + irho, where ReSum is the -normalized vertical field - sum of the polynomial terms in PolynomB. - -The kick is given by - - e L L delta L x -theta = - --- B + ------- - ----- , - x p y rho 2 - 0 rho - - e L -theta = --- B - y p x - 0 - - max_order - ---- - \ n - (B + iB )/ B rho = > (iA + B ) (x + iy) - y x / n n - ---- - n=0 - - ******************************************************************************/ - /* clang-format on */ - int i; +static void diff_kick(double *r6, double A0, double B0, double *A, double *B, int max_order, + double L, double irho, double rad_const, double diff_const, double *bdiff) { double ImSum = A[max_order]; double ReSum = B[max_order]; - double x, xpr, y, ypr, p_norm, dp_0, B2P, factor; double ReSumTemp; + double B2P, factor; + double p_norm = 1.0 / (1.0+r6[4]); + double x = r6[0]; + double y = r6[2]; + double dp_0 = r6[4]; /* save a copy of the initial value of dp/p */ + + /* calculate angles from momenta */ + double xpr = r6[1] * p_norm; + double ypr = r6[3] * p_norm; /* recursively calculate the local transverse magnetic field */ - for (i = max_order - 1; i >= 0; i--) { - ReSumTemp = ReSum * r6[0] - ImSum * r6[2] + B[i]; - ImSum = ImSum * r6[0] + ReSum * r6[2] + A[i]; + for (int i = max_order - 1; i >= 0; i--) { + ReSumTemp = ReSum * x - ImSum * y + B[i]; + ImSum = ImSum * x + ReSum * y + A[i]; ReSum = ReSumTemp; } - /* calculate angles from momenta */ - p_norm = 1.0 / (1.0+r6[4]); - x = r6[0]; - xpr = r6[1] * p_norm; - y = r6[2]; - ypr = r6[3] * p_norm; + ReSum += B0; + ImSum += A0; B2P = B2perp(ImSum, ReSum+irho, irho, x, xpr, y, ypr); factor = (1.0 + x*irho + (SQR(xpr) + SQR(ypr)) / 2.0) / SQR(p_norm) * L; @@ -159,8 +130,6 @@ theta = --- B thinkickB(r6, ReSum, ImSum, diff_const, B2P, factor, bdiff); } - dp_0 = r6[4]; /* save a copy of the initial value of dp/p */ - r6[4] -= rad_const * B2P * factor; /* recalculate momenta from angles after losing energy */ @@ -168,7 +137,7 @@ theta = --- B r6[1] = xpr / p_norm; r6[3] = ypr / p_norm; - r6[1] -= L * (ReSum - (dp_0 - r6[0] * irho) * irho); - r6[3] += L * ImSum; + r6[1] -= L * (ReSum - (dp_0-x*irho)*irho + irho*B[1]*x*y); + r6[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); r6[5] += L * irho * r6[0]; /* pathlength */ } \ No newline at end of file From 129602eff452c0da1a48e29013e636d9a97aa453 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 16:52:02 +0100 Subject: [PATCH 02/34] renaming --- atintegrators/ExactSectorBendPass.c | 19 +++++---- atintegrators/ExactSectorBendRadPass.c | 11 +++--- atintegrators/brad_bend_kick_bend.h | 2 + atintegrators/drift_exactbend.h | 39 +++++++++++++++++++ atintegrators/kick_k1h_kn.h | 20 ++++++++++ .../{exactkickrad.c => kickrad_k1h_kn.h} | 0 6 files changed, 75 insertions(+), 16 deletions(-) create mode 100644 atintegrators/brad_bend_kick_bend.h create mode 100644 atintegrators/drift_exactbend.h create mode 100644 atintegrators/kick_k1h_kn.h rename atintegrators/{exactkickrad.c => kickrad_k1h_kn.h} (100%) diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 3dfd642ae..c047b3109 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -1,8 +1,7 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "a_strthinkick.c" /* kick */ -#include "exactbend.c" +#include "b_bend_kick_bend.c" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -85,17 +84,17 @@ static void ExactSectorBend(double *r, double le, double bending_angle, bend_edge(r6, irho, -entrance_angle); if (num_int_steps == 0) { - exact_bend(r6, irho, le); + drift(r6, le, irho); } else { for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ - exact_bend(r6, irho, L1); - kick(r6, A0, B0, A, B, K1, max_order); - exact_bend(r6, irho, L2); - kick(r6, A0, B0, A, B, K2, max_order); - exact_bend(r6, irho, L2); - kick(r6, A0, B0, A, B, K1, max_order); - exact_bend(r6, irho, L1); + drift(r6, L1, irho); + kick(r6, A0, B0, A, B, max_order, K1, irho); + drift(r6, L2, irho); + kick(r6, A0, B0, A, B, max_order, K2, irho); + drift(r6, L2, irho); + kick(r6, A0, B0, A, B, max_order, K1, irho); + drift(r6, L1, irho); } } diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 9d5088bce..5beb88319 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -1,8 +1,7 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "exactbndkickrad.c" -#include "exactbend.c" +#include "brad_bend_kick_bend.h" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -88,13 +87,13 @@ static void ExactSectorBendRad(double *r, double le, double irho, bend_edge(r6, irho, -entrance_angle); for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ - exact_bend(r6, irho, L1); + drift(r6, L1, irho); kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, NULL); - exact_bend(r6, irho, L2); + drift(r6, L2, irho); kick(r6, A0, B0, A, B, max_order, K2, irho, rad_const, diff_const, NULL); - exact_bend(r6, irho, L2); + drift(r6, L2, irho); kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, NULL); - exact_bend(r6, irho, L1); + drift(r6, L1, irho); } /* Convert absolute path length to path lengthening */ diff --git a/atintegrators/brad_bend_kick_bend.h b/atintegrators/brad_bend_kick_bend.h new file mode 100644 index 000000000..63039de86 --- /dev/null +++ b/atintegrators/brad_bend_kick_bend.h @@ -0,0 +1,2 @@ +#include "drift_exactbend.h" +#include "kickrad_k1h_kn.h" diff --git a/atintegrators/drift_exactbend.h b/atintegrators/drift_exactbend.h new file mode 100644 index 000000000..079711c3d --- /dev/null +++ b/atintegrators/drift_exactbend.h @@ -0,0 +1,39 @@ +#include + +#ifndef PXYZ +#define PXYZ +static double pxyz(double dp1, double px, double py) +{ + return sqrt(dp1*dp1 - px*px - py*py); +} +#endif /*PXYZ*/ + +static void drift(double *r6, double irho, double L) +{ + /* Forest 12.18, bend-kick split, map W(L,irho) */ + + double dp1 = 1.0 + r6[delta_]; + double pz = pxyz(dp1, r6[px_], r6[py_]); + if (fabs(irho) < 1.e-6) { + double NormL = L / pz; + r6[x_] += r6[px_] * NormL; + r6[y_] += r6[py_] * NormL; + r6[ct_] += NormL * dp1; /* Absolute path length */ + } + else { + double pzmx = pz-(1.0+r6[x_]*irho); + double cs = cos(irho*L); + double sn = sin(irho*L); + double px = r6[px_]*cs + pzmx*sn; + double d2 = pxyz(dp1, 0.0, r6[py_]); + double dasin = L + (asin(r6[px_]/d2) - asin(px/d2))/irho; + double x = (pxyz(dp1,px,r6[py_]) - pzmx*cs + r6[px_]*sn - 1.0)/irho; + double dy = r6[py_]*dasin; + double dct = dp1*dasin; /* Absolute path length */ + + r6[x_] = x; + r6[px_] = px; + r6[y_] += dy; + r6[ct_] += dct; + } +} diff --git a/atintegrators/kick_k1h_kn.h b/atintegrators/kick_k1h_kn.h new file mode 100644 index 000000000..4940d577f --- /dev/null +++ b/atintegrators/kick_k1h_kn.h @@ -0,0 +1,20 @@ +static void kick(double* r, double A0, double B0, const double* A, const double* B, int max_order, double L,double irho) +{ + double ReSum = B[max_order]; + double ImSum = A[max_order]; + double ReSumTemp; + double x = r[0]; + double y = r[2]; + + /* recursively calculate the local transverse magnetic field */ + for (int i=max_order-1; i>=0; i--) { + ReSumTemp = ReSum*x - ImSum*y + B[i]; + ImSum = ImSum*x + ReSum*y + A[i]; + ReSum = ReSumTemp; + } + ReSum += B0; + ImSum += A0; + + r[1] -= L * (ReSum + irho*B[1]*x*y); + r[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); +} diff --git a/atintegrators/exactkickrad.c b/atintegrators/kickrad_k1h_kn.h similarity index 100% rename from atintegrators/exactkickrad.c rename to atintegrators/kickrad_k1h_kn.h From 9a2696453f53bed6961197c5c5dd73ef00548cbc Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 17:01:12 +0100 Subject: [PATCH 03/34] renaming --- atintegrators/b_drift_kick_drift_expanded.c | 4 ++-- atintegrators/{a_fastdrift.c => drift_fast.h} | 0 atintegrators/{a_h_k0h_k1h_kn.c => kick_h_k0h_k1h_kn.h} | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) rename atintegrators/{a_fastdrift.c => drift_fast.h} (100%) rename atintegrators/{a_h_k0h_k1h_kn.c => kick_h_k0h_k1h_kn.h} (92%) diff --git a/atintegrators/b_drift_kick_drift_expanded.c b/atintegrators/b_drift_kick_drift_expanded.c index e63bdf22e..a4e8bdfad 100644 --- a/atintegrators/b_drift_kick_drift_expanded.c +++ b/atintegrators/b_drift_kick_drift_expanded.c @@ -1,2 +1,2 @@ -#include "a_fastdrift.c" -#include "a_h_k0h_k1h_kn.c" \ No newline at end of file +#include "drift_fast.h" +#include "kick_h_k0h_k1h_kn.h" \ No newline at end of file diff --git a/atintegrators/a_fastdrift.c b/atintegrators/drift_fast.h similarity index 100% rename from atintegrators/a_fastdrift.c rename to atintegrators/drift_fast.h diff --git a/atintegrators/a_h_k0h_k1h_kn.c b/atintegrators/kick_h_k0h_k1h_kn.h similarity index 92% rename from atintegrators/a_h_k0h_k1h_kn.c rename to atintegrators/kick_h_k0h_k1h_kn.h index c463d2fe8..4ca53c241 100644 --- a/atintegrators/a_h_k0h_k1h_kn.c +++ b/atintegrators/kick_h_k0h_k1h_kn.h @@ -36,7 +36,7 @@ theta = --- B ReSum += B0; ImSum += A0; - r[1] -= L*(ReSum - (r[4]-x*irho)*irho + irho*B[1]*x*y); - r[3] += L*(ImSum - irho*B[1]*(x*x-y*y/2.0)); + r[1] -= L*(ReSum - (r[4]-x*irho)*irho); + r[3] += L*ImSum; r[5] += L*irho*x; /* pathlength */ } From 187984cd4a0ab24f2bab2a7847dc2eca824ba8e5 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 17:08:47 +0100 Subject: [PATCH 04/34] renaming --- atintegrators/b_bend_kick_bend.c | 2 + atintegrators/exactkickrad.c | 186 +++++++++++++++++++++++++++++++ atintegrators/kickrad_k1h_kn.h | 18 --- 3 files changed, 188 insertions(+), 18 deletions(-) create mode 100644 atintegrators/b_bend_kick_bend.c create mode 100644 atintegrators/exactkickrad.c diff --git a/atintegrators/b_bend_kick_bend.c b/atintegrators/b_bend_kick_bend.c new file mode 100644 index 000000000..d9bdd9bfc --- /dev/null +++ b/atintegrators/b_bend_kick_bend.c @@ -0,0 +1,2 @@ +#include "drift_exactbend.h" +#include "kick_k1h_kn.h" \ No newline at end of file diff --git a/atintegrators/exactkickrad.c b/atintegrators/exactkickrad.c new file mode 100644 index 000000000..726032bb7 --- /dev/null +++ b/atintegrators/exactkickrad.c @@ -0,0 +1,186 @@ +/*********************************************************************** + Note: in the US convention the transverse multipole field is written as: + + max_order+1 + --- + \ n-1 + (B + iB )/ B rho = > (ia + b ) (x + iy) + y x / n n + ---- + n=1 + is a polynomial in (x,y) with the highest order = MaxOrder + + + Using different index notation + + max_order + ---- + \ n + (B + iB )/ B rho = > (iA + B ) (x + iy) + y x / n n + ---- + n=0 + + A,B: i=0 ... max_order + [0] - dipole, [1] - quadrupole, [2] - sextupole ... + units for A,B[i] = 1/[m]^(i+1) + Coefficients are stored in the PolynomA, PolynomB field of the element + structure in MATLAB + + A[i] (C++,C) = PolynomA(i+1) (MATLAB) + B[i] (C++,C) = PolynomB(i+1) (MATLAB) + i = 0 .. MaxOrder + +*************************************************************************/ + +#define SQR(X) ((X)*(X)) + +static double StrB2perp(double bx, double by, + double x, double xpr, double y, double ypr) +/* Calculates sqr(|B x e|) , where e is a unit vector in the direction of velocity */ + +{ + /* components of the normalized velocity vector + double ex, ey, ez; + ex = xpr; + ey = ypr; + ez = sqrt(1 - xpr^2 - ypr^2); + + sqr(|B x e|) = sqr(|B|) * sqr(|e|) - sqr(B.e) + */ + + return SQR(bx) + SQR(by) - SQR(bx*xpr + by*ypr); +} + + +static double B2perp(double bx, double by, double irho, + double x, double xpr, double y, double ypr) + /* Calculates sqr(|e x B|) , where e is a unit vector in the direction of velocity */ + +{ + double nrm = SQR(1.0+x*irho); +// double v_norm2 = nrm + SQR(xpr) + SQR(ypr); + double v_norm2 = nrm + SQR(xpr)*(1.0-nrm) + SQR(ypr)*(1.0-nrm); + + /* components of the velocity vector + * double ex, ey, ez; + * ex = xpr; + * ey = ypr; + * ez = (1+x*irho) * sqrt(1 - xpr^2 - ypr^2); + */ + + return SQR(bx) + SQR(by) - SQR(bx*xpr + by*ypr)/v_norm2; +// return (SQR(by*(1+x*irho)) + SQR(bx*(1+x*irho)) + SQR(bx*ypr - by*xpr))/v_norm2 ; +} + +//static void ex_bndthinkickrad(double* r, double* A, double* B, double L, double irho, double E0, int max_order) +static void ex_bndthinkickrad(double* r, double* A, double* B, int max_order, + double L, double irho, double rad_const, double diff_const, double *bdiff) + +/***************************************************************************** +Calculate multipole kick in a curved element (bending magnet) +The reference coordinate system has the curvature given by the inverse +(design) radius irho. +IMPORTANT !!! +The magnetic field Bo that provides this curvature MUST NOT be included in the dipole term +PolynomB[1](MATLAB notation)(C: B[0] in this function) of the By field expansion +HOWEVER!!! to calculate the effect of classical radiation the full field must be +used in the square of the |v x B|. +When calling B2perp(Bx, By, ...), use the By = RESum + irho, where ImSum is the sum of +the polynomial terms in PolynomB. + + The kick is given by + + e L L delta L x + theta = - --- B + ------- - ----- , + x p y rho 2 + 0 rho + + e L + theta = --- B + y p x + 0 + + ******************************************************************************/ +{ + int i; + double ImSum = A[max_order]; + double ReSum = B[max_order]; + double ReSumTemp; + double x ,xpr, y, ypr, p_norm, B2P; + + for (i=max_order-1; i>=0; i--) { + ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i]; + ImSum = ImSum*r[0] + ReSum*r[2] + A[i]; + ReSum = ReSumTemp; + } + + /* calculate angles from momentums */ + p_norm = 1/(1+r[4]); + x = r[0]; + xpr = r[1]*p_norm; + y = r[2]; + ypr = r[3]*p_norm; + + B2P = B2perp(ImSum, ReSum+irho, irho, x , xpr, y ,ypr); + + /* Momentum loss */ + r[4] -= rad_const * SQR(1+r[4]) * B2P * (1.0+x*irho) * L / sqrt(1.0 - xpr*xpr - ypr*ypr); +// r[4] = r[4] - CRAD*SQR(1+r[4])*B2P*(1 + x*irho + (SQR(xpr)+SQR(ypr))/2 )*L; + + /* recalculate momentums from angles after losing energy for radiation */ + p_norm = 1/(1+r[4]); + r[1] = xpr/p_norm; + r[3] = ypr/p_norm; + + /* Multipole kick */ + r[1] -= L*ReSum; + r[3] += L*ImSum; +} + +//static void ex_strthinkickrad(double* r, const double* A, const double* B, double B0, double L, double E0, int max_order) +static void ex_strthinkickrad(double* r, const double* A, const double* B, int max_order, + double B0, double L, double rad_const, double diff_const, double *bdiff) +/***************************************************************************** + Calculate and apply a multipole kick to a 6-dimentional + phase space vector in a straight element ( quadrupole) + + IMPORTANT !!! + he reference coordinate system is straight but the field expansion may still + ontain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation, + [0], B[0] - C,C++ notation + + ******************************************************************************/ +{ + double ReSum = B[max_order]; + double ImSum = A[max_order]; + double ReSumTemp; + double x ,xpr, y, ypr, p_norm, B2P; + + for (int i=max_order-1; i>=0; i--) { + ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i]; + ImSum = ImSum*r[0] + ReSum*r[2] + A[i]; + ReSum = ReSumTemp; + } + + /* calculate angles from momentums */ + p_norm = 1/(1+r[4]); + x = r[0]; + xpr = r[1]*p_norm; + y = r[2]; + ypr = r[3]*p_norm; + + B2P = StrB2perp(ImSum, ReSum+B0 , x , xpr, y ,ypr); + + /* Momentum loss */ + r[4] -= rad_const * SQR(1+r[4]) * B2P * L / sqrt(1.0 - xpr*xpr - ypr*ypr); + + /* recalculate momentums from angles after losing energy for radiation */ + p_norm = 1/(1+r[4]); + r[1] = xpr/p_norm; + r[3] = ypr/p_norm; + + /* multipole kick */ + r[1] -= L*ReSum; + r[3] += L*ImSum; +} diff --git a/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h index 726032bb7..2e5e59bce 100644 --- a/atintegrators/kickrad_k1h_kn.h +++ b/atintegrators/kickrad_k1h_kn.h @@ -35,24 +35,6 @@ #define SQR(X) ((X)*(X)) -static double StrB2perp(double bx, double by, - double x, double xpr, double y, double ypr) -/* Calculates sqr(|B x e|) , where e is a unit vector in the direction of velocity */ - -{ - /* components of the normalized velocity vector - double ex, ey, ez; - ex = xpr; - ey = ypr; - ez = sqrt(1 - xpr^2 - ypr^2); - - sqr(|B x e|) = sqr(|B|) * sqr(|e|) - sqr(B.e) - */ - - return SQR(bx) + SQR(by) - SQR(bx*xpr + by*ypr); -} - - static double B2perp(double bx, double by, double irho, double x, double xpr, double y, double ypr) /* Calculates sqr(|e x B|) , where e is a unit vector in the direction of velocity */ From 7e1c8c8a6018b44245b752e6e013f510cfaade0e Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 17:27:41 +0100 Subject: [PATCH 05/34] renaming --- atintegrators/kick_h_k0h_k1h_kn.h | 6 +- atintegrators/kickrad_k1h_kn.h | 140 ++++++++++-------------------- 2 files changed, 49 insertions(+), 97 deletions(-) diff --git a/atintegrators/kick_h_k0h_k1h_kn.h b/atintegrators/kick_h_k0h_k1h_kn.h index 4ca53c241..93ceb1a77 100644 --- a/atintegrators/kick_h_k0h_k1h_kn.h +++ b/atintegrators/kick_h_k0h_k1h_kn.h @@ -36,7 +36,7 @@ theta = --- B ReSum += B0; ImSum += A0; - r[1] -= L*(ReSum - (r[4]-x*irho)*irho); - r[3] += L*ImSum; - r[5] += L*irho*x; /* pathlength */ + r[1] -= L * (ReSum - (r[4]-x*irho)*irho + irho*B[1]*x*y); + r[3] += L * (ImSum- irho*B[1]*(x*x-y*y/2.0)); + r[5] += L * irho*x; /* pathlength */ } diff --git a/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h index 2e5e59bce..a4a24747b 100644 --- a/atintegrators/kickrad_k1h_kn.h +++ b/atintegrators/kickrad_k1h_kn.h @@ -35,28 +35,25 @@ #define SQR(X) ((X)*(X)) -static double B2perp(double bx, double by, double irho, - double x, double xpr, double y, double ypr) - /* Calculates sqr(|e x B|) , where e is a unit vector in the direction of velocity */ - +static double B2perp(double bx, double by, double irho, double x, double xpr, double y, double ypr) +/* Calculates sqr(|e x B|), where e is a unit vector in the direction of velocity */ { + /* components of the velocity vector + double ex, ey, ez; + ex = xpr; + ey = ypr; + ez = (1+x*irho) * sqrt(1 - xpr^2 - ypr^2); + + sqr(|B x e|) = sqr(|B|) * sqr(|e|) - sqr(B.e) + */ double nrm = SQR(1.0+x*irho); -// double v_norm2 = nrm + SQR(xpr) + SQR(ypr); double v_norm2 = nrm + SQR(xpr)*(1.0-nrm) + SQR(ypr)*(1.0-nrm); - /* components of the velocity vector - * double ex, ey, ez; - * ex = xpr; - * ey = ypr; - * ez = (1+x*irho) * sqrt(1 - xpr^2 - ypr^2); - */ - return SQR(bx) + SQR(by) - SQR(bx*xpr + by*ypr)/v_norm2; -// return (SQR(by*(1+x*irho)) + SQR(bx*(1+x*irho)) + SQR(bx*ypr - by*xpr))/v_norm2 ; } //static void ex_bndthinkickrad(double* r, double* A, double* B, double L, double irho, double E0, int max_order) -static void ex_bndthinkickrad(double* r, double* A, double* B, int max_order, +static void kick(double* r, double A0, double B0, double* A, double* B, int max_order, double L, double irho, double rad_const, double diff_const, double *bdiff) /***************************************************************************** @@ -85,84 +82,39 @@ the polynomial terms in PolynomB. ******************************************************************************/ { - int i; - double ImSum = A[max_order]; - double ReSum = B[max_order]; - double ReSumTemp; - double x ,xpr, y, ypr, p_norm, B2P; - - for (i=max_order-1; i>=0; i--) { - ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i]; - ImSum = ImSum*r[0] + ReSum*r[2] + A[i]; - ReSum = ReSumTemp; - } - - /* calculate angles from momentums */ - p_norm = 1/(1+r[4]); - x = r[0]; - xpr = r[1]*p_norm; - y = r[2]; - ypr = r[3]*p_norm; - - B2P = B2perp(ImSum, ReSum+irho, irho, x , xpr, y ,ypr); - - /* Momentum loss */ - r[4] -= rad_const * SQR(1+r[4]) * B2P * (1.0+x*irho) * L / sqrt(1.0 - xpr*xpr - ypr*ypr); -// r[4] = r[4] - CRAD*SQR(1+r[4])*B2P*(1 + x*irho + (SQR(xpr)+SQR(ypr))/2 )*L; - - /* recalculate momentums from angles after losing energy for radiation */ - p_norm = 1/(1+r[4]); - r[1] = xpr/p_norm; - r[3] = ypr/p_norm; - - /* Multipole kick */ - r[1] -= L*ReSum; - r[3] += L*ImSum; -} - -//static void ex_strthinkickrad(double* r, const double* A, const double* B, double B0, double L, double E0, int max_order) -static void ex_strthinkickrad(double* r, const double* A, const double* B, int max_order, - double B0, double L, double rad_const, double diff_const, double *bdiff) -/***************************************************************************** - Calculate and apply a multipole kick to a 6-dimentional - phase space vector in a straight element ( quadrupole) - - IMPORTANT !!! - he reference coordinate system is straight but the field expansion may still - ontain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation, - [0], B[0] - C,C++ notation - - ******************************************************************************/ -{ - double ReSum = B[max_order]; - double ImSum = A[max_order]; - double ReSumTemp; - double x ,xpr, y, ypr, p_norm, B2P; - - for (int i=max_order-1; i>=0; i--) { - ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i]; - ImSum = ImSum*r[0] + ReSum*r[2] + A[i]; - ReSum = ReSumTemp; - } - - /* calculate angles from momentums */ - p_norm = 1/(1+r[4]); - x = r[0]; - xpr = r[1]*p_norm; - y = r[2]; - ypr = r[3]*p_norm; - - B2P = StrB2perp(ImSum, ReSum+B0 , x , xpr, y ,ypr); - - /* Momentum loss */ - r[4] -= rad_const * SQR(1+r[4]) * B2P * L / sqrt(1.0 - xpr*xpr - ypr*ypr); - - /* recalculate momentums from angles after losing energy for radiation */ - p_norm = 1/(1+r[4]); - r[1] = xpr/p_norm; - r[3] = ypr/p_norm; - - /* multipole kick */ - r[1] -= L*ReSum; - r[3] += L*ImSum; + double ReSum = B[max_order]; + double ImSum = A[max_order]; + double ReSumTemp; + double B2P; + double p_norm = 1.0 / (1.0+r6[4]); + double x = r6[0]; + double y = r6[2]; + double dp_0 = r6[4]; /* save a copy of the initial value of dp/p */ + + /* calculate angles from momenta */ + double xpr = r6[1] * p_norm; + double ypr = r6[3] * p_norm; + + /* recursively calculate the local transverse magnetic field */ + for (int i = max_order - 1; i >= 0; i--) { + ReSumTemp = ReSum * x - ImSum * y + B[i]; + ImSum = ImSum * x + ReSum * y + A[i]; + ReSum = ReSumTemp; + } + ReSum += B0; + ImSum += A0; + + B2P = B2perp(ImSum, ReSum+irho, irho, x , xpr, y ,ypr); + + /* Momentum loss */ + r[4] -= rad_const * SQR(1+r[4]) * B2P * (1.0+x*irho) * L / sqrt(1.0 - xpr*xpr - ypr*ypr); + + /* recalculate momentums from angles after losing energy for radiation */ + p_norm = 1/(1+r[4]); + r[1] = xpr/p_norm; + r[3] = ypr/p_norm; + + /* Multipole kick */ + r[1] -= L * (ReSum + irho*B[1]*x*y); + r[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); } From 9820e1efcb96fb4e9716a06ea3e6d0eff66bd118 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 17:35:16 +0100 Subject: [PATCH 06/34] renaming --- atintegrators/diff_h_k0h_k1h_kn.c | 4 ++-- atintegrators/kickrad_k1h_kn.h | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/atintegrators/diff_h_k0h_k1h_kn.c b/atintegrators/diff_h_k0h_k1h_kn.c index 702468aa3..5060dfffa 100644 --- a/atintegrators/diff_h_k0h_k1h_kn.c +++ b/atintegrators/diff_h_k0h_k1h_kn.c @@ -115,8 +115,8 @@ static void diff_kick(double *r6, double A0, double B0, double *A, double *B, in /* recursively calculate the local transverse magnetic field */ for (int i = max_order - 1; i >= 0; i--) { - ReSumTemp = ReSum * x - ImSum * y + B[i]; - ImSum = ImSum * x + ReSum * y + A[i]; + ReSumTemp = ReSum*x - ImSum*y + B[i]; + ImSum = ImSum*x + ReSum*y + A[i]; ReSum = ReSumTemp; } ReSum += B0; diff --git a/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h index a4a24747b..a9b2f6510 100644 --- a/atintegrators/kickrad_k1h_kn.h +++ b/atintegrators/kickrad_k1h_kn.h @@ -97,8 +97,8 @@ the polynomial terms in PolynomB. /* recursively calculate the local transverse magnetic field */ for (int i = max_order - 1; i >= 0; i--) { - ReSumTemp = ReSum * x - ImSum * y + B[i]; - ImSum = ImSum * x + ReSum * y + A[i]; + ReSumTemp = ReSum*x - ImSum*y + B[i]; + ImSum = ImSum*x + ReSum*y + A[i]; ReSum = ReSumTemp; } ReSum += B0; From dd232767338ecf49c4da12fee9ecc4cf751f2d73 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 17:42:30 +0100 Subject: [PATCH 07/34] renaming --- atintegrators/kick_h_k0h_k1h_kn.h | 12 ++++++------ atintegrators/kick_k1h_kn.h | 10 +++++----- atintegrators/kickrad_k1h_kn.h | 14 +++++++------- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/atintegrators/kick_h_k0h_k1h_kn.h b/atintegrators/kick_h_k0h_k1h_kn.h index 93ceb1a77..36318277b 100644 --- a/atintegrators/kick_h_k0h_k1h_kn.h +++ b/atintegrators/kick_h_k0h_k1h_kn.h @@ -1,4 +1,4 @@ -static void kick(double* r, double A0, double B0, double* A, double* B, double L, double irho, int max_order) +static void kick(double* r6, double A0, double B0, double* A, double* B, double L, double irho, int max_order) /***************************************************************************** Calculate multipole kick in a curved elemrnt (bending magnet) The reference coordinate system has the curvature given by the inverse @@ -24,8 +24,8 @@ theta = --- B double ReSum = B[max_order]; double ImSum = A[max_order]; double ReSumTemp; - double x = r[0]; - double y = r[2]; + double x = r6[0]; + double y = r6[2]; /* recursively calculate the local transverse magnetic field */ for (int i=max_order-1; i>=0; i--) { @@ -36,7 +36,7 @@ theta = --- B ReSum += B0; ImSum += A0; - r[1] -= L * (ReSum - (r[4]-x*irho)*irho + irho*B[1]*x*y); - r[3] += L * (ImSum- irho*B[1]*(x*x-y*y/2.0)); - r[5] += L * irho*x; /* pathlength */ + r6[1] -= L * (ReSum - (r6[4]-x*irho)*irho + irho*B[1]*x*y); + r6[3] += L * (ImSum- irho*B[1]*(x*x-y*y/2.0)); + r6[5] += L * irho*x; /* pathlength */ } diff --git a/atintegrators/kick_k1h_kn.h b/atintegrators/kick_k1h_kn.h index 4940d577f..58c2d8d1a 100644 --- a/atintegrators/kick_k1h_kn.h +++ b/atintegrators/kick_k1h_kn.h @@ -1,10 +1,10 @@ -static void kick(double* r, double A0, double B0, const double* A, const double* B, int max_order, double L,double irho) +static void kick(double* r6, double A0, double B0, const double* A, const double* B, int max_order, double L,double irho) { double ReSum = B[max_order]; double ImSum = A[max_order]; double ReSumTemp; - double x = r[0]; - double y = r[2]; + double x = r6[0]; + double y = r6[2]; /* recursively calculate the local transverse magnetic field */ for (int i=max_order-1; i>=0; i--) { @@ -15,6 +15,6 @@ static void kick(double* r, double A0, double B0, const double* A, const double* ReSum += B0; ImSum += A0; - r[1] -= L * (ReSum + irho*B[1]*x*y); - r[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); + r6[1] -= L * (ReSum + irho*B[1]*x*y); + r6[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); } diff --git a/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h index a9b2f6510..e990dbed4 100644 --- a/atintegrators/kickrad_k1h_kn.h +++ b/atintegrators/kickrad_k1h_kn.h @@ -53,7 +53,7 @@ static double B2perp(double bx, double by, double irho, double x, double xpr, do } //static void ex_bndthinkickrad(double* r, double* A, double* B, double L, double irho, double E0, int max_order) -static void kick(double* r, double A0, double B0, double* A, double* B, int max_order, +static void kick(double* r6, double A0, double B0, double* A, double* B, int max_order, double L, double irho, double rad_const, double diff_const, double *bdiff) /***************************************************************************** @@ -107,14 +107,14 @@ the polynomial terms in PolynomB. B2P = B2perp(ImSum, ReSum+irho, irho, x , xpr, y ,ypr); /* Momentum loss */ - r[4] -= rad_const * SQR(1+r[4]) * B2P * (1.0+x*irho) * L / sqrt(1.0 - xpr*xpr - ypr*ypr); + r6[4] -= rad_const * SQR(1+r6[4]) * B2P * (1.0+x*irho) * L / sqrt(1.0 - xpr*xpr - ypr*ypr); /* recalculate momentums from angles after losing energy for radiation */ - p_norm = 1/(1+r[4]); - r[1] = xpr/p_norm; - r[3] = ypr/p_norm; + p_norm = 1/(1+r6[4]); + r6[1] = xpr/p_norm; + r6[3] = ypr/p_norm; /* Multipole kick */ - r[1] -= L * (ReSum + irho*B[1]*x*y); - r[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); + r6[1] -= L * (ReSum + irho*B[1]*x*y); + r6[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); } From 5b7711e6ae44aa4b91d88380a9c94350413f66d1 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 25 Mar 2026 17:54:17 +0100 Subject: [PATCH 08/34] renaming --- atintegrators/kickrad_k1h_kn.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h index e990dbed4..4aaf85c00 100644 --- a/atintegrators/kickrad_k1h_kn.h +++ b/atintegrators/kickrad_k1h_kn.h @@ -110,7 +110,7 @@ the polynomial terms in PolynomB. r6[4] -= rad_const * SQR(1+r6[4]) * B2P * (1.0+x*irho) * L / sqrt(1.0 - xpr*xpr - ypr*ypr); /* recalculate momentums from angles after losing energy for radiation */ - p_norm = 1/(1+r6[4]); + p_norm = 1.0 / (1.0+r6[4]); r6[1] = xpr/p_norm; r6[3] = ypr/p_norm; From 7493422ab7ff7f52a819e1185b8ee36b25a70c87 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Fri, 27 Mar 2026 13:55:42 +0100 Subject: [PATCH 09/34] bug fix --- atintegrators/drift_exactbend.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atintegrators/drift_exactbend.h b/atintegrators/drift_exactbend.h index 079711c3d..336551383 100644 --- a/atintegrators/drift_exactbend.h +++ b/atintegrators/drift_exactbend.h @@ -8,7 +8,7 @@ static double pxyz(double dp1, double px, double py) } #endif /*PXYZ*/ -static void drift(double *r6, double irho, double L) +static void drift(double *r6, double L, double irho) { /* Forest 12.18, bend-kick split, map W(L,irho) */ From a471c14928f076acd60871846a8d34e9b50da9cc Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Sun, 29 Mar 2026 11:34:11 +0200 Subject: [PATCH 10/34] merged from master --- atintegrators/ExactMultipoleRadPass.c | 6 +----- atintegrators/ExactSectorBendRadPass.c | 2 +- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index e4a2f1a94..472a2a9fe 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/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; diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 5beb88319..f3fb43b39 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -185,7 +185,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; - gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error(); + gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);; check_error(); ExactSectorBendRad(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, From 81acad9360ef09ae4caa2a341a5c9a7cea6750e3 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Sun, 29 Mar 2026 13:58:48 +0200 Subject: [PATCH 11/34] typo --- atintegrators/ExactSectorBendRadPass.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index f3fb43b39..5beb88319 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -185,7 +185,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; - gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);; check_error(); + gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error(); ExactSectorBendRad(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, From 595f36393e22672e553f7bd490e7b98f94a83f3b Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 1 Apr 2026 11:55:50 +0200 Subject: [PATCH 12/34] bug fix --- atintegrators/diff_h_k0h_k1h_kn.c | 4 ++-- atintegrators/kick_h_k0h_k1h_kn.h | 4 ++-- atintegrators/kick_k1h_kn.h | 4 ++-- atintegrators/kickrad_k1h_kn.h | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/atintegrators/diff_h_k0h_k1h_kn.c b/atintegrators/diff_h_k0h_k1h_kn.c index 5060dfffa..eebef678f 100644 --- a/atintegrators/diff_h_k0h_k1h_kn.c +++ b/atintegrators/diff_h_k0h_k1h_kn.c @@ -137,7 +137,7 @@ static void diff_kick(double *r6, double A0, double B0, double *A, double *B, in r6[1] = xpr / p_norm; r6[3] = ypr / p_norm; - r6[1] -= L * (ReSum - (dp_0-x*irho)*irho + irho*B[1]*x*y); - r6[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); + r6[1] -= L * (ReSum - (dp_0-x*irho)*irho + irho*B[1]*(x*x-0.5*y*y)); + r6[3] += L * (ImSum + irho*B[1]*x*y); r6[5] += L * irho * r6[0]; /* pathlength */ } \ No newline at end of file diff --git a/atintegrators/kick_h_k0h_k1h_kn.h b/atintegrators/kick_h_k0h_k1h_kn.h index 36318277b..1e4545270 100644 --- a/atintegrators/kick_h_k0h_k1h_kn.h +++ b/atintegrators/kick_h_k0h_k1h_kn.h @@ -36,7 +36,7 @@ theta = --- B ReSum += B0; ImSum += A0; - r6[1] -= L * (ReSum - (r6[4]-x*irho)*irho + irho*B[1]*x*y); - r6[3] += L * (ImSum- irho*B[1]*(x*x-y*y/2.0)); + r6[1] -= L * (ReSum - (r6[4]-x*irho)*irho + irho*B[1]*(x*x-0.5*y*y)); + r6[3] += L * (ImSum + irho*B[1]*x*y); r6[5] += L * irho*x; /* pathlength */ } diff --git a/atintegrators/kick_k1h_kn.h b/atintegrators/kick_k1h_kn.h index 58c2d8d1a..11a1dca16 100644 --- a/atintegrators/kick_k1h_kn.h +++ b/atintegrators/kick_k1h_kn.h @@ -15,6 +15,6 @@ static void kick(double* r6, double A0, double B0, const double* A, const double ReSum += B0; ImSum += A0; - r6[1] -= L * (ReSum + irho*B[1]*x*y); - r6[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); + r6[1] -= L * (ReSum + irho*B[1]*(x*x-0.5*y*y)); + r6[3] += L * (ImSum + irho*B[1]*x*y); } diff --git a/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h index 4aaf85c00..833bec85d 100644 --- a/atintegrators/kickrad_k1h_kn.h +++ b/atintegrators/kickrad_k1h_kn.h @@ -115,6 +115,6 @@ the polynomial terms in PolynomB. r6[3] = ypr/p_norm; /* Multipole kick */ - r6[1] -= L * (ReSum + irho*B[1]*x*y); - r6[3] += L * (ImSum - irho*B[1]*(x*x-y*y/2.0)); + r6[1] -= L * (ReSum + irho*B[1]*(x*x-0.5*y*y)); + r6[3] += L * (ImSum + irho*B[1]*x*y); } From c1e1a850736a9a0623ff36943de43b5bdf873443 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 1 Apr 2026 15:10:57 +0200 Subject: [PATCH 13/34] add the last multipole to multipole fringe field --- atintegrators/ExactSectorBendPass.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index c047b3109..4a0b285d8 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -80,7 +80,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, if (FringeBendEntrance) bend_fringe(r6, irho, gK); if (FringeQuadEntrance) - multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + multipole_fringe(r6, le, A, B, max_order+1, 1.0, 1); bend_edge(r6, irho, -entrance_angle); if (num_int_steps == 0) { @@ -104,7 +104,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, -exit_angle); if (FringeQuadExit) - multipole_fringe(r6, le, A, B, max_order, -1.0, 1); + multipole_fringe(r6, le, A, B, max_order+1, -1.0, 1); if (FringeBendExit) bend_fringe(r6, -irho, gK); Yrot(r6, exit_angle); From 4317a0c292df744ff6c4427c2dab1d07d6602405 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 1 Apr 2026 15:11:59 +0200 Subject: [PATCH 14/34] add the last multipole to multipole fringe field --- atintegrators/ExactRectBendPass.c | 4 ++-- atintegrators/ExactRectangularBendPass.c | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index f8ab7a859..c5ef4a402 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -89,7 +89,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, if (FringeBendEntrance) bend_fringe(r6, irho, gK); if (FringeQuadEntrance) - multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + multipole_fringe(r6, le, A, B, max_order+1, 1.0, 1); bend_edge(r6, irho, phi2-entrance_angle); r6[0] += x0ref; @@ -115,7 +115,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); if (FringeQuadExit) - multipole_fringe(r6, le, A, B, max_order, -1.0, 1); + multipole_fringe(r6, le, A, B, max_order+1, -1.0, 1); if (FringeBendExit) bend_fringe(r6, -irho, gK); diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index a91fa447b..36d28d25f 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -90,7 +90,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, if (FringeBendEntrance) bend_fringe(r6, irho, gK); if (FringeQuadEntrance) - multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + multipole_fringe(r6, le, A, B, max_order+1, 1.0, 1); bend_edge(r6, irho, phi2-entrance_angle); /* integrator */ @@ -112,7 +112,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); if (FringeQuadExit) - multipole_fringe(r6, le, A, B, max_order, -1.0, 1); + multipole_fringe(r6, le, A, B, max_order+1, -1.0, 1); if (FringeBendExit) bend_fringe(r6, -irho, gK); From b9f729284c5fa5a6a6940f72c69b1ab7fe7f729f Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 1 Apr 2026 17:42:05 +0200 Subject: [PATCH 15/34] Updated tests --- atintegrators/BndMPoleSymplectic4Pass.c | 9 +++++---- atintegrators/BndMPoleSymplectic4RadPass.c | 3 ++- atintegrators/ExactRectangularBendPass.c | 4 ++-- atintegrators/ExactSectorBendPass.c | 7 ++++--- atintegrators/ExactSectorBendRadPass.c | 3 ++- atintegrators/b_bend_kick_bend.c | 2 -- atintegrators/b_drift_kick_drift_expanded.c | 2 -- atintegrators/kick_h_k0h_k1h_kn.h | 2 +- atintegrators/kick_k1h_kn.h | 2 +- pyat/test/test_physics.py | 21 +++++++++++---------- 10 files changed, 28 insertions(+), 27 deletions(-) delete mode 100644 atintegrators/b_bend_kick_bend.c delete mode 100644 atintegrators/b_drift_kick_drift_expanded.c diff --git a/atintegrators/BndMPoleSymplectic4Pass.c b/atintegrators/BndMPoleSymplectic4Pass.c index bb696e6d2..634cfefb8 100644 --- a/atintegrators/BndMPoleSymplectic4Pass.c +++ b/atintegrators/BndMPoleSymplectic4Pass.c @@ -2,7 +2,8 @@ #include "atelem.c" #include "atlalib.c" #include "atphyslib.c" -#include "b_drift_kick_drift_expanded.c" /* kick */ +#include "drift_fast.h" +#include "kick_h_k0h_k1h_kn.h" #include "quadfringe.c" /* QuadFringePassP, QuadFringePassN */ struct elem @@ -96,11 +97,11 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl /* integrator */ for (m=0; m < num_int_steps; m++) { /* Loop over slices */ drift(r6, NormL1); - kick(r6, A0, B0, A, B, K1, irho, max_order); + kick(r6, A0, B0, A, B, max_order, K1, irho); drift(r6, NormL2); - kick(r6, A0, B0, A, B, K2, irho, max_order); + kick(r6, A0, B0, A, B, max_order, K2, irho); drift(r6, NormL2); - kick(r6, A0, B0, A, B, K1, irho, max_order); + kick(r6, A0, B0, A, B, max_order, K1, irho); drift(r6, NormL1); } /* quadrupole gradient fringe */ diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index 4d361dc9e..1648c3472 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -1,7 +1,8 @@ #include "atelem.c" #include "atlalib.c" +#include "diff_drift.c" +#include "diff_h_k0h_k1h_kn.c" #include "diff_bend_fringe.c" -#include "diff_drift_kick_drift_expanded.c" #include "quadfringe.c" /* QuadFringePassP, QuadFringePassN */ struct elem diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 36d28d25f..a91fa447b 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -90,7 +90,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, if (FringeBendEntrance) bend_fringe(r6, irho, gK); if (FringeQuadEntrance) - multipole_fringe(r6, le, A, B, max_order+1, 1.0, 1); + multipole_fringe(r6, le, A, B, max_order, 1.0, 1); bend_edge(r6, irho, phi2-entrance_angle); /* integrator */ @@ -112,7 +112,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); if (FringeQuadExit) - multipole_fringe(r6, le, A, B, max_order+1, -1.0, 1); + multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) bend_fringe(r6, -irho, gK); diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 4a0b285d8..7b717bba6 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -1,7 +1,8 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "b_bend_kick_bend.c" +#include "drift_exactbend.h" +#include "kick_k1h_kn.h" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -80,7 +81,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, if (FringeBendEntrance) bend_fringe(r6, irho, gK); if (FringeQuadEntrance) - multipole_fringe(r6, le, A, B, max_order+1, 1.0, 1); + multipole_fringe(r6, le, A, B, max_order, 1.0, 1); bend_edge(r6, irho, -entrance_angle); if (num_int_steps == 0) { @@ -104,7 +105,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, -exit_angle); if (FringeQuadExit) - multipole_fringe(r6, le, A, B, max_order+1, -1.0, 1); + multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) bend_fringe(r6, -irho, gK); Yrot(r6, exit_angle); diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 5beb88319..248e86404 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -1,7 +1,8 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "brad_bend_kick_bend.h" +#include "drift_exactbend.h" +#include "kickrad_k1h_kn.h" #include "exactbendfringe.c" #include "exactmultipolefringe.c" diff --git a/atintegrators/b_bend_kick_bend.c b/atintegrators/b_bend_kick_bend.c deleted file mode 100644 index d9bdd9bfc..000000000 --- a/atintegrators/b_bend_kick_bend.c +++ /dev/null @@ -1,2 +0,0 @@ -#include "drift_exactbend.h" -#include "kick_k1h_kn.h" \ No newline at end of file diff --git a/atintegrators/b_drift_kick_drift_expanded.c b/atintegrators/b_drift_kick_drift_expanded.c deleted file mode 100644 index a4e8bdfad..000000000 --- a/atintegrators/b_drift_kick_drift_expanded.c +++ /dev/null @@ -1,2 +0,0 @@ -#include "drift_fast.h" -#include "kick_h_k0h_k1h_kn.h" \ No newline at end of file diff --git a/atintegrators/kick_h_k0h_k1h_kn.h b/atintegrators/kick_h_k0h_k1h_kn.h index 1e4545270..a699cb23a 100644 --- a/atintegrators/kick_h_k0h_k1h_kn.h +++ b/atintegrators/kick_h_k0h_k1h_kn.h @@ -1,4 +1,4 @@ -static void kick(double* r6, double A0, double B0, double* A, double* B, double L, double irho, int max_order) +static void kick(double* r6, double A0, double B0, const double* A, const double* B, int max_order, double L, double irho) /***************************************************************************** Calculate multipole kick in a curved elemrnt (bending magnet) The reference coordinate system has the curvature given by the inverse diff --git a/atintegrators/kick_k1h_kn.h b/atintegrators/kick_k1h_kn.h index 11a1dca16..b67c4727d 100644 --- a/atintegrators/kick_k1h_kn.h +++ b/atintegrators/kick_k1h_kn.h @@ -1,4 +1,4 @@ -static void kick(double* r6, double A0, double B0, const double* A, const double* B, int max_order, double L,double irho) +static void kick(double* r6, double A0, double B0, const double* A, const double* B, int max_order, double L, double irho) { double ReSum = B[max_order]; double ImSum = A[max_order]; diff --git a/pyat/test/test_physics.py b/pyat/test/test_physics.py index a2504ff6c..8d9c0b1c0 100644 --- a/pyat/test/test_physics.py +++ b/pyat/test/test_physics.py @@ -326,10 +326,10 @@ def test_get_tune_chrom(hmba_lattice): qharm = hmba_lattice.get_tune(method="interp_fft") qpharm = hmba_lattice.get_chrom(method="interp_fft") print(qlin, qharm) - assert_close(qlin, [0.2099983, 0.34001317], atol=1e-8) - assert_close(qharm, [0.20999833, 0.34001324], atol=1e-8) - assert_close(qplin, [5.734099, 3.917612], atol=1e-8) - assert_close(qpharm, [5.734123, 3.917639], atol=1e-8) + assert_close(qlin, [0.209998303584, 0.340013166682], atol=1e-8) + assert_close(qharm, [0.209998327527, 0.340013235807], atol=1e-8) + assert_close(qplin, [5.729114185134, 3.931703139652], atol=1e-8) + assert_close(qpharm, [5.729138523591, 3.931730165145], atol=1e-8) def test_nl_detuning_chromaticity(hmba_lattice): @@ -342,8 +342,8 @@ def test_nl_detuning_chromaticity(hmba_lattice): nlqplin, np.array( [ - [0.2101570, 5.730634, 151.87972, -18977.6808], - [0.3399707, 3.916998, 258.2324, -3529.81728], + [2.10157015e-01, 5.72564645e00, 1.51815354e02, -1.89809140e04], + [3.39970733e-01, 3.93108893e00, 2.58206212e02, -3.53120138e03], ] ), atol=1e-12, @@ -353,8 +353,8 @@ def test_nl_detuning_chromaticity(hmba_lattice): nlqpharm, np.array( [ - [0.2101570, 5.730630, 151.87968, -18977.7132], - [0.3399708, 3.916997, 258.23236, -3529.8072], + [2.10156985e-01, 5.72564318e00, 1.51815337e02, -1.89809436e04], + [3.39970750e-01, 3.93108797e00, 2.58206152e02, -3.53119072e03], ] ), atol=1e-12, @@ -362,13 +362,14 @@ def test_nl_detuning_chromaticity(hmba_lattice): ) assert_close( q0, - np.array([[0.210004, 0.340017], [0.210004, 0.340017]]), + np.array([[0.21000425, 0.34001688], [0.21000425, 0.34001688]]), atol=1e-12, rtol=1e-5, ) assert_close( q1, - np.array([[96183.925683, -104218.18371], [-104263.908197, 51684.400417]]), + np.array( + [[96172.00066329, -103658.3658934], [-103704.24788915, 51570.58307421]]), atol=1e-12, rtol=1e-5, ) From 3dad58f45f709f332b51619650bfa93cc9fc9bbc Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 1 Apr 2026 18:07:43 +0200 Subject: [PATCH 16/34] try with opposite sign --- atintegrators/exactmultipolefringe.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/atintegrators/exactmultipolefringe.c b/atintegrators/exactmultipolefringe.c index 1f5343fb8..93452af65 100644 --- a/atintegrators/exactmultipolefringe.c +++ b/atintegrators/exactmultipolefringe.c @@ -51,7 +51,8 @@ static void multipole_fringe(double *r6, double L, DU = B * DRX - A * DIX; DV = B * DIX + A * DRX; } - double f1 = -edge / 4.0 / (j + 1); +// double f1 = -edge / 4.0 / (j + 1); + double f1 = edge / 4.0 / (j + 1); U = U * f1; V = V * f1; From 370e7aa3325a6db53e81288c7ed5296ca67d320e Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 1 Apr 2026 18:20:16 +0200 Subject: [PATCH 17/34] relax Matlab test on orbit6 --- atintegrators/ExactRectBendPass.c | 4 ++-- atmat/attests/pytests.m | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index c5ef4a402..f8ab7a859 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -89,7 +89,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, if (FringeBendEntrance) bend_fringe(r6, irho, gK); if (FringeQuadEntrance) - multipole_fringe(r6, le, A, B, max_order+1, 1.0, 1); + multipole_fringe(r6, le, A, B, max_order, 1.0, 1); bend_edge(r6, irho, phi2-entrance_angle); r6[0] += x0ref; @@ -115,7 +115,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); if (FringeQuadExit) - multipole_fringe(r6, le, A, B, max_order+1, -1.0, 1); + multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) bend_fringe(r6, -irho, gK); diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m index fdc856fa8..450eca220 100644 --- a/atmat/attests/pytests.m +++ b/atmat/attests/pytests.m @@ -105,7 +105,7 @@ function orbit6(testCase,lat2,dp) porbit6=double(porbit6)'; % Matlab [~,morbit6]=findorbit6(lattice.m,dp=dp); - testCase.verifyEqual(morbit6,porbit6,AbsTol=2.E-12); + testCase.verifyEqual(morbit6,porbit6,AbsTol=3.E-12); end function m44(testCase,lat2,dp) From 8b3333bdc775e67bced488c68f617fd3a794ddf8 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 2 Apr 2026 11:36:30 +0200 Subject: [PATCH 18/34] test K1h on rectangularbend pass --- atintegrators/ExactRectangularBendPass.c | 33 +++++++++++++++++++++--- atintegrators/exactmultipolefringe.c | 3 +-- 2 files changed, 31 insertions(+), 5 deletions(-) diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index a91fa447b..b67445f32 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -34,6 +34,33 @@ struct elem double *KickAngle; }; +static void kick(double* r, const double* A, const double* B, double L, int max_order, double irho) +/***************************************************************************** + Calculate and apply a multipole kick to a 6-dimentional + phase space vector in a straight element (quadrupole) + + IMPORTANT !!! + The reference coordinate system is straight but the field expansion may still + contain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation, + A[0], B[0] - C,C++ notation + + ******************************************************************************/ +{ + int i; + double ReSum = B[max_order]; + double ImSum = A[max_order]; + double ReSumTemp; + double x = r[0]; + double y = r[2]; + for (i=max_order-1; i>=0; i--) { + ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i]; + ImSum = ImSum*r[0] + ReSum*r[2] + A[i]; + ReSum = ReSumTemp; + } + r[1] -= L * (ReSum + irho*B[1]*(x*x-0.5*y*y)); + r[3] += L * (ImSum + irho*B[1]*x*y); +} + static void ExactRectangularBend(double *r, double le, double bending_angle, double *A, double *B, int max_order, int num_int_steps, @@ -97,11 +124,11 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, r6[0] += x0ref; for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_drift(r6, L1); - strthinkick(r6, A, B, K1, max_order); + kick(r6, A, B, K1, max_order, irho); exact_drift(r6, L2); - strthinkick(r6, A, B, K2, max_order); + kick(r6, A, B, K2, max_order, irho); exact_drift(r6, L2); - strthinkick(r6, A, B, K1, max_order); + kick(r6, A, B, K1, max_order, irho); exact_drift(r6, L1); } r6[0] -= x0ref; diff --git a/atintegrators/exactmultipolefringe.c b/atintegrators/exactmultipolefringe.c index 93452af65..1f5343fb8 100644 --- a/atintegrators/exactmultipolefringe.c +++ b/atintegrators/exactmultipolefringe.c @@ -51,8 +51,7 @@ static void multipole_fringe(double *r6, double L, DU = B * DRX - A * DIX; DV = B * DIX + A * DRX; } -// double f1 = -edge / 4.0 / (j + 1); - double f1 = edge / 4.0 / (j + 1); + double f1 = -edge / 4.0 / (j + 1); U = U * f1; V = V * f1; From 6a50aab2b227fb2dc40ddf1265147487f35a6c36 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 2 Apr 2026 11:54:39 +0200 Subject: [PATCH 19/34] conversion of multipole fringe field --- pyat/at/load/xsuite.py | 27 ++++++++++----------------- 1 file changed, 10 insertions(+), 17 deletions(-) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 8e3373ef1..18fedc85b 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -115,11 +115,11 @@ |Straight magnet +------------------------+----------------+---------------------------+ | |*default* |"adaptive" |"drift-kick-drift-expanded"| +-----------------+------------------------+----------------+---------------------------+ - | |ExactSectorBendPass |"bend-kick-bend"| | + | |ExactSectorBendPass |"bend-kick-bend" | | +------------------------+----------------+---------------------------+ |Dipole |ExactRectangularBendPass|"drift-kick-drift-exact" | | +------------------------+----------------+---------------------------+ - | |*default* |"adaptive" |"rot-kick-rot" | + | |*default* |"adaptive" |"drift-kick-drift-expanded"| +-----------------+------------------------+----------------+---------------------------+ Longitudinal motion @@ -484,7 +484,12 @@ def xspoly(kmain: list[str], kerr: str) -> tuple[int, np.ndarray]: def _set_at_fringe(self) -> dict[str, Any]: """generate the AT fringe field description.""" - return {} + atparams = {} + if self.get("edge_entry_active", False): + atparams["FringeQuadEntrance"] = 1 + if self.get("edge_exit_active", False): + atparams["FringeQuadExit"] = 1 + return atparams def _set_xs_transforms(self, atparams: dict) -> None: """Generate Xsuite element displacements.""" @@ -540,7 +545,8 @@ def _set_xs_poly(self, atparams: dict) -> None: self["_isthick"] = length != 0.0 def _set_xs_fringe(self, atparams: dict) -> None: - """generate the Xsuite fringe field description.""" + self["edge_entry_active"] = bool(atparams.get("FringeQuadEntrance", 0)) + self["edge_exit_active"] = bool(atparams.get("FringeQuadExit", 0)) def _class_to_at(self, atparams: dict[str, Any]) -> type[elt.Element]: if atparams.get("Length", 0.0) == 0.0: @@ -571,19 +577,6 @@ class Quadrupole(Multipole): _atClass = elt.Quadrupole _mag_order: ClassVar[int] = 1 - def _set_at_fringe(self): - """generate the AT fringe field description.""" - atparams = {} - if self.get("edge_entry_active", 0): - atparams["FringeQuadEntrance"] = 1 - if self.get("edge_exit_active", 0): - atparams["FringeQuadExit"] = 1 - return atparams - - def _set_xs_fringe(self, atparams: dict): - self["edge_entry_active"] = atparams.get("FringeQuadEntrance", 0) - self["edge_exit_active"] = atparams.get("FringeQuadExit", 0) - class Sextupole(Multipole): """Xsuite Sextupole element.""" From 91472423fd0da3b743af41cbbbfe3d421756cf75 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 2 Apr 2026 15:05:43 +0200 Subject: [PATCH 20/34] Removed K1h correction in ExactRectangularBend --- atintegrators/ExactRectangularBendPass.c | 33 +++--------------------- 1 file changed, 3 insertions(+), 30 deletions(-) diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index b67445f32..a91fa447b 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -34,33 +34,6 @@ struct elem double *KickAngle; }; -static void kick(double* r, const double* A, const double* B, double L, int max_order, double irho) -/***************************************************************************** - Calculate and apply a multipole kick to a 6-dimentional - phase space vector in a straight element (quadrupole) - - IMPORTANT !!! - The reference coordinate system is straight but the field expansion may still - contain dipole terms: PolynomA(1), PolynomB(1) - in MATLAB notation, - A[0], B[0] - C,C++ notation - - ******************************************************************************/ -{ - int i; - double ReSum = B[max_order]; - double ImSum = A[max_order]; - double ReSumTemp; - double x = r[0]; - double y = r[2]; - for (i=max_order-1; i>=0; i--) { - ReSumTemp = ReSum*r[0] - ImSum*r[2] + B[i]; - ImSum = ImSum*r[0] + ReSum*r[2] + A[i]; - ReSum = ReSumTemp; - } - r[1] -= L * (ReSum + irho*B[1]*(x*x-0.5*y*y)); - r[3] += L * (ImSum + irho*B[1]*x*y); -} - static void ExactRectangularBend(double *r, double le, double bending_angle, double *A, double *B, int max_order, int num_int_steps, @@ -124,11 +97,11 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, r6[0] += x0ref; for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_drift(r6, L1); - kick(r6, A, B, K1, max_order, irho); + strthinkick(r6, A, B, K1, max_order); exact_drift(r6, L2); - kick(r6, A, B, K2, max_order, irho); + strthinkick(r6, A, B, K2, max_order); exact_drift(r6, L2); - kick(r6, A, B, K1, max_order, irho); + strthinkick(r6, A, B, K1, max_order); exact_drift(r6, L1); } r6[0] -= x0ref; From 27732850ed24eaa726e3f407a57aaf1cd5a177ac Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 2 Apr 2026 16:58:19 +0200 Subject: [PATCH 21/34] Added quadrupolar correction for wedge --- atintegrators/ExactRectangularBendPass.c | 10 ++++++++++ atintegrators/ExactSectorBendPass.c | 11 +++++++++++ atintegrators/exactbendfringe.c | 10 ++++++++++ 3 files changed, 31 insertions(+) diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index a91fa447b..0e590d564 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -56,12 +56,18 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = B[0]; double A0 = A[0]; + double k1_theta_entrance = 0.0; + double k1_theta_exit = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; A[0] += sin(KickAngle[1])/le; } B[0] += irho; + if (max_order >= 1) { + k1_theta_entrance = B[1] * entrance_angle; + k1_theta_exit = 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,\ @@ -91,6 +97,8 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, bend_fringe(r6, irho, gK); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + if (k1_theta_entrance != 0.0) + quad_wedge(r6, -k1_theta_entrance); bend_edge(r6, irho, phi2-entrance_angle); /* integrator */ @@ -111,6 +119,8 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); + if (k1_theta_exit != 0.0) + quad_wedge(r6, -k1_theta_exit); if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 7b717bba6..2cb979a27 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -52,11 +52,18 @@ static void ExactSectorBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = 0.0; double A0 = 0.0; + double k1_theta_entrance = 0.0; + double k1_theta_exit = 0.0; + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B0 = -sin(KickAngle[0]) / le; A0 = sin(KickAngle[1]) / le; } + if (max_order >= 1) { + k1_theta_entrance = B[1] * entrance_angle; + k1_theta_exit = 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,\ @@ -82,6 +89,8 @@ static void ExactSectorBend(double *r, double le, double bending_angle, bend_fringe(r6, irho, gK); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + if (k1_theta_entrance != 0.0) + quad_wedge(r6, -k1_theta_entrance); bend_edge(r6, irho, -entrance_angle); if (num_int_steps == 0) { @@ -104,6 +113,8 @@ static void ExactSectorBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, -exit_angle); + if (k1_theta_exit != 0.0) + quad_wedge(r6, -k1_theta_exit); if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) diff --git a/atintegrators/exactbendfringe.c b/atintegrators/exactbendfringe.c index 3c8e8ac52..b4469efd6 100644 --- a/atintegrators/exactbendfringe.c +++ b/atintegrators/exactbendfringe.c @@ -118,6 +118,16 @@ static void bend_edge(double *r6, double rhoinv, double theta) } } +static void quad_wedge(double *r6, double k1_theta) +{ + double x = r6[x_]; + double y = r6[y_]; + double dpx = k1_theta * (x*x - 0.5*y*y); + double dpy = k1_theta * x * y; + r6[px_] -= dpx; + r6[py_] += dpy; +} + static void bend_fringe_test(double *r6, double irho, double gK) { /* Forest 13.13, bend fringe in the hard-edge limit */ From 1061a18545cef9dcea59b38d98058dc13918a127 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Fri, 3 Apr 2026 20:23:59 +0200 Subject: [PATCH 22/34] num_multipople_kicks = 3 * NumIntSteps --- atintegrators/kick_k1h_kn.h | 15 ++++++++++----- pyat/at/load/xsuite.py | 27 ++++++++++++++------------- 2 files changed, 24 insertions(+), 18 deletions(-) diff --git a/atintegrators/kick_k1h_kn.h b/atintegrators/kick_k1h_kn.h index b67c4727d..65a1b1fb6 100644 --- a/atintegrators/kick_k1h_kn.h +++ b/atintegrators/kick_k1h_kn.h @@ -8,13 +8,18 @@ static void kick(double* r6, double A0, double B0, const double* A, const double /* recursively calculate the local transverse magnetic field */ for (int i=max_order-1; i>=0; i--) { - ReSumTemp = ReSum*x - ImSum*y + B[i]; - ImSum = ImSum*x + ReSum*y + A[i]; - ReSum = ReSumTemp; + ReSumTemp = ReSum*x - ImSum*y + B[i]; + ImSum = ImSum*x + ReSum*y + A[i]; + ReSum = ReSumTemp; } ReSum += B0; ImSum += A0; - r6[1] -= L * (ReSum + irho*B[1]*(x*x-0.5*y*y)); - r6[3] += L * (ImSum + irho*B[1]*x*y); + r6[1] -= L * ReSum; + r6[3] += L * ImSum; + + if (max_order >= 1) { /* K1h correction */ + r6[1] -= L * irho * B[1] * (x*x-0.5*y*y); + r6[3] += L * irho * B[1] * x * y; + } } diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 18fedc85b..302188639 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -329,8 +329,6 @@ def from_at(cls, match_model: bool = False, **atparams) -> XsElement: xs_model = cls._at2xsuite_model.get(None, None) if (integrator := cls._at_integrator) is not None: xsparams["integrator"] = integrator - else: - xsparams.pop("num_multipole_kicks", None) if xs_model is not None: xsparams["model"] = xs_model # Set the Xsuite class @@ -407,7 +405,6 @@ class Multipole(XsElement): _at_integrator = "yoshida4" _xsuite2at_attr = XsElement._xsuite2at_attr | { "order": "MaxOrder", - "num_multipole_kicks": "NumIntSteps", } def _set_at_transforms(self) -> dict: @@ -478,6 +475,8 @@ def xspoly(kmain: list[str], kerr: str) -> tuple[int, np.ndarray]: "PolynomB": polyb[: maxorder + 1], "PolynomA": polya[: maxorder + 1], } + if (nk := self.get("num_multipole_kicks", 0)) != 0: + atparams["NumIntSteps"] = (nk-1) // 3 + 1 if (taper := self.get("delta_taper")) is not None: atparams["FieldScaling"] = 1.0 + taper return atparams @@ -521,7 +520,7 @@ def _set_xs_transforms(self, atparams: dict) -> None: self.update(misalign) - def _set_xs_poly(self, atparams: dict) -> None: + def _set_xs_poly(self, atparams: dict, match_model: bool = False) -> None: """Generate the AT field expansion.""" pata = atparams.get("PolynomA", np.zeros(4)) pola = np.fromiter(poly_to_mad(pata), dtype=float, count=pata.size) @@ -540,6 +539,8 @@ def _set_xs_poly(self, atparams: dict) -> None: if np.any(pola) or np.any(polb): self["knl"] = list(polb) self["ksl"] = list(pola) + if match_model: + self["num_multipole_kicks"] = 3 * atparams["NumIntSteps"] if (scaling := atparams.get("FieldScaling")) is not None: self["delta_taper"] = scaling - 1.0 self["_isthick"] = length != 0.0 @@ -562,9 +563,9 @@ def _params_to_at(self, **atparams) -> dict[str, Any]: return atparams @classmethod - def from_at(cls, **atparams): - elem = super().from_at(**atparams) - elem._set_xs_poly(atparams) + def from_at(cls, match_model: bool = False, **atparams): + elem = super().from_at(match_model=match_model, **atparams) + elem._set_xs_poly(atparams, match_model=match_model) elem._set_xs_fringe(atparams) elem._set_xs_transforms(atparams) return elem @@ -649,8 +650,8 @@ def _params_to_at(self, **atparams) -> dict[str, Any]: return atparams @classmethod - def from_at(cls, **atparams): - elem = super().from_at(**atparams) + def from_at(cls, match_model: bool = False, **atparams): + elem = super().from_at(match_model=match_model, **atparams) elem["k0_from_h"] = True return elem @@ -698,8 +699,8 @@ def _params_to_at(self, **atparams) -> dict[str, Any]: return atparams @classmethod - def from_at(cls, **atparams): - elem = super().from_at(**atparams) + def from_at(cls, match_model: bool = False, **atparams): + elem = super().from_at(match_model=match_model, **atparams) elem["rbend_model"] = "straight-body" hangle = 0.5 * elem["angle"] elem["edge_entry_angle"] -= hangle @@ -739,8 +740,8 @@ def _params_to_at(self, **atparams) -> dict[str, Any]: return atparams @classmethod - def from_at(cls, **atparams): - elem = super().from_at(**atparams) + def from_at(cls, match_model: bool = False, **atparams): + elem = super().from_at(match_model=match_model, **atparams) elem._set_xs_lag(atparams) return elem From e840ef6032ddd66be8238c793843c5a39f233d4f Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Tue, 21 Apr 2026 16:56:23 +0200 Subject: [PATCH 23/34] corrected dipole export to xsuite, quad wedge correction applied in all "exact" integrators --- atintegrators/ExactRectBendPass.c | 10 +++++++ atintegrators/ExactRectBendRadPass.c | 10 +++++++ atintegrators/ExactRectangularBendRadPass.c | 10 +++++++ atintegrators/ExactSectorBendRadPass.c | 10 +++++++ pyat/at/load/xsuite.py | 29 ++++++++++++++------- 5 files changed, 60 insertions(+), 9 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index f8ab7a859..b9ba3894c 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -56,11 +56,17 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = B[0]; double A0 = A[0]; + double k1_theta_entrance = 0.0; + double k1_theta_exit = 0.0; 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_theta_entrance = B[1] * entrance_angle; + k1_theta_exit = 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,\ @@ -90,6 +96,8 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, bend_fringe(r6, irho, gK); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + if (k1_theta_entrance != 0.0) + quad_wedge(r6, -k1_theta_entrance); bend_edge(r6, irho, phi2-entrance_angle); r6[0] += x0ref; @@ -114,6 +122,8 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); + if (k1_theta_exit != 0.0) + quad_wedge(r6, -k1_theta_exit); if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 203385f0d..e47d1bae0 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -59,11 +59,17 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double A0 = A[0]; double rad_const = RAD_CONST*pow(gamma, 3); double diff_const = DIF_CONST*pow(gamma, 5); + double k1_theta_entrance = 0.0; + double k1_theta_exit = 0.0; 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_theta_entrance = B[1] * entrance_angle; + k1_theta_exit = 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,\ @@ -93,6 +99,8 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, bend_fringe(r6, irho, gK); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + if (k1_theta_entrance != 0.0) + quad_wedge(r6, -k1_theta_entrance); bend_edge(r6, irho, phi2-entrance_angle); r6[0] += x0ref; @@ -112,6 +120,8 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); + if (k1_theta_exit != 0.0) + quad_wedge(r6, -k1_theta_exit); if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index 66f2d9d61..e5dfbf76c 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -59,12 +59,18 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double A0 = A[0]; double rad_const = RAD_CONST*pow(gamma, 3); double diff_const = DIF_CONST*pow(gamma, 5); + double k1_theta_entrance = 0.0; + double k1_theta_exit = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; A[0] += sin(KickAngle[1])/le; } B[0] += irho; + if (max_order >= 1) { + k1_theta_entrance = B[1] * entrance_angle; + k1_theta_exit = 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,\ @@ -94,6 +100,8 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, bend_fringe(r6, irho, gK); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + if (k1_theta_entrance != 0.0) + quad_wedge(r6, -k1_theta_entrance); bend_edge(r6, irho, phi2-entrance_angle); /* integrator */ @@ -114,6 +122,8 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* edge focus */ bend_edge(r6, irho, phi2-exit_angle); + if (k1_theta_exit != 0.0) + quad_wedge(r6, -k1_theta_exit); if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 248e86404..8366a509b 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -54,11 +54,17 @@ static void ExactSectorBendRad(double *r, double le, double irho, double diff_const = DIF_CONST*pow(gamma, 5); double B0 = 0.0; double A0 = 0.0; + double k1_theta_entrance = 0.0; + double k1_theta_exit = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B0 = -sin(KickAngle[0]) / le; A0 = sin(KickAngle[1]) / le; } + if (max_order >= 1) { + k1_theta_entrance = B[1] * entrance_angle; + k1_theta_exit = 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,\ @@ -85,6 +91,8 @@ static void ExactSectorBendRad(double *r, double le, double irho, bend_fringe(r6, irho, gK); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + if (k1_theta_entrance != 0.0) + quad_wedge(r6, -k1_theta_entrance); bend_edge(r6, irho, -entrance_angle); for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ @@ -102,6 +110,8 @@ static void ExactSectorBendRad(double *r, double le, double irho, /* edge focus */ bend_edge(r6, irho, -exit_angle); + if (k1_theta_exit != 0.0) + quad_wedge(r6, -k1_theta_exit); if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 302188639..4ce97d0d8 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -476,7 +476,7 @@ def xspoly(kmain: list[str], kerr: str) -> tuple[int, np.ndarray]: "PolynomA": polya[: maxorder + 1], } if (nk := self.get("num_multipole_kicks", 0)) != 0: - atparams["NumIntSteps"] = (nk-1) // 3 + 1 + atparams["NumIntSteps"] = (nk - 1) // 3 + 1 if (taper := self.get("delta_taper")) is not None: atparams["FieldScaling"] = 1.0 + taper return atparams @@ -522,25 +522,36 @@ def _set_xs_transforms(self, atparams: dict) -> None: def _set_xs_poly(self, atparams: dict, match_model: bool = False) -> None: """Generate the AT field expansion.""" + + def extract(poly, ord): + try: + v = poly[ord] + except IndexError: + v = 0.0 + else: + poly[ord] = 0.0 + return v + pata = atparams.get("PolynomA", np.zeros(4)) pola = np.fromiter(poly_to_mad(pata), dtype=float, count=pata.size) patb = atparams.get("PolynomB", np.zeros(4)) polb = np.fromiter(poly_to_mad(patb), dtype=float, count=patb.size) length = atparams.get("Length") korder = getattr(self, "_mag_order", None) - if korder is not None: - self["k" + str(korder)] = polb[korder] - self["k" + str(korder) + "s"] = pola[korder] - pola[korder] = 0.0 - polb[korder] = 0.0 + if korder == 0: # dipole + self["k1"] = extract(polb, 1) + self["k2"] = extract(polb, 2) + elif korder is not None: # quadrupole, sextupole, octupole + self["k" + str(korder)] = extract(polb, korder) + self["k" + str(korder) + "s"] = extract(pola, korder) if length > 0.0: polb *= length pola *= length if np.any(pola) or np.any(polb): self["knl"] = list(polb) self["ksl"] = list(pola) - if match_model: - self["num_multipole_kicks"] = 3 * atparams["NumIntSteps"] + if match_model and (numintsteps := atparams.get("NumIntSteps")) is not None: + self["num_multipole_kicks"] = 3 * numintsteps if (scaling := atparams.get("FieldScaling")) is not None: self["delta_taper"] = scaling - 1.0 self["_isthick"] = length != 0.0 @@ -613,7 +624,7 @@ class Bend(Multipole): "edge_entry_angle": "EntranceAngle", "edge_exit_angle": "ExitAngle", } - _mag_order = 1 + _mag_order = 0 _edge_to_xs: ClassVar[dict[bool, dict[bool, str]]] = { True: {True: "full", False: "dipole-only"}, False: {True: "linear", False: "linear"}, From 243ce90aee5a6f54ab6e3fb64ce6e4b69a80a9b3 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 22 Apr 2026 10:44:01 +0200 Subject: [PATCH 24/34] Processing of the fringe field integrals fint1 and fint2 --- atintegrators/ExactRectBendPass.c | 37 +++++++++++++-------- atintegrators/ExactRectBendRadPass.c | 37 +++++++++++++-------- atintegrators/ExactRectangularBendPass.c | 37 +++++++++++++-------- atintegrators/ExactRectangularBendRadPass.c | 37 +++++++++++++-------- atintegrators/ExactSectorBendPass.c | 33 ++++++++++-------- atintegrators/ExactSectorBendRadPass.c | 32 +++++++++++------- 6 files changed, 136 insertions(+), 77 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index b9ba3894c..f630b8329 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -22,7 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gKEntrance; + double gKExit; double x0ref; double refdz; double *R1; @@ -40,7 +41,8 @@ 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 gKEntrance, double gKExit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -59,7 +61,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double k1_theta_entrance = 0.0; double k1_theta_exit = 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; } @@ -70,7 +72,7 @@ static void ExactRectangularBend(double *r, double le, double bending_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,\ + irho,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ LR,le,phi2) @@ -93,7 +95,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ if (FringeBendEntrance) - bend_fringe(r6, irho, gK); + bend_fringe(r6, irho, gKEntrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (k1_theta_entrance != 0.0) @@ -127,7 +129,7 @@ static void ExactRectangularBend(double *r, double le, double bending_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, gKExit); /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -168,7 +170,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(); @@ -202,7 +206,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->FringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gKEntrance=FullGap * FringeInt1; + Elem->gKExit=FullGap * FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -217,7 +222,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->gKEntrance, Elem->gKExit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -251,7 +257,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(); @@ -277,7 +285,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, + FullGap*FringeInt1, FullGap*FringeInt1, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { @@ -294,12 +303,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")); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index e47d1bae0..f9f4430a2 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -23,7 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gKEntrance; + double gKExit; double x0ref; double refdz; double *R1; @@ -41,7 +42,8 @@ static void ExactRectangularBendRad(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 gKEntrance, double gKExit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -62,7 +64,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double k1_theta_entrance = 0.0; double k1_theta_exit = 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; } @@ -73,7 +75,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_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,\ + irho,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ le,phi2,rad_const,diff_const) @@ -96,7 +98,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* edge focus */ if (FringeBendEntrance) - bend_fringe(r6, irho, gK); + bend_fringe(r6, irho, gKEntrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (k1_theta_entrance != 0.0) @@ -125,7 +127,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_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, gKExit); /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -168,7 +170,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(); @@ -202,7 +206,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->FringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gKEntrance=FullGap * FringeInt1; + Elem->gKExit=FullGap * FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -220,7 +225,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->gKEntrance, Elem->gKExit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, gamma, num_particles); @@ -258,7 +264,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(); @@ -283,7 +291,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, + FullGap*FringeInt1, FullGap*FringeInt1, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); } else if (nrhs == 0) { @@ -300,13 +309,15 @@ 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(16, 1); + plhs[1] = mxCreateCellMatrix(18, 1); mxSetCell(plhs[1], i1++, mxCreateString("Energy")); 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")); diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 0e590d564..6ccf24ab2 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -22,7 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gKEntrance; + double gKExit; double x0ref; double refdz; double *R1; @@ -40,7 +41,8 @@ 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 gKEntrance, double gKExit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -59,7 +61,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double k1_theta_entrance = 0.0; double k1_theta_exit = 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; } @@ -71,7 +73,7 @@ static void ExactRectangularBend(double *r, double le, double bending_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,\ + irho,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ LR,le,phi2) @@ -94,7 +96,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* edge focus */ if (FringeBendEntrance) - bend_fringe(r6, irho, gK); + bend_fringe(r6, irho, gKEntrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (k1_theta_entrance != 0.0) @@ -124,7 +126,7 @@ static void ExactRectangularBend(double *r, double le, double bending_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, gKExit); /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -165,7 +167,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(); @@ -195,7 +199,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->FringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gKEntrance=FullGap * FringeInt1; + Elem->gKExit=FullGap * FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -210,7 +215,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->gKEntrance, Elem->gKExit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -244,7 +250,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(); @@ -266,7 +274,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, + FullGap*FringeInt1, FullGap*FringeInt1, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { @@ -283,12 +292,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")); diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index e5dfbf76c..91795c421 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -23,7 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gKEntrance; + double gKExit; double x0ref; double refdz; double *R1; @@ -41,7 +42,8 @@ static void ExactRectangularBendRad(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 gKEntrance, double gKExit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -62,7 +64,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double k1_theta_entrance = 0.0; double k1_theta_exit = 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; } @@ -74,7 +76,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_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,\ + irho,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ LR,le,phi2,rad_const,diff_const) @@ -97,7 +99,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* edge focus */ if (FringeBendEntrance) - bend_fringe(r6, irho, gK); + bend_fringe(r6, irho, gKEntrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (k1_theta_entrance != 0.0) @@ -127,7 +129,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_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, gKExit); /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -170,7 +172,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(); @@ -204,7 +208,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->FringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gKEntrance=FullGap * FringeInt1; + Elem->gKExit=FullGap * FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -221,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->gKEntrance, Elem->gKExit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, gamma, num_particles); @@ -259,7 +265,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(); @@ -284,7 +292,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, + FullGap*FringeInt1, FullGap*FringeInt1, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); } else if (nrhs == 0) { @@ -301,13 +310,15 @@ 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(16, 1); + plhs[1] = mxCreateCellMatrix(18, 1); mxSetCell(plhs[1], i1++, mxCreateString("Energy")); 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")); diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 2cb979a27..ab938b96c 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -22,7 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gKEntrance; + double gKExit; double *R1; double *R2; double *T1; @@ -38,7 +39,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, - double gK, + double gKEntrance, double gKExit, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -55,7 +56,6 @@ static void ExactSectorBend(double *r, double le, double bending_angle, double k1_theta_entrance = 0.0; double k1_theta_exit = 0.0; - if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B0 = -sin(KickAngle[0]) / le; A0 = sin(KickAngle[1]) / le; @@ -67,7 +67,7 @@ static void ExactSectorBend(double *r, double le, double bending_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,\ + irho,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ entrance_angle,exit_angle,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,le) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gKEntrance=FullGap * FringeInt1; + Elem->gKExit=FullGap * FringeInt2; Elem->R1=R1; Elem->R2=R2; Elem->T1=T1; @@ -200,7 +203,7 @@ 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->gKEntrance, Elem->gKExit, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -234,7 +237,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 *R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error(); double *R2=atGetOptionalDoubleArray(ElemData,"R2"); check_error(); double *T1=atGetOptionalDoubleArray(ElemData,"T1"); check_error(); @@ -258,7 +263,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, + FullGap*FringeInt1, FullGap*FringeInt1, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { @@ -275,12 +280,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(13, 1); + plhs[1] = mxCreateCellMatrix(15, 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("T1")); mxSetCell(plhs[1], i1++, mxCreateString("T2")); mxSetCell(plhs[1], i1++, mxCreateString("R1")); diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 8366a509b..5c7615bec 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -23,7 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gKEntrance; + double gKExit; double *R1; double *R2; double *T1; @@ -39,7 +40,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, - double gK, + double gKEntrance, double gKExit, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -68,7 +69,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, #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,\ + irho,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ entrance_angle,exit_angle,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,le,\ rad_const, diff_const) @@ -88,7 +89,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, Yrot(r6, entrance_angle); if (FringeBendEntrance) - bend_fringe(r6, irho, gK); + bend_fringe(r6, irho, gKEntrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (k1_theta_entrance != 0.0) @@ -115,7 +116,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); if (FringeBendExit) - bend_fringe(r6, -irho, gK); + bend_fringe(r6, -irho, gKExit); Yrot(r6, exit_angle); /* Check physical apertures at the exit of the magnet */ @@ -154,7 +155,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 *R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error(); double *R2=atGetOptionalDoubleArray(ElemData,"R2"); check_error(); double *T1=atGetOptionalDoubleArray(ElemData,"T1"); check_error(); @@ -186,7 +189,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->FringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gKEntrance=FullGap * FringeInt1; + Elem->gKExit=FullGap * FringeInt2; Elem->R1=R1; Elem->R2=R2; Elem->T1=T1; @@ -203,7 +207,7 @@ 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->gKEntrance, Elem->gKExit, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, gamma, num_particles); @@ -242,7 +246,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 *R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error(); double *R2=atGetOptionalDoubleArray(ElemData,"R2"); check_error(); double *T1=atGetOptionalDoubleArray(ElemData,"T1"); check_error(); @@ -266,7 +272,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, + FullGap*FringeInt1, FullGap*FringeInt1, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); } else if (nrhs == 0) { @@ -283,13 +289,15 @@ 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(14, 1); + plhs[1] = mxCreateCellMatrix(16, 1); mxSetCell(plhs[1], i1++, mxCreateString("Energy")); 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("T1")); mxSetCell(plhs[1], i1++, mxCreateString("T2")); mxSetCell(plhs[1], i1++, mxCreateString("R1")); From 517770f53d020a132989eb36c3494a15e39aba27 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 23 Apr 2026 15:29:04 +0200 Subject: [PATCH 25/34] Processing of the fringe field integrals fint1 and fint2 --- atintegrators/ExactRectBendPass.c | 51 ++++++++++--------- atintegrators/ExactRectBendRadPass.c | 50 ++++++++++--------- atintegrators/ExactRectangularBendPass.c | 54 ++++++++++++--------- atintegrators/ExactRectangularBendRadPass.c | 50 ++++++++++--------- atintegrators/ExactSectorBendPass.c | 45 +++++++++-------- atintegrators/ExactSectorBendRadPass.c | 46 +++++++++--------- 6 files changed, 160 insertions(+), 136 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index f630b8329..6f69ec1de 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -22,8 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gKEntrance; - double gKExit; + double gK_entrance; + double gK_exit; double x0ref; double refdz; double *R1; @@ -41,7 +41,7 @@ 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 gKEntrance, double gKExit, + double gK_entrance, double gK_exit, double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, @@ -50,7 +50,9 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, { 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; @@ -58,24 +60,24 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = B[0]; double A0 = A[0]; - double k1_theta_entrance = 0.0; - double k1_theta_exit = 0.0; + double k1_entrance_angle = 0.0; + double k1_exit_angle = 0.0; 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_theta_entrance = B[1] * entrance_angle; - k1_theta_exit = B[1] * exit_angle; + 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,gKEntrance,gKExit,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; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gKEntrance=FullGap * FringeInt1; - Elem->gKExit=FullGap * FringeInt2; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -222,7 +226,7 @@ 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->gKEntrance, Elem->gKExit, + Elem->gK_entrance, Elem->gK_exit, Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, @@ -281,11 +285,12 @@ 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, - FullGap*FringeInt1, FullGap*FringeInt1, + FullGap*FringeInt1, FullGap*FringeInt2, x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index f9f4430a2..bfc81b3a5 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -23,8 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gKEntrance; - double gKExit; + double gK_entrance; + double gK_exit; double x0ref; double refdz; double *R1; @@ -42,7 +42,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, - double gKEntrance, double gKExit, + double gK_entrance, double gK_exit, double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, @@ -51,7 +51,9 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, { 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; @@ -61,24 +63,24 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double A0 = A[0]; double rad_const = RAD_CONST*pow(gamma, 3); double diff_const = DIF_CONST*pow(gamma, 5); - double k1_theta_entrance = 0.0; - double k1_theta_exit = 0.0; + double k1_entrance_angle = 0.0; + double k1_exit_angle = 0.0; 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_theta_entrance = B[1] * entrance_angle; - k1_theta_exit = B[1] * exit_angle; + 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,gKEntrance,gKExit,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,\ - le,phi2,rad_const,diff_const) + le,phi_entrance,phi_exit,rad_const,diff_const) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gKEntrance=FullGap * FringeInt1; - Elem->gKExit=FullGap * FringeInt2; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -225,7 +229,7 @@ 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->gKEntrance, Elem->gKExit, + Elem->gK_entrance, Elem->gK_exit, Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, @@ -291,7 +295,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - FullGap*FringeInt1, FullGap*FringeInt1, + FullGap*FringeInt1, FullGap*FringeInt2, x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 6ccf24ab2..59f56cfac 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -22,8 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gKEntrance; - double gKExit; + double gK_entrance; + double gK_exit; double x0ref; double refdz; double *R1; @@ -41,7 +41,7 @@ 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 gKEntrance, double gKExit, + double gK_entrance, double gK_exit, double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, @@ -50,7 +50,9 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, { 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; @@ -58,8 +60,8 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = B[0]; double A0 = A[0]; - double k1_theta_entrance = 0.0; - double k1_theta_exit = 0.0; + double k1_entrance_angle = 0.0; + double k1_exit_angle = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; @@ -67,16 +69,16 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, } B[0] += irho; if (max_order >= 1) { - k1_theta_entrance = B[1] * entrance_angle; - k1_theta_exit = B[1] * exit_angle; + 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,gKEntrance,gKExit,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; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gKEntrance=FullGap * FringeInt1; - Elem->gKExit=FullGap * FringeInt2; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -211,11 +215,12 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } - ExactRectangularBend(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, + ExactRectangularBend(r_in, Elem->Length, Elem->BendingAngle, + Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, Elem->FringeQuadEntrance, Elem->FringeQuadExit, - Elem->gKEntrance, Elem->gKExit, + Elem->gK_entrance, Elem->gK_exit, Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, @@ -270,11 +275,12 @@ 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, - FullGap*FringeInt1, FullGap*FringeInt1, + FullGap*FringeInt1, FullGap*FringeInt2, x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index 91795c421..02b95ecbc 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -23,8 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gKEntrance; - double gKExit; + double gK_entrance; + double gK_exit; double x0ref; double refdz; double *R1; @@ -42,7 +42,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, - double gKEntrance, double gKExit, + double gK_entrance, double gK_exit, double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, @@ -51,7 +51,9 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, { 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; @@ -61,8 +63,8 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, double A0 = A[0]; double rad_const = RAD_CONST*pow(gamma, 3); double diff_const = DIF_CONST*pow(gamma, 5); - double k1_theta_entrance = 0.0; - double k1_theta_exit = 0.0; + double k1_entrance_angle = 0.0; + double k1_exit_angle = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; @@ -70,16 +72,16 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, } B[0] += irho; if (max_order >= 1) { - k1_theta_entrance = B[1] * entrance_angle; - k1_theta_exit = B[1] * exit_angle; + 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,gKEntrance,gKExit,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,rad_const,diff_const) + le,phi_entrance, phi_exit,rad_const,diff_const) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gKEntrance=FullGap * FringeInt1; - Elem->gKExit=FullGap * FringeInt2; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -226,7 +230,7 @@ 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->gKEntrance, Elem->gKExit, + Elem->gK_entrance, Elem->gK_exit, Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, @@ -292,7 +296,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - FullGap*FringeInt1, FullGap*FringeInt1, + FullGap*FringeInt1, FullGap*FringeInt2, x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index ab938b96c..ea2759195 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -22,8 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gKEntrance; - double gKExit; + double gK_entrance; + double gK_exit; double *R1; double *R2; double *T1; @@ -39,7 +39,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, - double gKEntrance, double gKExit, + double gK_entrance, double gK_exit, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -53,22 +53,22 @@ static void ExactSectorBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = 0.0; double A0 = 0.0; - double k1_theta_entrance = 0.0; - double k1_theta_exit = 0.0; + double k1_entrance_angle = 0.0; + double k1_exit_angle = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B0 = -sin(KickAngle[0]) / le; A0 = sin(KickAngle[1]) / le; } if (max_order >= 1) { - k1_theta_entrance = B[1] * entrance_angle; - k1_theta_exit = B[1] * exit_angle; + 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,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ - entrance_angle,exit_angle,\ + 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,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,le) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gKEntrance=FullGap * FringeInt1; - Elem->gKExit=FullGap * FringeInt2; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->R1=R1; Elem->R2=R2; Elem->T1=T1; @@ -203,7 +205,7 @@ 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->gKEntrance, Elem->gKExit, + Elem->gK_entrance, Elem->gK_exit, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -259,11 +261,12 @@ 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]); + ExactSectorBend(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - FullGap*FringeInt1, FullGap*FringeInt1, + FullGap*FringeInt1, FullGap*FringeInt2, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 5c7615bec..e0ff6dfa6 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -23,8 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gKEntrance; - double gKExit; + double gK_entrance; + double gK_exit; double *R1; double *R2; double *T1; @@ -40,7 +40,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, - double gKEntrance, double gKExit, + double gK_entrance, double gK_exit, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -55,22 +55,22 @@ static void ExactSectorBendRad(double *r, double le, double irho, double diff_const = DIF_CONST*pow(gamma, 5); double B0 = 0.0; double A0 = 0.0; - double k1_theta_entrance = 0.0; - double k1_theta_exit = 0.0; + double k1_entrance_angle = 0.0; + double k1_exit_angle = 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B0 = -sin(KickAngle[0]) / le; A0 = sin(KickAngle[1]) / le; } if (max_order >= 1) { - k1_theta_entrance = B[1] * entrance_angle; - k1_theta_exit = B[1] * exit_angle; + 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,gKEntrance,gKExit,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ - entrance_angle,exit_angle,\ + 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,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,le,\ rad_const, diff_const) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gKEntrance=FullGap * FringeInt1; - Elem->gKExit=FullGap * FringeInt2; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->R1=R1; Elem->R2=R2; Elem->T1=T1; @@ -207,7 +209,7 @@ 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->gKEntrance, Elem->gKExit, + Elem->gK_entrance, Elem->gK_exit, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, gamma, num_particles); @@ -272,13 +274,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - FullGap*FringeInt1, FullGap*FringeInt1, + FullGap*FringeInt1, FullGap*FringeInt2, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); } else if (nrhs == 0) { /* list of required fields */ int i0 = 0; - plhs[0] = mxCreateCellMatrix(9, 1); + plhs[0] = mxCreateCellMatrix(8, 1); mxSetCell(plhs[0], i0++, mxCreateString("Length")); mxSetCell(plhs[0], i0++, mxCreateString("PolynomA")); mxSetCell(plhs[0], i0++, mxCreateString("PolynomB")); From c27cb90ed723ff2300d1c6da31de38a1930e62a0 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 29 Apr 2026 09:50:40 +0200 Subject: [PATCH 26/34] fix conversion of gaps --- pyat/at/load/xsuite.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 4ce97d0d8..17ce3fa8f 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -638,7 +638,7 @@ def _set_at_fringe(self) -> dict[str, Any]: msg = "Entry and Exit gaps for dipole are different, use entry" warnings.warn(AtWarning(msg), stacklevel=2) if entry_hgap is not None: - atparams["FullGap"] = entry_hgap + atparams["FullGap"] = 2.0 * entry_hgap if self.get("edge_entry_model", "linear") in ["linear", "full"]: atparams["FringeQuadEntrance"] = 1 @@ -648,8 +648,8 @@ def _set_at_fringe(self) -> dict[str, Any]: def _set_xs_fringe(self, atparams: dict): if (gap := atparams.get("FullGap")) is not None: - self["edge_entry_gap"] = gap - self["edge_exit_gap"] = gap + self["edge_entry_hgap"] = 0.5 * gap + self["edge_exit_hgap"] = 0.5 * gap exact = atparams.get("PassMethod", "").startswith("Exact") qentry = atparams.get("FringeQuadEntrance", 0) qexit = atparams.get("FringeQuadExit", 0) @@ -1051,9 +1051,7 @@ def from_at(cls, ring: Lattice, match_model: bool = False, **kwargs) -> XsLine: def refpart(rng): prt = rng.particle if prt.name == "relativistic": - with warnings.catch_warnings(): - warnings.simplefilter("ignore", category=UserWarning) - prt = Particle("electron") + prt = Particle("electron") gamma0 = rng.energy / prt.rest_energy beta0 = sqrt(1.0 - 1.0 / gamma0 / gamma0) return { From ed70746025a1a56b26f58dde7433acb68443a08f Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 6 May 2026 10:36:11 +0200 Subject: [PATCH 27/34] Handle AT Corrector, handle "suppress" edge model --- pyat/at/load/xsuite.py | 88 +++++++++++++++++++++++++++++++++--------- 1 file changed, 69 insertions(+), 19 deletions(-) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 17ce3fa8f..f13ad0021 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -157,10 +157,11 @@ import json import warnings +from abc import abstractmethod from math import sqrt from pathlib import Path from collections.abc import Callable -from typing import Any, ClassVar +from typing import Any, ClassVar, Protocol import contextlib import numpy as np @@ -225,7 +226,15 @@ def default(self, obj): return super().default(obj) -class XsElement(dict): +class _XsFactory(Protocol): + """Base class for Xsuite element factories.""" + + @classmethod + @abstractmethod + def from_at(cls, **atparams) -> XsElement: ... + + +class XsElement(dict, _XsFactory): """Base class for Xsuite elements.""" # Class attributes @@ -465,8 +474,8 @@ def xspoly(kmain: list[str], kerr: str) -> tuple[int, np.ndarray]: return porder, np.fromiter(poly_from_mad(poly), dtype=float, count=lpoly) length = self.get("length", 0.0) - xsorder = self.get("order", 0) - atorder = getattr(self._atClass, "DefaultOrder", xsorder) + xsorder: int = self.get("order", 0) + atorder: int = getattr(self._atClass, "DefaultOrder", xsorder) aorder, polya = xspoly(["k0s", "k1s", "k2s", "k3s"], "ksl") border, polyb = xspoly(["k0", "k1", "k2", "k3"], "knl") maxorder = max(aorder, border) @@ -625,12 +634,32 @@ class Bend(Multipole): "edge_exit_angle": "ExitAngle", } _mag_order = 0 - _edge_to_xs: ClassVar[dict[bool, dict[bool, str]]] = { - True: {True: "full", False: "dipole-only"}, - False: {True: "linear", False: "linear"}, + _default_bend_fringe: ClassVar[dict[bool, int]] = { + True: 4, + False: 1, + } + _at2xsuite_edge: ClassVar[dict[int, str]] = { + 0: "suppress", + 1: "linear", + 2: "linear", + 3: "linear", + 4: "dipole-only", + } + _xsuite2at_edge: ClassVar[dict[str, tuple]] = { + "suppress": (0, None), + "linear": (None, None), + "dipole-only": (None, None), + "full": (None, 1), } def _set_at_fringe(self) -> dict[str, Any]: + def edge_model(xskey, bendkey, quadkey): + [bend, quad] = self._at2xsuite_edge[xskey] + if bend is not None: + atparams[bendkey] = bend + if quad is not None: + atparams[quadkey] = quad + atparams = {} entry_hgap = self.get("edge_entry_hgap") exit_hgap = self.get("edge_exit_hgap") @@ -640,21 +669,27 @@ def _set_at_fringe(self) -> dict[str, Any]: if entry_hgap is not None: atparams["FullGap"] = 2.0 * entry_hgap - if self.get("edge_entry_model", "linear") in ["linear", "full"]: - atparams["FringeQuadEntrance"] = 1 - if self.get("edge_exit_model", "linear") in ["linear", "full"]: - atparams["FringeQuadExit"] = 1 + edge_model("edge_entry_model", "FringeBendEntrance", "FringeQuadEntrance") + edge_model("edge_exit_model", "FringeBendExit", "FringeQuadExit") return atparams def _set_xs_fringe(self, atparams: dict): + def edge_model(bendkey, quadkey): + at_quad_fringe = atparams.get(quadkey, 0) + if at_quad_fringe > 0: + return "full" + else: + return self._at2xsuite_edge[atparams.get(bendkey, default_bend_fringe)] + if (gap := atparams.get("FullGap")) is not None: self["edge_entry_hgap"] = 0.5 * gap self["edge_exit_hgap"] = 0.5 * gap exact = atparams.get("PassMethod", "").startswith("Exact") - qentry = atparams.get("FringeQuadEntrance", 0) - qexit = atparams.get("FringeQuadExit", 0) - self["edge_entry_model"] = self._edge_to_xs[exact][qentry > 0] - self["edge_exit_model"] = self._edge_to_xs[exact][qexit > 0] + default_bend_fringe = self._default_bend_fringe[exact] + self["edge_entry_model"] = edge_model( + "FringeBendEntrance", "FringeQuadEntrance" + ) + self["edge_exit_model"] = edge_model("FringeBendExit", "FringeQuadExit") def _params_to_at(self, **atparams) -> dict[str, Any]: atparams = super()._params_to_at(EntranceAngle=0.0, ExitAngle=0.0, **atparams) @@ -757,7 +792,7 @@ def from_at(cls, match_model: bool = False, **atparams): return elem -class NotInAT: +class NotInAT(_XsFactory): """Class for Xsuite elements without AT equivalent.""" @classmethod @@ -776,7 +811,7 @@ def from_dict( return xsclass.from_dict(xsparams, name=name, warn=warn) -class NotInXsuite: +class NotInXsuite(_XsFactory): """Class for AT elements without Xsuite equivalent.""" @classmethod @@ -788,7 +823,7 @@ def from_at(cls, **atparams) -> XsElement: return xsclass.from_at(**atparams) -class Dipole: +class Dipole(_XsFactory): """Class for handling AT dipoles.""" @classmethod @@ -802,6 +837,20 @@ def from_at(cls, **atparams) -> XsElement: return Bend.from_at(**atparams) +class Corrector(_XsFactory): + """Class for handling AT correctors.""" + + # noinspection PyPep8Naming + @classmethod + def from_at(cls, KickAngle=(0.0, 0.0), **atparams) -> XsElement: + pola = np.array([KickAngle[1]]) + polb = np.array([-KickAngle[0]]) + if (length := atparams["Length"]) != 0.0: + pola /= length + polb /= length + return Multipole.from_at(PolynomA=pola, PolynomB=polb, MaxOrder=0, **atparams) + + _xsclass: dict[str, type[XsElement]] = { "Marker": Marker, "Drift": Drift, @@ -815,7 +864,7 @@ def from_at(cls, **atparams) -> XsElement: } -_at2xsclass: dict[type[elt.Element], type[XsElement]] = { +_at2xsclass: dict[type[elt.Element], type[_XsFactory]] = { elt.Marker: Marker, elt.Monitor: Marker, elt.Drift: Drift, @@ -826,6 +875,7 @@ def from_at(cls, **atparams) -> XsElement: elt.Multipole: Multipole, elt.ThinMultipole: Multipole, elt.Dipole: Dipole, + elt.Corrector: Corrector, } From a591a33af1a3ea061bbfa49ee279fe1cd9101c3d Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 6 May 2026 11:08:22 +0200 Subject: [PATCH 28/34] bug fix --- pyat/at/load/xsuite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index f13ad0021..8efdb35e2 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -654,7 +654,7 @@ class Bend(Multipole): def _set_at_fringe(self) -> dict[str, Any]: def edge_model(xskey, bendkey, quadkey): - [bend, quad] = self._at2xsuite_edge[xskey] + [bend, quad] = self._xsuite2at_edge[self.get(xskey, "linear")] if bend is not None: atparams[bendkey] = bend if quad is not None: From 3285e736068cae900a9a4da98e09c5becc6a421f Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 6 May 2026 11:38:42 +0200 Subject: [PATCH 29/34] Apply x0ref to the whole magnet including edge effects --- atintegrators/ExactRectBendPass.c | 4 ++-- atintegrators/ExactRectBendRadPass.c | 4 ++-- atintegrators/ExactRectangularBendPass.c | 4 ++-- atintegrators/ExactRectangularBendRadPass.c | 4 ++-- pyat/at/load/xsuite.py | 2 +- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index 6f69ec1de..5e41c414a 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -87,6 +87,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); + r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -105,7 +106,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, bend_edge(r6, irho, phi_entrance); } - r6[0] += x0ref; if (num_int_steps == 0) { exact_straight_bend(r6, irho, LR); } @@ -120,7 +120,6 @@ 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); @@ -143,6 +142,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ + r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index bfc81b3a5..076b359e2 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -90,6 +90,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); + r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -108,7 +109,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, bend_edge(r6, irho, phi_entrance); } - r6[0] += x0ref; for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_straight_bend(r6, irho, L1); ex_strthinkickrad(r6, A, B, max_order, irho, K1, rad_const, diff_const, NULL); @@ -118,7 +118,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, ex_strthinkickrad(r6, A, B, max_order, irho, K1, rad_const, diff_const, NULL); exact_straight_bend(r6, irho, L1); } - r6[0] -= x0ref; /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); @@ -141,6 +140,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ + r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 59f56cfac..e9e247c08 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -88,6 +88,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); + r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -107,7 +108,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, } /* integrator */ - r6[0] += x0ref; for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_drift(r6, L1); strthinkick(r6, A, B, K1, max_order); @@ -117,7 +117,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, strthinkick(r6, A, B, K1, max_order); exact_drift(r6, L1); } - r6[0] -= x0ref; /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); @@ -140,6 +139,7 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ + r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index 02b95ecbc..bed02ffc4 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -91,6 +91,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); + r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -110,7 +111,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, } /* integrator */ - r6[0] += x0ref; for (int m = 0; m < num_int_steps; m++) { /* Loop over slices */ exact_drift(r6, L1); ex_strthinkickrad(r6, A, B, max_order, 0.0, K1, rad_const, diff_const, NULL); @@ -120,7 +120,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, ex_strthinkickrad(r6, A, B, max_order, 0.0, K1, rad_const, diff_const, NULL); exact_drift(r6, L1); } - r6[0] -= x0ref; /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); @@ -143,6 +142,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ + r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 8efdb35e2..30e7e1650 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -647,7 +647,7 @@ class Bend(Multipole): } _xsuite2at_edge: ClassVar[dict[str, tuple]] = { "suppress": (0, None), - "linear": (None, None), + "linear": (1, None), "dipole-only": (None, None), "full": (None, 1), } From 9e523fb188d4a6a81e4b5f832a2b459fc49a8ca2 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 6 May 2026 14:47:56 +0200 Subject: [PATCH 30/34] correct edge_model "suppressed" --- pyat/at/load/xsuite.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 30e7e1650..b35c08487 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -639,14 +639,14 @@ class Bend(Multipole): False: 1, } _at2xsuite_edge: ClassVar[dict[int, str]] = { - 0: "suppress", + 0: "suppressed", 1: "linear", 2: "linear", 3: "linear", 4: "dipole-only", } _xsuite2at_edge: ClassVar[dict[str, tuple]] = { - "suppress": (0, None), + "suppressed": (0, None), "linear": (1, None), "dipole-only": (None, None), "full": (None, 1), From e1aadcf06f539960ede6ddd578499a5006ca6978 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 6 May 2026 21:11:43 +0200 Subject: [PATCH 31/34] change rectangular magnet shift --- atintegrators/ExactRectBendPass.c | 12 +- atintegrators/ExactRectBendRadPass.c | 12 +- atintegrators/ExactRectangularBendPass.c | 12 +- atintegrators/ExactRectangularBendRadPass.c | 12 +- atintegrators/ExactSectorBendPass.c | 4 +- atintegrators/ExactSectorBendRadPass.c | 4 +- atintegrators/exactbendfringe.c | 155 ++++++++++++-------- 7 files changed, 124 insertions(+), 87 deletions(-) diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index 5e41c414a..2ae9a17f0 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -87,7 +87,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); - r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -96,13 +95,14 @@ 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_entrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (phi_entrance != 0.0) { - if (k1_entrance_angle != 0.0) quad_wedge(r6, -k1_entrance_angle); + if (k1_entrance_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_entrance_angle); bend_edge(r6, irho, phi_entrance); } @@ -124,15 +124,16 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); - /* edge focus */ + /* Exit face */ if (phi_exit != 0.0) { bend_edge(r6, irho, phi_exit); - if (k1_exit_angle != 0.0) quad_wedge(r6, -k1_exit_angle); + 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_exit); + r6[0] -= x0ref; /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -142,7 +143,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ - r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 076b359e2..7ffa14c1e 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -90,7 +90,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); - r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -99,12 +98,13 @@ static void ExactRectangularBendRad(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_entrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); - if (phi_entrance != 0.0) { + if (phi_entrance != 0.0 && FringeQuadEntrance) { if (k1_entrance_angle != 0.0) quad_wedge(r6, -k1_entrance_angle); bend_edge(r6, irho, phi_entrance); } @@ -122,15 +122,16 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); - /* edge focus */ + /* Exit face */ if (phi_exit != 0.0) { bend_edge(r6, irho, phi_exit); - if (k1_exit_angle != 0.0) quad_wedge(r6, -k1_exit_angle); + 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_exit); + r6[0] -= x0ref; /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -140,7 +141,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ - r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index e9e247c08..31b1a07c5 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -88,7 +88,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); - r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -97,13 +96,14 @@ 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_entrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (phi_entrance != 0.0) { - if (k1_entrance_angle != 0.0) quad_wedge(r6, -k1_entrance_angle); + if (k1_entrance_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_entrance_angle); bend_edge(r6, irho, phi_entrance); } @@ -121,15 +121,16 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); - /* edge focus */ + /* Exit face */ if (phi_exit != 0.0) { bend_edge(r6, irho, phi_exit); - if (k1_exit_angle != 0.0) quad_wedge(r6, -k1_exit_angle); + 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_exit); + r6[0] -= x0ref; /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -139,7 +140,6 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ - r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index bed02ffc4..2e8e680d0 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -91,7 +91,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* misalignment at entrance */ if (T1) ATaddvv(r6,T1); if (R1) ATmultmv(r6,R1); - r6[0] += x0ref; /* Change to the magnet referential */ Yrot(r6, entrance_angle); @@ -100,13 +99,14 @@ static void ExactRectangularBendRad(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_entrance); if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (phi_entrance != 0.0) { - if (k1_entrance_angle != 0.0) quad_wedge(r6, -k1_entrance_angle); + if (k1_entrance_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_entrance_angle); bend_edge(r6, irho, phi_entrance); } @@ -124,15 +124,16 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, /* Convert absolute path length to path lengthening */ r6[5] -= (le+refdz); - /* edge focus */ + /* Exit face */ if (phi_exit != 0.0) { bend_edge(r6, irho, phi_exit); - if (k1_exit_angle != 0.0) quad_wedge(r6, -k1_exit_angle); + 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_exit); + r6[0] -= x0ref; /* Check physical apertures at the exit of the magnet */ if (RApertures) checkiflostRectangularAp(r6, RApertures); @@ -142,7 +143,6 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, Yrot(r6, exit_angle); /* Misalignment at exit */ - r6[0] -= x0ref; if (R2) ATmultmv(r6,R2); if (T2) ATaddvv(r6,T2); diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index ea2759195..c076d7728 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -90,7 +90,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (entrance_angle != 0.0) { - if (k1_entrance_angle != 0.0) quad_wedge(r6, -k1_entrance_angle); + if (k1_entrance_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_entrance_angle); bend_edge(r6, irho, -entrance_angle); } @@ -115,7 +115,7 @@ static void ExactSectorBend(double *r, double le, double bending_angle, /* edge focus */ if (exit_angle != 0.0) { bend_edge(r6, irho, -exit_angle); - if (k1_exit_angle != 0.0) quad_wedge(r6, -k1_exit_angle); + if (k1_exit_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_exit_angle); } if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 1); diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index e0ff6dfa6..c16a3de42 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -93,7 +93,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 1); if (entrance_angle != 0.0) { - if (k1_entrance_angle != 0.0) quad_wedge(r6, -k1_entrance_angle); + if (k1_entrance_angle != 0.0 && FringeQuadEntrance) quad_wedge(r6, -k1_entrance_angle); bend_edge(r6, irho, -entrance_angle); } @@ -113,7 +113,7 @@ static void ExactSectorBendRad(double *r, double le, double irho, /* edge focus */ if (exit_angle != 0.0) { bend_edge(r6, irho, -exit_angle); - if (k1_exit_angle != 0.0) quad_wedge(r6, -k1_exit_angle); + 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); diff --git a/atintegrators/exactbendfringe.c b/atintegrators/exactbendfringe.c index b4469efd6..f0c33bb95 100644 --- a/atintegrators/exactbendfringe.c +++ b/atintegrators/exactbendfringe.c @@ -18,22 +18,25 @@ static void Yrot(double *r6, double phi) if (phi != 0.0) { double dp1 = 1.0 + r6[delta_]; + double x = r6[x_]; + double px = r6[px_]; + double py = r6[py_]; double c = cos(phi); double s = sin(phi); - double pz = pxyz(dp1, r6[px_], r6[py_]); - double p = c*pz - s*r6[px_]; - double px = s*pz + c*r6[px_]; - double x = r6[x_]*pz/p; - double dy = r6[x_]*r6[py_]*s/p; - double dct = dp1*r6[x_]*s/p; - r6[x_] = x; - r6[px_] = px; + double pz = pxyz(dp1, px, py); + double p = c*pz - s*px; + double new_px = s*pz + c*px; + double new_x = x*pz/p; + double dy = x*py*s/p; + double dct = dp1*x*s/p; + r6[x_] = new_x; + r6[px_] = new_px; r6[y_] += dy; r6[ct_] += dct; } } -static void bend_fringe(double *r6, double irho, double gK) +static void bend_fringe_orig(double *r6, double irho, double gK) { /* Forest 13.13, bend fringe in the hard-edge limit */ @@ -101,21 +104,27 @@ static void bend_edge(double *r6, double rhoinv, double theta) double dp1 = 1.0 + r6[4]; double c = cos(theta); double s = sin(theta); - double pz = pxyz(dp1, r6[px_], r6[py_]); - double d2 = pxyz(dp1, 0.0, r6[py_]); - double px = r6[px_]*c + (pz - rhoinv*r6[x_])*s; - double dasin = asin(r6[px_]/d2) - asin(px/d2); - double num = r6[x_]*(r6[px_]*sin(2.0*theta) + s*s*(2.0*pz - rhoinv*r6[x_])); - double den = pxyz(dp1, px, r6[py_]) + pxyz(dp1, r6[px_], r6[py_])*c - r6[px_]*s; - double x = r6[x_]*c + num/den; - double dy = r6[py_]*theta/rhoinv + r6[py_]/rhoinv*dasin; + double x = r6[x_]; + double px = r6[px_]; + double py = r6[py_]; + double pz = pxyz(dp1, px, py); + double d2 = pxyz(dp1, 0.0, py); + double new_px = px*c + (pz - rhoinv*x)*s; + double dasin = asin(px/d2) - asin(new_px/d2); + double num = x*s*(2.0*px*c + s*(2.0*pz - rhoinv*x)); + double den = pxyz(dp1, new_px, py) + pz*c - px*s; + double new_x = x*c + num/den; + double dy = py*theta/rhoinv + py/rhoinv*dasin; double dct = dp1/rhoinv*(theta + dasin); - r6[x_] = x; - r6[px_] = px; + r6[x_] = new_x; + r6[px_] = new_px; r6[y_] += dy; r6[ct_] += dct; } + else { + Yrot(r6, theta); + } } static void quad_wedge(double *r6, double k1_theta) @@ -128,45 +137,73 @@ static void quad_wedge(double *r6, double k1_theta) r6[py_] += dpy; } -static void bend_fringe_test(double *r6, double irho, double gK) +void bend_fringe(double *r6, double irho, double gK) { - /* Forest 13.13, bend fringe in the hard-edge limit */ - - double b0 = irho; - - double pz = pxyz(1.0+r6[delta_], r6[px_], r6[py_]); - double px = r6[px_]; - double py = r6[py_]; - double d = r6[delta_]; - double xp = px / pz; - double yp = py / pz; - - double psi = b0*xp/(1+yp*yp); - - /* these are the partial derivatives of psi with respect to px, py and delta - total horror from Mathematica. This could benefit from some mini-TPSA */ - - double px2 = px*px; - double py2 = py*py; - double pz2 = pz*pz; - double pz4 = pz2*pz2; - double denom = (pz2 + py2); - - double dpx = b0*(pz4 + pz2*denom - px2*py2)/pz/denom/denom; - double dpy = -b0*yp*px/denom; - double dd = b0*xp*(1.0+d)*(py2-pz2)/denom/denom; - - /* solve quadratic equation in yf (Forest fringe_part_I.pdf) */ - - double yf = (2 * r6[y_]) / (1 + sqrt(1 - 2 * dpy * r6[y_])); - double dxf = 0.5 * dpx * yf * yf; - double dct = 0.5 * dd * yf * yf; - double dpyf = psi * yf; - // atPrintf("Fringe dx, dpy, dct: %g, %g, %g\n", dxf, dpyf, dct); - - - r6[x_] += dxf; - r6[y_] = yf; - r6[py_] -= dpyf; - r6[ct_] -= dct; + double factor = 0.0; + if (fabs(gK) > 1.0e-6) { + factor = SQR(irho) / 9.0 / gK; + } + const double irho_g_fint = irho * gK; + const double y = r6[y_]; + const double px = r6[px_]; + const double py = r6[py_]; + const double dp1 = r6[delta_] + 1.0; + + const double pz2 = SQR(dp1) - SQR(px) - SQR(py); + const double pz = sqrt(pz2); + const double xp = px / pz; + const double yp = py / pz; + + const double dpz_dpx = -xp; + const double dpz_dpy = -yp; + const double dpz_ddelta = dp1 / pz; + + const double dxp_dpx = -px/pz2 * dpz_dpx + 1/pz; + const double dxp_dpy = -px/pz2 * dpz_dpy; + const double dxp_ddelta = -px/pz2 * dpz_ddelta; + + const double dyp_dpx = -py/pz2 * dpz_dpx; + const double dyp_dpy = -py/pz2 * dpz_dpy + 1/pz; + const double dyp_ddelta = -py/pz2 * dpz_ddelta; + + const double phi0 = xp / (1.0 + SQR(yp)); + const double dphi0_dxp = 1.0 / (1.0 + SQR(yp)); + const double dphi0_dyp = -2 * xp * yp / SQR(1.0 + SQR(yp)); + + const double dphi0_dpx = dphi0_dxp * dxp_dpx + dphi0_dyp * dyp_dpx; + const double dphi0_dpy = dphi0_dxp * dxp_dpy + dphi0_dyp * dyp_dpy; + const double dphi0_ddelta = dphi0_dxp * dxp_ddelta + dphi0_dyp * dyp_ddelta; + + const double phi1 = 1.0 + 2.0 * SQR(xp) + SQR(xp) * SQR(yp); + const double dphi1_dxp = 4.0 * xp + 2.0 * SQR(yp) * xp; + const double dphi1_dyp = 2.0 * SQR(xp) * yp; + + const double dphi1_dpx = dphi1_dxp * dxp_dpx + dphi1_dyp * dyp_dpx; + const double dphi1_dpy = dphi1_dxp * dxp_dpy + dphi1_dyp * dyp_dpy; + const double dphi1_ddelta = dphi1_dxp * dxp_ddelta + dphi1_dyp * dyp_ddelta; + + const double phi2 = atan(phi0) - irho_g_fint * pz * phi1; + const double dphi2_dpx = dphi0_dpx / (1.0 + SQR(phi0)) + - irho_g_fint * (pz * dphi1_dpx + phi1 * dpz_dpx); + const double dphi2_dpy = dphi0_dpy / (1.0 + SQR(phi0)) + - irho_g_fint * (pz * dphi1_dpy + phi1 * dpz_dpy); + const double dphi2_ddelta = dphi0_ddelta / (1.0 + SQR(phi0)) + - irho_g_fint * (pz * dphi1_ddelta + phi1 * dpz_ddelta); + + const double Phi0 = irho * tan(phi2); + const double cphi2 = cos(phi2); + const double irho_c2 = irho / cphi2 / cphi2; + const double dPhi0_dpx = irho_c2 * dphi2_dpx; + const double dPhi0_dpy = irho_c2 * dphi2_dpy; + const double dPhi0_ddelta = irho_c2 * dphi2_ddelta; + + const double new_y = 2.0 * y / (1.0 + sqrt(1.0 - 2.0 * dPhi0_dpy * y)); + const double delta_x = dPhi0_dpx * SQR(new_y) / 2; + const double delta_py = -Phi0 * new_y - factor / dp1 * SQR(new_y) * new_y; + const double delta_l = -dPhi0_ddelta * SQR(new_y) / 2.0; + + r6[x_] += delta_x; + r6[y_] = new_y; + r6[py_] += delta_py; + r6[ct_] += delta_l; } From 151ef3a68fe6a9ad5316996a523355344c8c085e Mon Sep 17 00:00:00 2001 From: AmarildoTopalli Date: Wed, 13 May 2026 10:57:40 +0200 Subject: [PATCH 32/34] fix(xsuite): remove duplicate keyword args in Corrector.from_at (#1078) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit atelem.to_file() serializes all element attributes including PolynomA, PolynomB and MaxOrder into atparams. These were then passed twice to Multipole.from_at() — once explicitly and once via **atparams — causing a TypeError. Fix by popping PolynomA, PolynomB and MaxOrder from atparams before the call, so only the values recomputed from KickAngle are used. --- pyat/at/load/xsuite.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index b35c08487..5c3a1bc2f 100644 --- a/pyat/at/load/xsuite.py +++ b/pyat/at/load/xsuite.py @@ -848,6 +848,13 @@ def from_at(cls, KickAngle=(0.0, 0.0), **atparams) -> XsElement: if (length := atparams["Length"]) != 0.0: pola /= length polb /= length + + #################################### + atparams.pop("PolynomA", None) + atparams.pop("PolynomB", None) + atparams.pop("MaxOrder", None) + ##################################### + return Multipole.from_at(PolynomA=pola, PolynomB=polb, MaxOrder=0, **atparams) From ce52a966412b7d68b7248395eddb32898224e3da Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Fri, 15 May 2026 15:25:53 +0200 Subject: [PATCH 33/34] rbend_tune inactive if no multipole --- atintegrators/ExactRectangularBendPass.c | 15 +++++---------- pyat/at/lattice/elements/rectangular_bend.py | 7 +++---- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 31b1a07c5..1316259c6 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -60,23 +60,18 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = B[0]; double A0 = A[0]; - double k1_entrance_angle = 0.0; - double k1_exit_angle = 0.0; + double B1 = (max_order >= 1) ? B[1] : 0.0; if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; A[0] += sin(KickAngle[1])/le; } B[0] += irho; - 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_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,\ + irho,gK_entrance,gK_exit,A0,B0,B1,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ + entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ LR,le,phi_entrance,phi_exit) for (int c = 0; c Date: Fri, 15 May 2026 15:58:21 +0200 Subject: [PATCH 34/34] rbend_tune inactive if no multipole --- pyat/at/lattice/elements/rectangular_bend.py | 37 +++++++++----------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/pyat/at/lattice/elements/rectangular_bend.py b/pyat/at/lattice/elements/rectangular_bend.py index ed23d22c0..82f8258e7 100644 --- a/pyat/at/lattice/elements/rectangular_bend.py +++ b/pyat/at/lattice/elements/rectangular_bend.py @@ -1,4 +1,4 @@ -"""Additional method for rectangular bending magnets""" +"""Additional method for rectangular bending magnets.""" from math import sin, cos @@ -12,7 +12,7 @@ def rbendtune(self: Dipole) -> None: # noinspection PyUnresolvedReferences - r"""Set *X0ref* and *RefDZ* for rectangular bending magnets + r"""Set *X0ref* and *RefDZ* for rectangular bending magnets. This method must be called after creating a rectangular bending magnet or after setting its *PolynomA/B* attributes. It will set the correct *X0ref* @@ -31,33 +31,30 @@ def rbendtune(self: Dipole) -> None: """ def cross(x0r: float): - """Return the horizontal exit angle of the reference particle""" + """Return the horizontal exit angle of the reference particle.""" elem.X0ref = x0r out = elem.track(np.zeros(6)) return out[1] - def checkmul(el): - """Check if there are multipoles""" - for order in range(el.MaxOrder + 1): - if el.PolynomB[order] != 0.0: - return True - return False - passmethod = self.PassMethod.replace("RadPass", "Pass") if passmethod in { "BndStrMPoleSymplectic4Pass", "ExactRectangularBendPass", "ExactRectBendPass", - } and checkmul(self): - elem = self.copy() - elem.PassMethod = passmethod - theta = elem.BendingAngle - - # Analytical estimate - x0ref = elem.Length * ((cos(0.5 * theta) - 1.0) / theta + sin(0.5 * theta) / 12) - - # cancel output angle - x0ref = float(fsolve(cross, x0ref)) + }: + # check if there are multipoles + if any(self.PolynomB[order] != 0.0 for order in range(self.MaxOrder + 1)): + elem = self.copy() + elem.PassMethod = passmethod + tta = elem.BendingAngle + + # Analytical estimate + x0ref = elem.Length * ((cos(0.5 * tta) - 1.0) / tta + sin(0.5 * tta) / 12) + + # cancel output angle + x0ref = float(fsolve(cross, x0ref)) + else: + x0ref = 0.0 self.X0ref = x0ref rout = self.track(np.zeros(6))