diff --git a/atintegrators/BndMPoleSymplectic4Pass.c b/atintegrators/BndMPoleSymplectic4Pass.c index 1aef9d8e2..634cfefb8 100644 --- a/atintegrators/BndMPoleSymplectic4Pass.c +++ b/atintegrators/BndMPoleSymplectic4Pass.c @@ -2,7 +2,8 @@ #include "atelem.c" #include "atlalib.c" #include "atphyslib.c" -#include "driftkick.c" /* fastdrift and bndthinkick */ +#include "drift_fast.h" +#include "kick_h_k0h_k1h_kn.h" #include "quadfringe.c" /* QuadFringePassP, QuadFringePassN */ struct elem @@ -55,12 +56,12 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl double K2 = SL*KICK2; bool useLinFrEleEntrance = (fringeIntM0 != NULL && fringeIntP0 != NULL && FringeQuadEntrance==2); bool useLinFrEleExit = (fringeIntM0 != NULL && fringeIntP0 != NULL && FringeQuadExit==2); - double B0 = B[0]; - double A0 = A[0]; + double B0 = 0.0; + double A0 = 0.0; - if (KickAngle) { /* Convert corrector component to polynomial coefficients */ - B[0] -= sin(KickAngle[0])/le; - A[0] += sin(KickAngle[1])/le; + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + B0 = -sin(KickAngle[0]) / le; + A0 = sin(KickAngle[1]) / le; } #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \ @@ -95,13 +96,13 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl } /* integrator */ for (m=0; m < num_int_steps; m++) { /* Loop over slices */ - fastdrift(r6, NormL1); - bndthinkick(r6, A, B, K1, irho, max_order); - fastdrift(r6, NormL2); - bndthinkick(r6, A, B, K2, irho, max_order); - fastdrift(r6, NormL2); - bndthinkick(r6, A, B, K1, irho, max_order); - fastdrift(r6, NormL1); + drift(r6, NormL1); + kick(r6, A0, B0, A, B, max_order, K1, irho); + drift(r6, NormL2); + kick(r6, A0, B0, A, B, max_order, K2, irho); + drift(r6, NormL2); + kick(r6, A0, B0, A, B, max_order, K1, irho); + drift(r6, NormL1); } /* quadrupole gradient fringe */ if (FringeQuadExit && B[1]!=0) { @@ -122,10 +123,6 @@ void BndMPoleSymplectic4Pass(double *r, double le, double irho, double *A, doubl if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling); } } - if (KickAngle) { /* Remove corrector component in polynomial coefficients */ - B[0] = B0; - A[0] = A0; - } } #if defined(MATLAB_MEX_FILE) || defined(PYAT) diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index 897906ff9..1648c3472 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -1,8 +1,8 @@ #include "atelem.c" #include "atlalib.c" -#include "diff_bend_fringe.c" -#include "diff_bnd_kick.c" #include "diff_drift.c" +#include "diff_h_k0h_k1h_kn.c" +#include "diff_bend_fringe.c" #include "quadfringe.c" /* QuadFringePassP, QuadFringePassN */ struct elem @@ -57,15 +57,16 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do double K2 = SL*KICK2; bool useLinFrEleEntrance = (fringeIntM0 != NULL && fringeIntP0 != NULL && FringeQuadEntrance==2); bool useLinFrEleExit = (fringeIntM0 != NULL && fringeIntP0 != NULL && FringeQuadExit==2); - double B0 = B[0]; - double A0 = A[0]; double rad_const = RAD_CONST*pow(gamma, 3); double diff_const = DIF_CONST*pow(gamma, 5); + double B0 = 0.0; + double A0 = 0.0; - if (KickAngle) { /* Convert corrector component to polynomial coefficients */ - B[0] -= sin(KickAngle[0])/le; - A[0] += sin(KickAngle[1])/le; + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + B0 = -sin(KickAngle[0]) / le; + A0 = sin(KickAngle[1]) / le; } + #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \ shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,bdiff,\ irho,gap,A,B,L1,L2,K1,K2,max_order,num_int_steps,rad_const, diff_const,scaling,\ @@ -95,11 +96,11 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do /* integrator */ for (m=0; m < num_int_steps; m++) { /* Loop over slices */ diff_drift(r6, L1, bdiff); - diff_bnd_kick(r6, A, B, max_order, K1, irho, rad_const, diff_const, bdiff); + diff_kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, bdiff); diff_drift(r6, L2, bdiff); - diff_bnd_kick(r6, A, B, max_order, K2, irho, rad_const, diff_const, bdiff); + diff_kick(r6, A0, B0, A, B, max_order, K2, irho, rad_const, diff_const, bdiff); diff_drift(r6, L2, bdiff); - diff_bnd_kick(r6, A, B, max_order, K1, irho, rad_const, diff_const, bdiff); + diff_kick(r6, A0, B0, A, B, max_order, K1, irho, rad_const, diff_const, bdiff); diff_drift(r6, L1, bdiff); } /* quadrupole gradient fringe */ @@ -121,10 +122,6 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling); } } - if (KickAngle) { /* Remove corrector component in polynomial coefficients */ - B[0] = B0; - A[0] = A0; - } } #if defined(MATLAB_MEX_FILE) || defined(PYAT) 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/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index f8ab7a859..2ae9a17f0 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -22,7 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gK_entrance; + double gK_exit; double x0ref; double refdz; double *R1; @@ -40,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 gK_entrance, double gK_exit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -48,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; @@ -56,18 +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_entrance_angle = 0.0; + double k1_exit_angle = 0.0; - if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; A[0] += sin(KickAngle[1])/le; } + if (max_order >= 1) { + k1_entrance_angle = B[1] * entrance_angle; + k1_exit_angle = B[1] * exit_angle; + } #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \ shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,\ - irho,gK,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ - entrance_angle,exit_angle,x0ref,refdz,\ + irho,gK_entrance,gK_exit,A0,B0,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ + k1_entrance_angle,k1_exit_angle,entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ - LR,le,phi2) + LR,le,phi_entrance,phi_exit) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -207,7 +226,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, Elem->FringeQuadEntrance, Elem->FringeQuadExit, - Elem->gK,Elem->x0ref,Elem->refdz, + Elem->gK_entrance, Elem->gK_exit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -241,7 +261,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) int FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); - double gK=atGetOptionalDouble(ElemData,"gK", 0.0); check_error(); + double FullGap=atGetOptionalDouble(ElemData,"FullGap",0.0); check_error(); + double FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1",0.0); check_error(); + double FringeInt2=atGetOptionalDouble(ElemData,"FringeInt2",0.0); check_error(); double x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error(); double refdz=atGetOptionalDouble(ElemData,"RefDZ", 0.0); check_error(); double *R1=atGetOptionalDoubleArray(ElemData,"R1"); check_error(); @@ -263,11 +285,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + ExactRectangularBend(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, + FullGap*FringeInt1, FullGap*FringeInt2, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { @@ -284,12 +308,14 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mxSetCell(plhs[0], i0++, mxCreateString("ExitAngle")); if (nlhs>1) { /* list of optional fields */ int i1 = 0; - plhs[1] = mxCreateCellMatrix(15, 1); + plhs[1] = mxCreateCellMatrix(17, 1); mxSetCell(plhs[1], i1++, mxCreateString("FringeBendEntrance")); mxSetCell(plhs[1], i1++, mxCreateString("FringeBendExit")); mxSetCell(plhs[1], i1++, mxCreateString("FringeQuadEntrance")); mxSetCell(plhs[1], i1++, mxCreateString("FringeQuadExit")); - mxSetCell(plhs[1], i1++, mxCreateString("gK")); + mxSetCell(plhs[1], i1++, mxCreateString("FullGap")); + mxSetCell(plhs[1], i1++, mxCreateString("FringeInt1")); + mxSetCell(plhs[1], i1++, mxCreateString("FringeInt2")); mxSetCell(plhs[1], i1++, mxCreateString("X0ref")); mxSetCell(plhs[1], i1++, mxCreateString("RefDZ")); mxSetCell(plhs[1], i1++, mxCreateString("T1")); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 203385f0d..7ffa14c1e 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -23,7 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gK_entrance; + double gK_exit; 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 gK_entrance, double gK_exit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -49,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; @@ -59,18 +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_entrance_angle = 0.0; + double k1_exit_angle = 0.0; - if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; A[0] += sin(KickAngle[1])/le; } + if (max_order >= 1) { + k1_entrance_angle = B[1] * entrance_angle; + k1_exit_angle = B[1] * exit_angle; + } #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \ shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,\ - irho,gK,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ - entrance_angle,exit_angle,x0ref,refdz,\ + irho,gK_entrance,gK_exit,A0,B0,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ + k1_entrance_angle,k1_exit_angle,entrance_angle,exit_angle,x0ref,refdz,\ FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ - 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->gK=gK; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -210,7 +229,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, Elem->FringeQuadEntrance, Elem->FringeQuadExit, - Elem->gK,Elem->x0ref,Elem->refdz, + Elem->gK_entrance, Elem->gK_exit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, gamma, num_particles); @@ -248,7 +268,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(); @@ -273,7 +295,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*FringeInt2, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); } else if (nrhs == 0) { @@ -290,13 +313,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 a91fa447b..1316259c6 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -22,7 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gK_entrance; + double gK_exit; double x0ref; double refdz; double *R1; @@ -40,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 gK_entrance, double gK_exit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -48,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; @@ -56,8 +60,9 @@ static void ExactRectangularBend(double *r, double le, double bending_angle, double K2 = SL*KICK2; double B0 = B[0]; double A0 = A[0]; + double B1 = (max_order >= 1) ? B[1] : 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; } @@ -65,10 +70,10 @@ 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,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,phi2) + LR,le,phi_entrance,phi_exit) for (int c = 0; cFringeBendExit=FringeBendExit; Elem->FringeQuadEntrance=FringeQuadEntrance; Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -196,11 +210,13 @@ 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->gK,Elem->x0ref,Elem->refdz, + Elem->gK_entrance, Elem->gK_exit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -234,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(); @@ -252,11 +270,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + ExactRectangularBend(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, + FullGap*FringeInt1, FullGap*FringeInt2, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { @@ -273,12 +293,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 66f2d9d61..2e8e680d0 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -23,7 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gK_entrance; + double gK_exit; 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 gK_entrance, double gK_exit, + double x0ref, double refdz, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -49,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; @@ -59,19 +63,25 @@ 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_entrance_angle = 0.0; + double k1_exit_angle = 0.0; - if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ B[0] -= sin(KickAngle[0])/le; A[0] += sin(KickAngle[1])/le; } 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,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->gK=gK; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->x0ref=x0ref; Elem->refdz=refdz; Elem->R1=R1; @@ -211,7 +230,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, Elem->FringeQuadEntrance, Elem->FringeQuadExit, - Elem->gK,Elem->x0ref,Elem->refdz, + Elem->gK_entrance, Elem->gK_exit, + Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, gamma, num_particles); @@ -249,7 +269,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(); @@ -274,7 +296,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*FringeInt2, + x0ref, refdz, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Gamma, num_particles); } else if (nrhs == 0) { @@ -291,13 +314,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 44b910834..c076d7728 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -1,8 +1,8 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "driftkick.c" /* strthinkick.c */ -#include "exactbend.c" +#include "drift_exactbend.h" +#include "kick_k1h_kn.h" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -22,7 +22,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gK_entrance; + double gK_exit; 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 gK_entrance, double gK_exit, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -50,18 +51,24 @@ 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; + 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 (KickAngle) { /* Convert corrector component to polynomial coefficients */ + B0 = -sin(KickAngle[0]) / le; + A0 = sin(KickAngle[1]) / le; + } + if (max_order >= 1) { + k1_entrance_angle = B[1] * entrance_angle; + k1_exit_angle = B[1] * exit_angle; } #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \ shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,\ - irho,gK,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ - entrance_angle,exit_angle,\ + 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->gK=gK; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->R1=R1; Elem->R2=R2; Elem->T1=T1; @@ -192,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->gK, + Elem->gK_entrance, Elem->gK_exit, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, num_particles); @@ -226,7 +239,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(); @@ -246,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, - gK, + FullGap*FringeInt1, FullGap*FringeInt2, T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, num_particles); } else if (nrhs == 0) { @@ -267,12 +283,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 d4ed0b06b..c16a3de42 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -1,8 +1,8 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "exactkickrad.c" -#include "exactbend.c" +#include "drift_exactbend.h" +#include "kickrad_k1h_kn.h" #include "exactbendfringe.c" #include "exactmultipolefringe.c" @@ -23,7 +23,8 @@ struct elem int FringeBendExit; int FringeQuadEntrance; int FringeQuadExit; - double gK; + double gK_entrance; + double gK_exit; 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 gK_entrance, double gK_exit, double *T1, double *T2, double *R1, double *R2, double *RApertures, double *EApertures, @@ -50,20 +51,26 @@ 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; + 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 (KickAngle) { /* Convert corrector component to polynomial coefficients */ + B0 = -sin(KickAngle[0]) / le; + A0 = sin(KickAngle[1]) / le; + } + if (max_order >= 1) { + k1_entrance_angle = B[1] * entrance_angle; + k1_exit_angle = B[1] * exit_angle; } #pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD) default(none) \ shared(r,num_particles,R1,T1,R2,T2,RApertures,EApertures,\ - irho,gK,A,B,L1,L2,K1,K2,max_order,num_int_steps,scaling,\ - entrance_angle,exit_angle,\ + 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->gK=gK; + Elem->gK_entrance=FullGap*FringeInt1; + Elem->gK_exit=FullGap*FringeInt2; Elem->R1=R1; Elem->R2=R2; Elem->T1=T1; @@ -196,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->gK, + 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); @@ -235,7 +248,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(); @@ -259,13 +274,13 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit, - gK, + 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")); @@ -276,13 +291,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")); 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/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..eebef678f 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*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/drift_exactbend.h b/atintegrators/drift_exactbend.h new file mode 100644 index 000000000..336551383 --- /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 L, double irho) +{ + /* 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/drift_fast.h b/atintegrators/drift_fast.h new file mode 100644 index 000000000..6ebad8752 --- /dev/null +++ b/atintegrators/drift_fast.h @@ -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/exactbendfringe.c b/atintegrators/exactbendfringe.c index 3c8e8ac52..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,62 +104,106 @@ 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 bend_fringe_test(double *r6, double irho, double gK) +static void quad_wedge(double *r6, double k1_theta) { - /* 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); - + 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; +} - r6[x_] += dxf; - r6[y_] = yf; - r6[py_] -= dpyf; - r6[ct_] -= dct; +void bend_fringe(double *r6, double irho, double gK) +{ + 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; } diff --git a/atintegrators/kick_h_k0h_k1h_kn.h b/atintegrators/kick_h_k0h_k1h_kn.h new file mode 100644 index 000000000..a699cb23a --- /dev/null +++ b/atintegrators/kick_h_k0h_k1h_kn.h @@ -0,0 +1,42 @@ +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 +(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 = r6[0]; + double y = r6[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; + + 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 new file mode 100644 index 000000000..65a1b1fb6 --- /dev/null +++ b/atintegrators/kick_k1h_kn.h @@ -0,0 +1,25 @@ +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 = r6[0]; + double y = r6[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; + + 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/atintegrators/kickrad_k1h_kn.h b/atintegrators/kickrad_k1h_kn.h new file mode 100644 index 000000000..833bec85d --- /dev/null +++ b/atintegrators/kickrad_k1h_kn.h @@ -0,0 +1,120 @@ +/*********************************************************************** + 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 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)*(1.0-nrm) + SQR(ypr)*(1.0-nrm); + + return SQR(bx) + SQR(by) - SQR(bx*xpr + by*ypr)/v_norm2; +} + +//static void ex_bndthinkickrad(double* r, double* A, double* B, double L, double irho, double E0, 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) + +/***************************************************************************** +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 + + ******************************************************************************/ +{ + 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 */ + 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.0 / (1.0+r6[4]); + r6[1] = xpr/p_norm; + r6[3] = ypr/p_norm; + + /* Multipole kick */ + 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/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) diff --git a/pyat/at/lattice/elements/rectangular_bend.py b/pyat/at/lattice/elements/rectangular_bend.py index 62545f3f6..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,34 +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", }: - elem = self.copy() - elem.PassMethod = passmethod - theta = elem.BendingAngle + # 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 * theta) - 1.0) / theta + sin(0.5 * theta) / 12) + # Analytical estimate + x0ref = elem.Length * ((cos(0.5 * tta) - 1.0) / tta + sin(0.5 * tta) / 12) - # Search if there are multipoles - if checkmul(self): + # cancel output angle x0ref = float(fsolve(cross, x0ref)) + else: + x0ref = 0.0 self.X0ref = x0ref rout = self.track(np.zeros(6)) diff --git a/pyat/at/load/xsuite.py b/pyat/at/load/xsuite.py index 8e3373ef1..5c3a1bc2f 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 @@ -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 @@ -329,8 +338,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 +414,6 @@ class Multipole(XsElement): _at_integrator = "yoshida4" _xsuite2at_attr = XsElement._xsuite2at_attr | { "order": "MaxOrder", - "num_multipole_kicks": "NumIntSteps", } def _set_at_transforms(self) -> dict: @@ -468,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) @@ -478,13 +484,20 @@ 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 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.""" @@ -516,31 +529,45 @@ 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.""" + + 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 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 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: @@ -556,9 +583,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 @@ -571,19 +598,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.""" @@ -619,13 +633,33 @@ class Bend(Multipole): "edge_entry_angle": "EntranceAngle", "edge_exit_angle": "ExitAngle", } - _mag_order = 1 - _edge_to_xs: ClassVar[dict[bool, dict[bool, str]]] = { - True: {True: "full", False: "dipole-only"}, - False: {True: "linear", False: "linear"}, + _mag_order = 0 + _default_bend_fringe: ClassVar[dict[bool, int]] = { + True: 4, + False: 1, + } + _at2xsuite_edge: ClassVar[dict[int, str]] = { + 0: "suppressed", + 1: "linear", + 2: "linear", + 3: "linear", + 4: "dipole-only", + } + _xsuite2at_edge: ClassVar[dict[str, tuple]] = { + "suppressed": (0, None), + "linear": (1, 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._xsuite2at_edge[self.get(xskey, "linear")] + 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") @@ -633,31 +667,37 @@ 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 - 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_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) - 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) 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 @@ -705,8 +745,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 @@ -746,13 +786,13 @@ 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 -class NotInAT: +class NotInAT(_XsFactory): """Class for Xsuite elements without AT equivalent.""" @classmethod @@ -771,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 @@ -783,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 @@ -797,6 +837,27 @@ 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 + + #################################### + atparams.pop("PolynomA", None) + atparams.pop("PolynomB", None) + atparams.pop("MaxOrder", None) + ##################################### + + return Multipole.from_at(PolynomA=pola, PolynomB=polb, MaxOrder=0, **atparams) + + _xsclass: dict[str, type[XsElement]] = { "Marker": Marker, "Drift": Drift, @@ -810,7 +871,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, @@ -821,6 +882,7 @@ def from_at(cls, **atparams) -> XsElement: elt.Multipole: Multipole, elt.ThinMultipole: Multipole, elt.Dipole: Dipole, + elt.Corrector: Corrector, } @@ -1046,9 +1108,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 { 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, )