diff --git a/atintegrators/AperturePass.c b/atintegrators/AperturePass.c index c8101784d2..9d03f4d2b1 100644 --- a/atintegrators/AperturePass.c +++ b/atintegrators/AperturePass.c @@ -46,7 +46,7 @@ MODULE_DEF(AperturePass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index 4c95aabae3..ed239a2987 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -51,7 +51,8 @@ void write_buffer(double *data, double *buffer, int datasize, int buffersize){ void BeamLoadingCavityPass(double *r_in,int num_particles,int nbunch, double *bunch_spos,double *bunch_currents, - double circumference,int nturn,struct elem *Elem) { + double circumference,int nturn,double energy, + struct elem *Elem) { /* * r_in - 6-by-N matrix of initial conditions reshaped into * 1-d array of 6*N elements @@ -63,7 +64,6 @@ void BeamLoadingCavityPass(double *r_in,int num_particles,int nbunch, long buffersize = Elem->buffersize; double normfact = Elem->normfact; double le = Elem->Length; - double energy = Elem->Energy; double rffreq = Elem->Frequency; double harmn = Elem->HarmNumber; double tlag = Elem->TimeLag; @@ -156,6 +156,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { double rl = Param->RingLength; + double energy; int nturn=Param->nturn; if (!Elem) { long nslice,nturns,blmode,cavitymode, buffersize; @@ -176,7 +177,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, /*attributes for RF cavity*/ Length=atGetDouble(ElemData,"Length"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); /*attributes for resonator*/ @@ -202,7 +202,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, vbeam_buffer=atGetDoubleArray(ElemData,"_vbeam_buffer"); check_error(); vbunch_buffer=atGetDoubleArray(ElemData,"_vbunch_buffer"); check_error(); /*optional attributes*/ - + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); z_cuts=atGetOptionalDoubleArray(ElemData,"ZCuts"); check_error(); int dimsth[] = {Param->nbunch*nslice*nturns, 4}; @@ -239,6 +239,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->vbeam_buffer = vbeam_buffer; Elem->vbunch_buffer = vbunch_buffer; } + energy = atEnergy(Param->energy, Elem->Energy); + if(num_particlesnbunch){ atError("Number of particles has to be greater or equal to the number of bunches."); }else if (num_particles%Param->nbunch!=0){ @@ -257,7 +259,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } #endif BeamLoadingCavityPass(r_in,num_particles,Param->nbunch,Param->bunch_spos, - Param->bunch_currents,rl,nturn,Elem); + Param->bunch_currents,rl,nturn,energy,Elem); return Elem; } @@ -267,10 +269,10 @@ MODULE_DEF(BeamLoadingCavityPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -293,7 +295,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double *vbunch_buffer; /*attributes for RF cavity*/ Length=atGetDouble(ElemData,"Length"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); /*attributes for resonator*/ @@ -319,6 +320,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) vbeam_buffer=atGetDoubleArray(ElemData,"_vbeam_buffer"); check_error(); vbunch_buffer=atGetDoubleArray(ElemData,"_vbunch_buffer"); check_error(); /*optional attributes*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); z_cuts=atGetOptionalDoubleArray(ElemData,"ZCuts"); check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); @@ -348,6 +350,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) Elem->vgen_buffer = vgen_buffer; Elem->vbeam_buffer = vbeam_buffer; Elem->vbunch_buffer = vbunch_buffer; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); /* ALLOCATE memory for the output array of the same size as the input */ @@ -356,7 +359,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double bspos = 0.0; double bcurr = 0.0; - BeamLoadingCavityPass(r_in,num_particles,1,&bspos,&bcurr,1,0,Elem); + BeamLoadingCavityPass(r_in,num_particles,1,&bspos,&bcurr,1,0,Energy,Elem); } else if (nrhs == 0) { /* return list of required fields */ diff --git a/atintegrators/BeamMomentsPass.c b/atintegrators/BeamMomentsPass.c index 3af8ab98b3..a8974dfe41 100644 --- a/atintegrators/BeamMomentsPass.c +++ b/atintegrators/BeamMomentsPass.c @@ -98,7 +98,7 @@ MODULE_DEF(BeamMomentsPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; diff --git a/atintegrators/BendLinearPass.c b/atintegrators/BendLinearPass.c index 5a96beb3da..bcf4aa44db 100644 --- a/atintegrators/BendLinearPass.c +++ b/atintegrators/BendLinearPass.c @@ -220,7 +220,7 @@ MODULE_DEF(BendLinearPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2 ) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, K, ByError, FringeInt1, FringeInt2, FullGap; double *R1, *R2, *T1, *T2; double *r_in; diff --git a/atintegrators/BndMPoleSymplectic4E2Pass.c b/atintegrators/BndMPoleSymplectic4E2Pass.c index f2a8c82e71..28329151d4 100644 --- a/atintegrators/BndMPoleSymplectic4E2Pass.c +++ b/atintegrators/BndMPoleSymplectic4E2Pass.c @@ -232,7 +232,7 @@ MODULE_DEF(BndMPoleSymplectic4E2Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double irho; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2; int MaxOrder, NumIntSteps; diff --git a/atintegrators/BndMPoleSymplectic4E2RadPass.c b/atintegrators/BndMPoleSymplectic4E2RadPass.c index 0756450e88..2c5ae5ec4b 100644 --- a/atintegrators/BndMPoleSymplectic4E2RadPass.c +++ b/atintegrators/BndMPoleSymplectic4E2RadPass.c @@ -202,7 +202,7 @@ void BndMPoleSymplectic4E2RadPass(double *r, double le, double irho, double *A, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { - double irho; + double irho, energy; if (!Elem) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, Energy, FringeInt1, FringeInt2; @@ -216,8 +216,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1",0); check_error(); @@ -258,13 +258,15 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + BndMPoleSymplectic4E2RadPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeInt1, Elem->FringeInt2, Elem->FullGap, Elem->h1, Elem->h2, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -275,7 +277,9 @@ MODULE_DEF(BndMPoleSymplectic4E2RadPass) /* Dummy module initialisation * #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; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps; @@ -293,8 +297,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap", 0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1", 0); check_error(); @@ -309,6 +313,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); irho = BendingAngle/Length; + 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]); diff --git a/atintegrators/BndMPoleSymplectic4Pass.c b/atintegrators/BndMPoleSymplectic4Pass.c index 91818201eb..1aef9d8e26 100644 --- a/atintegrators/BndMPoleSymplectic4Pass.c +++ b/atintegrators/BndMPoleSymplectic4Pass.c @@ -214,7 +214,7 @@ MODULE_DEF(BndMPoleSymplectic4Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/BndMPoleSymplectic4QuantPass.c b/atintegrators/BndMPoleSymplectic4QuantPass.c index b542a70f94..b4bb1304fc 100644 --- a/atintegrators/BndMPoleSymplectic4QuantPass.c +++ b/atintegrators/BndMPoleSymplectic4QuantPass.c @@ -180,7 +180,7 @@ void BndMPoleSymplectic4QuantPass(double *r, double le, double irho, double *A, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { - double irho; + double irho, energy; if (!Elem) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; @@ -195,8 +195,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -245,6 +245,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + BndMPoleSymplectic4QuantPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, @@ -253,7 +255,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->fringeIntM0, Elem->fringeIntP0, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, + Elem->KickAngle, Elem->Scaling, energy, Param->thread_rng, num_particles); return Elem; } @@ -265,12 +267,14 @@ MODULE_DEF(BndMPoleSymplectic4QuantPass) /* Dummy module initialisation * #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; + double rest_energy = 0.0; + double charge = -1.0; double irho; double *r_in; const mxArray *ElemData = prhs[0]; @@ -285,8 +289,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -305,10 +309,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); irho = BendingAngle/Length; + 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]); r_in = mxGetDoubles(plhs[0]); + BndMPoleSymplectic4QuantPass(r_in, Length, irho, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index acb47b6603..7767416483 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -128,7 +128,7 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { - double irho; + double irho, energy; if (!Elem) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; @@ -143,8 +143,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -193,6 +193,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + BndMPoleSymplectic4RadPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, @@ -201,7 +203,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->fringeIntM0, Elem->fringeIntP0, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -212,7 +214,9 @@ MODULE_DEF(BndMPoleSymplectic4RadPass) /* Dummy module initialisation */ #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 Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, @@ -232,8 +236,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -252,10 +256,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); irho = BendingAngle/Length; + 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]); r_in = mxGetDoubles(plhs[0]); + BndMPoleSymplectic4RadPass(r_in, Length, irho, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/BndStrMPoleSymplectic4Pass.c b/atintegrators/BndStrMPoleSymplectic4Pass.c index 6f79d9e548..b626574b83 100644 --- a/atintegrators/BndStrMPoleSymplectic4Pass.c +++ b/atintegrators/BndStrMPoleSymplectic4Pass.c @@ -271,7 +271,7 @@ MODULE_DEF(BndStrMPoleSymplectic4Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, X0ref, ByError, RefDZ; int MaxOrder, NumIntSteps; diff --git a/atintegrators/CavityPass.c b/atintegrators/CavityPass.c index fe3e202f63..bce190f57a 100644 --- a/atintegrators/CavityPass.c +++ b/atintegrators/CavityPass.c @@ -63,12 +63,14 @@ void CavityPass(double *r_in, double le, double nv, double freq, double lag, dou ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy = Param->energy; if (!Elem) { double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag; Length=atGetDouble(ElemData,"Length"); check_error(); Voltage=atGetDouble(ElemData,"Voltage"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); + /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); @@ -79,7 +81,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->TimeLag=TimeLag; Elem->PhaseLag=PhaseLag; } - CavityPass(r_in,Elem->Length,Elem->Voltage/Elem->Energy,Elem->Frequency, + if (energy == 0.0) energy = Elem->Energy; + + CavityPass(r_in,Elem->Length,Elem->Voltage/energy,Elem->Frequency, Elem->TimeLag,Elem->PhaseLag,num_particles); return Elem; } @@ -92,20 +96,25 @@ MODULE_DEF(CavityPass) /* Dummy module initialisation */ 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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); double Length=atGetDouble(ElemData,"Length"); double Voltage=atGetDouble(ElemData,"Voltage"); - double Energy=atGetDouble(ElemData,"Energy"); double Frequency=atGetDouble(ElemData,"Frequency"); + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); double TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); double PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + CavityPass(r_in,Length,Voltage/Energy,Frequency,TimeLag,PhaseLag,num_particles); } else if (nrhs == 0) { /* return list of required fields */ diff --git a/atintegrators/ChangePRefPass.c b/atintegrators/ChangePRefPass.c index 1810a8379b..0161f9f66e 100644 --- a/atintegrators/ChangePRefPass.c +++ b/atintegrators/ChangePRefPass.c @@ -43,7 +43,7 @@ MODULE_DEF(ChangePRefPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/CorrectorPass.c b/atintegrators/CorrectorPass.c index 5e58cb46fe..f5470c986c 100644 --- a/atintegrators/CorrectorPass.c +++ b/atintegrators/CorrectorPass.c @@ -132,7 +132,7 @@ MODULE_DEF(CorrectorPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/DeltaQPass.c b/atintegrators/DeltaQPass.c index e92f340820..ef4c53dca4 100644 --- a/atintegrators/DeltaQPass.c +++ b/atintegrators/DeltaQPass.c @@ -141,7 +141,7 @@ MODULE_DEF(DeltaQPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/DriftPass.c b/atintegrators/DriftPass.c index bdbeb9fbbd..4935dc5984 100644 --- a/atintegrators/DriftPass.c +++ b/atintegrators/DriftPass.c @@ -91,7 +91,7 @@ MODULE_DEF(DriftPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/EAperturePass.c b/atintegrators/EAperturePass.c index 3277961cc2..66ca6e60fd 100644 --- a/atintegrators/EAperturePass.c +++ b/atintegrators/EAperturePass.c @@ -44,7 +44,7 @@ MODULE_DEF(EAperturePass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/EnergyLossRadPass.c b/atintegrators/EnergyLossRadPass.c index 9d639584b3..9dcfa840c9 100644 --- a/atintegrators/EnergyLossRadPass.c +++ b/atintegrators/EnergyLossRadPass.c @@ -62,7 +62,7 @@ MODULE_DEF(EnergyLossRadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Check if the number of input arguments is correct. */ - if (nrhs == 2) { + if (nrhs >= 2) { /* Get the input arguments. */ double *r_in; const mxArray *ElemData = prhs[0]; diff --git a/atintegrators/ExactDriftPass.c b/atintegrators/ExactDriftPass.c index 2c04e14a58..5f6df8af0d 100644 --- a/atintegrators/ExactDriftPass.c +++ b/atintegrators/ExactDriftPass.c @@ -83,7 +83,7 @@ MODULE_DEF(ExactDriftPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactMultipolePass.c b/atintegrators/ExactMultipolePass.c index de1ea59795..f0e3f8030a 100644 --- a/atintegrators/ExactMultipolePass.c +++ b/atintegrators/ExactMultipolePass.c @@ -160,7 +160,7 @@ MODULE_DEF(ExactMultipolePass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index 248d4c8e90..ce3094005d 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -106,6 +106,7 @@ static void multipole_pass( 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(); @@ -148,12 +149,14 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, 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->Energy, Elem->Scaling, num_particles); + Elem->KickAngle, energy, Elem->Scaling, num_particles); return Elem; } @@ -163,7 +166,9 @@ MODULE_DEF(ExactMultipoleRadPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -174,7 +179,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int MaxOrder = atGetLong(ElemData, "MaxOrder"); check_error(); int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); 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(); @@ -189,6 +194,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (NumIntSteps <= 0) { atError("NumIntSteps must be positive"); 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]); r_in = mxGetDoubles(plhs[0]); diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index d1460dd90e..f8ab7a8596 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -221,7 +221,7 @@ MODULE_DEF(ExactRectBendPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 26b6d0bb5e..3272c16652 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -139,6 +139,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, 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(); @@ -197,6 +198,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, @@ -205,7 +208,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->gK,Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -216,7 +219,9 @@ MODULE_DEF(ExactRectBendRadPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -231,7 +236,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); 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(); @@ -251,10 +256,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) 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]); r_in = mxGetDoubles(plhs[0]); + ExactRectangularBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 82bc6d6eb5..a91fa447b8 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -214,7 +214,7 @@ MODULE_DEF(ExactRectangularBendPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index 99e353c931..f409888775 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -141,6 +141,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, 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(); @@ -199,6 +200,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, @@ -206,7 +209,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->gK,Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -217,7 +220,9 @@ MODULE_DEF(ExactRectangularBendRadPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -232,7 +237,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); 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(); @@ -252,6 +257,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (NumIntSteps <= 0) { atError("NumIntSteps must be positive"); 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]); diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 64254cec73..44b9108349 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -206,7 +206,7 @@ MODULE_DEF(ExactSectorBendPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index df7a1c3f9b..42469e4f63 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -127,6 +127,7 @@ static void ExactSectorBendRad(double *r, double le, double bending_angle, 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(); @@ -181,6 +182,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + ExactSectorBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, @@ -189,7 +192,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->gK, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -200,7 +203,9 @@ MODULE_DEF(ExactSectorBendRadPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -215,7 +220,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); 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(); @@ -233,10 +238,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) 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]); r_in = mxGetDoubles(plhs[0]); + ExactSectorBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 26b5e945c6..9d783f56c0 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -219,7 +219,7 @@ MODULE_DEF(GWigSymplecticPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index 802aab5844..af306484a3 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -178,6 +178,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double *R1, *R2, *T1, *T2; double *By, *Bx; @@ -185,7 +186,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, int Nstep, Nmeth; int NHharm, NVharm; - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); Ltot = atGetDouble(ElemData, "Length"); check_error(); Lw = atGetDouble(ElemData, "Lw"); check_error(); Bmax = atGetDouble(ElemData, "Bmax"); check_error(); @@ -196,6 +196,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, By = atGetDoubleArray(ElemData, "By"); check_error(); Bx = atGetDoubleArray(ElemData, "Bx"); check_error(); /* Optional fields */ + Energy = atGetOptionalDouble(ElemData, "Energy",Param->energy); check_error(); R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); @@ -218,7 +219,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T1=T1; Elem->T2=T2; } - GWigSymplecticRadPass(r_in, Elem->Energy, Elem->Length, Elem->Lw, + energy = atEnergy(Param->energy, Elem->Energy); + + GWigSymplecticRadPass(r_in, energy, Elem->Length, Elem->Lw, Elem->Bmax, Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm, Elem->By, Elem->Bx, Elem->T1, Elem->T2, Elem->R1, Elem->R2, num_particles); @@ -236,7 +239,9 @@ MODULE_DEF(GWigSymplecticRadPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -245,8 +250,8 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs double Ltot, Lw, Bmax, Energy; int Nstep, Nmeth; int NHharm, NVharm; + if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); - Energy = atGetDouble(ElemData, "Energy"); check_error(); Ltot = atGetDouble(ElemData, "Length"); check_error(); Lw = atGetDouble(ElemData, "Lw"); check_error(); Bmax = atGetDouble(ElemData, "Bmax"); check_error(); @@ -257,11 +262,13 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs By = atGetDoubleArray(ElemData, "By"); check_error(); Bx = atGetDoubleArray(ElemData, "Bx"); check_error(); /* Optional fields */ + Energy = atGetOptionalDouble(ElemData, "Energy",0.0); check_error(); R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); - if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); + 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]); r_in = mxGetDoubles(plhs[0]); diff --git a/atintegrators/IdTablePass.c b/atintegrators/IdTablePass.c index 1739e85c9b..c697f38a71 100644 --- a/atintegrators/IdTablePass.c +++ b/atintegrators/IdTablePass.c @@ -150,7 +150,7 @@ MODULE_DEF(IdTablePass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/IdentityPass.c b/atintegrators/IdentityPass.c index 1ba8085d4a..791da90fee 100644 --- a/atintegrators/IdentityPass.c +++ b/atintegrators/IdentityPass.c @@ -74,7 +74,7 @@ MODULE_DEF(IdentityPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ImpedanceTablePass.c b/atintegrators/ImpedanceTablePass.c index d2bc29b2de..183509c1ba 100755 --- a/atintegrators/ImpedanceTablePass.c +++ b/atintegrators/ImpedanceTablePass.c @@ -302,7 +302,7 @@ MODULE_DEF(ImpedanceTablePass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/Matrix66Pass.c b/atintegrators/Matrix66Pass.c index 11e126bfba..8dd7f89efe 100644 --- a/atintegrators/Matrix66Pass.c +++ b/atintegrators/Matrix66Pass.c @@ -60,7 +60,7 @@ MODULE_DEF(Matrix66Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/MatrixTijkPass.c b/atintegrators/MatrixTijkPass.c index 81789fbbbe..e7a10aa9a9 100644 --- a/atintegrators/MatrixTijkPass.c +++ b/atintegrators/MatrixTijkPass.c @@ -96,7 +96,7 @@ MODULE_DEF(MatrixTijkPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/QuantDiffPass.c b/atintegrators/QuantDiffPass.c index 3920dde9df..0399c3b28d 100644 --- a/atintegrators/QuantDiffPass.c +++ b/atintegrators/QuantDiffPass.c @@ -68,7 +68,7 @@ MODULE_DEF(QuantDiffPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double* r_in; const mxArray* ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index c2ac702f6f..60bf4c72a2 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -37,12 +37,14 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, { int nturn=Param->nturn; double T0=Param->T0; + double energy = Param->energy; if (!Elem) { double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag; Length=atGetDouble(ElemData,"Length"); check_error(); Voltage=atGetDouble(ElemData,"Voltage"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); + /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); @@ -54,7 +56,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->TimeLag=TimeLag; Elem->PhaseLag=PhaseLag; } - RFCavityPass(r_in, Elem->Length, Elem->Voltage/Elem->Energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, + if (energy == 0.0) energy = Elem->Energy; + + RFCavityPass(r_in, Elem->Length, Elem->Voltage/energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, Elem->PhaseLag, nturn, T0, num_particles); return Elem; } @@ -67,23 +71,28 @@ MODULE_DEF(RFCavityPass) /* Dummy module initialisation */ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if(nrhs == 2) - { + if (nrhs >= 2) { + double Energy = 0.0; + 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"); double Voltage=atGetDouble(ElemData,"Voltage"); - double Energy=atGetDouble(ElemData,"Energy"); + Energy=atGetOptionalDouble(ElemData,"Energy",Energy); double Frequency=atGetDouble(ElemData,"Frequency"); double TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); double PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); double T0=1.0/Frequency; /* Does not matter since nturns == 0 */ double HarmNumber=round(Frequency*T0); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); } diff --git a/atintegrators/SimpleQuantDiffPass.c b/atintegrators/SimpleQuantDiffPass.c index 988fad0be8..6887f80363 100644 --- a/atintegrators/SimpleQuantDiffPass.c +++ b/atintegrators/SimpleQuantDiffPass.c @@ -86,7 +86,7 @@ MODULE_DEF(SimpleQuantDiffPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/SimpleRadiationPass.c b/atintegrators/SimpleRadiationPass.c index 2bed3b5155..3443aff029 100644 --- a/atintegrators/SimpleRadiationPass.c +++ b/atintegrators/SimpleRadiationPass.c @@ -110,7 +110,7 @@ MODULE_DEF(SimpleRadiationPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/SliceMomentsPass.c b/atintegrators/SliceMomentsPass.c index 8ba718a788..b388ce2dec 100644 --- a/atintegrators/SliceMomentsPass.c +++ b/atintegrators/SliceMomentsPass.c @@ -186,7 +186,7 @@ MODULE_DEF(SliceMomentsPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; diff --git a/atintegrators/SolenoidLinearPass.c b/atintegrators/SolenoidLinearPass.c index 5cd1698cbf..95fc10547d 100644 --- a/atintegrators/SolenoidLinearPass.c +++ b/atintegrators/SolenoidLinearPass.c @@ -100,7 +100,7 @@ MODULE_DEF(SolenoidLinearPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2 ) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/StrMPoleSymplectic4Pass.c b/atintegrators/StrMPoleSymplectic4Pass.c index 6105eee811..9ac88ad4f8 100644 --- a/atintegrators/StrMPoleSymplectic4Pass.c +++ b/atintegrators/StrMPoleSymplectic4Pass.c @@ -172,7 +172,7 @@ MODULE_DEF(StrMPoleSymplectic4Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/StrMPoleSymplectic4QuantPass.c b/atintegrators/StrMPoleSymplectic4QuantPass.c index 83fa35448b..e641f151ff 100644 --- a/atintegrators/StrMPoleSymplectic4QuantPass.c +++ b/atintegrators/StrMPoleSymplectic4QuantPass.c @@ -159,6 +159,7 @@ void StrMPoleSymplectic4QuantPass(double *r, double le, double *A, double *B, 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, Energy, Scaling; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; @@ -168,8 +169,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -204,6 +205,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + StrMPoleSymplectic4QuantPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->FringeQuadEntrance, Elem->FringeQuadExit, @@ -211,7 +214,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, - Elem->Energy, Param->thread_rng, num_particles); + energy, Param->thread_rng, num_particles); return Elem; } @@ -222,7 +225,9 @@ MODULE_DEF(StrMPoleSymplectic4QuantPass) /* Dummy module initialisation * #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -236,8 +241,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -250,10 +255,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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]); r_in = mxGetDoubles(plhs[0]); + StrMPoleSymplectic4QuantPass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit, diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index 88fb00470f..6b9971355a 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/atintegrators/StrMPoleSymplectic4RadPass.c @@ -109,6 +109,7 @@ void StrMPoleSymplectic4RadPass(double *r, double le, double *A, double *B, 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, Energy, Scaling; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; @@ -118,8 +119,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -154,13 +155,15 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + StrMPoleSymplectic4RadPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->FringeQuadEntrance, Elem->FringeQuadExit, Elem->fringeIntM0, Elem->fringeIntP0, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -171,11 +174,13 @@ MODULE_DEF(StrMPoleSymplectic4RadPass) /* Dummy module initialisation */ #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 *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); - double Length, Energy, Scaling; + double Length, Scaling, Energy; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; if (mxGetM(prhs[1]) != 6) mexErrMsgTxt("Second argument must be a 6 x N matrix"); @@ -185,8 +190,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -199,10 +204,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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]); r_in = mxGetDoubles(plhs[0]); + StrMPoleSymplectic4RadPass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit, diff --git a/atintegrators/TestRandomPass.c b/atintegrators/TestRandomPass.c index 8585784686..40dea1c2cc 100644 --- a/atintegrators/TestRandomPass.c +++ b/atintegrators/TestRandomPass.c @@ -64,7 +64,7 @@ MODULE_DEF(TestRandomPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ThinMPolePass.c b/atintegrators/ThinMPolePass.c index e83ec304c9..af8e7dcfd9 100644 --- a/atintegrators/ThinMPolePass.c +++ b/atintegrators/ThinMPolePass.c @@ -129,7 +129,7 @@ MODULE_DEF(ThinMPolePass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/VariableThinMPolePass.c b/atintegrators/VariableThinMPolePass.c index 0c25e9afac..8ba2f88fa7 100644 --- a/atintegrators/VariableThinMPolePass.c +++ b/atintegrators/VariableThinMPolePass.c @@ -186,7 +186,7 @@ MODULE_DEF(VariableThinMPolePass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double* r_in; const mxArray* ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/WakeFieldPass.c b/atintegrators/WakeFieldPass.c index 0e705f1af2..703f2fac34 100755 --- a/atintegrators/WakeFieldPass.c +++ b/atintegrators/WakeFieldPass.c @@ -165,7 +165,7 @@ MODULE_DEF(WakeFieldPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/WiggLinearPass.c b/atintegrators/WiggLinearPass.c index 20976281da..beeb13659c 100644 --- a/atintegrators/WiggLinearPass.c +++ b/atintegrators/WiggLinearPass.c @@ -125,7 +125,7 @@ MODULE_DEF(WiggLinearPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2 ) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 6fde043485..04af7ccca5 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -6,6 +6,8 @@ #define ATELEM_C #include "atcommon.h" +#include "attypes.h" +#include "atconstants.h" /*----------------------------------------------------*/ /* For the integrator code */ @@ -70,6 +72,29 @@ typedef mxArray atElem; #define atError(...) mexErrMsgIdAndTxt("AT:PassError", __VA_ARGS__) #define atWarning(...) mexWarnMsgIdAndTxt("AT:PassWarning", __VA_ARGS__) #define atPrintf(...) mexPrintf(__VA_ARGS__) +#include "ringproperties.c" + +double atEnergy(double ringenergy, double elemenergy) +{ + if (ringenergy!=0.0) + return ringenergy; + else + if (elemenergy!=0.0) + return elemenergy; + else { + atError("Energy not defined."); + return 0.0; /* Never reached but makes the compiler happy */ + } +} + +double atGamma(double ringenergy, double elemenergy, double rest_energy) +{ + double energy = atEnergy(ringenergy, elemenergy); + if (rest_energy == 0.0) + return 1.0E-9 * energy / __E0; + else + return energy / rest_energy; +} static mxArray *get_field(const mxArray *pm, const char *fieldname) { @@ -180,6 +205,8 @@ typedef PyObject atElem; #define atError(...) return (struct elem *) PyErr_Format(PyExc_ValueError, __VA_ARGS__) #define atWarning(...) if (PyErr_WarnFormat(PyExc_RuntimeWarning, 0, __VA_ARGS__) != 0) return NULL #define atPrintf(...) PySys_WriteStdout(__VA_ARGS__) +#define atEnergy(ringenergy,elemenergy) (ringenergy) +#define atGamma(ringenergy,elemenergy,rest_energy) ((rest_energy) == 0.0 ? 1.0E-9*(ringenergy)/__E0 : (ringenergy)/(rest_energy)) static int array_imported = 0; @@ -327,7 +354,6 @@ static double *atGetOptionalDoubleArray(const PyObject *element, char *name) #endif /* defined(PYAT) */ #if defined(PYAT) || defined(MATLAB_MEX_FILE) -#include "attypes.h" #ifdef __cplusplus #define C_LINK extern "C" diff --git a/atintegrators/ringproperties.c b/atintegrators/ringproperties.c new file mode 100644 index 0000000000..baabef3dbf --- /dev/null +++ b/atintegrators/ringproperties.c @@ -0,0 +1,34 @@ + +static double atGetOptionalDoubleProp(const mxArray *obj, const char *fieldname, double default_value) +{ + mxArray *field=mxGetProperty(obj, 0, fieldname); + return (field) ? mxGetScalar(field) : default_value; +} + +static void atParticle(const mxArray *opts, double *rest_energy, double *charge) +{ + const mxArray *part = mxGetField(opts, 0, "Particle"); + if (part) { + if (mxIsClass(part, "atparticle")) { /* OK */ + *rest_energy = atGetOptionalDoubleProp(part, "rest_energy", 0.0); + *charge = atGetOptionalDoubleProp(part, "charge", -1.0); + } + else { /* particle is not a Particle object */ + mexErrMsgIdAndTxt("Atpass:WrongParameter","Particle must be an 'atparticle' object"); + } + } +} + +static void atProperties(const mxArray *opts, double *energy, double *rest_energy, double *charge) +{ + mxArray *field; + if (!mxIsStruct(opts)) { + mexErrMsgIdAndTxt("Atpass:WrongParameter","ring properties must be a struct"); + } + field = mxGetField(opts, 0, "Energy"); + if (field) { + double ener = mxGetScalar(field); + if (ener != 0.0) *energy = ener; + atParticle(opts, rest_energy, charge); + } +} diff --git a/atmat/atphysics/LinearOptics/findelemm44.m b/atmat/atphysics/LinearOptics/findelemm44.m index 3e5976d038..d67278714e 100644 --- a/atmat/atphysics/LinearOptics/findelemm44.m +++ b/atmat/atphysics/LinearOptics/findelemm44.m @@ -12,10 +12,19 @@ % The transverse matrix is momentum-dependent, % the 5-th component of ORBITIN is used as the DP value % +% M66=FINDELEMM44(...,'Energy',ENERGY) +% Use ENERGY and ignore the 'Energy' field of elements +% +% M66=FINDELEMM44(...,'Particle',PARTICLE) +% Use PARTICLE (default is relativistic) +% +% % See also FINDELEMM66 [XYStep,varargs]=getoption(varargin,'XYStep'); [R0,varargs]=getoption(varargs,'orbit',zeros(6,1)); +[energy,varargs]=getoption(varargs,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',[]); [MethodName,R0]=getargs(varargs,ELEM.PassMethod,R0); % Build a diagonal matrix of initial conditions @@ -25,6 +34,6 @@ % Add to the orbit_in RIN = R0 + [D4 -D4]; % Propagate through the element -ROUT = feval(MethodName,ELEM,RIN); +ROUT=elempass(ELEM,RIN,'PassMethod',MethodName,'Energy',energy,'Particle',particle); % Calculate numerical derivative M44 = ((ROUT(1:4,1:4)-ROUT(1:4,5:8))./scaling); diff --git a/atmat/atphysics/LinearOptics/findelemm66.m b/atmat/atphysics/LinearOptics/findelemm66.m index 7fb7b361f6..4cbf32e588 100644 --- a/atmat/atphysics/LinearOptics/findelemm66.m +++ b/atmat/atphysics/LinearOptics/findelemm66.m @@ -9,11 +9,19 @@ % M66=FINDELEMM66(...,'orbit',ORBITIN) % ORBITIN - 6-by-1 phase space coordinates at the entrance % (default: zeros(6,1)) +% +% M66=FINDELEMM66(...,'Energy',ENERGY) +% Use ENERGY and ignore the 'Energy' field of elements +% +% M66=FINDELEMM66(...,'Particle',PARTICLE) +% Use PARTICLE (default is relativistic) % % See also FINDELEMM44 [XYStep,varargs]=getoption(varargin,'XYStep'); [R0,varargs]=getoption(varargs,'orbit',zeros(6,1)); +[energy,varargs]=getoption(varargs,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',[]); [MethodName,R0]=getargs(varargs,ELEM.PassMethod,R0); % Build a diagonal matrix of initial conditions @@ -23,6 +31,6 @@ % Add to the orbit_in RIN = R0 + [D6, -D6]; % Propagate through the element -ROUT = feval(MethodName,ELEM,RIN); +ROUT=elempass(ELEM,RIN,'PassMethod',MethodName,'Energy',energy,'Particle',particle); % Calculate numerical derivative M66 = (ROUT(:,1:6)-ROUT(:,7:12))./scaling; diff --git a/atmat/atphysics/ParameterSummaryFunctions/atx.m b/atmat/atphysics/ParameterSummaryFunctions/atx.m index 1278ce643d..d6f2eb9682 100644 --- a/atmat/atphysics/ParameterSummaryFunctions/atx.m +++ b/atmat/atphysics/ParameterSummaryFunctions/atx.m @@ -192,7 +192,7 @@ if has_cavity try - [envelope,espread,blength,m,T]=ohmienvelope(ron,radindex,refpts); + [envelope,espread,blength,m,T]=ohmienvelope(ron,radindex,refpts,energy); [tns,chi]=atdampingrates(m); fs=abs(tns(3))*cell_revfreq; alpha=chi*cell_revfreq; diff --git a/atmat/atphysics/Radiation/atdiffmat.m b/atmat/atphysics/Radiation/atdiffmat.m new file mode 100644 index 0000000000..0dc61611cb --- /dev/null +++ b/atmat/atphysics/Radiation/atdiffmat.m @@ -0,0 +1,51 @@ +function [BCUM,Batbeg] = atdiffmat(ring,energy,varargin) +%quantumDiff Compute the radiation-diffusion matrix +% +%[BCUM,BS]=ATDIFFMAT(RING,ENERGY) +% RING: Closed ring AT structure, containing radiative elements and +% RF cavity. Radiative elements are identified by a +% PassMethod ending with 'RadPass'. +% ENERGY: Lattice energy +% +% BCUM: Cumulated diffusion matrix +% BS: Cumulative diffusion matrix at the beginning of each element +% +%[BCUM,BS]=ATDIFFMAT(...,'orbit',ORBITIN) +% ORBITIN: Initial 6-D closed orbit. +% In this mode, RING may be a section of a ring. + +NumElements=length(ring); + +[orb0,varargs]=getoption(varargin, 'orbit', []); %#ok +if isempty(orb0) + orbit=num2cell(findorbit6(ring, 1:NumElements),1)'; +else + orbit=num2cell(linepass(ring, orb0, 1:NumElements),1)'; +end + +zr=zeros(6,6); +% Calculate cumulative Radiation-Diffusion matrix for the ring +BCUM = zeros(6,6); +% Batbeg{i} is the cumulative diffusion matrix from +% 0 to the beginning of the i-th element +Batbeg=[zr;cellfun(@cumulb,ring,orbit,'UniformOutput',false)]; + + function btx=cumulb(elem,orbit) + if endsWith(elem.PassMethod,'RadPass') + if isfield(elem,'Bmax') + b=FDW(elem,orbit,energy); % For wigglers + else + b=findmpoleraddiffmatrix(elem, orbit, energy); % For other elements + end + else + b=zr; + end + % Calculate 6-by-6 linear transfer matrix in each element + % near the equilibrium orbit + m=findelemm66(elem,elem.PassMethod,orbit); + % Cumulative diffusion matrix of the entire ring + BCUM = m*BCUM*m' + b; + btx=BCUM; + end +end + diff --git a/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c b/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c index d3f203d35a..92e69de50a 100644 --- a/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c +++ b/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c @@ -522,7 +522,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxCreateDoubleMatrix(6,6,mxREAL); bdiff = mxGetDoubles(plhs[0]); for (i=0; i<36; i++) bdiff[i]=0.0; - + diffmatrix(mxElem, orb, energy, bdiff); } #endif /*MATLAB_MEX_FILE*/ diff --git a/atmat/atphysics/Radiation/ohmienvelope.m b/atmat/atphysics/Radiation/ohmienvelope.m index fb72a4245a..e1d5f4e183 100644 --- a/atmat/atphysics/Radiation/ohmienvelope.m +++ b/atmat/atphysics/Radiation/ohmienvelope.m @@ -1,4 +1,4 @@ -function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts) +function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts,energy) %#ok %OHMIENVELOPE calculates equilibrium beam envelope in a % circular accelerator using Ohmi's beam envelope formalism [1]. % [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994) @@ -33,37 +33,18 @@ check_6d(ring,true,'strict',0); NumElements = length(ring); +if nargin<4, energy=atGetRingProperties(ring, 'Energy'); end if nargin<3, refpts=1; end -% Erase wigglers from the radiative element list. -% Diffusion matrix to be computed with separate FDW function. -Wig=atgetcells(ring,'Bmax'); -radindex = radindex & ~Wig; - -[mring, ms, orbit] = findm66(ring,1:NumElements+1); -mt=squeeze(num2cell(ms,[1 2])); -orb=num2cell(orbit,1)'; - -zr={zeros(6,6)}; -B=zr(ones(NumElements,1)); % B{i} is the diffusion matrix of the i-th element - -% calculate Radiation-Diffusion matrix B for elements with radiation -B(radindex)=cellfun(@findmpoleraddiffmatrix,... - ring(radindex),orb(radindex),'UniformOutput',false); -B(Wig)=cellfun(@FDW,... - ring(Wig),orb(Wig),'UniformOutput',false); - -% Calculate cumulative Radiation-Diffusion matrix for the ring -BCUM = zeros(6,6); -% Batbeg{i} is the cumulative diffusion matrix from -% 0 to the beginning of the i-th element -Batbeg=[zr;cellfun(@cumulb,ring,orb(1:end-1),B,'UniformOutput',false)]; +orb0=findorbit(ring); +[BCUM,Batbeg]=atdiffmat(ring,energy,'orbit',orb0); +[mring, ms, orbit] = findm66(ring,1:NumElements+1,'orbit',orb0); % ------------------------------------------------------------------------ % Equation for the moment matrix R is % R = MRING*R*MRING' + BCUM; % We rewrite it in the form of Sylvester-Lyapunov equation -% to use MATLAB's SYLVERTER function: +% to use MATLAB's SYLVESTER function: % AA*R + R*BB = CC % where % AA = inv(MRING) @@ -77,8 +58,9 @@ R = sylvester(AA,BB,CC); % Envelope matrix at the ring entrance rmsdp = sqrt(R(5,5)); % R.M.S. energy spread -rmsbl = sqrt(R(6,6)); % R.M.S. bunch length +rmsbl = sqrt(R(6,6)); % R.M.S. bunch lenght +mt=squeeze(num2cell(ms,[1 2])); [rr,tt,ss]=cellfun(@propag,mt(refpts),Batbeg(refpts),'UniformOutput',false); envelope=struct('R',rr,'Sigma',ss,'Tilt',tt); @@ -88,15 +70,6 @@ if nout>=2, varargout{2}=ms(:,:,refpts); end if nout>=1, varargout{1}=mring; end - function btx=cumulb(elem,orbit,b) - % Calculate 6-by-6 linear transfer matrix in each element - % near the equilibrium orbit - m=findelemm66(elem,elem.PassMethod,orbit); - % Cumulative diffusion matrix of the entire ring - BCUM = m*BCUM*m' + b; - btx=BCUM; - end - function [r,tilt,sigma]=propag(m,cumb) r=m*R*m'+cumb; [u,dr] = eig(r([1 3],[1 3])); diff --git a/atmat/atphysics/Radiation/quantumDiff.m b/atmat/atphysics/Radiation/quantumDiff.m index 23aa0233e9..3b37f4ebc1 100644 --- a/atmat/atphysics/Radiation/quantumDiff.m +++ b/atmat/atphysics/Radiation/quantumDiff.m @@ -1,65 +1,32 @@ function DiffMat = quantumDiff(elems, varargin) -%quantumDiff Compute the radiation-diffusion matrix +%QUANTUMDIFF Compute the radiation-diffusion matrix % %DIFFMAT=QUANTUMDIFF(RING) % RING: Closed ring AT structure, containing radiative elements and % RF cavity. Radiative elements are identified by a % PassMethod ending with 'RadPass'. % -%DIFFMAT=QUANTUMDIFF(RING,RADINDEX) -% RADINDEX: Indices of elements where diffusion occurs, typically dipoles -% and quadrupoles. -% %DIFFMAT=QUANTUMDIFF(LINE,RADINDEX,ORBITIN) (Deprecated syntax) %DIFFMAT=QUANTUMDIFF(...,'orbit',ORBITIN) +% RADINDEX: Ignored % ORBITIN: Initial 6-D closed orbit. % In this mode, LINE may be a section of a ring. +% +%DIFFMAT=QUANTUMDIFF(...,'Energy',ENERGY) +% ENERGY: Lattice energy -NumElements=length(elems); - -[orb0,varargs]=getoption(varargin, 'orbit', []); -if length(varargs) >= 2 % QUANTUMDIFF(RING,RADINDEX,ORBITIN) - orb0 = varargs{2}; -end -if length(varargs) >= 1 % QUANTUMDIFF(RING,RADINDEX,...) - radindex=varargs{1}; -else % QUANTUMDIFF(RING) - radindex=atgetcells(elems,'PassMethod',@(elem, pass) endsWith(pass,'RadPass')); -end +[orb0,varargs]=getoption(varargin, 'orbit',[]); +[energy,varargs]=getoption(varargs,'Energy',[]); +[radindex,orb0,varargs]=getargs(varargs,[],orb0); %#ok -%[mring, ms, orbit] = findm66(ring,1:NumElements+1); if isempty(orb0) - orb=num2cell(findorbit6(elems, 1:NumElements),1)'; -else - orb=num2cell(linepass(elems, orb0, 1:NumElements),1)'; + orb0=findorbit6(elems); +end +if isempty(energy) + energy=atGetRingProperties(elems,'Energy'); end -zr={zeros(6,6)}; -B=zr(ones(NumElements,1)); % B{i} is the diffusion matrix of the i-th element - -% calculate Radiation-Diffusion matrix B for elements with radiation -B(radindex)=cellfun(@findmpoleraddiffmatrix,... - elems(radindex),orb(radindex),'UniformOutput',false); - -% Calculate cumulative Radiation-Diffusion matrix for the ring -BCUM = zeros(6,6); -% Batbeg{i} is the cumulative diffusion matrix from -% 0 to the beginning of the i-th element -Batbeg=[zr;cellfun(@cumulb,elems,orb,B,'UniformOutput',false)]; %#ok - -DiffCum = BCUM; - -DiffMat=(DiffCum+DiffCum')/2; - -%Lmat=chol((DiffCum+DiffCum')/2); +BCUM=atdiffmat(elems,energy,'orbit',orb0); - function btx=cumulb(elem,orbit,b) - % Calculate 6-by-6 linear transfer matrix in each element - % near the equilibrium orbit - m=findelemm66(elem,elem.PassMethod,orbit); - % Cumulative diffusion matrix of the entire ring - BCUM = m*BCUM*m' + b; - btx=BCUM; - end +DiffMat=(BCUM+BCUM')/2; end - diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m index ed8c9d547b..567aa3471c 100644 --- a/atmat/attests/pytests.m +++ b/atmat/attests/pytests.m @@ -184,7 +184,7 @@ function tunechrom6(testCase,lat2,dp) ptune=double(plat.get_tune()); pchrom=double(plat.get_chrom()); testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=1.e-9); - testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=2.e-4); + testCase.verifyEqual(mchrom*periodicity,pchrom,RelTol=1.e-4,AbsTol=3.e-4); end function linopt1(testCase,dp) @@ -209,6 +209,7 @@ function linopt1(testCase,dp) end function avlin1(testCase, lat, dp) + % Check average beta, dispersion, phase advance lattice=testCase.ring4.(lat); mrefs=true(1,length(lattice.m)); prefs=py.numpy.array(mrefs); @@ -272,5 +273,27 @@ function checkattr(rm, rp) end end end + + function emittances(testCase, lat2) + % Check emittances, tunes and damping rates + lattice=testCase.ring6.(lat2); + % python + pdata=cell(lattice.p.ohmi_envelope()); + pemit=double(py.getattr(pdata{2},'mode_emittances')); + ptunes=double(py.getattr(pdata{2},'tunes')); + pdamprate=double(py.getattr(pdata{2},'damping_rates')); + %matlab + [mdata,~,~,m]=ohmienvelope(lattice.m); + jmt=jmat(3); + aa=amat(m); + nn=-aa'*jmt*mdata.R*jmt*aa; + memit=0.5*[nn(1,1)+nn(2,2) nn(3,3)+nn(4,4) nn(5,5)+nn(6,6)]; + [mtunes,mdamprate]=atdampingrates(m); + %check + testCase.verifyEqual(memit,pemit,AbsTol=1.e-30,RelTol=1.e-6); + testCase.verifyEqual(mtunes,ptunes,AbsTol=1.e-10); + testCase.verifyEqual(mdamprate,pdamprate,AbsTol=1.e-10); + end + end end diff --git a/atmat/attrack/atpass.c b/atmat/attrack/atpass.c index c46b5362c5..6057991a52 100644 --- a/atmat/attrack/atpass.c +++ b/atmat/attrack/atpass.c @@ -11,6 +11,7 @@ #include "attypes.h" #include "elempass.h" #include "atrandom.c" +#include "ringproperties.c" /* Get ready for R2018a C matrix API */ #ifndef mxGetDoubles @@ -235,42 +236,7 @@ static mxDouble *passhook(mxArray *mxPassArg[], mxArray *mxElem, mxArray *func) return tempdoubleptr; } -static double getoptionaldoubleprop(const mxArray *obj, const char *fieldname, double default_value) -{ - mxArray *field=mxGetProperty(obj, 0, fieldname); - return (field) ? mxGetScalar(field) : default_value; -} - -static void getparticle(const mxArray *opts, double *rest_energy, double *charge) -{ - const mxArray *part = mxGetField(opts, 0, "Particle"); - if (part) { - if (mxIsClass(part, "atparticle")) { /* OK */ - *rest_energy = getoptionaldoubleprop(part, "rest_energy", 0.0); - *charge = getoptionaldoubleprop(part, "charge", -1.0); - } - else { /* particle is not a Particle object */ - mexErrMsgIdAndTxt("Atpass:WrongParameter","Particle must be an 'atparticle' object"); - } - } -} - -static void getproperties(const mxArray *opts, double *energy, double *rest_energy, double *charge) -{ - mxArray *field; - if (!mxIsStruct(opts)) { - mexErrMsgIdAndTxt("Atpass:WrongParameter","ring properties must be a struct"); - } - field = mxGetField(opts, 0, "Energy"); - if (field) { - *energy = mxGetScalar(field); - getparticle(opts, rest_energy, charge); - } -} - -/*! - * - +/* @param[in] [0] LATTICE @param[in,out] [1] INITCONDITIONS @param[in] [2] NEWLATTICE @@ -345,7 +311,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) param.nturn = counter; if (nrhs >= 10) { - getproperties(RINGPROPERTIES, ¶m.energy, ¶m.rest_energy, ¶m.charge); + atProperties(RINGPROPERTIES, ¶m.energy, ¶m.rest_energy, ¶m.charge); } if (nlhs >= 2) { diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m new file mode 100644 index 0000000000..8d567a71ca --- /dev/null +++ b/atmat/attrack/elempass.m @@ -0,0 +1,31 @@ +function rout = elempass(elem,rin,varargin) +%ELEMPASS Tracks particles through a single element +% +%ROUT=ELEMPASS(ELEM,RIN) Tracks particles through ELEM +% +% ELEM: lattice element +% RIN 6xN matrix: input coordinates of N particles +% +% ROUT 6xN matrix: output coordinates of N particles +% +% ROUT=ELEMPASS(...,'PassMethod',PASSMETHOD,...) +% Use PASSMETHOD (default: ELEM.PassMethod) +% +% ROUT=ELEMPASS(...,'Energy',ENERGY,...) +% Use ENERGY and ignore the 'Energy' field of elements +% +% ROUT=ELEMPASS(...,'Particle',PARTICLE,...) +% Use PARTICLE (default: relativistic) +% +% See also: RINGPASS, LINEPASS + +[props.Energy,varargs]=getoption(varargin,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',[]); +[methodname,varargs]=getoption(varargs,'PassMethod',elem.PassMethod); %#ok +if ~isempty(particle) + props.Particle=particle; +end + +rout = feval(methodname,elem,rin,props); + +end \ No newline at end of file diff --git a/machine_data/noenergy.mat b/machine_data/noenergy.mat new file mode 100644 index 0000000000..fa41926121 Binary files /dev/null and b/machine_data/noenergy.mat differ diff --git a/machine_data/noringparam.mat b/machine_data/noringparam.mat new file mode 100644 index 0000000000..40185d4700 Binary files /dev/null and b/machine_data/noringparam.mat differ diff --git a/pyat/at/load/matfile.py b/pyat/at/load/matfile.py index 66acea2e3f..721c3af5d3 100644 --- a/pyat/at/load/matfile.py +++ b/pyat/at/load/matfile.py @@ -11,6 +11,7 @@ from os.path import abspath, basename, splitext from typing import Any from collections.abc import Sequence, Generator +from math import isfinite from warnings import warn import numpy as np @@ -147,7 +148,8 @@ def ringparam_filter( for k, v in elem.items(): k2 = _m2p.get(k, k) if k2 is not None: - params.setdefault(k2, v) + if k2 != "cell_harmnumber" or isfinite(v): + params.setdefault(k2, v) if keep_all: pars = vars(elem).copy() name = pars.pop("FamName") diff --git a/pyat/at/load/utils.py b/pyat/at/load/utils.py index 8d4872ab97..7aa6eab428 100644 --- a/pyat/at/load/utils.py +++ b/pyat/at/load/utils.py @@ -86,7 +86,11 @@ class RingParam(elt.Element): "Periodicity", ] _conversions = dict( - elt.Element._conversions, Energy=float, Periodicity=int, Particle=_particle + elt.Element._conversions, + Energy=float, + Periodicity=int, + Particle=_particle, + cell_harmnumber=float, ) # noinspection PyPep8Naming diff --git a/pyat/machine_data/noenergy.mat b/pyat/machine_data/noenergy.mat new file mode 100644 index 0000000000..fa41926121 Binary files /dev/null and b/pyat/machine_data/noenergy.mat differ diff --git a/pyat/machine_data/noringparam.mat b/pyat/machine_data/noringparam.mat new file mode 100644 index 0000000000..40185d4700 Binary files /dev/null and b/pyat/machine_data/noringparam.mat differ diff --git a/pyat/test/conftest.py b/pyat/test/conftest.py index c8364d24bf..fa1afc6a3b 100644 --- a/pyat/test/conftest.py +++ b/pyat/test/conftest.py @@ -58,3 +58,17 @@ def hmba_lattice(): with as_file(files(machine_data) / 'hmba.mat') as path: ring = at.load_lattice(path) return ring + + +@pytest.fixture(scope='session') +def noenergy_lattice(): + with as_file(files(machine_data) / 'noenergy.mat') as path: + ring = at.load_lattice(path) + return ring + + +@pytest.fixture(scope='session') +def noringparam_lattice(): + with as_file(files(machine_data) / 'noringparam.mat') as path: + ring = at.load_lattice(path) + return ring diff --git a/pyat/test/test_linopt6.py b/pyat/test/test_linopt6.py index ba3c60862e..f22bf1fab0 100644 --- a/pyat/test/test_linopt6.py +++ b/pyat/test/test_linopt6.py @@ -4,7 +4,10 @@ from at import linopt2, linopt4, linopt6, get_optics -@pytest.mark.parametrize("lattice", ["dba_lattice", "hmba_lattice"]) +@pytest.mark.parametrize( + "lattice", + ["dba_lattice", "noenergy_lattice", "hmba_lattice", "noringparam_lattice"], +) def test_linopt6_norad(request, lattice): """Compare the results of linopt2 and linopt6 in 4d""" lattice = request.getfixturevalue(lattice) @@ -19,13 +22,16 @@ def test_linopt6_norad(request, lattice): assert_close(ld2.W, ld6.W, atol=1e-6, rtol=0) -@pytest.mark.parametrize("lattice", ["hmba_lattice"]) +@pytest.mark.parametrize( + "lattice", + ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) def test_linopt6_rad(request, lattice): """Compare the results with and without radiation""" lattice = request.getfixturevalue(lattice) refpts = range(len(lattice) + 1) # Turn cavity ON, without radiation - radlattice = lattice.radiation_on(dipole_pass=None, copy=True) + radlattice = lattice.enable_6d(dipole_pass=None, copy=True) ld04, rd4, ld4 = linopt6(lattice, refpts, get_w=True) ld06, rd6, ld6 = linopt6(radlattice, refpts, get_w=True) @@ -47,12 +53,12 @@ def test_linopt6_line(request, lattice, dp, method): refpts = lattice.uint32_refpts(range(len(lattice) + 1)) ld04, rd4, ld4 = get_optics(lattice, refpts, dp=dp, method=method) - twin = dict( - alpha=ld04.alpha, - beta=ld04.beta, - dispersion=ld04.dispersion, - closed_orbit=ld04.closed_orbit, - ) + twin = { + "alpha": ld04.alpha, + "beta": ld04.beta, + "dispersion": ld04.dispersion, + "closed_orbit": ld04.closed_orbit, + } # twin = ld04 tr04, bd4, tr4 = get_optics(lattice, refpts, dp=dp, twiss_in=twin, method=method) for field in ["s_pos", "closed_orbit", "dispersion", "alpha", "beta", "mu"]: diff --git a/pyat/test/test_physics.py b/pyat/test/test_physics.py index cd2519acfd..65af3e268c 100644 --- a/pyat/test/test_physics.py +++ b/pyat/test/test_physics.py @@ -1,12 +1,12 @@ -import at import numpy +import pytest from numpy.testing import assert_allclose as assert_close from numpy.testing import assert_equal -import pytest + +import at from at import AtWarning, physics -from at import lattice_track from at import lattice_pass, internal_lpass - +from at import lattice_track DP = 1e-5 DP2 = 0.005 @@ -101,10 +101,13 @@ def test_find_m44_returns_same_answer_as_matlab(dba_lattice, refpts): assert mstack.shape == (len(refpts), 4, 4) +@pytest.mark.parametrize( + "lattice", ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) @pytest.mark.parametrize('refpts', ([145], [20], [1, 2, 3])) -def test_find_m66(hmba_lattice, refpts): - hmba_lattice = hmba_lattice.radiation_on(copy=True) - m66, mstack = physics.find_m66(hmba_lattice, refpts=refpts) +def test_find_m66(request, lattice, refpts): + lattice = request.getfixturevalue(lattice).enable_6d(copy=True) + m66, mstack = lattice.find_m66(refpts=refpts) assert_close(m66, M66_MATLAB, rtol=0, atol=1e-8) stack_size = 0 if refpts is None else len(refpts) assert mstack.shape == (stack_size, 6, 6) @@ -139,10 +142,13 @@ def test_find_sync_orbit_finds_zeros(dba_lattice): numpy.testing.assert_equal(sync_orbit, numpy.zeros(6)) -def test_find_orbit6(hmba_lattice): - hmba_lattice = hmba_lattice.radiation_on(copy=True) - refpts = numpy.ones(len(hmba_lattice), dtype=bool) - orbit6, all_points = physics.find_orbit6(hmba_lattice, refpts) +@pytest.mark.parametrize( + "lattice", ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) +def test_find_orbit6(request, lattice): + lattice = request.getfixturevalue(lattice).enable_6d(copy=True) + refpts = numpy.ones(len(lattice), dtype=bool) + orbit6, all_points = lattice.find_orbit6(refpts) assert_close(orbit6, orbit6_MATLAB, rtol=0, atol=1e-12) @@ -350,10 +356,13 @@ def test_simple_ring(): assert_close(ring.get_tune(), [0.1, 0.2], atol=1e-10) +@pytest.mark.parametrize( + "lattice", ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) @pytest.mark.parametrize('refpts', ([121], [0, 40, 121])) -def test_ohmi_envelope(hmba_lattice, refpts): - hmba_lattice = hmba_lattice.radiation_on(copy=True) - emit0, beamdata, emit = hmba_lattice.ohmi_envelope(refpts) +def test_ohmi_envelope(request, lattice, refpts): + lattice = request.getfixturevalue(lattice).enable_6d(copy=True) + emit0, beamdata, emit = lattice.ohmi_envelope(refpts) obs = emit[-1] # All expected values are Matlab results