Skip to content

Forbid negative lattice or element energies #1066

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
swhite2401 merged 3 commits into from
Mar 28, 2026
Merged

Forbid negative lattice or element energies #1066

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
316e709
Forbid negative lattice or element energies
swhite2401 Mar 25, 2026
10e6cb0
bugfix: matlab error
swhite2401 Mar 25, 2026
ab88feb
python specific energy check
swhite2401 Mar 27, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 2 additions & 6 deletions atintegrators/BeamLoadingCavityPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,11 +234,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
feedback_angle_offset=atGetOptionalDouble(ElemData,"feedback_angle_offset", 0.0); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();
Copy link
Copy Markdown
Contributor

@lfarv lfarv Mar 26, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test as it is was implemented by @oscarxblanco in #1016 to solve the same issue as you do.
What you write here is exactly what I suggested in #1016 (comment) but was not implemented. So it's ok for me


int dimsth[] = {Param->nbunch*nslice*nturns, 4};
atCheckArrayDims(ElemData,"_turnhistory", 2, dimsth); check_error();
Expand Down Expand Up @@ -279,7 +275,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->ts = ts;
Elem->system_harmonic = system_harmonic;
}
energy = atEnergy(Param->energy, Elem->Energy);
energy = atEnergy(Param->energy, Elem->Energy); check_error();

if(num_particles<Param->nbunch){
atError("Number of particles has to be greater or equal to the number of bunches."); check_error();
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/BndMPoleSymplectic4E2RadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -233,11 +233,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand Down Expand Up @@ -265,7 +261,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->KickAngle=KickAngle;
}
irho = Elem->BendingAngle/Elem->Length;
energy = atEnergy(Param->energy, Elem->Energy);
energy = atEnergy(Param->energy, Elem->Energy); check_error();

BndMPoleSymplectic4E2RadPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/BndMPoleSymplectic4QuantPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,11 +216,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand Down Expand Up @@ -252,7 +248,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->KickAngle=KickAngle;
}
irho = Elem->BendingAngle/Elem->Length;
energy = atEnergy(Param->energy, Elem->Energy);
energy = atEnergy(Param->energy, Elem->Energy); check_error();

BndMPoleSymplectic4QuantPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/BndMPoleSymplectic4RadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -169,11 +169,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand Down Expand Up @@ -205,7 +201,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->KickAngle=KickAngle;
}
irho = Elem->BendingAngle/Elem->Length;
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

BndMPoleSymplectic4RadPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/CrabCavityPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,11 +102,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
SigVV = atGetOptionalDouble(ElemData,"SigVV",0.0); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand All @@ -125,7 +121,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
SigPhi=Elem->SigPhi;
SigVV=Elem->SigVV;
}
energy = atEnergy(Param->energy, Elem->Energy);
energy = atEnergy(Param->energy, Elem->Energy); check_error(); check_error();
nvx = Elem->Vx;
nvy = Elem->Vy;

Expand Down
2 changes: 1 addition & 1 deletion atintegrators/ExactMultipoleRadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem,
Elem->RApertures = RApertures;
Elem->KickAngle = KickAngle;
}
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

multipole_pass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/ExactRectBendRadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,11 +174,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
}

/* Check energy */
Energy = atEnergy(Param->energy, Energy);
if (Energy == 0) {
atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
check_error();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand Down Expand Up @@ -207,7 +203,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->RApertures=RApertures;
Elem->KickAngle=KickAngle;
}
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle,
Elem->PolynomA, Elem->PolynomB,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/ExactRectangularBendRadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,11 +176,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
}

/* Check energy */
Energy = atEnergy(Param->energy, Energy);
if (Energy == 0) {
atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
check_error();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand Down Expand Up @@ -209,7 +205,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->RApertures=RApertures;
Elem->KickAngle=KickAngle;
}
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/ExactSectorBendRadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,11 +161,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
}

/* Check energy */
Energy = atEnergy(Param->energy, Energy);
if (Energy == 0) {
atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
check_error();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand Down Expand Up @@ -193,7 +189,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->KickAngle=KickAngle;
}
irho = Elem->BendingAngle/Elem->Length;
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

ExactSectorBendRad(r_in, Elem->Length, irho,
Elem->PolynomA, Elem->PolynomB,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/GWigSymplecticPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,11 +145,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
T2 = atGetOptionalDoubleArray(ElemData, "T2"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Energy=Energy;
Expand All @@ -168,7 +164,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->T1=T1;
Elem->T2=T2;
}
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

GWigSymplecticPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax,
Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/GWigSymplecticRadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,11 +271,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
T2 = atGetOptionalDoubleArray(ElemData, "T2"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Energy=Energy;
Expand All @@ -294,7 +290,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->T1=T1;
Elem->T2=T2;
}
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

GWigSymplecticRadPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax,
Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/RFCavityPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand All @@ -64,7 +60,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->TimeLag=TimeLag;
Elem->PhaseLag=PhaseLag;
}
energy = atEnergy(energy, Elem->Energy);
energy = atEnergy(energy, Elem->Energy); check_error();

RFCavityPass(r_in, Elem->Length, Elem->Voltage/energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag,
Elem->PhaseLag, nturn, T0, num_particles);
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/StrMPoleSymplectic4QuantPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -185,11 +185,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand All @@ -212,7 +208,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->RApertures=RApertures;
Elem->KickAngle=KickAngle;
}
energy = atEnergy(Param->energy, Elem->Energy);
energy = atEnergy(Param->energy, Elem->Energy); check_error();

StrMPoleSymplectic4QuantPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps,
Expand Down
8 changes: 2 additions & 6 deletions atintegrators/StrMPoleSymplectic4RadPass.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -145,11 +145,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); 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();
}
Energy = atEnergy(Param->energy, Energy); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->Length=Length;
Expand All @@ -172,7 +168,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->RApertures=RApertures;
Elem->KickAngle=KickAngle;
}
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); check_error();

StrMPoleSymplectic4RadPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB,
Elem->MaxOrder, Elem->NumIntSteps,
Expand Down
36 changes: 22 additions & 14 deletions atintegrators/atelem.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,26 +76,17 @@ typedef mxArray atElem;

double atEnergy(double ringenergy, double elemenergy)
{
if (ringenergy!=0.0)
if (ringenergy > 0.0)
return ringenergy;
else
if (elemenergy!=0.0)
if (elemenergy > 0.0)
return elemenergy;
else {
atError("Energy not defined.");
atError("Energy must be positive. Check lattice or passmethod parameters.");
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)
{
mxArray *field;
Expand Down Expand Up @@ -206,11 +197,19 @@ typedef PyObject atElem;
#define atError(...) PyErr_Format(PyExc_ValueError, __VA_ARGS__)
#define atWarning(...) PyErr_WarnFormat(PyExc_RuntimeWarning, 0, __VA_ARGS__)
#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;

double atEnergy(double ringenergy, double elemenergy)
{
if (ringenergy > 0.0)
return ringenergy;
else {
atError("Energy must be positive. Check lattice or set the track() keyword argument.");
return 0.0; /* Never reached but makes the compiler happy */
}
}

static NUMPY_IMPORT_ARRAY_TYPE init_numpy(void)
{
import_array();
Expand Down Expand Up @@ -362,6 +361,15 @@ static double *atGetOptionalDoubleArray(const PyObject *element, char *name)
#define C_LINK
#endif

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;
}

Comment on lines +364 to +372
Copy link
Copy Markdown
Contributor

@lfarv lfarv Mar 26, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For python, the macro as it is is enough (and more efficient).

Copy link
Copy Markdown
Contributor Author

@swhite2401 swhite2401 Mar 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not so sure about this, how about a simulation of energy ramp where the energy is set to negative during the ramp? I think we still need the energy check in atGamma()

C_LINK ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, double *r_in,
int num_particles, struct parameters *Param);

Expand Down
6 changes: 5 additions & 1 deletion atintegrators/ringproperties.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,11 @@ static void atProperties(const mxArray *opts, double *energy, double *rest_energ
field = mxGetField(opts, 0, "Energy");
if (field) {
double ener = mxGetScalar(field);
if (ener != 0.0) *energy = ener;
if (ener > 0.0) {
*energy = ener;
}else{
mexErrMsgIdAndTxt("Atpass:WrongParameter", "Energy must be positive. Check lattice or passmethod parameters.");
}
atParticle(opts, rest_energy, charge);
}
}
Loading