From da61a44e78de54e2ca2dcce776f43c62857db43e Mon Sep 17 00:00:00 2001 From: kevin Date: Wed, 18 Feb 2026 11:35:31 +0100 Subject: [PATCH 1/4] Make ExactQuantPass methods from ExactRadPass methods and SymplecticQuantPass methods --- atintegrators/ExactMultipoleQuantPass.c | 284 ++++++++++++++ atintegrators/ExactRectangularBendQuantPass.c | 362 ++++++++++++++++++ 2 files changed, 646 insertions(+) create mode 100644 atintegrators/ExactMultipoleQuantPass.c create mode 100644 atintegrators/ExactRectangularBendQuantPass.c diff --git a/atintegrators/ExactMultipoleQuantPass.c b/atintegrators/ExactMultipoleQuantPass.c new file mode 100644 index 0000000000..2962df7f4a --- /dev/null +++ b/atintegrators/ExactMultipoleQuantPass.c @@ -0,0 +1,284 @@ +#include "atconstants.h" +#include "atelem.c" +#include "atlalib.c" +#include "atquantlib.c" +#include "driftkick.c" /* fastdrift.c, strthinkick.c */ +#include "exactdrift.c" +#include "exactmultipolefringe.c" + +struct elem { + double Length; + double *PolynomA; + double *PolynomB; + int MaxOrder; + int NumIntSteps; + double Energy; + /* Optional fields */ + double Scaling; + int FringeQuadEntrance; + int FringeQuadExit; + double *R1; + double *R2; + double *T1; + double *T2; + double *RApertures; + double *EApertures; + double *KickAngle; +}; + +static void multipole_pass( + double *r, double le, double *A, double *B, int max_order, int num_int_steps, + int FringeQuadEntrance, int FringeQuadExit, /* 0 (no fringe), else */ + double *T1, double *T2, double *R1, double *R2, double *RApertures, + double *EApertures, double *KickAngle, double scaling, double E0, pcg32_random_t* rng, int num_particles) +{ + double SL = le / num_int_steps; + double L1 = SL * DRIFT1; + double L2 = SL * DRIFT2; + double K1 = SL * KICK1; + double K2 = SL * KICK2; + double qe = 1.60217733e-19; + double epsilon0 = 8.854187817e-12; + double clight = 2.99792458e8; + double emass = 510998.9461; /* electron mass in eV */ /* 9.10938188e-31; in kg*/ + double hbar = 1.054571726e-34; + double pi = 3.14159265358979; + double alpha0 = qe * qe / (4 * pi * epsilon0 * hbar * clight); + double B0 = B[0]; + double A0 = A[0]; + + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + B[0] -= sin(KickAngle[0]) / le; + A[0] += 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, A, B, L1, L2, K1, K2, max_order, \ + FringeQuadEntrance, FringeQuadExit, \ + num_int_steps, scaling, le) + */ + for (int c = 0; c < num_particles; c++) { /*Loop over particles */ + double *r6 = r + c * 6; + if (!atIsNaN(r6[0])) { + int m; + + /* Check for change of reference momentum */ + if (scaling != 1.0) ATChangePRef(r6, scaling); + + /* misalignment at entrance */ + if (T1) ATaddvv(r6, T1); + if (R1) ATmultmv(r6, R1); + + /* Check physical apertures at the entrance of the magnet */ + if (RApertures) checkiflostRectangularAp(r6, RApertures); + if (EApertures) checkiflostEllipticalAp(r6, EApertures); + + /* Fringe field effect */ + if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 0); + + /* integrator */ + for (m = 0; m < num_int_steps; m++) { /* Loop over slices */ + int i; + double ng, ec, de, energy, gamma, cstec, cstng; + double ds, rho, dxp, dyp; + int nph; + double p_norm = 1.0 / (1.0 + r6[4]); + double NormL1 = L1 * p_norm; + double NormL2 = L2 * p_norm; + double dpp0 = r6[4]; + double xp0 = r6[1] * p_norm; + double yp0 = r6[3] * p_norm; + double s0 = r6[5]; + + exact_drift(r6, L1); + strthinkick(r6, A, B, K1, max_order); + exact_drift(r6, L2); + strthinkick(r6, A, B, K2, max_order); + exact_drift(r6, L2); + strthinkick(r6, A, B, K1, max_order); + exact_drift(r6, L1); + + energy = dpp0 * E0 + E0; + + gamma = energy / emass; /* emass in eV */ + cstec = 3.0 * gamma * gamma * gamma * clight / (2.0) * hbar / qe; + cstng = 5.0 * sqrt(3.0) * alpha0 * gamma / (6.0); + + dxp = r6[1] * p_norm - xp0; + dyp = r6[3] * p_norm - yp0; + ds = r6[5] - s0; + + rho = (SL + ds) / sqrt(dxp * dxp + dyp * dyp); + + ng = cstng / rho * (SL + ds); + ec = cstec / rho; + + nph = atrandp_r(rng, ng); + + de = 0.0; + for (i = 0; i < nph; i++) { + de = de + getEnergy(rng, ec); + }; + r6[4] = r6[4] - de / E0; + r6[1] = r6[1] * p_norm * (1 + r6[4]); + r6[3] = r6[3] * p_norm * (1 + r6[4]); + } + + /* Convert absolute path length to path lengthening */ + r6[5] -= le; + + /* Fringe field effect */ + if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 0); + + /* Check physical apertures at the exit of the magnet */ + if (RApertures) checkiflostRectangularAp(r6, RApertures); + if (EApertures) checkiflostEllipticalAp(r6, EApertures); + + /* Misalignment at exit */ + if (R2) ATmultmv(r6, R2); + if (T2) ATaddvv(r6, T2); + + /* Check for change of reference momentum */ + 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) +ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, + double *r_in, int num_particles, + struct parameters *Param) { + double energy; + if (!Elem) { + double Length = atGetDouble(ElemData, "Length"); check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); check_error(); + /*optional fields*/ + double Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); + double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); + int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); + int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); + double *R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); + double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); + double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); check_error(); + + if (NumIntSteps <= 0) { + atError("NumIntSteps must be positive"); check_error(); + } + + Elem = (struct elem *)atMalloc(sizeof(struct elem)); + Elem->Length = Length; + Elem->PolynomA = PolynomA; + Elem->PolynomB = PolynomB; + Elem->MaxOrder = MaxOrder; + Elem->NumIntSteps = NumIntSteps; + /*optional fields*/ + Elem->Energy=Energy; + Elem->Scaling=Scaling; + Elem->FringeQuadEntrance=FringeQuadEntrance; + Elem->FringeQuadExit=FringeQuadExit; + Elem->R1 = R1; + Elem->R2 = R2; + Elem->T1 = T1; + Elem->T2 = T2; + Elem->EApertures = EApertures; + Elem->RApertures = RApertures; + Elem->KickAngle = KickAngle; + } + energy = atEnergy(Param->energy, Elem->Energy); + + multipole_pass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, + Elem->MaxOrder, Elem->NumIntSteps, + Elem->FringeQuadEntrance, Elem->FringeQuadExit, + Elem->T1, Elem->T2, Elem->R1, Elem->R2, + Elem->RApertures, Elem->EApertures, + Elem->KickAngle, Elem->Scaling, energy, Param->thread_rng, num_particles); + return Elem; +} + +MODULE_DEF(ExactMultipoleQuantPass) /* Dummy module initialisation */ + +#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ + +#if defined(MATLAB_MEX_FILE) +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { + if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; + double *r_in; + const mxArray *ElemData = prhs[0]; + int num_particles = mxGetN(prhs[1]); + + double Length = atGetDouble(ElemData, "Length"); check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); check_error(); + /*optional fields*/ + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); + double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); + int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); + int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); + double *R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); + double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); + double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); check_error(); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + + if (NumIntSteps <= 0) { + atError("NumIntSteps must be positive"); check_error(); + } + /* ALLOCATE memory for the output array of the same size as the input */ + plhs[0] = mxDuplicateArray(prhs[1]); + r_in = mxGetDoubles(plhs[0]); + + multipole_pass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, + FringeQuadEntrance, FringeQuadExit, + T1, T2, R1, R2, + RApertures, EApertures, + KickAngle, Scaling, Energy, &pcg32_global, num_particles); + } else if (nrhs == 0) { + /* list of required fields */ + int i0 = 0; + plhs[0] = mxCreateCellMatrix(5, 1); + mxSetCell(plhs[0], i0++, mxCreateString("Length")); + mxSetCell(plhs[0], i0++, mxCreateString("PolynomA")); + mxSetCell(plhs[0], i0++, mxCreateString("PolynomB")); + mxSetCell(plhs[0], i0++, mxCreateString("MaxOrder")); + mxSetCell(plhs[0], i0++, mxCreateString("NumIntSteps")); + + if (nlhs > 1) { + /* list of optional fields */ + int i1 = 0; + plhs[1] = mxCreateCellMatrix(11, 1); + mxSetCell(plhs[1], i1++, mxCreateString("Energy")); + mxSetCell(plhs[1], i1++, mxCreateString("FieldScaling")); + mxSetCell(plhs[1], i1++, mxCreateString("FringeQuadEntrance")); + mxSetCell(plhs[1], i1++, mxCreateString("FringeQuadExit")); + mxSetCell(plhs[1], i1++, mxCreateString("T1")); + mxSetCell(plhs[1], i1++, mxCreateString("T2")); + mxSetCell(plhs[1], i1++, mxCreateString("R1")); + mxSetCell(plhs[1], i1++, mxCreateString("R2")); + mxSetCell(plhs[1], i1++, mxCreateString("RApertures")); + mxSetCell(plhs[1], i1++, mxCreateString("EApertures")); + mxSetCell(plhs[1], i1++, mxCreateString("KickAngle")); + } + } else { + mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments"); + } +} +#endif /*MATLAB_MEX_FILE*/ diff --git a/atintegrators/ExactRectangularBendQuantPass.c b/atintegrators/ExactRectangularBendQuantPass.c new file mode 100644 index 0000000000..f1734a67d0 --- /dev/null +++ b/atintegrators/ExactRectangularBendQuantPass.c @@ -0,0 +1,362 @@ +#include "atconstants.h" +#include "atelem.c" +#include "atlalib.c" +#include "atphyslib.c" +#include "atquantlib.c" +#include "exactdrift.c" +#include "driftkick.c" /* strthinkick.c */ +#include "exactbendfringe.c" +#include "exactmultipolefringe.c" +#include + +struct elem +{ + double Length; + double *PolynomA; + double *PolynomB; + int MaxOrder; + int NumIntSteps; + double BendingAngle; + double EntranceAngle; + double ExitAngle; + double Energy; + /* Optional fields */ + double Scaling; + int FringeBendEntrance; + int FringeBendExit; + int FringeQuadEntrance; + int FringeQuadExit; + double gK; + double x0ref; + double refdz; + double *R1; + double *R2; + double *T1; + double *T2; + double *RApertures; + double *EApertures; + double *KickAngle; +}; + +static void ExactRectangularBendQuant(double *r, double le, double bending_angle, + double *A, double *B, + int max_order, int num_int_steps, + double entrance_angle, double exit_angle, + int FringeBendEntrance, int FringeBendExit, + int FringeQuadEntrance, int FringeQuadExit, + double gK, double x0ref, double refdz, + double *T1, double *T2, + double *R1, double *R2, + double *RApertures, double *EApertures, + double *KickAngle, double scaling, double E0, + pcg32_random_t* rng, int num_particles) +{ + double irho = bending_angle / le; + double phi2 = 0.5 * bending_angle; + double LR = phi2 < 1.e-10 ? le : le *sin(phi2) / phi2; + double SL = LR/num_int_steps; + double L1 = SL*DRIFT1; + double L2 = SL*DRIFT2; + double K1 = SL*KICK1; + double K2 = SL*KICK2; + double qe = 1.60217733e-19; + double epsilon0 = 8.854187817e-12; + double clight = 2.99792458e8; + double emass = 510998.9461; /* electron mass in eV */ /* 9.10938188e-31; in kg*/ + double hbar = 1.054571726e-34; + double pi = 3.14159265358979; + double alpha0 = qe * qe / (4 * pi * epsilon0 * hbar * clight); + double B0 = B[0]; + double A0 = A[0]; + + if (KickAngle) { /* Convert corrector component to polynomial coefficients */ + B[0] -= sin(KickAngle[0])/le; + A[0] += sin(KickAngle[1])/le; + } + B[0] += 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,\ + entrance_angle,exit_angle,x0ref,refdz,\ + FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ + LR,le,phi2) + */ + for (int c = 0; cenergy); check_error(); + double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); + int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); + 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 x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error(); + double refdz=atGetOptionalDouble(ElemData,"RefDZ", 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(); + double *T2=atGetOptionalDoubleArray(ElemData,"T2"); check_error(); + double *EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); + double *RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); + double *KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); + + if (NumIntSteps <= 0) { + atError("NumIntSteps must be positive"); check_error(); + } + + Elem = (struct elem*)atMalloc(sizeof(struct elem)); + Elem->Length=Length; + Elem->PolynomA=PolynomA; + Elem->PolynomB=PolynomB; + Elem->MaxOrder=MaxOrder; + Elem->NumIntSteps=NumIntSteps; + Elem->BendingAngle=BendingAngle; + Elem->EntranceAngle=EntranceAngle; + Elem->ExitAngle=ExitAngle; + /*optional fields*/ + Elem->Energy=Energy; + Elem->Scaling=Scaling; + Elem->FringeBendEntrance=FringeBendEntrance; + Elem->FringeBendExit=FringeBendExit; + Elem->FringeQuadEntrance=FringeQuadEntrance; + Elem->FringeQuadExit=FringeQuadExit; + Elem->gK=gK; + Elem->x0ref=x0ref; + Elem->refdz=refdz; + Elem->R1=R1; + Elem->R2=R2; + Elem->T1=T1; + Elem->T2=T2; + Elem->EApertures=EApertures; + Elem->RApertures=RApertures; + Elem->KickAngle=KickAngle; + } + irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + + ExactRectangularBendQuant(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->T1, Elem->T2, Elem->R1, Elem->R2, + Elem->RApertures, Elem->EApertures, + Elem->KickAngle, Elem->Scaling, energy, + Param->thread_rng, num_particles); + return Elem; +} + +MODULE_DEF(ExactRectangularBendQuant) /* Dummy module initialisation */ + +#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ + +#if defined(MATLAB_MEX_FILE) +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; + double irho, Energy; + double *r_in; + const mxArray *ElemData = prhs[0]; + int num_particles = mxGetN(prhs[1]); + if (mxGetM(prhs[1]) != 6) mexErrMsgTxt("Second argument must be a 6 x N matrix"); + + double Length=atGetDouble(ElemData,"Length"); check_error(); + double *PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); + double *PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); + int MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); + int NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); + double BendingAngle=atGetOptionalDouble(ElemData,"BendingAngle",0.0); check_error(); + double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); + double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); + /*optional fields*/ + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); + double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); + int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); + 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 x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error(); + double refdz=atGetOptionalDouble(ElemData,"RefDZ", 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(); + double *T2=atGetOptionalDoubleArray(ElemData,"T2"); check_error(); + double *EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); + double *RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); + double *KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); + double irho = BendingAngle/Length; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + + if (NumIntSteps <= 0) { + atError("NumIntSteps must be positive"); check_error(); + } + + /* ALLOCATE memory for the output array of the same size as the input */ + plhs[0] = mxDuplicateArray(prhs[1]); + r_in = mxGetDoubles(plhs[0]); + ExactRectangularBendQuant(r_in, Length, BendingAngle, PolynomA, PolynomB, + MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, + FringeBendEntrance, FringeBendExit, + FringeQuadEntrance, FringeQuadExit, + gK, x0ref, refdz, + T1, T2, R1, R2, RApertures, EApertures, + KickAngle, Scaling, Energy, &pcg32_global, num_particles); + } else if (nrhs == 0) { + /* list of required fields */ + int i0 = 0; + plhs[0] = mxCreateCellMatrix(8, 1); + mxSetCell(plhs[0], i0++, mxCreateString("Length")); + mxSetCell(plhs[0], i0++, mxCreateString("PolynomA")); + mxSetCell(plhs[0], i0++, mxCreateString("PolynomB")); + mxSetCell(plhs[0], i0++, mxCreateString("MaxOrder")); + mxSetCell(plhs[0], i0++, mxCreateString("NumIntSteps")); + mxSetCell(plhs[0], i0++, mxCreateString("BendingAngle")); + mxSetCell(plhs[0], i0++, mxCreateString("EntranceAngle")); + mxSetCell(plhs[0], i0++, mxCreateString("ExitAngle")); + + if (nlhs>1) { /* list of optional fields */ + int i1 = 0; + 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("X0ref")); + mxSetCell(plhs[1], i1++, mxCreateString("RefDZ")); + mxSetCell(plhs[1], i1++, mxCreateString("T1")); + mxSetCell(plhs[1], i1++, mxCreateString("T2")); + mxSetCell(plhs[1], i1++, mxCreateString("R1")); + mxSetCell(plhs[1], i1++, mxCreateString("R2")); + mxSetCell(plhs[1], i1++, mxCreateString("RApertures")); + mxSetCell(plhs[1], i1++, mxCreateString("EApertures")); + mxSetCell(plhs[1], i1++, mxCreateString("KickAngle")); + mxSetCell(plhs[1], i1++, mxCreateString("FieldScaling")); + } + } + else { + mexErrMsgIdAndTxt("AT:WrongArg","Needs 0 or 2 arguments"); + } +} +#endif /* MATLAB_MEX_FILE */ From 37373358ee7c3e1edbe6384017a14b4ce6a3fdbb Mon Sep 17 00:00:00 2001 From: kevin Date: Wed, 18 Feb 2026 16:16:57 +0100 Subject: [PATCH 2/4] fix build and created ExactSectorBendQuantPass --- atintegrators/ExactMultipoleQuantPass.c | 217 +++++---- atintegrators/ExactRectangularBendQuantPass.c | 378 +++++++++------- atintegrators/ExactSectorBendQuantPass.c | 418 ++++++++++++++++++ 3 files changed, 788 insertions(+), 225 deletions(-) create mode 100644 atintegrators/ExactSectorBendQuantPass.c diff --git a/atintegrators/ExactMultipoleQuantPass.c b/atintegrators/ExactMultipoleQuantPass.c index 2962df7f4a..2a28d2ec1e 100644 --- a/atintegrators/ExactMultipoleQuantPass.c +++ b/atintegrators/ExactMultipoleQuantPass.c @@ -2,11 +2,12 @@ #include "atelem.c" #include "atlalib.c" #include "atquantlib.c" -#include "driftkick.c" /* fastdrift.c, strthinkick.c */ +#include "driftkick.c" /* fastdrift.c, strthinkick.c */ #include "exactdrift.c" #include "exactmultipolefringe.c" -struct elem { +struct elem +{ double Length; double *PolynomA; double *PolynomB; @@ -27,10 +28,10 @@ struct elem { }; static void multipole_pass( - double *r, double le, double *A, double *B, int max_order, int num_int_steps, - int FringeQuadEntrance, int FringeQuadExit, /* 0 (no fringe), else */ - double *T1, double *T2, double *R1, double *R2, double *RApertures, - double *EApertures, double *KickAngle, double scaling, double E0, pcg32_random_t* rng, int num_particles) + double *r, double le, double *A, double *B, int max_order, int num_int_steps, + int FringeQuadEntrance, int FringeQuadExit, /* 0 (no fringe), else */ + double *T1, double *T2, double *R1, double *R2, double *RApertures, + double *EApertures, double *KickAngle, double scaling, double E0, pcg32_random_t *rng, int num_particles) { double SL = le / num_int_steps; double L1 = SL * DRIFT1; @@ -47,7 +48,8 @@ static void multipole_pass( double B0 = B[0]; double A0 = A[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; } @@ -59,27 +61,36 @@ static void multipole_pass( FringeQuadEntrance, FringeQuadExit, \ num_int_steps, scaling, le) */ - for (int c = 0; c < num_particles; c++) { /*Loop over particles */ + for (int c = 0; c < num_particles; c++) + { /*Loop over particles */ double *r6 = r + c * 6; - if (!atIsNaN(r6[0])) { + if (!atIsNaN(r6[0])) + { int m; /* Check for change of reference momentum */ - if (scaling != 1.0) ATChangePRef(r6, scaling); + if (scaling != 1.0) + ATChangePRef(r6, scaling); /* misalignment at entrance */ - if (T1) ATaddvv(r6, T1); - if (R1) ATmultmv(r6, R1); + if (T1) + ATaddvv(r6, T1); + if (R1) + ATmultmv(r6, R1); /* Check physical apertures at the entrance of the magnet */ - if (RApertures) checkiflostRectangularAp(r6, RApertures); - if (EApertures) checkiflostEllipticalAp(r6, EApertures); + if (RApertures) + checkiflostRectangularAp(r6, RApertures); + if (EApertures) + checkiflostEllipticalAp(r6, EApertures); /* Fringe field effect */ - if (FringeQuadEntrance) multipole_fringe(r6, le, A, B, max_order, 1.0, 0); + if (FringeQuadEntrance) + multipole_fringe(r6, le, A, B, max_order, 1.0, 0); /* integrator */ - for (m = 0; m < num_int_steps; m++) { /* Loop over slices */ + for (m = 0; m < num_int_steps; m++) + { /* Loop over slices */ int i; double ng, ec, de, energy, gamma, cstec, cstng; double ds, rho, dxp, dyp; @@ -110,16 +121,17 @@ static void multipole_pass( dyp = r6[3] * p_norm - yp0; ds = r6[5] - s0; - rho = (SL + ds) / sqrt(dxp * dxp + dyp * dyp); + rho = (ds) / sqrt(dxp * dxp + dyp * dyp); - ng = cstng / rho * (SL + ds); + ng = cstng / rho * (ds); ec = cstec / rho; nph = atrandp_r(rng, ng); de = 0.0; - for (i = 0; i < nph; i++) { - de = de + getEnergy(rng, ec); + for (i = 0; i < nph; i++) + { + de = de + getEnergy(rng, ec); }; r6[4] = r6[4] - de / E0; r6[1] = r6[1] * p_norm * (1 + r6[4]); @@ -130,18 +142,24 @@ static void multipole_pass( r6[5] -= le; /* Fringe field effect */ - if (FringeQuadExit) multipole_fringe(r6, le, A, B, max_order, -1.0, 0); + if (FringeQuadExit) + multipole_fringe(r6, le, A, B, max_order, -1.0, 0); /* Check physical apertures at the exit of the magnet */ - if (RApertures) checkiflostRectangularAp(r6, RApertures); - if (EApertures) checkiflostEllipticalAp(r6, EApertures); + if (RApertures) + checkiflostRectangularAp(r6, RApertures); + if (EApertures) + checkiflostEllipticalAp(r6, EApertures); /* Misalignment at exit */ - if (R2) ATmultmv(r6, R2); - if (T2) ATaddvv(r6, T2); + if (R2) + ATmultmv(r6, R2); + if (T2) + ATaddvv(r6, T2); /* Check for change of reference momentum */ - if (scaling != 1.0) ATChangePRef(r6, 1.0/scaling); + if (scaling != 1.0) + ATChangePRef(r6, 1.0 / scaling); } } /* Remove corrector component in polynomial coefficients */ @@ -151,30 +169,49 @@ static void multipole_pass( #if defined(MATLAB_MEX_FILE) || defined(PYAT) ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, - double *r_in, int num_particles, - struct parameters *Param) { + double *r_in, int num_particles, struct parameters *Param) +{ double energy; - if (!Elem) { - double Length = atGetDouble(ElemData, "Length"); check_error(); - double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); check_error(); - double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); check_error(); - int MaxOrder = atGetLong(ElemData, "MaxOrder"); check_error(); - int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); check_error(); + if (!Elem) + { + double Length = atGetDouble(ElemData, "Length"); + check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); + check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); + check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); + check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); + check_error(); /*optional fields*/ - double Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); - double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); - int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); - int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); - double *R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); - double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); - double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); - double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); - double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); check_error(); - double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); check_error(); - double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); check_error(); - - if (NumIntSteps <= 0) { - atError("NumIntSteps must be positive"); check_error(); + double Energy = atGetOptionalDouble(ElemData, "Energy", Param->energy); + check_error(); + double Scaling = atGetOptionalDouble(ElemData, "FieldScaling", 1.0); + check_error(); + int FringeQuadEntrance = atGetOptionalLong(ElemData, "FringeQuadEntrance", 0); + check_error(); + int FringeQuadExit = atGetOptionalLong(ElemData, "FringeQuadExit", 0); + check_error(); + double *R1 = atGetOptionalDoubleArray(ElemData, "R1"); + check_error(); + double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); + check_error(); + double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); + check_error(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); + check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); + check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); + check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); + check_error(); + + if (NumIntSteps <= 0) + { + atError("NumIntSteps must be positive"); + check_error(); } Elem = (struct elem *)atMalloc(sizeof(struct elem)); @@ -184,10 +221,10 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, Elem->MaxOrder = MaxOrder; Elem->NumIntSteps = NumIntSteps; /*optional fields*/ - Elem->Energy=Energy; - Elem->Scaling=Scaling; - Elem->FringeQuadEntrance=FringeQuadEntrance; - Elem->FringeQuadExit=FringeQuadExit; + Elem->Energy = Energy; + Elem->Scaling = Scaling; + Elem->FringeQuadEntrance = FringeQuadEntrance; + Elem->FringeQuadExit = FringeQuadExit; Elem->R1 = R1; Elem->R2 = R2; Elem->T1 = T1; @@ -212,35 +249,56 @@ MODULE_DEF(ExactMultipoleQuantPass) /* Dummy module initialisation */ #endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ #if defined(MATLAB_MEX_FILE) -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs >= 2) { +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + if (nrhs >= 2) + { double rest_energy = 0.0; double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); - double Length = atGetDouble(ElemData, "Length"); check_error(); - double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); check_error(); - double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); check_error(); - int MaxOrder = atGetLong(ElemData, "MaxOrder"); check_error(); - int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); check_error(); + double Length = atGetDouble(ElemData, "Length"); + check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); + check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); + check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); + check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); + check_error(); /*optional fields*/ - double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); - double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); - int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); - int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); - double *R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); - double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); - double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); - double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); - double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); check_error(); - double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); check_error(); - double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); check_error(); - if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); - - if (NumIntSteps <= 0) { - atError("NumIntSteps must be positive"); check_error(); + double Energy = atGetOptionalDouble(ElemData, "Energy", 0.0); + check_error(); + double Scaling = atGetOptionalDouble(ElemData, "FieldScaling", 1.0); + check_error(); + int FringeQuadEntrance = atGetOptionalLong(ElemData, "FringeQuadEntrance", 0); + check_error(); + int FringeQuadExit = atGetOptionalLong(ElemData, "FringeQuadExit", 0); + check_error(); + double *R1 = atGetOptionalDoubleArray(ElemData, "R1"); + check_error(); + double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); + check_error(); + double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); + check_error(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); + check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); + check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); + check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); + check_error(); + if (nrhs > 2) + atProperties(prhs[2], &Energy, &rest_energy, &charge); + + if (NumIntSteps <= 0) + { + atError("NumIntSteps must be positive"); + check_error(); } /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); @@ -251,7 +309,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { T1, T2, R1, R2, RApertures, EApertures, KickAngle, Scaling, Energy, &pcg32_global, num_particles); - } else if (nrhs == 0) { + } + else if (nrhs == 0) + { /* list of required fields */ int i0 = 0; plhs[0] = mxCreateCellMatrix(5, 1); @@ -261,7 +321,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { mxSetCell(plhs[0], i0++, mxCreateString("MaxOrder")); mxSetCell(plhs[0], i0++, mxCreateString("NumIntSteps")); - if (nlhs > 1) { + if (nlhs > 1) + { /* list of optional fields */ int i1 = 0; plhs[1] = mxCreateCellMatrix(11, 1); @@ -277,7 +338,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { mxSetCell(plhs[1], i1++, mxCreateString("EApertures")); mxSetCell(plhs[1], i1++, mxCreateString("KickAngle")); } - } else { + } + else + { mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments"); } } diff --git a/atintegrators/ExactRectangularBendQuantPass.c b/atintegrators/ExactRectangularBendQuantPass.c index f1734a67d0..33b320f085 100644 --- a/atintegrators/ExactRectangularBendQuantPass.c +++ b/atintegrators/ExactRectangularBendQuantPass.c @@ -1,13 +1,11 @@ #include "atconstants.h" #include "atelem.c" #include "atlalib.c" -#include "atphyslib.c" #include "atquantlib.c" #include "exactdrift.c" -#include "driftkick.c" /* strthinkick.c */ +#include "driftkick.c" /* strthinkick.c */ #include "exactbendfringe.c" #include "exactmultipolefringe.c" -#include struct elem { @@ -39,26 +37,26 @@ struct elem }; static void ExactRectangularBendQuant(double *r, double le, double bending_angle, - double *A, double *B, - int max_order, int num_int_steps, - double entrance_angle, double exit_angle, - int FringeBendEntrance, int FringeBendExit, - int FringeQuadEntrance, int FringeQuadExit, - double gK, double x0ref, double refdz, - double *T1, double *T2, - double *R1, double *R2, - double *RApertures, double *EApertures, - double *KickAngle, double scaling, double E0, - pcg32_random_t* rng, int num_particles) + double *A, double *B, + int max_order, int num_int_steps, + double entrance_angle, double exit_angle, + int FringeBendEntrance, int FringeBendExit, + int FringeQuadEntrance, int FringeQuadExit, + double gK, double x0ref, double refdz, + double *T1, double *T2, + double *R1, double *R2, + double *RApertures, double *EApertures, + double *KickAngle, double scaling, double E0, + pcg32_random_t *rng, int num_particles) { double irho = bending_angle / le; double phi2 = 0.5 * bending_angle; - double LR = phi2 < 1.e-10 ? le : le *sin(phi2) / phi2; - double SL = LR/num_int_steps; - double L1 = SL*DRIFT1; - double L2 = SL*DRIFT2; - double K1 = SL*KICK1; - double K2 = SL*KICK2; + double LR = phi2 < 1.e-10 ? le : le * sin(phi2) / phi2; + double SL = LR / num_int_steps; + double L1 = SL * DRIFT1; + double L2 = SL * DRIFT2; + double K1 = SL * KICK1; + double K2 = SL * KICK2; double qe = 1.60217733e-19; double epsilon0 = 8.854187817e-12; double clight = 2.99792458e8; @@ -69,9 +67,10 @@ static void ExactRectangularBendQuant(double *r, double le, double bending_angle double B0 = B[0]; double A0 = A[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 */ + B[0] -= sin(KickAngle[0]) / le; + A[0] += sin(KickAngle[1]) / le; } B[0] += irho; /* @@ -82,33 +81,41 @@ static void ExactRectangularBendQuant(double *r, double le, double bending_angle FringeBendEntrance,FringeBendExit,FringeQuadEntrance,FringeQuadExit,\ LR,le,phi2) */ - for (int c = 0; cenergy); check_error(); - double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); - int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); - 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 x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error(); - double refdz=atGetOptionalDouble(ElemData,"RefDZ", 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(); - double *T2=atGetOptionalDoubleArray(ElemData,"T2"); check_error(); - double *EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); - double *RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); - double *KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); - - if (NumIntSteps <= 0) { - atError("NumIntSteps must be positive"); check_error(); + double Energy = atGetOptionalDouble(ElemData, "Energy", Param->energy); + check_error(); + double Scaling = atGetOptionalDouble(ElemData, "FieldScaling", 1.0); + check_error(); + int FringeBendEntrance = atGetOptionalLong(ElemData, "FringeBendEntrance", 1); + check_error(); + 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 x0ref = atGetOptionalDouble(ElemData, "X0ref", 0.0); + check_error(); + double refdz = atGetOptionalDouble(ElemData, "RefDZ", 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(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); + check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); + check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); + check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); + check_error(); + + if (NumIntSteps <= 0) + { + atError("NumIntSteps must be positive"); + check_error(); + } + + /* 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(); } - Elem = (struct elem*)atMalloc(sizeof(struct elem)); - Elem->Length=Length; - Elem->PolynomA=PolynomA; - Elem->PolynomB=PolynomB; - Elem->MaxOrder=MaxOrder; - Elem->NumIntSteps=NumIntSteps; - Elem->BendingAngle=BendingAngle; - Elem->EntranceAngle=EntranceAngle; - Elem->ExitAngle=ExitAngle; + Elem = (struct elem *)atMalloc(sizeof(struct elem)); + Elem->Length = Length; + Elem->PolynomA = PolynomA; + Elem->PolynomB = PolynomB; + Elem->MaxOrder = MaxOrder; + Elem->NumIntSteps = NumIntSteps; + Elem->BendingAngle = BendingAngle; + Elem->EntranceAngle = EntranceAngle; + Elem->ExitAngle = ExitAngle; /*optional fields*/ - Elem->Energy=Energy; - Elem->Scaling=Scaling; - Elem->FringeBendEntrance=FringeBendEntrance; - Elem->FringeBendExit=FringeBendExit; - Elem->FringeQuadEntrance=FringeQuadEntrance; - Elem->FringeQuadExit=FringeQuadExit; - Elem->gK=gK; - Elem->x0ref=x0ref; - Elem->refdz=refdz; - Elem->R1=R1; - Elem->R2=R2; - Elem->T1=T1; - Elem->T2=T2; - Elem->EApertures=EApertures; - Elem->RApertures=RApertures; - Elem->KickAngle=KickAngle; + Elem->Energy = Energy; + Elem->Scaling = Scaling; + Elem->FringeBendEntrance = FringeBendEntrance; + Elem->FringeBendExit = FringeBendExit; + Elem->FringeQuadEntrance = FringeQuadEntrance; + Elem->FringeQuadExit = FringeQuadExit; + Elem->gK = gK; + Elem->x0ref = x0ref; + Elem->refdz = refdz; + Elem->R1 = R1; + Elem->R2 = R2; + Elem->T1 = T1; + Elem->T2 = T2; + Elem->EApertures = EApertures; + Elem->RApertures = RApertures; + Elem->KickAngle = KickAngle; } - irho = Elem->BendingAngle/Elem->Length; + irho = Elem->BendingAngle / Elem->Length; energy = atEnergy(Param->energy, Elem->Energy); ExactRectangularBendQuant(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->T1, Elem->T2, Elem->R1, Elem->R2, - Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, energy, - Param->thread_rng, num_particles); + Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, + Elem->FringeBendEntrance, Elem->FringeBendExit, + Elem->FringeQuadEntrance, Elem->FringeQuadExit, + Elem->gK, Elem->x0ref, Elem->refdz, + Elem->T1, Elem->T2, Elem->R1, Elem->R2, + Elem->RApertures, Elem->EApertures, + Elem->KickAngle, Elem->Scaling, energy, + Param->thread_rng, num_particles); return Elem; } -MODULE_DEF(ExactRectangularBendQuant) /* Dummy module initialisation */ +MODULE_DEF(ExactRectangularBendQuantPass) /* Dummy module initialisation */ #endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs >= 2) { + if (nrhs >= 2) + { double rest_energy = 0.0; double charge = -1.0; double irho, Energy; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); - if (mxGetM(prhs[1]) != 6) mexErrMsgTxt("Second argument must be a 6 x N matrix"); - - double Length=atGetDouble(ElemData,"Length"); check_error(); - double *PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); - double *PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); - int MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); - int NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - double BendingAngle=atGetOptionalDouble(ElemData,"BendingAngle",0.0); check_error(); - double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); - double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); + if (mxGetM(prhs[1]) != 6) + mexErrMsgTxt("Second argument must be a 6 x N matrix"); + + double Length = atGetDouble(ElemData, "Length"); + check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); + check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); + check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); + check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); + check_error(); + double BendingAngle = atGetOptionalDouble(ElemData, "BendingAngle", 0.0); + check_error(); + double EntranceAngle = atGetDouble(ElemData, "EntranceAngle"); + check_error(); + double ExitAngle = atGetDouble(ElemData, "ExitAngle"); + check_error(); /*optional fields*/ - double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); - double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); - int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); - 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 x0ref=atGetOptionalDouble(ElemData,"X0ref", 0.0); check_error(); - double refdz=atGetOptionalDouble(ElemData,"RefDZ", 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(); - double *T2=atGetOptionalDoubleArray(ElemData,"T2"); check_error(); - double *EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); - double *RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); - double *KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); - double irho = BendingAngle/Length; - if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); - - if (NumIntSteps <= 0) { - atError("NumIntSteps must be positive"); check_error(); + double Energy = atGetOptionalDouble(ElemData, "Energy", 0.0); + check_error(); + double Scaling = atGetOptionalDouble(ElemData, "FieldScaling", 1.0); + check_error(); + int FringeBendEntrance = atGetOptionalLong(ElemData, "FringeBendEntrance", 1); + check_error(); + 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 x0ref = atGetOptionalDouble(ElemData, "X0ref", 0.0); + check_error(); + double refdz = atGetOptionalDouble(ElemData, "RefDZ", 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(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); + check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); + check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); + check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); + check_error(); + double irho = BendingAngle / Length; + if (nrhs > 2) + atProperties(prhs[2], &Energy, &rest_energy, &charge); + + if (NumIntSteps <= 0) + { + atError("NumIntSteps must be positive"); + check_error(); } /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); ExactRectangularBendQuant(r_in, Length, BendingAngle, PolynomA, PolynomB, - MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, - FringeBendEntrance, FringeBendExit, - FringeQuadEntrance, FringeQuadExit, - gK, x0ref, refdz, - T1, T2, R1, R2, RApertures, EApertures, - KickAngle, Scaling, Energy, &pcg32_global, num_particles); - } else if (nrhs == 0) { + MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, + FringeBendEntrance, FringeBendExit, + FringeQuadEntrance, FringeQuadExit, + gK, x0ref, refdz, + T1, T2, R1, R2, RApertures, EApertures, + KickAngle, Scaling, Energy, &pcg32_global, num_particles); + } + else if (nrhs == 0) + { /* list of required fields */ int i0 = 0; plhs[0] = mxCreateCellMatrix(8, 1); @@ -334,7 +414,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mxSetCell(plhs[0], i0++, mxCreateString("EntranceAngle")); mxSetCell(plhs[0], i0++, mxCreateString("ExitAngle")); - if (nlhs>1) { /* list of optional fields */ + if (nlhs > 1) + { /* list of optional fields */ int i1 = 0; plhs[1] = mxCreateCellMatrix(16, 1); mxSetCell(plhs[1], i1++, mxCreateString("Energy")); @@ -355,8 +436,9 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) mxSetCell(plhs[1], i1++, mxCreateString("FieldScaling")); } } - else { - mexErrMsgIdAndTxt("AT:WrongArg","Needs 0 or 2 arguments"); + else + { + mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments"); } } #endif /* MATLAB_MEX_FILE */ diff --git a/atintegrators/ExactSectorBendQuantPass.c b/atintegrators/ExactSectorBendQuantPass.c new file mode 100644 index 0000000000..a033a45ea7 --- /dev/null +++ b/atintegrators/ExactSectorBendQuantPass.c @@ -0,0 +1,418 @@ +#include "atconstants.h" +#include "atelem.c" +#include "atlalib.c" +#include "atquantlib.c" +#include "driftkick.c" /* strthinkick.c */ +#include "exactbend.c" +#include "exactbendfringe.c" +#include "exactmultipolefringe.c" + +struct elem +{ + double Length; + double *PolynomA; + double *PolynomB; + int MaxOrder; + int NumIntSteps; + double BendingAngle; + double EntranceAngle; + double ExitAngle; + double Energy; + /* Optional fields */ + double Scaling; + int FringeBendEntrance; + int FringeBendExit; + int FringeQuadEntrance; + int FringeQuadExit; + double gK; + double *R1; + double *R2; + double *T1; + double *T2; + double *RApertures; + double *EApertures; + double *KickAngle; +}; + +static void ExactSectorBendQuant( + double *r, double le, double irho, double *A, double *B, int max_order, int num_int_steps, + double entrance_angle, double exit_angle, + int FringeBendEntrance, int FringeBendExit, + int FringeQuadEntrance, int FringeQuadExit, + double gK, + double *T1, double *T2, double *R1, double *R2, + double *RApertures, double *EApertures, + double *KickAngle, double scaling, double E0, + pcg32_random_t *rng, int num_particles) +{ + double SL = le / num_int_steps; + double L1 = SL * DRIFT1; + double L2 = SL * DRIFT2; + double K1 = SL * KICK1; + double K2 = SL * KICK2; + double B0 = B[0]; + double A0 = A[0]; + double qe = 1.60217733e-19; + double epsilon0 = 8.854187817e-12; + double clight = 2.99792458e8; + double emass = 510998.9461; /* electron mass in eV */ /* 9.10938188e-31; in kg*/ + double hbar = 1.054571726e-34; + double pi = 3.14159265358979; + double alpha0 = qe * qe / (4 * pi * epsilon0 * hbar * clight); + + if (KickAngle) + { /* Convert corrector component to polynomial coefficients */ + B[0] -= sin(KickAngle[0]) / le; + A[0] += 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, \ + irho, gK, 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) + */ + for (int c = 0; c < num_particles; c++) + { /* Loop over particles */ + double *r6 = r + 6 * c; + if (!atIsNaN(r6[0])) + { + /* Check for change of reference momentum */ + if (scaling != 1.0) + ATChangePRef(r6, scaling); + + /* misalignment at entrance */ + if (T1) + ATaddvv(r6, T1); + if (R1) + ATmultmv(r6, R1); + + /* Check physical apertures at the entrance of the magnet */ + if (RApertures) + checkiflostRectangularAp(r6, RApertures); + if (EApertures) + checkiflostEllipticalAp(r6, EApertures); + + Yrot(r6, entrance_angle); + if (FringeBendEntrance) + bend_fringe(r6, irho, gK); + if (FringeQuadEntrance) + multipole_fringe(r6, le, A, B, max_order, 1.0, 1); + bend_edge(r6, irho, -entrance_angle); + + for (int m = 0; m < num_int_steps; m++) + { /* Loop over slices */ + int i; + double ng, ec, de, energy, gamma, cstec, cstng; + double ds, rho, dxp, dyp; + int nph; + double p_norm = 1.0 / (1.0 + r6[4]); + double NormL1 = L1 * p_norm; + double NormL2 = L2 * p_norm; + double dpp0 = r6[4]; + double xp0 = r6[1] * p_norm; + double yp0 = r6[3] * p_norm; + double s0 = r6[5]; + + exact_bend(r6, irho, L1); + strthinkick(r6, A, B, K1, max_order); + exact_bend(r6, irho, L2); + strthinkick(r6, A, B, K2, max_order); + exact_bend(r6, irho, L2); + strthinkick(r6, A, B, K1, max_order); + exact_bend(r6, irho, L1); + + energy = dpp0 * E0 + E0; + + gamma = energy / emass; /* emass in eV */ + cstec = 3.0 * gamma * gamma * gamma * clight / (2.0) * hbar / qe; + cstng = 5.0 * sqrt(3.0) * alpha0 * gamma / (6.0); + + dxp = r6[1] * p_norm - xp0; + dyp = r6[3] * p_norm - yp0; + ds = r6[5] - s0; + + rho = (SL + ds) / sqrt(dxp * dxp + dyp * dyp); + + ng = cstng / rho * (SL + ds); + ec = cstec / rho; + + nph = atrandp_r(rng, ng); + + de = 0.0; + for (i = 0; i < nph; i++) + { + de = de + getEnergy(rng, ec); + }; + r6[4] = r6[4] - de / E0; + r6[1] = r6[1] * p_norm * (1 + r6[4]); + r6[3] = r6[3] * p_norm * (1 + r6[4]); + } + + /* Convert absolute path length to path lengthening */ + r6[5] -= le; + + /* edge focus */ + bend_edge(r6, irho, -exit_angle); + if (FringeQuadExit) + multipole_fringe(r6, le, A, B, max_order, -1.0, 1); + if (FringeBendExit) + bend_fringe(r6, -irho, gK); + Yrot(r6, exit_angle); + + /* Check physical apertures at the exit of the magnet */ + if (RApertures) + checkiflostRectangularAp(r6, RApertures); + if (EApertures) + checkiflostEllipticalAp(r6, EApertures); + + /* Misalignment at exit */ + if (R2) + ATmultmv(r6, R2); + if (T2) + ATaddvv(r6, T2); + + /* Check for change of reference momentum */ + 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) +ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, + double *r_in, int num_particles, struct parameters *Param) +{ + double irho, energy; + + if (!Elem) + { + double Length = atGetDouble(ElemData, "Length"); + check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); + check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); + check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); + check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); + check_error(); + double BendingAngle = atGetOptionalDouble(ElemData, "BendingAngle", 0.0); + check_error(); + double EntranceAngle = atGetDouble(ElemData, "EntranceAngle"); + check_error(); + double ExitAngle = atGetDouble(ElemData, "ExitAngle"); + check_error(); + /*optional fields*/ + double Energy = atGetOptionalDouble(ElemData, "Energy", Param->energy); + check_error(); + double Scaling = atGetOptionalDouble(ElemData, "FieldScaling", 1.0); + check_error(); + int FringeBendEntrance = atGetOptionalLong(ElemData, "FringeBendEntrance", 1); + check_error(); + 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 *R1 = atGetOptionalDoubleArray(ElemData, "R1"); + check_error(); + double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); + check_error(); + double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); + check_error(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); + check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); + check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); + check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); + check_error(); + + if (NumIntSteps == 0) + { + atError("NumIntSteps == 0 not allowed with radiation"); + check_error(); + } + + /* 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(); + } + + Elem = (struct elem *)atMalloc(sizeof(struct elem)); + Elem->Length = Length; + Elem->PolynomA = PolynomA; + Elem->PolynomB = PolynomB; + Elem->MaxOrder = MaxOrder; + Elem->NumIntSteps = NumIntSteps; + Elem->BendingAngle = BendingAngle; + Elem->EntranceAngle = EntranceAngle; + Elem->ExitAngle = ExitAngle; + /*optional fields*/ + Elem->Energy = Energy; + Elem->Scaling = Scaling; + Elem->FringeBendEntrance = FringeBendEntrance; + Elem->FringeBendExit = FringeBendExit; + Elem->FringeQuadEntrance = FringeQuadEntrance; + Elem->FringeQuadExit = FringeQuadExit; + Elem->gK = gK; + Elem->R1 = R1; + Elem->R2 = R2; + Elem->T1 = T1; + Elem->T2 = T2; + Elem->EApertures = EApertures; + Elem->RApertures = RApertures; + Elem->KickAngle = KickAngle; + } + irho = Elem->BendingAngle / Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + + ExactSectorBendQuant(r_in, Elem->Length, irho, + Elem->PolynomA, Elem->PolynomB, + Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, + Elem->FringeBendEntrance, Elem->FringeBendExit, + Elem->FringeQuadEntrance, Elem->FringeQuadExit, + Elem->gK, + Elem->T1, Elem->T2, Elem->R1, Elem->R2, + Elem->RApertures, Elem->EApertures, + Elem->KickAngle, Elem->Scaling, energy, + Param->thread_rng, num_particles); + return Elem; +} + +MODULE_DEF(ExactSectorBendQuantPass) /* Dummy module initialisation */ + +#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ + +#if defined(MATLAB_MEX_FILE) +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) +{ + if (nrhs >= 2) + { + double rest_energy = 0.0; + double charge = -1.0; + double irho; + double *r_in; + const mxArray *ElemData = prhs[0]; + int num_particles = mxGetN(prhs[1]); + if (mxGetM(prhs[1]) != 6) + mexErrMsgTxt("Second argument must be a 6 x N matrix"); + + double Length = atGetDouble(ElemData, "Length"); + check_error(); + double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); + check_error(); + double *PolynomB = atGetDoubleArray(ElemData, "PolynomB"); + check_error(); + int MaxOrder = atGetLong(ElemData, "MaxOrder"); + check_error(); + int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); + check_error(); + double BendingAngle = atGetOptionalDouble(ElemData, "BendingAngle", 0.0); + check_error(); + double EntranceAngle = atGetDouble(ElemData, "EntranceAngle"); + check_error(); + double ExitAngle = atGetDouble(ElemData, "ExitAngle"); + check_error(); + /*optional fields*/ + double Energy = atGetOptionalDouble(ElemData, "Energy", 0.0); + check_error(); + double Scaling = atGetOptionalDouble(ElemData, "FieldScaling", 1.0); + check_error(); + int FringeBendEntrance = atGetOptionalLong(ElemData, "FringeBendEntrance", 1); + check_error(); + 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 *R1 = atGetOptionalDoubleArray(ElemData, "R1"); + check_error(); + double *R2 = atGetOptionalDoubleArray(ElemData, "R2"); + check_error(); + double *T1 = atGetOptionalDoubleArray(ElemData, "T1"); + check_error(); + double *T2 = atGetOptionalDoubleArray(ElemData, "T2"); + check_error(); + double *EApertures = atGetOptionalDoubleArray(ElemData, "EApertures"); + check_error(); + double *RApertures = atGetOptionalDoubleArray(ElemData, "RApertures"); + check_error(); + double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); + check_error(); + + if (NumIntSteps == 0) + { + atError("NumIntSteps == 0 not allowed with radiation"); + check_error(); + } + if (nrhs > 2) + atProperties(prhs[2], &Energy, &rest_energy, &charge); + + /* ALLOCATE memory for the output array of the same size as the input */ + plhs[0] = mxDuplicateArray(prhs[1]); + irho = BendingAngle / Length; + r_in = mxGetDoubles(plhs[0]); + + ExactSectorBendQuant(r_in, Length, irho, PolynomA, PolynomB, + MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, + FringeBendEntrance, FringeBendExit, + FringeQuadEntrance, FringeQuadExit, + gK, + T1, T2, R1, R2, RApertures, EApertures, + KickAngle, Scaling, Energy, &pcg32_global, num_particles); + } + else if (nrhs == 0) + { + /* list of required fields */ + int i0 = 0; + plhs[0] = mxCreateCellMatrix(9, 1); + mxSetCell(plhs[0], i0++, mxCreateString("Length")); + mxSetCell(plhs[0], i0++, mxCreateString("PolynomA")); + mxSetCell(plhs[0], i0++, mxCreateString("PolynomB")); + mxSetCell(plhs[0], i0++, mxCreateString("MaxOrder")); + mxSetCell(plhs[0], i0++, mxCreateString("NumIntSteps")); + mxSetCell(plhs[0], i0++, mxCreateString("BendingAngle")); + mxSetCell(plhs[0], i0++, mxCreateString("EntranceAngle")); + mxSetCell(plhs[0], i0++, mxCreateString("ExitAngle")); + if (nlhs > 1) + { /* list of optional fields */ + int i1 = 0; + plhs[1] = mxCreateCellMatrix(14, 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("T1")); + mxSetCell(plhs[1], i1++, mxCreateString("T2")); + mxSetCell(plhs[1], i1++, mxCreateString("R1")); + mxSetCell(plhs[1], i1++, mxCreateString("R2")); + mxSetCell(plhs[1], i1++, mxCreateString("RApertures")); + mxSetCell(plhs[1], i1++, mxCreateString("EApertures")); + mxSetCell(plhs[1], i1++, mxCreateString("KickAngle")); + mxSetCell(plhs[1], i1++, mxCreateString("FieldScaling")); + } + } + else + { + mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments"); + } +} +#endif /* MATLAB_MEX_FILE */ From 752d25f178c977122cce5f2531be6e9f9a821fd9 Mon Sep 17 00:00:00 2001 From: kevin Date: Fri, 20 Feb 2026 09:04:35 +0100 Subject: [PATCH 3/4] correct SR emissions in ExactBendQuantPass, fix MATLAB build --- atintegrators/ExactRectangularBendQuantPass.c | 32 +++++++++---------- atintegrators/ExactSectorBendQuantPass.c | 10 +++--- 2 files changed, 22 insertions(+), 20 deletions(-) diff --git a/atintegrators/ExactRectangularBendQuantPass.c b/atintegrators/ExactRectangularBendQuantPass.c index 33b320f085..a554ec1195 100644 --- a/atintegrators/ExactRectangularBendQuantPass.c +++ b/atintegrators/ExactRectangularBendQuantPass.c @@ -36,18 +36,18 @@ struct elem double *KickAngle; }; -static void ExactRectangularBendQuant(double *r, double le, double bending_angle, - double *A, double *B, - int max_order, int num_int_steps, - double entrance_angle, double exit_angle, - int FringeBendEntrance, int FringeBendExit, - int FringeQuadEntrance, int FringeQuadExit, - double gK, double x0ref, double refdz, - double *T1, double *T2, - double *R1, double *R2, - double *RApertures, double *EApertures, - double *KickAngle, double scaling, double E0, - pcg32_random_t *rng, int num_particles) +static void ExactRectangularBendQuant( + double *r, double le, double bending_angle, + double *A, double *B, + int max_order, int num_int_steps, + double entrance_angle, double exit_angle, + int FringeBendEntrance, int FringeBendExit, + int FringeQuadEntrance, int FringeQuadExit, + double gK, double x0ref, double refdz, + double *T1, double *T2, double *R1, double *R2, + double *RApertures, double *EApertures, + double *KickAngle, double scaling, double E0, + pcg32_random_t *rng, int num_particles) { double irho = bending_angle / le; double phi2 = 0.5 * bending_angle; @@ -146,9 +146,9 @@ static void ExactRectangularBendQuant(double *r, double le, double bending_angle dyp = r6[3] * p_norm - yp0; ds = r6[5] - s0; - rho = (SL + ds) / sqrt(dxp * dxp + dyp * dyp); + rho = (ds + SL) / sqrt(dxp * dxp + dyp * dyp); - ng = cstng / rho * (SL + ds); + ng = cstng / rho * (ds); ec = cstec / rho; nph = atrandp_r(rng, ng); @@ -165,7 +165,7 @@ static void ExactRectangularBendQuant(double *r, double le, double bending_angle r6[0] -= x0ref; /* Convert absolute path length to path lengthening */ - r6[5] -= LR; /*(le + refdz);*/ + r6[5] -= (le + refdz); /* edge focus */ bend_edge(r6, irho, phi2 - exit_angle); @@ -323,7 +323,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double rest_energy = 0.0; double charge = -1.0; - double irho, Energy; + double irho; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactSectorBendQuantPass.c b/atintegrators/ExactSectorBendQuantPass.c index a033a45ea7..3167f35a6c 100644 --- a/atintegrators/ExactSectorBendQuantPass.c +++ b/atintegrators/ExactSectorBendQuantPass.c @@ -35,7 +35,9 @@ struct elem }; static void ExactSectorBendQuant( - double *r, double le, double irho, double *A, double *B, int max_order, int num_int_steps, + double *r, double le, double irho, + double *A, double *B, + int max_order, int num_int_steps, double entrance_angle, double exit_angle, int FringeBendEntrance, int FringeBendExit, int FringeQuadEntrance, int FringeQuadExit, @@ -129,13 +131,13 @@ static void ExactSectorBendQuant( cstec = 3.0 * gamma * gamma * gamma * clight / (2.0) * hbar / qe; cstng = 5.0 * sqrt(3.0) * alpha0 * gamma / (6.0); - dxp = r6[1] * p_norm - xp0; + dxp = r6[1] * p_norm - xp0 - irho * SL; dyp = r6[3] * p_norm - yp0; ds = r6[5] - s0; - rho = (SL + ds) / sqrt(dxp * dxp + dyp * dyp); + rho = (ds) / sqrt(dxp * dxp + dyp * dyp); - ng = cstng / rho * (SL + ds); + ng = cstng / rho * (ds); ec = cstec / rho; nph = atrandp_r(rng, ng); From e42d7397134ef313d362799579f1f964bd6eebd8 Mon Sep 17 00:00:00 2001 From: kevin Date: Fri, 20 Feb 2026 09:12:16 +0100 Subject: [PATCH 4/4] fix re-declaration error in MATLAB --- atintegrators/ExactRectangularBendQuantPass.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/atintegrators/ExactRectangularBendQuantPass.c b/atintegrators/ExactRectangularBendQuantPass.c index a554ec1195..e5d8fe4a32 100644 --- a/atintegrators/ExactRectangularBendQuantPass.c +++ b/atintegrators/ExactRectangularBendQuantPass.c @@ -379,7 +379,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) check_error(); double *KickAngle = atGetOptionalDoubleArray(ElemData, "KickAngle"); check_error(); - double irho = BendingAngle / Length; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); @@ -391,6 +391,7 @@ 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]); + irho = BendingAngle / Length; r_in = mxGetDoubles(plhs[0]); ExactRectangularBendQuant(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle,