From 252a317d21d67f832ef1f577168dea6faa66363e Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 10 Nov 2025 13:25:40 +0100 Subject: [PATCH 01/56] check energy is pass for wiggler methods --- atmat/attrack/elempass.m | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m index 8d567a71c..dcbf52e4c 100644 --- a/atmat/attrack/elempass.m +++ b/atmat/attrack/elempass.m @@ -19,13 +19,22 @@ % % See also: RINGPASS, LINEPASS -[props.Energy,varargs]=getoption(varargin,'Energy',0.0); +props = struct; +[varenergy,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 +% check pass methods that require energy +needs_energy = {'GWigSymplecticPass', 'GWigSymplecticRadPass'}; +if any(strcmp(needs_energy,methodname)) && varenergy == 0 + error('Energy needs to be defined'); +else + props.Energy = varenergy; +end + rout = feval(methodname,elem,rin,props); end \ No newline at end of file From c74340ccb2191fd52fb6b8543b80411e587cd6ef Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 10 Nov 2025 13:49:59 +0100 Subject: [PATCH 02/56] black, isor, flake8, numpy as np.; and energy check --- pyat/at/tracking/track.py | 80 ++++++++++++++++++++++++++------------- 1 file changed, 53 insertions(+), 27 deletions(-) diff --git a/pyat/at/tracking/track.py b/pyat/at/tracking/track.py index d39d2c2c9..292875b7f 100755 --- a/pyat/at/tracking/track.py +++ b/pyat/at/tracking/track.py @@ -13,20 +13,35 @@ import multiprocessing from collections.abc import Iterable +from enum import Enum from functools import partial from warnings import warn -from enum import Enum - -import numpy -from .atpass import atpass as _atpass, elempass as _elempass, reset_rng -from .utils import fortran_align, has_collective, format_results -from .utils import initialize_lpass, disable_varelem, variable_refs -from ..lattice import AtError, AtWarning, DConstant, random -from ..lattice import Lattice, Element, Refpts, End -from ..lattice import get_uint32_index -from ..cconfig import iscuda -from ..cconfig import isopencl +import numpy as np + +from ..cconfig import iscuda, isopencl +from ..lattice import ( + AtError, + AtWarning, + DConstant, + Element, + End, + Lattice, + Refpts, + get_uint32_index, + random, +) +from .atpass import atpass as _atpass +from .atpass import elempass as _elempass +from .atpass import reset_rng +from .utils import ( + disable_varelem, + format_results, + fortran_align, + has_collective, + initialize_lpass, + variable_refs, +) class MPMode(Enum): @@ -47,10 +62,10 @@ class MPMode(Enum): } if iscuda() or isopencl(): - from .gpu import gpupass as _gpupass from .gpu import gpuinfo as _gpuinfo + from .gpu import gpupass as _gpupass -_imax = numpy.iinfo(int).max +_imax = np.iinfo(int).max _globring: list[Element] | None = None @@ -71,7 +86,7 @@ def _atpass_spawn(ring, seed, rank, rin, **kwargs): def _pass(ring, r_in, pool_size, start_method, seed, **kwargs): ctx = multiprocessing.get_context(start_method) # Split input in as many slices as processes - args = enumerate(numpy.array_split(r_in, pool_size, axis=1)) + args = enumerate(np.array_split(r_in, pool_size, axis=1)) # Generate a new starting point for C RNGs global _globring _globring = ring @@ -325,7 +340,7 @@ def lattice_track( trackparam = {} part_kw = ["energy", "particle"] try: - npart = numpy.shape(r_in)[1] + npart = np.shape(r_in)[1] except IndexError: npart = 1 @@ -337,12 +352,12 @@ def lattice_track( lattice = initialize_lpass(lattice, nturns, kwargs) ldtype = [ - ("islost", numpy.bool_), - ("turn", numpy.uint32), - ("elem", numpy.uint32), - ("coord", numpy.float64, (6,)), + ("islost", np.bool_), + ("turn", np.uint32), + ("elem", np.uint32), + ("coord", np.float64, (6,)), ] - loss_map = numpy.recarray((npart,), ldtype) + loss_map = np.recarray((npart,), ldtype) lat_kw = ["turn"] [trackparam.update((kw, kwargs.get(kw))) for kw in kwargs if kw in lat_kw] trackparam.update({"refpts": get_uint32_index(lattice, refpts), "nturns": nturns}) @@ -375,15 +390,19 @@ def lattice_track( return rout, trackparam, trackdata -def element_track(element: Element, r_in, in_place: bool = False, **kwargs): - """ - :py:func:`element_track` tracks particles through one element of a - calling the element-specific tracking function specified in the - Element's *PassMethod* field +def element_track( + element: Element, r_in: np.ndarray, in_place: bool = False, **kwargs: dict[str, any] +) -> np.ndarray: + """:py:func:`element_track` tracks particles through one element. + + It calls the element-specific tracking function specified in the + Element's *PassMethod* field. Usage: >>> element_track(element, r_in) >>> element.track(r_in) + >>> element_track(element, r_in, energy=1e9) + Parameters: element: element to track through @@ -405,12 +424,19 @@ def element_track(element: Element, r_in, in_place: bool = False, **kwargs): Returns: r_out: (6, N, R, T) array containing output coordinates of N particles at R reference points for T turns + + Raises: + AtError: when energy is required by the pass method. """ if not in_place: r_in = r_in.copy() - rout = _element_pass(element, r_in, **kwargs) - return rout + # check if Energy is required + needs_energy = ("GWigSymplecticPass", "GWigSymplecticRadPass") + if (element.PassMethod in needs_energy) and ("energy" not in kwargs): + raise AtError("Energy is required.") + + return _element_pass(element, r_in, **kwargs) def gpu_core_count(gpuPool: list[int] | None): From 7218a927d9de8d208d682ea1b970628609a7564d Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 10 Nov 2025 14:01:32 +0100 Subject: [PATCH 03/56] energy=1e9; and black --- pyat/test/test_integrators.py | 64 +++++++++++++++++------------------ 1 file changed, 31 insertions(+), 33 deletions(-) diff --git a/pyat/test/test_integrators.py b/pyat/test/test_integrators.py index 2a896cd99..c754c2189 100644 --- a/pyat/test/test_integrators.py +++ b/pyat/test/test_integrators.py @@ -10,57 +10,54 @@ from at import element_pass, internal_epass -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_exact_hamiltonian_pass(rin, func): - drift = elements.Multipole('m1', 1, [0, 0, 0, 0], [0, 0, 0, 0]) + drift = elements.Multipole("m1", 1, [0, 0, 0, 0], [0, 0, 0, 0]) drift.Type = 0 - drift.PassMethod = 'ExactHamiltonianPass' + drift.PassMethod = "ExactHamiltonianPass" drift.BendingAngle = 0 func(drift, rin) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_exact_hamiltonian_pass_with_dls_dipole(rin, func): - bend = elements.Multipole('rb', 0.15, [0, 0, 0, 0], - [-0.0116333, 3.786786, 0, 0]) + bend = elements.Multipole("rb", 0.15, [0, 0, 0, 0], [-0.0116333, 3.786786, 0, 0]) bend.Type = 1 - bend.PassMethod = 'ExactHamiltonianPass' + bend.PassMethod = "ExactHamiltonianPass" bend.BendingAngle = -0.001745 bend.Energy = 3.5e9 bend.MaxOrder = 3 - if func==element_track: + if func == element_track: func(bend, rin, in_place=True) else: func(bend, rin) # Results from Matlab - expected = numpy.array([9.23965e-9, 1.22319e-5, 0, - 0, 0, -4.8100e-10]).reshape(6, 1) + expected = numpy.array([9.23965e-9, 1.22319e-5, 0, 0, 0, -4.8100e-10]).reshape(6, 1) numpy.testing.assert_allclose(rin, expected, rtol=1e-5, atol=1e-6) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) -@pytest.mark.parametrize('passmethod', - ('GWigSymplecticPass', 'GWigSymplecticRadPass')) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("passmethod", ("GWigSymplecticPass", "GWigSymplecticRadPass")) def test_gwig_symplectic_pass(rin, passmethod, func): # Parameters copied from one of the Diamond wigglers. - wiggler = elements.Wiggler('w', 1.15, 0.05, 0.8) + wiggler = elements.Wiggler("w", 1.15, 0.05, 0.8) wiggler.PassMethod = passmethod - func(wiggler, rin) + func(wiggler, rin, energy=1e9) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_bndstrmpole_symplectic_4_pass(rin, func): - bend = elements.Dipole('b', 1.0) - bend.PassMethod = 'BndStrMPoleSymplectic4Pass' + bend = elements.Dipole("b", 1.0) + bend.PassMethod = "BndStrMPoleSymplectic4Pass" func(bend, rin) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_pydrift(func): - pydrift = elements.Drift('drift', 1.0, PassMethod='pyDriftPass') - cdrift = elements.Drift('drift', 1.0, PassMethod='DriftPass') - pyout, *_ = func(pydrift, numpy.zeros(6)+1.0e-6) - cout, *_ = func(cdrift, numpy.zeros(6)+1.0e-6) + pydrift = elements.Drift("drift", 1.0, PassMethod="pyDriftPass") + cdrift = elements.Drift("drift", 1.0, PassMethod="DriftPass") + pyout, *_ = func(pydrift, numpy.zeros(6) + 1.0e-6) + cout, *_ = func(cdrift, numpy.zeros(6) + 1.0e-6) numpy.testing.assert_equal(pyout, cout) shift_elem(pydrift, 1.0e-3, 1.0e-3) @@ -71,23 +68,24 @@ def test_pydrift(func): tilt_elem(pydrift, 1.0e-3, 1.0e-3) tilt_elem(cdrift, 1.0e-3, 1.0e-3) - pyout, *_ = func(pydrift, numpy.zeros(6)+1.0e-6) - cout, *_ = func(cdrift, numpy.zeros(6)+1.0e-6) + pyout, *_ = func(pydrift, numpy.zeros(6) + 1.0e-6) + cout, *_ = func(cdrift, numpy.zeros(6) + 1.0e-6) numpy.testing.assert_equal(pyout, cout) # Multiple particles - pyout, *_ = func(pydrift, numpy.zeros((6, 2))+1.0e-6) - cout, *_ = func(cdrift, numpy.zeros((6, 2))+1.0e-6) + pyout, *_ = func(pydrift, numpy.zeros((6, 2)) + 1.0e-6) + cout, *_ = func(cdrift, numpy.zeros((6, 2)) + 1.0e-6) numpy.testing.assert_equal(pyout, cout) -@pytest.mark.parametrize('func', (lattice_track, lattice_pass, internal_lpass)) +@pytest.mark.parametrize("func", (lattice_track, lattice_pass, internal_lpass)) def test_pyintegrator(hmba_lattice, func): - params = {'Length': 0, - 'PassMethod': 'pyIdentityPass', - } - id_elem = Element('py_id', **params) - pin = numpy.zeros((6, 2))+1.0e-6 + params = { + "Length": 0, + "PassMethod": "pyIdentityPass", + } + id_elem = Element("py_id", **params) + pin = numpy.zeros((6, 2)) + 1.0e-6 pout1, *_ = func(hmba_lattice, pin.copy(), nturns=1) hmba_lattice = hmba_lattice + [id_elem] pout2, *_ = func(hmba_lattice, pin.copy(), nturns=1) From a943191e1f5204b9aab41aabc06cd4e262c0d334 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 1 Dec 2025 17:02:28 +0100 Subject: [PATCH 04/56] add RFCavityPass --- pyat/at/tracking/track.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyat/at/tracking/track.py b/pyat/at/tracking/track.py index 292875b7f..70e0ed29e 100755 --- a/pyat/at/tracking/track.py +++ b/pyat/at/tracking/track.py @@ -432,7 +432,7 @@ def element_track( r_in = r_in.copy() # check if Energy is required - needs_energy = ("GWigSymplecticPass", "GWigSymplecticRadPass") + needs_energy = ("GWigSymplecticPass", "GWigSymplecticRadPass", "RFCavityPass") if (element.PassMethod in needs_energy) and ("energy" not in kwargs): raise AtError("Energy is required.") From 09e05bc61db834a541045c74a65aa47ae48a95c0 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 1 Dec 2025 17:53:12 +0100 Subject: [PATCH 05/56] add RFCavityPass --- atmat/attrack/elempass.m | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m index dcbf52e4c..5c95d41a1 100644 --- a/atmat/attrack/elempass.m +++ b/atmat/attrack/elempass.m @@ -28,7 +28,8 @@ end % check pass methods that require energy -needs_energy = {'GWigSymplecticPass', 'GWigSymplecticRadPass'}; +needs_energy = {'GWigSymplecticPass', 'GWigSymplecticRadPass', ... + 'RFCavityPass'}; if any(strcmp(needs_energy,methodname)) && varenergy == 0 error('Energy needs to be defined'); else From 7724c5f47d555a3d63e02448c8aeed6f494b6219 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 1 Dec 2025 18:32:13 +0100 Subject: [PATCH 06/56] check if Energy exists --- pyat/at/tracking/track.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/pyat/at/tracking/track.py b/pyat/at/tracking/track.py index 70e0ed29e..ae2b45cbe 100755 --- a/pyat/at/tracking/track.py +++ b/pyat/at/tracking/track.py @@ -433,7 +433,11 @@ def element_track( # check if Energy is required needs_energy = ("GWigSymplecticPass", "GWigSymplecticRadPass", "RFCavityPass") - if (element.PassMethod in needs_energy) and ("energy" not in kwargs): + if ( + element.PassMethod in needs_energy + and not hasattr(element, "Energy") + and ("energy" not in kwargs) + ): raise AtError("Energy is required.") return _element_pass(element, r_in, **kwargs) From bb8a4229d8b33d811a70ee69135670506b395ff8 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 1 Dec 2025 18:32:52 +0100 Subject: [PATCH 07/56] check if Energy exists --- atmat/attrack/elempass.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m index 5c95d41a1..03991b64a 100644 --- a/atmat/attrack/elempass.m +++ b/atmat/attrack/elempass.m @@ -30,8 +30,10 @@ % check pass methods that require energy needs_energy = {'GWigSymplecticPass', 'GWigSymplecticRadPass', ... 'RFCavityPass'}; -if any(strcmp(needs_energy,methodname)) && varenergy == 0 - error('Energy needs to be defined'); +if any(strcmp(needs_energy,methodname)) + if ~isfield(elem,'Energy') && varenergy == 0 + error('Energy needs to be defined'); + end else props.Energy = varenergy; end From e550f11b750b91de01cc92179708de9b70a8aab5 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 3 Dec 2025 16:54:40 +0100 Subject: [PATCH 08/56] restore from master --- atmat/attrack/elempass.m | 14 +------ pyat/at/tracking/track.py | 84 +++++++++++++-------------------------- 2 files changed, 28 insertions(+), 70 deletions(-) diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m index 03991b64a..8d567a71c 100644 --- a/atmat/attrack/elempass.m +++ b/atmat/attrack/elempass.m @@ -19,25 +19,13 @@ % % See also: RINGPASS, LINEPASS -props = struct; -[varenergy,varargs]=getoption(varargin,'Energy',0.0); +[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 -% check pass methods that require energy -needs_energy = {'GWigSymplecticPass', 'GWigSymplecticRadPass', ... - 'RFCavityPass'}; -if any(strcmp(needs_energy,methodname)) - if ~isfield(elem,'Energy') && varenergy == 0 - error('Energy needs to be defined'); - end -else - props.Energy = varenergy; -end - rout = feval(methodname,elem,rin,props); end \ No newline at end of file diff --git a/pyat/at/tracking/track.py b/pyat/at/tracking/track.py index ae2b45cbe..d39d2c2c9 100755 --- a/pyat/at/tracking/track.py +++ b/pyat/at/tracking/track.py @@ -13,35 +13,20 @@ import multiprocessing from collections.abc import Iterable -from enum import Enum from functools import partial from warnings import warn +from enum import Enum + +import numpy -import numpy as np - -from ..cconfig import iscuda, isopencl -from ..lattice import ( - AtError, - AtWarning, - DConstant, - Element, - End, - Lattice, - Refpts, - get_uint32_index, - random, -) -from .atpass import atpass as _atpass -from .atpass import elempass as _elempass -from .atpass import reset_rng -from .utils import ( - disable_varelem, - format_results, - fortran_align, - has_collective, - initialize_lpass, - variable_refs, -) +from .atpass import atpass as _atpass, elempass as _elempass, reset_rng +from .utils import fortran_align, has_collective, format_results +from .utils import initialize_lpass, disable_varelem, variable_refs +from ..lattice import AtError, AtWarning, DConstant, random +from ..lattice import Lattice, Element, Refpts, End +from ..lattice import get_uint32_index +from ..cconfig import iscuda +from ..cconfig import isopencl class MPMode(Enum): @@ -62,10 +47,10 @@ class MPMode(Enum): } if iscuda() or isopencl(): - from .gpu import gpuinfo as _gpuinfo from .gpu import gpupass as _gpupass + from .gpu import gpuinfo as _gpuinfo -_imax = np.iinfo(int).max +_imax = numpy.iinfo(int).max _globring: list[Element] | None = None @@ -86,7 +71,7 @@ def _atpass_spawn(ring, seed, rank, rin, **kwargs): def _pass(ring, r_in, pool_size, start_method, seed, **kwargs): ctx = multiprocessing.get_context(start_method) # Split input in as many slices as processes - args = enumerate(np.array_split(r_in, pool_size, axis=1)) + args = enumerate(numpy.array_split(r_in, pool_size, axis=1)) # Generate a new starting point for C RNGs global _globring _globring = ring @@ -340,7 +325,7 @@ def lattice_track( trackparam = {} part_kw = ["energy", "particle"] try: - npart = np.shape(r_in)[1] + npart = numpy.shape(r_in)[1] except IndexError: npart = 1 @@ -352,12 +337,12 @@ def lattice_track( lattice = initialize_lpass(lattice, nturns, kwargs) ldtype = [ - ("islost", np.bool_), - ("turn", np.uint32), - ("elem", np.uint32), - ("coord", np.float64, (6,)), + ("islost", numpy.bool_), + ("turn", numpy.uint32), + ("elem", numpy.uint32), + ("coord", numpy.float64, (6,)), ] - loss_map = np.recarray((npart,), ldtype) + loss_map = numpy.recarray((npart,), ldtype) lat_kw = ["turn"] [trackparam.update((kw, kwargs.get(kw))) for kw in kwargs if kw in lat_kw] trackparam.update({"refpts": get_uint32_index(lattice, refpts), "nturns": nturns}) @@ -390,19 +375,15 @@ def lattice_track( return rout, trackparam, trackdata -def element_track( - element: Element, r_in: np.ndarray, in_place: bool = False, **kwargs: dict[str, any] -) -> np.ndarray: - """:py:func:`element_track` tracks particles through one element. - - It calls the element-specific tracking function specified in the - Element's *PassMethod* field. +def element_track(element: Element, r_in, in_place: bool = False, **kwargs): + """ + :py:func:`element_track` tracks particles through one element of a + calling the element-specific tracking function specified in the + Element's *PassMethod* field Usage: >>> element_track(element, r_in) >>> element.track(r_in) - >>> element_track(element, r_in, energy=1e9) - Parameters: element: element to track through @@ -424,23 +405,12 @@ def element_track( Returns: r_out: (6, N, R, T) array containing output coordinates of N particles at R reference points for T turns - - Raises: - AtError: when energy is required by the pass method. """ if not in_place: r_in = r_in.copy() - # check if Energy is required - needs_energy = ("GWigSymplecticPass", "GWigSymplecticRadPass", "RFCavityPass") - if ( - element.PassMethod in needs_energy - and not hasattr(element, "Energy") - and ("energy" not in kwargs) - ): - raise AtError("Energy is required.") - - return _element_pass(element, r_in, **kwargs) + rout = _element_pass(element, r_in, **kwargs) + return rout def gpu_core_count(gpuPool: list[int] | None): From 0bd4453564b2e241656c7999254fd7848ee7fc54 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 3 Dec 2025 18:06:43 +0100 Subject: [PATCH 09/56] check Param->energy --- atintegrators/GWigSymplecticPass.c | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 58bc8e7f8..128324f73 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -121,6 +121,11 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, { double gamma; + if (Param->energy == 0){ + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + if (!Elem) { double *R1, *R2, *T1, *T2; double *By, *Bx; From de1759065dc1b27115d6baf9ba8d7fa912ea2792 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 3 Dec 2025 18:44:32 +0100 Subject: [PATCH 10/56] check gamma --- atintegrators/GWigSymplecticPass.c | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 128324f73..7100e4d02 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -121,11 +121,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, { double gamma; - if (Param->energy == 0){ - atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); - check_error(); - } - if (!Elem) { double *R1, *R2, *T1, *T2; double *By, *Bx; @@ -149,6 +144,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); + /* Check energy */ + gamma = atGamma(Param->energy, Energy, Param->rest_energy); + if (gamma == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Energy=Energy; Elem->Length=Ltot; @@ -165,8 +167,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->R2=R2; Elem->T1=T1; Elem->T2=T2; - } - gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); + }else{ + gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); + }; GWigSymplecticPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm, From b8f0876c7b3a6ba653bd61578196007c17b9fd98 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 3 Dec 2025 19:46:28 +0100 Subject: [PATCH 11/56] check RFCavityPass --- atintegrators/RFCavityPass.c | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index 60bf4c72a..b737ed707 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -47,6 +47,14 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); + + /* Check energy */ + energy = atEnergy(energy, Energy); + if (energy == 0){ + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + }; + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->Voltage=Voltage; @@ -56,7 +64,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->TimeLag=TimeLag; Elem->PhaseLag=PhaseLag; } - if (energy == 0.0) energy = Elem->Energy; + energy = atEnergy(energy, Elem->Energy); RFCavityPass(r_in, Elem->Length, Elem->Voltage/energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, Elem->PhaseLag, nturn, T0, num_particles); @@ -93,6 +101,9 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + /* Check energy */ + Energy = atEnergy(Energy, 0); + RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); } From b1bcbcbe4f9ba5062048829d7ee2f5a6542e28f0 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 17 Dec 2025 18:20:35 +0100 Subject: [PATCH 12/56] add energy parameter --- pyat/test/test_basic_elements.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyat/test/test_basic_elements.py b/pyat/test/test_basic_elements.py index 5cab9dbbb..f62b61b55 100644 --- a/pyat/test/test_basic_elements.py +++ b/pyat/test/test_basic_elements.py @@ -507,9 +507,9 @@ def test_rfcavity(rin, func): rin[4, 0] = 1e-6 rin[5, 0] = 1e-6 if func == lattice_track: - func(lattice, rin, in_place=True) + func(lattice, rin, energy=6.e+9, in_place=True) else: - func(lattice, rin) + func(lattice, rin, energy=6.e+9) expected = numpy.array([0., 0., 0., 0., 9.990769e-7, 1.e-6]).reshape(6, 1) numpy.testing.assert_allclose(rin, expected, atol=1e-12) From 8819a2862a3c265feb12ad346c1cef963c64c483 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 17 Dec 2025 18:24:52 +0100 Subject: [PATCH 13/56] remove energy from lattice_track --- pyat/test/test_basic_elements.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyat/test/test_basic_elements.py b/pyat/test/test_basic_elements.py index f62b61b55..6188eefed 100644 --- a/pyat/test/test_basic_elements.py +++ b/pyat/test/test_basic_elements.py @@ -507,7 +507,7 @@ def test_rfcavity(rin, func): rin[4, 0] = 1e-6 rin[5, 0] = 1e-6 if func == lattice_track: - func(lattice, rin, energy=6.e+9, in_place=True) + func(lattice, rin, in_place=True) else: func(lattice, rin, energy=6.e+9) expected = numpy.array([0., 0., 0., 0., 9.990769e-7, 1.e-6]).reshape(6, 1) From 2494d4f8173f30304d7d88c3a242ce79d68e5502 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Wed, 17 Dec 2025 18:51:14 +0100 Subject: [PATCH 14/56] rename lattice to rflattice --- pyat/test/test_basic_elements.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/pyat/test/test_basic_elements.py b/pyat/test/test_basic_elements.py index 6188eefed..c5bad7975 100644 --- a/pyat/test/test_basic_elements.py +++ b/pyat/test/test_basic_elements.py @@ -503,13 +503,13 @@ def test_quad(rin, func): @pytest.mark.parametrize('func', (lattice_track, lattice_pass, internal_lpass)) def test_rfcavity(rin, func): rf = elements.RFCavity('rfcavity', 0.0, 187500, 3.5237e+8, 31, 6.e+9) - lattice = [rf, rf, rf, rf] + rflattice = [rf, rf, rf, rf] rin[4, 0] = 1e-6 rin[5, 0] = 1e-6 if func == lattice_track: - func(lattice, rin, in_place=True) + func(rflattice, rin, in_place=True) else: - func(lattice, rin, energy=6.e+9) + func(rflattice, rin) expected = numpy.array([0., 0., 0., 0., 9.990769e-7, 1.e-6]).reshape(6, 1) numpy.testing.assert_allclose(rin, expected, atol=1e-12) From 4e6639326576434152aa9bf774d660ee66c1391b Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 17:09:51 +0100 Subject: [PATCH 15/56] pyat redefine atEnergy, atGamma --- atintegrators/atelem.c | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 8b597af37..e8a60183a 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -206,8 +206,6 @@ 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; @@ -352,6 +350,29 @@ static double *atGetOptionalDoubleArray(const PyObject *element, char *name) return atGetOptionalDoubleArraySz(element, name, &msz, &nsz); } +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; +} + + #endif /* defined(PYAT) */ #if defined(PYAT) || defined(MATLAB_MEX_FILE) From 2e997fff7074567495e228d92b7583e3546d437b Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 17:31:32 +0100 Subject: [PATCH 16/56] check error --- atintegrators/RFCavityPass.c | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index b737ed707..8f0d0ccae 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -50,10 +50,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, /* Check energy */ energy = atEnergy(energy, Energy); - if (energy == 0){ - atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); - check_error(); - }; + check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; @@ -65,6 +62,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->PhaseLag=PhaseLag; } 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); From a926212fa262a698c0aaecc16462affed71421ac Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 17:31:53 +0100 Subject: [PATCH 17/56] check error --- atintegrators/GWigSymplecticPass.c | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 7100e4d02..eb6314528 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -146,10 +146,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, /* Check energy */ gamma = atGamma(Param->energy, Energy, Param->rest_energy); - if (gamma == 0) { - atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); - check_error(); - } + check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Energy=Energy; @@ -168,7 +165,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T1=T1; Elem->T2=T2; }else{ - 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, @@ -224,6 +222,7 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); GWigSymplecticPass(r_in, Gamma, Ltot, Lw, Bmax, Nstep, Nmeth, NHharm, From 4f765282e8d81e40852240419470053dfa6b274f Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 17:33:09 +0100 Subject: [PATCH 18/56] check error --- atintegrators/GWigSymplecticRadPass.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index 3b508915b..5afba37f5 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -288,6 +288,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T2=T2; } 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, @@ -341,6 +342,7 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); GWigSymplecticRadPass(r_in, Gamma, Ltot, Lw, Bmax, Nstep, Nmeth, From a98e1a0133da5246c19f840ed9378f01afa63cdf Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 17:46:20 +0100 Subject: [PATCH 19/56] add check_error --- atintegrators/BeamLoadingCavityPass.c | 1 + atintegrators/BndMPoleSymplectic4E2RadPass.c | 1 + atintegrators/BndMPoleSymplectic4QuantPass.c | 1 + atintegrators/BndMPoleSymplectic4RadPass.c | 2 + atintegrators/CrabCavityPass.c | 1 + atintegrators/ExactMultipoleRadPass.c | 2 + atintegrators/ExactRectBendRadPass.c | 2 + atintegrators/ExactRectangularBendRadPass.c | 2 + atintegrators/ExactSectorBendRadPass.c | 2 + atintegrators/RFCavityPass.c | 1 + atintegrators/StrMPoleSymplectic4QuantPass.c | 1 + atintegrators/StrMPoleSymplectic4RadPass.c | 2 + atintegrators/VariableThinMPolePass.asv | 293 +++++++++++++++++++ 13 files changed, 311 insertions(+) create mode 100644 atintegrators/VariableThinMPolePass.asv diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index 0cee79530..63afbc0e3 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -262,6 +262,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->fbmode = fbmode; } energy = atEnergy(Param->energy, Elem->Energy); + check_error(); if(num_particlesnbunch){ atError("Number of particles has to be greater or equal to the number of bunches."); check_error(); diff --git a/atintegrators/BndMPoleSymplectic4E2RadPass.c b/atintegrators/BndMPoleSymplectic4E2RadPass.c index 2c5ae5ec4..9ff5154ea 100644 --- a/atintegrators/BndMPoleSymplectic4E2RadPass.c +++ b/atintegrators/BndMPoleSymplectic4E2RadPass.c @@ -259,6 +259,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; 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, diff --git a/atintegrators/BndMPoleSymplectic4QuantPass.c b/atintegrators/BndMPoleSymplectic4QuantPass.c index b4bb1304f..0b6922102 100644 --- a/atintegrators/BndMPoleSymplectic4QuantPass.c +++ b/atintegrators/BndMPoleSymplectic4QuantPass.c @@ -246,6 +246,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; 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, diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index e4c35e230..94a71f3fd 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -199,6 +199,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; 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, @@ -267,6 +268,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); BndMPoleSymplectic4RadPass(r_in, Length, irho, PolynomA, PolynomB, diff --git a/atintegrators/CrabCavityPass.c b/atintegrators/CrabCavityPass.c index f61ee1554..6e8474742 100644 --- a/atintegrators/CrabCavityPass.c +++ b/atintegrators/CrabCavityPass.c @@ -119,6 +119,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, SigVV=Elem->SigVV; } energy = atEnergy(Param->energy, Elem->Energy); + check_error(); nvx = Elem->Vx; nvy = Elem->Vy; diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index 1140ca0df..a39127d20 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -158,6 +158,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, Elem->KickAngle = KickAngle; } 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, @@ -208,6 +209,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); multipole_pass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 86d080724..a520584da 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -201,6 +201,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); + check_error(); ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, @@ -264,6 +265,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); ExactRectangularBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index f1a957c28..8ddb7885e 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -203,6 +203,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } 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, @@ -265,6 +266,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); ExactRectangularBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 9c72907bc..f44e450c7 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -187,6 +187,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); + check_error(); ExactSectorBendRad(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, @@ -250,6 +251,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxDuplicateArray(prhs[1]); irho = BendingAngle/Length; Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); ExactSectorBendRad(r_in, Length, irho, PolynomA, PolynomB, diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index 8f0d0ccae..3e63ec3dd 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -101,6 +101,7 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* Check energy */ Energy = atEnergy(Energy, 0); + check_error(); RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); diff --git a/atintegrators/StrMPoleSymplectic4QuantPass.c b/atintegrators/StrMPoleSymplectic4QuantPass.c index e641f151f..32d5b37b4 100644 --- a/atintegrators/StrMPoleSymplectic4QuantPass.c +++ b/atintegrators/StrMPoleSymplectic4QuantPass.c @@ -206,6 +206,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } energy = atEnergy(Param->energy, Elem->Energy); + check_error(); StrMPoleSymplectic4QuantPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index b2daf7a78..02c6e69fd 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/atintegrators/StrMPoleSymplectic4RadPass.c @@ -166,6 +166,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); + check_error(); StrMPoleSymplectic4RadPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, @@ -220,6 +221,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); + check_error(); r_in = mxGetDoubles(plhs[0]); StrMPoleSymplectic4RadPass(r_in, Length, PolynomA, PolynomB, diff --git a/atintegrators/VariableThinMPolePass.asv b/atintegrators/VariableThinMPolePass.asv new file mode 100644 index 000000000..8bae47a01 --- /dev/null +++ b/atintegrators/VariableThinMPolePass.asv @@ -0,0 +1,293 @@ +/* VariableThinMPolePass + Accelerator Toolbox + S.White +*/ + +#include "atconstants.h" +#include "atelem.c" +#include "atlalib.c" +#include "atrandom.c" +#include "driftkick.c" + +struct elemab { + double* Amplitude; + double Frequency; + double Phase; + int NSamples; + double* Func; + double* Func1deriv; + double* Func2deriv; + double* Func3deriv; +}; + +struct elem { + double* PolynomA; + double* PolynomB; + struct elemab* ElemA; + struct elemab* ElemB; + int Seed; + int Mode; + int MaxOrder; + double* Ramps; + int Periodic; +}; + +double get_amp(double amp, double* ramps, double t) +{ + double ampt = amp; + if (ramps) { + if (t <= ramps[0]) { + ampt = 0.0; + } else if (t <= ramps[1]) { + ampt = amp * (t - ramps[0]) / (ramps[1] - ramps[0]); + } else if (t <= ramps[2]) { + ampt = amp; + } else if (t <= ramps[3]) { + ampt = amp - amp * (t - ramps[2]) / (ramps[3] - ramps[2]); + } else { + ampt = 0.0; + } + } + return ampt; +} + +double get_pol(struct elemab* elem, double* ramps, int mode, + double t, int turn, int seed, int order, int periodic) +{ + int idx; + double ampt, freq, ph, val; + double* func, func1deriv, func2deriv, func3deriv; + double* amp = elem->Amplitude; + if (!amp) { + return 0.0; + } + tmod = t; + ampt = get_amp(amp[order], ramps, turn); + switch (mode) { + case 0: + freq = elem->Frequency; + ph = elem->Phase; + ampt *= sin(TWOPI * freq * t + ph); + return ampt; + case 1: + val = atrandn(0.0, 1.0); + ampt *= val; + return ampt; + case 2: + if (periodic || turn < elem->NSamples) { + func = elem->Func; + func1deriv = elem->Func1deriv; + func2deriv = elem->Func2deriv; + func3deriv = elem->Func3deriv; + idx = turn % elem->NSamples; + t2 = t*t; + oneoversix = 0.166666666666667; + ampt *= func[idx] + func1deriv[idx]*t + 0.5*func2deriv[idx]*t2 + + oneoversix*func3deriv[idx]*t2*t; + return ampt; + } else { + return 0.0; + } + default: + return 0.0; + } +} + +void VariableThinMPolePass(double* r, struct elem* Elem, double t0, int turn, int num_particles) +{ + + int i, c; + double* r6; + double t = t0 * turn; + + int maxorder = Elem->MaxOrder; + int periodic = Elem->Periodic; + double* pola = Elem->PolynomA; + double* polb = Elem->PolynomB; + int seed = Elem->Seed; + int mode = Elem->Mode; + struct elemab* ElemA = Elem->ElemA; + struct elemab* ElemB = Elem->ElemB; + double* ramps = Elem->Ramps; + + /*if (mode != 0) { + for (i = 0; i < maxorder + 1; i++) { + pola[i] = get_pol(ElemA, ramps, mode, t, turn, seed, i, periodic); + polb[i] = get_pol(ElemB, ramps, mode, t, turn, seed, i, periodic); + }; + };*/ + + for (c = 0; c < num_particles; c++) { + r6 = r + c * 6; + if (!atIsNaN(r6[0])) { +// if (mode == 0) { + double tpart = t + r6[5] / C0; + for (i = 0; i < maxorder + 1; i++) { + pola[i] = get_pol(ElemA, ramps, mode, tpart, turn, seed, i, periodic); + polb[i] = get_pol(ElemB, ramps, mode, tpart, turn, seed, i, periodic); + }; +// }; + strthinkick(r6, pola, polb, 1.0, maxorder); + } + } +} + +#if defined(MATLAB_MEX_FILE) || defined(PYAT) +ExportMode struct elem* trackFunction(const atElem* ElemData, struct elem* Elem, + double* r_in, int num_particles, struct parameters* Param) +{ + if (!Elem) { + int MaxOrder, Mode, Seed, NSamplesA, NSamplesB, Periodic; + double *PolynomA, *PolynomB, *AmplitudeA, *AmplitudeB; + double *Ramps, *FuncA, *FuncB; + double *FuncA1deriv, *FuncB1deriv; + double *FuncA2deriv, *FuncB2deriv; + double *FuncA3deriv, *FuncB3deriv; + double FrequencyA, FrequencyB; + double PhaseA, PhaseB; + struct elemab *ElemA, *ElemB; + MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); + Mode=atGetLong(ElemData,"Mode"); check_error(); + PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); + PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); + AmplitudeA=atGetOptionalDoubleArray(ElemData,"AmplitudeA"); check_error(); + AmplitudeB=atGetOptionalDoubleArray(ElemData,"AmplitudeB"); check_error(); + FrequencyA=atGetOptionalDouble(ElemData,"FrequencyA", 0); check_error(); + FrequencyB=atGetOptionalDouble(ElemData,"FrequencyB", 0); check_error(); + PhaseA=atGetOptionalDouble(ElemData,"PhaseA", 0); check_error(); + PhaseB=atGetOptionalDouble(ElemData,"PhaseB", 0); check_error(); + Ramps=atGetOptionalDoubleArray(ElemData, "Ramps"); check_error(); + Seed=atGetOptionalLong(ElemData, "Seed", 0); check_error(); + NSamplesA=atGetOptionalLong(ElemData, "NSamplesA", 1); check_error(); + NSamplesB=atGetOptionalLong(ElemData, "NSamplesB", 1); check_error(); + FuncA=atGetOptionalDoubleArray(ElemData,"FuncA"); check_error(); + FuncB=atGetOptionalDoubleArray(ElemData,"FuncB"); check_error(); + FuncA1deriv=atGetOptionalDoubleArray(ElemData,"FuncA1deriv"); check_error(); + FuncB1deriv=atGetOptionalDoubleArray(ElemData,"FuncB1deriv"); check_error(); + FuncA2deriv=atGetOptionalDoubleArray(ElemData,"FuncA2deriv"); check_error(); + FuncB2deriv=atGetOptionalDoubleArray(ElemData,"FuncB2deriv"); check_error(); + FuncA3deriv=atGetOptionalDoubleArray(ElemData,"FuncA3deriv"); check_error(); + FuncB3deriv=atGetOptionalDoubleArray(ElemData,"FuncB3deriv"); check_error(); + Periodic=atGetOptionalLong(ElemData,"Periodic", 1); check_error(); + Elem = (struct elem*)atMalloc(sizeof(struct elem)); + ElemA = (struct elemab*)atMalloc(sizeof(struct elemab)); + ElemB = (struct elemab*)atMalloc(sizeof(struct elemab)); + Elem->PolynomA = PolynomA; + Elem->PolynomB = PolynomB; + Elem->Ramps = Ramps; + Elem->Seed = Seed; + Elem->Mode = Mode; + Elem->MaxOrder = MaxOrder; + Elem->Periodic = Periodic; + ElemA->Amplitude = AmplitudeA; + ElemB->Amplitude = AmplitudeB; + ElemA->Frequency = FrequencyA; + ElemB->Frequency = FrequencyB; + ElemA->Phase = PhaseA; + ElemB->Phase = PhaseB; + ElemA->NSamples = NSamplesA; + ElemB->NSamples = NSamplesB; + ElemA->Func = FuncA; + ElemB->Func = FuncB; + ElemA->Func1deriv = FuncA1deriv; + ElemB->Func1deriv = FuncB1deriv; + ElemA->Func1deriv = FuncA1deriv; + ElemB->Func1deriv = FuncB1deriv; + + Elem->ElemA = ElemA; + Elem->ElemB = ElemB; + } + double t0 = Param->T0; + int turn = Param->nturn; + VariableThinMPolePass(r_in, Elem, t0, turn, num_particles); + return Elem; +} + +MODULE_DEF(VariableThinMPolePass) /* Dummy module initialisation */ + +#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ + +#ifdef MATLAB_MEX_FILE +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (nrhs == 2) { + double* r_in; + const mxArray* ElemData = prhs[0]; + int num_particles = mxGetN(prhs[1]); + int MaxOrder, Mode, Seed, NSamplesA, NSamplesB, Periodic; + double *PolynomA, *PolynomB, *AmplitudeA, *AmplitudeB; + double *Ramps, *FuncA, *FuncB; + double FrequencyA, FrequencyB; + double PhaseA, PhaseB; + struct elemab ElA, *ElemA = &ElA; + struct elemab ElB, *ElemB = &ElB; + struct elem El, *Elem = &El; + MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); + Mode=atGetLong(ElemData,"Mode"); check_error(); + PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); + PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); + AmplitudeA=atGetOptionalDoubleArray(ElemData,"AmplitudeA"); check_error(); + AmplitudeB=atGetOptionalDoubleArray(ElemData,"AmplitudeB"); check_error(); + FrequencyA=atGetOptionalDouble(ElemData,"FrequencyA", 0); check_error(); + FrequencyB=atGetOptionalDouble(ElemData,"FrequencyB", 0); check_error(); + PhaseA=atGetOptionalDouble(ElemData,"PhaseA", 0); check_error(); + PhaseB=atGetOptionalDouble(ElemData,"PhaseB", 0); check_error(); + Ramps=atGetOptionalDoubleArray(ElemData, "Ramps"); check_error(); + Seed=atGetOptionalLong(ElemData, "Seed", 0); check_error(); + NSamplesA=atGetOptionalLong(ElemData, "NSamplesA", 0); check_error(); + NSamplesB=atGetOptionalLong(ElemData, "NSamplesB", 0); check_error(); + FuncA=atGetOptionalDoubleArray(ElemData,"FuncA"); check_error(); + FuncB=atGetOptionalDoubleArray(ElemData,"FuncB"); check_error(); + Periodic=atGetOptionalLong(ElemData,"Periodic", 1); check_error(); + Elem->PolynomA = PolynomA; + Elem->PolynomB = PolynomB; + Elem->Ramps = Ramps; + Elem->Seed = Seed; + Elem->Mode = Mode; + Elem->MaxOrder = MaxOrder; + Elem->Periodic = Periodic; + ElemA->Amplitude = AmplitudeA; + ElemB->Amplitude = AmplitudeB; + ElemA->Frequency = FrequencyA; + ElemB->Frequency = FrequencyB; + ElemA->Phase = PhaseA; + ElemB->Phase = PhaseB; + ElemA->NSamples = NSamplesA; + ElemB->NSamples = NSamplesB; + ElemA->Func = FuncA; + ElemB->Func = FuncB; + Elem->ElemA = ElemA; + Elem->ElemB = ElemB; + /* ALLOCATE memory for the output array of the same size as the input */ + plhs[0] = mxDuplicateArray(prhs[1]); + r_in = mxGetDoubles(plhs[0]); + VariableThinMPolePass(r_in, Elem, 0, 0, num_particles); + } else if (nrhs == 0) { + /* list of required fields */ + plhs[0] = mxCreateCellMatrix(4, 1); + mxSetCell(plhs[0], 0, mxCreateString("MaxOrder")); + mxSetCell(plhs[0], 1, mxCreateString("Mode")); + mxSetCell(plhs[0], 2, mxCreateString("PolynomA")); + mxSetCell(plhs[0], 3, mxCreateString("PolynomB")); + if (nlhs > 1) { + /* list of optional fields */ + plhs[1] = mxCreateCellMatrix(13, 1); + mxSetCell(plhs[1], 0, mxCreateString("AmplitudeA")); + mxSetCell(plhs[1], 1, mxCreateString("AmplitudeB")); + mxSetCell(plhs[1], 2, mxCreateString("FrequencyA")); + mxSetCell(plhs[1], 3, mxCreateString("FrequencyB")); + mxSetCell(plhs[1], 4, mxCreateString("PhaseA")); + mxSetCell(plhs[1], 5, mxCreateString("PhaseB")); + mxSetCell(plhs[1], 6, mxCreateString("Ramps")); + mxSetCell(plhs[1], 7, mxCreateString("Seed")); + mxSetCell(plhs[1], 8, mxCreateString("FuncA")); + mxSetCell(plhs[1], 9, mxCreateString("FuncB")); + mxSetCell(plhs[1], 10, mxCreateString("NSamplesA")); + mxSetCell(plhs[1], 11, mxCreateString("NSamplesB")); + mxSetCell(plhs[1], 12, mxCreateString("Periodic")); + } + } else { + mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments"); + } +} +#endif /* MATLAB_MEX_FILE */ From 6da92a7a2c2e68ccfb2332e4214dcdedb62d208a Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 17:51:05 +0100 Subject: [PATCH 20/56] rm asv --- atintegrators/VariableThinMPolePass.asv | 293 ------------------------ 1 file changed, 293 deletions(-) delete mode 100644 atintegrators/VariableThinMPolePass.asv diff --git a/atintegrators/VariableThinMPolePass.asv b/atintegrators/VariableThinMPolePass.asv deleted file mode 100644 index 8bae47a01..000000000 --- a/atintegrators/VariableThinMPolePass.asv +++ /dev/null @@ -1,293 +0,0 @@ -/* VariableThinMPolePass - Accelerator Toolbox - S.White -*/ - -#include "atconstants.h" -#include "atelem.c" -#include "atlalib.c" -#include "atrandom.c" -#include "driftkick.c" - -struct elemab { - double* Amplitude; - double Frequency; - double Phase; - int NSamples; - double* Func; - double* Func1deriv; - double* Func2deriv; - double* Func3deriv; -}; - -struct elem { - double* PolynomA; - double* PolynomB; - struct elemab* ElemA; - struct elemab* ElemB; - int Seed; - int Mode; - int MaxOrder; - double* Ramps; - int Periodic; -}; - -double get_amp(double amp, double* ramps, double t) -{ - double ampt = amp; - if (ramps) { - if (t <= ramps[0]) { - ampt = 0.0; - } else if (t <= ramps[1]) { - ampt = amp * (t - ramps[0]) / (ramps[1] - ramps[0]); - } else if (t <= ramps[2]) { - ampt = amp; - } else if (t <= ramps[3]) { - ampt = amp - amp * (t - ramps[2]) / (ramps[3] - ramps[2]); - } else { - ampt = 0.0; - } - } - return ampt; -} - -double get_pol(struct elemab* elem, double* ramps, int mode, - double t, int turn, int seed, int order, int periodic) -{ - int idx; - double ampt, freq, ph, val; - double* func, func1deriv, func2deriv, func3deriv; - double* amp = elem->Amplitude; - if (!amp) { - return 0.0; - } - tmod = t; - ampt = get_amp(amp[order], ramps, turn); - switch (mode) { - case 0: - freq = elem->Frequency; - ph = elem->Phase; - ampt *= sin(TWOPI * freq * t + ph); - return ampt; - case 1: - val = atrandn(0.0, 1.0); - ampt *= val; - return ampt; - case 2: - if (periodic || turn < elem->NSamples) { - func = elem->Func; - func1deriv = elem->Func1deriv; - func2deriv = elem->Func2deriv; - func3deriv = elem->Func3deriv; - idx = turn % elem->NSamples; - t2 = t*t; - oneoversix = 0.166666666666667; - ampt *= func[idx] + func1deriv[idx]*t + 0.5*func2deriv[idx]*t2 - + oneoversix*func3deriv[idx]*t2*t; - return ampt; - } else { - return 0.0; - } - default: - return 0.0; - } -} - -void VariableThinMPolePass(double* r, struct elem* Elem, double t0, int turn, int num_particles) -{ - - int i, c; - double* r6; - double t = t0 * turn; - - int maxorder = Elem->MaxOrder; - int periodic = Elem->Periodic; - double* pola = Elem->PolynomA; - double* polb = Elem->PolynomB; - int seed = Elem->Seed; - int mode = Elem->Mode; - struct elemab* ElemA = Elem->ElemA; - struct elemab* ElemB = Elem->ElemB; - double* ramps = Elem->Ramps; - - /*if (mode != 0) { - for (i = 0; i < maxorder + 1; i++) { - pola[i] = get_pol(ElemA, ramps, mode, t, turn, seed, i, periodic); - polb[i] = get_pol(ElemB, ramps, mode, t, turn, seed, i, periodic); - }; - };*/ - - for (c = 0; c < num_particles; c++) { - r6 = r + c * 6; - if (!atIsNaN(r6[0])) { -// if (mode == 0) { - double tpart = t + r6[5] / C0; - for (i = 0; i < maxorder + 1; i++) { - pola[i] = get_pol(ElemA, ramps, mode, tpart, turn, seed, i, periodic); - polb[i] = get_pol(ElemB, ramps, mode, tpart, turn, seed, i, periodic); - }; -// }; - strthinkick(r6, pola, polb, 1.0, maxorder); - } - } -} - -#if defined(MATLAB_MEX_FILE) || defined(PYAT) -ExportMode struct elem* trackFunction(const atElem* ElemData, struct elem* Elem, - double* r_in, int num_particles, struct parameters* Param) -{ - if (!Elem) { - int MaxOrder, Mode, Seed, NSamplesA, NSamplesB, Periodic; - double *PolynomA, *PolynomB, *AmplitudeA, *AmplitudeB; - double *Ramps, *FuncA, *FuncB; - double *FuncA1deriv, *FuncB1deriv; - double *FuncA2deriv, *FuncB2deriv; - double *FuncA3deriv, *FuncB3deriv; - double FrequencyA, FrequencyB; - double PhaseA, PhaseB; - struct elemab *ElemA, *ElemB; - MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); - Mode=atGetLong(ElemData,"Mode"); check_error(); - PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); - PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); - AmplitudeA=atGetOptionalDoubleArray(ElemData,"AmplitudeA"); check_error(); - AmplitudeB=atGetOptionalDoubleArray(ElemData,"AmplitudeB"); check_error(); - FrequencyA=atGetOptionalDouble(ElemData,"FrequencyA", 0); check_error(); - FrequencyB=atGetOptionalDouble(ElemData,"FrequencyB", 0); check_error(); - PhaseA=atGetOptionalDouble(ElemData,"PhaseA", 0); check_error(); - PhaseB=atGetOptionalDouble(ElemData,"PhaseB", 0); check_error(); - Ramps=atGetOptionalDoubleArray(ElemData, "Ramps"); check_error(); - Seed=atGetOptionalLong(ElemData, "Seed", 0); check_error(); - NSamplesA=atGetOptionalLong(ElemData, "NSamplesA", 1); check_error(); - NSamplesB=atGetOptionalLong(ElemData, "NSamplesB", 1); check_error(); - FuncA=atGetOptionalDoubleArray(ElemData,"FuncA"); check_error(); - FuncB=atGetOptionalDoubleArray(ElemData,"FuncB"); check_error(); - FuncA1deriv=atGetOptionalDoubleArray(ElemData,"FuncA1deriv"); check_error(); - FuncB1deriv=atGetOptionalDoubleArray(ElemData,"FuncB1deriv"); check_error(); - FuncA2deriv=atGetOptionalDoubleArray(ElemData,"FuncA2deriv"); check_error(); - FuncB2deriv=atGetOptionalDoubleArray(ElemData,"FuncB2deriv"); check_error(); - FuncA3deriv=atGetOptionalDoubleArray(ElemData,"FuncA3deriv"); check_error(); - FuncB3deriv=atGetOptionalDoubleArray(ElemData,"FuncB3deriv"); check_error(); - Periodic=atGetOptionalLong(ElemData,"Periodic", 1); check_error(); - Elem = (struct elem*)atMalloc(sizeof(struct elem)); - ElemA = (struct elemab*)atMalloc(sizeof(struct elemab)); - ElemB = (struct elemab*)atMalloc(sizeof(struct elemab)); - Elem->PolynomA = PolynomA; - Elem->PolynomB = PolynomB; - Elem->Ramps = Ramps; - Elem->Seed = Seed; - Elem->Mode = Mode; - Elem->MaxOrder = MaxOrder; - Elem->Periodic = Periodic; - ElemA->Amplitude = AmplitudeA; - ElemB->Amplitude = AmplitudeB; - ElemA->Frequency = FrequencyA; - ElemB->Frequency = FrequencyB; - ElemA->Phase = PhaseA; - ElemB->Phase = PhaseB; - ElemA->NSamples = NSamplesA; - ElemB->NSamples = NSamplesB; - ElemA->Func = FuncA; - ElemB->Func = FuncB; - ElemA->Func1deriv = FuncA1deriv; - ElemB->Func1deriv = FuncB1deriv; - ElemA->Func1deriv = FuncA1deriv; - ElemB->Func1deriv = FuncB1deriv; - - Elem->ElemA = ElemA; - Elem->ElemB = ElemB; - } - double t0 = Param->T0; - int turn = Param->nturn; - VariableThinMPolePass(r_in, Elem, t0, turn, num_particles); - return Elem; -} - -MODULE_DEF(VariableThinMPolePass) /* Dummy module initialisation */ - -#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/ - -#ifdef MATLAB_MEX_FILE -void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) -{ - if (nrhs == 2) { - double* r_in; - const mxArray* ElemData = prhs[0]; - int num_particles = mxGetN(prhs[1]); - int MaxOrder, Mode, Seed, NSamplesA, NSamplesB, Periodic; - double *PolynomA, *PolynomB, *AmplitudeA, *AmplitudeB; - double *Ramps, *FuncA, *FuncB; - double FrequencyA, FrequencyB; - double PhaseA, PhaseB; - struct elemab ElA, *ElemA = &ElA; - struct elemab ElB, *ElemB = &ElB; - struct elem El, *Elem = &El; - MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); - Mode=atGetLong(ElemData,"Mode"); check_error(); - PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); - PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); - AmplitudeA=atGetOptionalDoubleArray(ElemData,"AmplitudeA"); check_error(); - AmplitudeB=atGetOptionalDoubleArray(ElemData,"AmplitudeB"); check_error(); - FrequencyA=atGetOptionalDouble(ElemData,"FrequencyA", 0); check_error(); - FrequencyB=atGetOptionalDouble(ElemData,"FrequencyB", 0); check_error(); - PhaseA=atGetOptionalDouble(ElemData,"PhaseA", 0); check_error(); - PhaseB=atGetOptionalDouble(ElemData,"PhaseB", 0); check_error(); - Ramps=atGetOptionalDoubleArray(ElemData, "Ramps"); check_error(); - Seed=atGetOptionalLong(ElemData, "Seed", 0); check_error(); - NSamplesA=atGetOptionalLong(ElemData, "NSamplesA", 0); check_error(); - NSamplesB=atGetOptionalLong(ElemData, "NSamplesB", 0); check_error(); - FuncA=atGetOptionalDoubleArray(ElemData,"FuncA"); check_error(); - FuncB=atGetOptionalDoubleArray(ElemData,"FuncB"); check_error(); - Periodic=atGetOptionalLong(ElemData,"Periodic", 1); check_error(); - Elem->PolynomA = PolynomA; - Elem->PolynomB = PolynomB; - Elem->Ramps = Ramps; - Elem->Seed = Seed; - Elem->Mode = Mode; - Elem->MaxOrder = MaxOrder; - Elem->Periodic = Periodic; - ElemA->Amplitude = AmplitudeA; - ElemB->Amplitude = AmplitudeB; - ElemA->Frequency = FrequencyA; - ElemB->Frequency = FrequencyB; - ElemA->Phase = PhaseA; - ElemB->Phase = PhaseB; - ElemA->NSamples = NSamplesA; - ElemB->NSamples = NSamplesB; - ElemA->Func = FuncA; - ElemB->Func = FuncB; - Elem->ElemA = ElemA; - Elem->ElemB = ElemB; - /* ALLOCATE memory for the output array of the same size as the input */ - plhs[0] = mxDuplicateArray(prhs[1]); - r_in = mxGetDoubles(plhs[0]); - VariableThinMPolePass(r_in, Elem, 0, 0, num_particles); - } else if (nrhs == 0) { - /* list of required fields */ - plhs[0] = mxCreateCellMatrix(4, 1); - mxSetCell(plhs[0], 0, mxCreateString("MaxOrder")); - mxSetCell(plhs[0], 1, mxCreateString("Mode")); - mxSetCell(plhs[0], 2, mxCreateString("PolynomA")); - mxSetCell(plhs[0], 3, mxCreateString("PolynomB")); - if (nlhs > 1) { - /* list of optional fields */ - plhs[1] = mxCreateCellMatrix(13, 1); - mxSetCell(plhs[1], 0, mxCreateString("AmplitudeA")); - mxSetCell(plhs[1], 1, mxCreateString("AmplitudeB")); - mxSetCell(plhs[1], 2, mxCreateString("FrequencyA")); - mxSetCell(plhs[1], 3, mxCreateString("FrequencyB")); - mxSetCell(plhs[1], 4, mxCreateString("PhaseA")); - mxSetCell(plhs[1], 5, mxCreateString("PhaseB")); - mxSetCell(plhs[1], 6, mxCreateString("Ramps")); - mxSetCell(plhs[1], 7, mxCreateString("Seed")); - mxSetCell(plhs[1], 8, mxCreateString("FuncA")); - mxSetCell(plhs[1], 9, mxCreateString("FuncB")); - mxSetCell(plhs[1], 10, mxCreateString("NSamplesA")); - mxSetCell(plhs[1], 11, mxCreateString("NSamplesB")); - mxSetCell(plhs[1], 12, mxCreateString("Periodic")); - } - } else { - mexErrMsgIdAndTxt("AT:WrongArg", "Needs 0 or 2 arguments"); - } -} -#endif /* MATLAB_MEX_FILE */ From 052f68c28a7233b7d35049d78988358174003d89 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 18:15:57 +0100 Subject: [PATCH 21/56] black isort flake8; numpy as np --- pyat/test/test_basic_elements.py | 250 ++++++++++++++++--------------- 1 file changed, 132 insertions(+), 118 deletions(-) diff --git a/pyat/test/test_basic_elements.py b/pyat/test/test_basic_elements.py index c5bad7975..be8dae96f 100644 --- a/pyat/test/test_basic_elements.py +++ b/pyat/test/test_basic_elements.py @@ -1,43 +1,51 @@ -import pytest -import numpy import warnings -from at import element_track, lattice_track -from at import lattice_pass, internal_lpass -from at import element_pass, internal_epass + +import numpy as np +import pytest +from at import ( + element_pass, + element_track, + elements, + internal_epass, + internal_lpass, + lattice, + lattice_pass, + lattice_track, +) from at.lattice.elements.conversions import _array, _array66 -from at import elements, lattice -from numpy.testing import assert_equal def test_data_checks(): - val = numpy.zeros([6, 6]) + val = np.zeros([6, 6]) assert _array(val).shape == (36,) assert _array66(val).shape == (6, 6) def test_element_string_ordering(): - d = elements.Drift('D0', 1, attr=numpy.array(0)) - assert d.__str__() == ("Drift:\n FamName: D0\n Length: 1.0\n" - " PassMethod: DriftPass\n attr: 0") + d = elements.Drift("D0", 1, attr=np.array(0)) + assert d.__str__() == ( + "Drift:\n FamName: D0\n Length: 1.0\n" + " PassMethod: DriftPass\n attr: 0" + ) assert d.__repr__() == "Drift('D0', 1.0, attr=array(0))" def test_element_creation_raises_exception(): with pytest.raises(ValueError): - elements.Element('family_name', R1='not_an_array') + elements.Element("family_name", R1="not_an_array") def test_base_element_methods(): - e = elements.Element('family_name') + e = elements.Element("family_name") assert e.divide([0.2, 0.5, 0.3]) == [e] assert id(e.copy()) != id(e) def test_argument_checks(): - q = elements.Quadrupole('quad', 1.0, 0.5) + q = elements.Quadrupole("quad", 1.0, 0.5) # Test type with pytest.raises(ValueError): - q.Length = 'a' + q.Length = "a" # Test shape with pytest.raises(ValueError): q.T1 = [0.0, 0.0] @@ -51,94 +59,93 @@ def test_argument_checks(): def test_dipole(): - d = elements.Dipole('dipole', 1.0, 0.01) + d = elements.Dipole("dipole", 1.0, 0.01) assert d.MaxOrder == 0 assert len(d.PolynomA) == 2 assert d.K == 0.0 - d = elements.Dipole('dipole', 1.0, 0.01, -0.5) + d = elements.Dipole("dipole", 1.0, 0.01, -0.5) assert d.MaxOrder == 1 assert len(d.PolynomA) == 2 assert d.K == -0.5 - d = elements.Dipole('dipole', 1.0, 0.01, PolynomB=[0.0, 0.1, 0.0]) + d = elements.Dipole("dipole", 1.0, 0.01, PolynomB=[0.0, 0.1, 0.0]) assert d.MaxOrder == 1 assert len(d.PolynomA) == 3 assert d.K == 0.1 - d = elements.Dipole('dipole', 1.0, 0.01, PolynomB=[0.0, 0.0, 0.005]) + d = elements.Dipole("dipole", 1.0, 0.01, PolynomB=[0.0, 0.0, 0.005]) assert d.MaxOrder == 2 assert len(d.PolynomA) == 3 assert d.K == 0.0 - d = elements.Dipole('dipole', 1.0, 0.01, PolynomB=[0.0, 0.0, 0.005], - MaxOrder=0) + d = elements.Dipole("dipole", 1.0, 0.01, PolynomB=[0.0, 0.0, 0.005], MaxOrder=0) assert d.MaxOrder == 0 assert len(d.PolynomA) == 3 assert d.K == 0.0 def test_quadrupole(): - q = elements.Quadrupole('quadrupole', 1.0) + q = elements.Quadrupole("quadrupole", 1.0) assert q.MaxOrder == 1 assert len(q.PolynomA) == 2 assert q.K == 0.0 - q = elements.Quadrupole('quadrupole', 1.0, -0.5) + q = elements.Quadrupole("quadrupole", 1.0, -0.5) assert q.MaxOrder == 1 assert len(q.PolynomA) == 2 assert q.K == -0.5 - q = elements.Quadrupole('quadrupole', 1.0, PolynomB=[0.0, 0.0, 0.005]) + q = elements.Quadrupole("quadrupole", 1.0, PolynomB=[0.0, 0.0, 0.005]) assert q.MaxOrder == 2 assert len(q.PolynomA) == 3 assert q.K == 0.0 - q = elements.Quadrupole('quadrupole', 1.0, PolynomB=[0.0, 0.5, 0.005], - MaxOrder=1) + q = elements.Quadrupole("quadrupole", 1.0, PolynomB=[0.0, 0.5, 0.005], MaxOrder=1) assert q.MaxOrder == 1 assert len(q.PolynomA) == 3 assert q.K == 0.5 def test_sextupole(): - s = elements.Sextupole('sextupole', 1.0) + s = elements.Sextupole("sextupole", 1.0) assert s.MaxOrder == 2 assert len(s.PolynomA) == 3 assert s.H == 0.0 - s = elements.Sextupole('sextupole', 1.0, -0.5) + s = elements.Sextupole("sextupole", 1.0, -0.5) assert s.MaxOrder == 2 assert len(s.PolynomA) == 3 assert s.H == -0.5 - s = elements.Sextupole('sextupole', 1.0, PolynomB=[0.0, 0.0, 0.005, 0.0]) + s = elements.Sextupole("sextupole", 1.0, PolynomB=[0.0, 0.0, 0.005, 0.0]) assert s.MaxOrder == 2 assert len(s.PolynomA) == 4 assert s.H == 0.005 - s = elements.Sextupole('sextupole', 1.0, PolynomB=[0.0, 0.0, 0.005, 0.001]) + s = elements.Sextupole("sextupole", 1.0, PolynomB=[0.0, 0.0, 0.005, 0.001]) assert s.MaxOrder == 3 assert len(s.PolynomA) == 4 assert s.H == 0.005 - s = elements.Sextupole('sextupole', 1.0, PolynomB=[0.0, 0.5, 0.005, 0.001], - MaxOrder=2) + s = elements.Sextupole( + "sextupole", 1.0, PolynomB=[0.0, 0.5, 0.005, 0.001], MaxOrder=2 + ) assert s.MaxOrder == 2 assert len(s.PolynomA) == 4 assert s.H == 0.005 def test_octupole(): - o = elements.Octupole('octupole', 1.0, [], [0.0, 0.0, 0.0, 0.0]) + o = elements.Octupole("octupole", 1.0, [], [0.0, 0.0, 0.0, 0.0]) assert o.MaxOrder == 3 assert len(o.PolynomA) == 4 def test_thinmultipole(): - m = elements.ThinMultipole('thin', [], [0.0, 0.0, 0.0, 0.0]) + m = elements.ThinMultipole("thin", [], [0.0, 0.0, 0.0, 0.0]) assert m.MaxOrder == 0 assert len(m.PolynomA) == 4 - m = elements.ThinMultipole('thin', [], [0.0, 0.0, 1.0, 0.0]) + m = elements.ThinMultipole("thin", [], [0.0, 0.0, 1.0, 0.0]) assert m.MaxOrder == 2 assert len(m.PolynomA) == 4 def test_multipole(): - m = elements.Multipole('multi', 1.0, [], [0.0, 0.0, 0.0, 0.0]) + m = elements.Multipole("multi", 1.0, [], [0.0, 0.0, 0.0, 0.0]) assert m.Length == 1.0 assert m.MaxOrder == 0 assert m.NumIntSteps == 10 - assert m.PassMethod == 'StrMPoleSymplectic4Pass' + assert m.PassMethod == "StrMPoleSymplectic4Pass" @pytest.mark.parametrize( @@ -306,19 +313,20 @@ def test_sextupolar_strength_prioritisation(): def test_divide_splits_attributes_correctly(): - pre = elements.Drift('drift', 1) + pre = elements.Drift("drift", 1) post = pre.divide([0.2, 0.5, 0.3]) assert len(post) == 3 assert sum([e.Length for e in post]) == pre.Length - pre = elements.Dipole('dipole', 1, KickAngle=[0.5, -0.5], BendingAngle=0.2) + pre = elements.Dipole("dipole", 1, KickAngle=[0.5, -0.5], BendingAngle=0.2) post = pre.divide([0.2, 0.5, 0.3]) assert len(post) == 3 assert sum([e.Length for e in post]) == pre.Length assert sum([e.KickAngle[0] for e in post]) == pre.KickAngle[0] assert sum([e.KickAngle[1] for e in post]) == pre.KickAngle[1] assert sum([e.BendingAngle for e in post]) == pre.BendingAngle - pre = elements.RFCavity('rfc', 1, voltage=187500, frequency=3.5237e+8, - harmonic_number=31, energy=6.e+9) + pre = elements.RFCavity( + "rfc", 1, voltage=187500, frequency=3.5237e8, harmonic_number=31, energy=6.0e9 + ) post = pre.divide([0.2, 0.5, 0.3]) assert len(post) == 3 assert sum([e.Length for e in post]) == pre.Length @@ -327,46 +335,47 @@ def test_divide_splits_attributes_correctly(): def test_insert_into_drift(): # Create elements - drift = elements.Drift('drift', 1) - monitor = elements.Monitor('bpm') - quad = elements.Quadrupole('quad', 0.3) + drift = elements.Drift("drift", 1) + monitor = elements.Monitor("bpm") + quad = elements.Quadrupole("quad", 0.3) # Test None splitting behaviour - el_list = drift.insert([(0., None), (0.3, None), (0.7, None), (1., None)]) + el_list = drift.insert([(0.0, None), (0.3, None), (0.7, None), (1.0, None)]) assert len(el_list) == 3 - numpy.testing.assert_almost_equal([e.Length for e in el_list], - [0.3, 0.4, 0.3]) + np.testing.assert_almost_equal([e.Length for e in el_list], [0.3, 0.4, 0.3]) # Test normal insertion el_list = drift.insert([(0.3, monitor), (0.7, quad)]) assert len(el_list) == 5 - numpy.testing.assert_almost_equal([e.Length for e in el_list], - [0.3, 0.0, 0.25, 0.3, 0.15]) + np.testing.assert_almost_equal( + [e.Length for e in el_list], [0.3, 0.0, 0.25, 0.3, 0.15] + ) # Test insertion at either end produces -ve length drifts el_list = drift.insert([(0.0, quad), (1.0, quad)]) assert len(el_list) == 5 - numpy.testing.assert_almost_equal([e.Length for e in el_list], - [-0.15, 0.3, 0.7, 0.3, -0.15]) + np.testing.assert_almost_equal( + [e.Length for e in el_list], [-0.15, 0.3, 0.7, 0.3, -0.15] + ) -@pytest.mark.parametrize('func', (lattice_track, lattice_pass, internal_lpass)) +@pytest.mark.parametrize("func", (lattice_track, lattice_pass, internal_lpass)) def test_correct_dimensions_does_not_raise_error(rin, func): func([], rin, 1) - rin = numpy.zeros((6,)) + rin = np.zeros((6,)) func([], rin, 1) - rin = numpy.array(numpy.zeros((6, 2), order='F')) + rin = np.array(np.zeros((6, 2), order="F")) func([], rin, 1) @pytest.mark.parametrize("dipole_class", (elements.Dipole, elements.Bend)) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_dipole_bend_synonym(rin, dipole_class, func): - b = dipole_class('dipole', 1.0, 0.1, EntranceAngle=0.05, ExitAngle=0.05) + b = dipole_class("dipole", 1.0, 0.1, EntranceAngle=0.05, ExitAngle=0.05) rin[0, 0] = 1e-6 if func == element_track: func(b, rin, in_place=True) else: func(b, rin) - rin_expected = numpy.array([1e-6, 0, 0, 0, 0, 1e-7]).reshape((6, 1)) - numpy.testing.assert_almost_equal(rin, rin_expected) + rin_expected = np.array([1e-6, 0, 0, 0, 0, 1e-7]).reshape((6, 1)) + np.testing.assert_almost_equal(rin, rin_expected) assert b.K == 0.0 b.PolynomB[1] = 0.2 assert b.K == 0.2 @@ -374,35 +383,35 @@ def test_dipole_bend_synonym(rin, dipole_class, func): assert b.PolynomB[1] == 0.1 -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_marker(rin, func): - m = elements.Marker('marker') + m = elements.Marker("marker") assert m.Length == 0 - rin = numpy.array(numpy.random.rand(*rin.shape), order='F') - rin_orig = numpy.array(rin, copy=True, order='F') + rin = np.array(np.random.rand(*rin.shape), order="F") + rin_orig = np.array(rin, copy=True, order="F") if func == element_track: func(m, rin, in_place=True) else: func(m, rin) - numpy.testing.assert_equal(rin, rin_orig) + np.testing.assert_equal(rin, rin_orig) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_monitor(rin, func): - mon = elements.Monitor('monitor') + mon = elements.Monitor("monitor") assert mon.Length == 0 - rin = numpy.array(numpy.random.rand(*rin.shape), order='F') + rin = np.array(np.random.rand(*rin.shape), order="F") rin_orig = rin.copy() if func == element_track: func(mon, rin, in_place=True) else: func(mon, rin) - numpy.testing.assert_equal(rin, rin_orig) + np.testing.assert_equal(rin, rin_orig) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_aperture_inside_limits(rin, func): - a = elements.Aperture('aperture', [-1e-3, 1e-3, -1e-4, 1e-4]) + a = elements.Aperture("aperture", [-1e-3, 1e-3, -1e-4, 1e-4]) assert a.Length == 0 rin[0, 0] = 1e-5 rin[2, 0] = -1e-5 @@ -411,12 +420,12 @@ def test_aperture_inside_limits(rin, func): func(a, rin, in_place=True) else: func(a, rin) - numpy.testing.assert_equal(rin, rin_orig) + np.testing.assert_equal(rin, rin_orig) -@pytest.mark.parametrize('func', (lattice_track, lattice_pass, internal_lpass)) +@pytest.mark.parametrize("func", (lattice_track, lattice_pass, internal_lpass)) def test_aperture_outside_limits(rin, func): - a = elements.Aperture('aperture', [-1e-3, 1e-3, -1e-4, 1e-4]) + a = elements.Aperture("aperture", [-1e-3, 1e-3, -1e-4, 1e-4]) assert a.Length == 0 lattice = [a] rin[0, 0] = 1e-2 @@ -425,13 +434,13 @@ def test_aperture_outside_limits(rin, func): func(lattice, rin, in_place=True) else: func(lattice, rin) - assert numpy.isnan(rin[0, 0]) + assert np.isnan(rin[0, 0]) assert rin[2, 0] == 0.0 # Only the 1st coordinate is nan, the rest is zero -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_drift_offset(rin, func): - d = elements.Drift('drift', 1) + d = elements.Drift("drift", 1) rin[0, 0] = 1e-6 rin[2, 0] = 2e-6 rin_orig = rin.copy() @@ -439,12 +448,12 @@ def test_drift_offset(rin, func): func(d, rin, in_place=True) else: func(d, rin) - numpy.testing.assert_equal(rin, rin_orig) + np.testing.assert_equal(rin, rin_orig) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_drift_divergence(rin, func): - d = elements.Drift('drift', 1.0) + d = elements.Drift("drift", 1.0) assert d.Length == 1 rin[1, 0] = 1e-6 rin[3, 0] = -2e-6 @@ -453,16 +462,15 @@ def test_drift_divergence(rin, func): else: func(d, rin) # results from Matlab - rin_expected = numpy.array([1e-6, 1e-6, -2e-6, -2e-6, 0, - 2.5e-12]).reshape(6, 1) - numpy.testing.assert_equal(rin, rin_expected) + rin_expected = np.array([1e-6, 1e-6, -2e-6, -2e-6, 0, 2.5e-12]).reshape(6, 1) + np.testing.assert_equal(rin, rin_expected) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_drift_two_particles(rin, func): - d = elements.Drift('drift', 1.0) + d = elements.Drift("drift", 1.0) assert d.Length == 1 - two_rin = numpy.array(numpy.concatenate((rin, rin), axis=1), order='F') + two_rin = np.array(np.concatenate((rin, rin), axis=1), order="F") # particle one is offset two_rin[0, 0] = 1e-6 two_rin[2, 0] = 2e-6 @@ -475,24 +483,27 @@ def test_drift_two_particles(rin, func): else: func(d, two_rin) # results from Matlab - p1_expected = numpy.array(two_rin_orig[:, 0]).reshape(6, 1) - p2_expected = numpy.array([1e-6, 1e-6, -2e-6, -2e-6, 0, - 2.5e-12]).reshape(6, 1) - two_rin_expected = numpy.concatenate((p1_expected, p2_expected), axis=1) - numpy.testing.assert_equal(two_rin, two_rin_expected) + p1_expected = np.array(two_rin_orig[:, 0]).reshape(6, 1) + p2_expected = np.array([1e-6, 1e-6, -2e-6, -2e-6, 0, 2.5e-12]).reshape(6, 1) + two_rin_expected = np.concatenate((p1_expected, p2_expected), axis=1) + np.testing.assert_equal(two_rin, two_rin_expected) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_quad(rin, func): - q = elements.Quadrupole('quad', 0.4, k=1) + q = elements.Quadrupole("quad", 0.4, k=1) rin[0, 0] = 1e-6 if func == element_track: func(q, rin, in_place=True) else: func(q, rin) - expected = numpy.array([0.9210610203854122, -0.3894182419439, 0, - 0, 0, 0.0000000103303797478]).reshape(6, 1) * 1e-6 - numpy.testing.assert_allclose(rin, expected) + expected = ( + np.array( + [0.9210610203854122, -0.3894182419439, 0, 0, 0, 0.0000000103303797478] + ).reshape(6, 1) + * 1e-6 + ) + np.testing.assert_allclose(rin, expected) assert q.K == 1 q.PolynomB[1] = 0.2 assert q.K == 0.2 @@ -500,40 +511,41 @@ def test_quad(rin, func): assert q.PolynomB[1] == 0.1 -@pytest.mark.parametrize('func', (lattice_track, lattice_pass, internal_lpass)) -def test_rfcavity(rin, func): - rf = elements.RFCavity('rfcavity', 0.0, 187500, 3.5237e+8, 31, 6.e+9) - rflattice = [rf, rf, rf, rf] +@pytest.mark.parametrize("func", [lattice_track, lattice_pass, internal_lpass]) +def test_rfcavity(rin: np.ndarray, func: any) -> None: + rf0 = elements.RFCavity("rfcavity", 0.0, 187500, 3.5237e8, 31, 6.0e9) + rflattice = [rf0, rf0, rf0, rf0] rin[4, 0] = 1e-6 rin[5, 0] = 1e-6 if func == lattice_track: func(rflattice, rin, in_place=True) else: func(rflattice, rin) - expected = numpy.array([0., 0., 0., 0., 9.990769e-7, 1.e-6]).reshape(6, 1) - numpy.testing.assert_allclose(rin, expected, atol=1e-12) + result = np.array([0.0, 0.0, 0.0, 0.0, 9.990769e-7, 1.0e-6]).reshape(6, 1) + np.testing.assert_allclose(rin, result, atol=1e-12) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) @pytest.mark.parametrize("n", (0, 1, 2, 3, 4, 5)) def test_m66(rin, n, func): - m = numpy.random.rand(6, 6) - m66 = elements.M66('m66', m) + m = np.random.rand(6, 6) + m66 = elements.M66("m66", m) assert m66.Length == 0 rin[n, 0] = 1e-6 if func == element_track: func(m66, rin, in_place=True) else: func(m66, rin) - expected = numpy.array([m[0, n], m[1, n], m[2, n], m[3, n], m[4, n], - m[5, n]]).reshape(6, 1) * 1e-6 - numpy.testing.assert_equal(rin, expected) + expected = ( + np.array([m[0, n], m[1, n], m[2, n], m[3, n], m[4, n], m[5, n]]).reshape(6, 1) + * 1e-6 + ) + np.testing.assert_equal(rin, expected) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_corrector(rin, func): - c = elements.Corrector('corrector', 0.0, numpy.array([0.9, 0.5], - dtype=numpy.float64)) + c = elements.Corrector("corrector", 0.0, np.array([0.9, 0.5], dtype=np.float64)) assert c.Length == 0 rin[0, 0] = 1e-6 rin_orig = rin.copy() @@ -543,29 +555,31 @@ def test_corrector(rin, func): func(c, rin, in_place=True) else: func(c, rin) - numpy.testing.assert_equal(rin, rin_orig) + np.testing.assert_equal(rin, rin_orig) -@pytest.mark.parametrize('func', (element_track, element_pass, internal_epass)) +@pytest.mark.parametrize("func", (element_track, element_pass, internal_epass)) def test_wiggler(rin, func): period = 0.05 periods = 23 bmax = 1 - by = numpy.array([1, 1, 0, 1, 1, 0], dtype=numpy.float64) - c = elements.Wiggler('wiggler', period * periods, period, bmax, By=by) + by = np.array([1, 1, 0, 1, 1, 0], dtype=np.float64) + c = elements.Wiggler("wiggler", period * periods, period, bmax, By=by) assert abs(c.Length - 1.15) < 1e-10 # Expected value from Matlab AT. - expected = numpy.array(rin, copy=True) + expected = np.array(rin, copy=True) expected[5] = 0.000000181809691064259 if func == element_track: func(c, rin, energy=3e9, in_place=True) else: func(c, rin, energy=3e9) - numpy.testing.assert_allclose(rin, expected, atol=1e-12) + np.testing.assert_allclose(rin, expected, atol=1e-12) def test_exit_entrance(): - q = elements.Quadrupole('quad', 0.4, k=1) + q = elements.Quadrupole("quad", 0.4, k=1) for kin, kout in zip(q._entrance_fields, q._exit_fields): - assert_equal(kin.replace('Entrance', ''). replace('1', ''), - kout.replace('Exit', '').replace('2', '')) + np.testing.assert_equal( + kin.replace("Entrance", "").replace("1", ""), + kout.replace("Exit", "").replace("2", ""), + ) From 1c06e640175dfd6e268eed6b69bd160a6d6578a1 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 22 Dec 2025 18:32:16 +0100 Subject: [PATCH 22/56] remove double definition atEnergy atGamma --- atintegrators/atelem.c | 31 +++++-------------------------- 1 file changed, 5 insertions(+), 26 deletions(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index e8a60183a..f984cfaf7 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -74,27 +74,6 @@ typedef mxArray atElem; #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) { @@ -350,6 +329,11 @@ static double *atGetOptionalDoubleArray(const PyObject *element, char *name) return atGetOptionalDoubleArraySz(element, name, &msz, &nsz); } + +#endif /* defined(PYAT) */ + +#if defined(PYAT) || defined(MATLAB_MEX_FILE) + double atEnergy(double ringenergy, double elemenergy) { if (ringenergy!=0.0) @@ -372,11 +356,6 @@ double atGamma(double ringenergy, double elemenergy, double rest_energy) return energy / rest_energy; } - -#endif /* defined(PYAT) */ - -#if defined(PYAT) || defined(MATLAB_MEX_FILE) - #ifdef __cplusplus #define C_LINK extern "C" #else From ba8c5545f45d7aa7586522e9cb01e92364ac4767 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:12:50 +0100 Subject: [PATCH 23/56] restore atelem --- atintegrators/atelem.c | 46 +++++++++++++++++++++--------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index f984cfaf7..8b597af37 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -74,6 +74,27 @@ typedef mxArray atElem; #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) { @@ -185,6 +206,8 @@ 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; @@ -329,33 +352,10 @@ static double *atGetOptionalDoubleArray(const PyObject *element, char *name) return atGetOptionalDoubleArraySz(element, name, &msz, &nsz); } - #endif /* defined(PYAT) */ #if defined(PYAT) || defined(MATLAB_MEX_FILE) -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; -} - #ifdef __cplusplus #define C_LINK extern "C" #else From 616152c3d3ce66885a37f8580b92152bca67dd60 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:27:36 +0100 Subject: [PATCH 24/56] add check_error --- atintegrators/atelem.c | 1 + 1 file changed, 1 insertion(+) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 8b597af37..7298699b2 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -83,6 +83,7 @@ double atEnergy(double ringenergy, double elemenergy) return elemenergy; else { atError("Energy not defined."); + check_error(); return 0.0; /* Never reached but makes the compiler happy */ } } From 9234171b9977ff9350452ba2914389979df9a771 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:27:58 +0100 Subject: [PATCH 25/56] check gamma == 0 --- atintegrators/GWigSymplecticPass.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index eb6314528..9299cf9d6 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -146,7 +146,10 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, /* Check energy */ gamma = atGamma(Param->energy, Energy, Param->rest_energy); - check_error(); + if (gamma == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Energy=Energy; @@ -166,7 +169,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T2=T2; }else{ gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); - check_error(); }; GWigSymplecticPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, From 86712c13e67652df1c41ef0ddc90bbd1c4ddf767 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:32:25 +0100 Subject: [PATCH 26/56] check gamma == 0 --- atintegrators/GWigSymplecticRadPass.c | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index 5afba37f5..b525f6fb0 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -270,6 +270,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); + /* Check energy */ + gamma = atGamma(Param->energy, Energy, Param->rest_energy); + if (gamma == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Energy=Energy; Elem->Length=Ltot; @@ -286,9 +293,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->R2=R2; Elem->T1=T1; Elem->T2=T2; - } - gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); - check_error(); + }else{ + gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); + }; GWigSymplecticRadPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm, From dbe20a5fd45421b74778a26e30bec69fdeb8dc88 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:44:47 +0100 Subject: [PATCH 27/56] check Energy --- atintegrators/BeamLoadingCavityPass.c | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index 63afbc0e3..5cbaa3c25 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -224,6 +224,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, z_cuts=atGetOptionalDoubleArray(ElemData,"ZCuts"); check_error(); 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(); + } + int dimsth[] = {Param->nbunch*nslice*nturns, 4}; atCheckArrayDims(ElemData,"_turnhistory", 2, dimsth); check_error(); int dimsvb[] = {Param->nbunch, 2}; @@ -260,9 +267,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->vbunch_buffer = vbunch_buffer; Elem->feedback_angle_offset = feedback_angle_offset; Elem->fbmode = fbmode; - } - energy = atEnergy(Param->energy, Elem->Energy); - check_error(); + }else{ + energy = atEnergy(Param->energy, Elem->Energy); + }; if(num_particlesnbunch){ atError("Number of particles has to be greater or equal to the number of bunches."); check_error(); From 2c3033443876e68e950e862ab9722168fba1bb53 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:51:00 +0100 Subject: [PATCH 28/56] check Energy --- atintegrators/BndMPoleSymplectic4E2RadPass.c | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/atintegrators/BndMPoleSymplectic4E2RadPass.c b/atintegrators/BndMPoleSymplectic4E2RadPass.c index 9ff5154ea..d446a4330 100644 --- a/atintegrators/BndMPoleSymplectic4E2RadPass.c +++ b/atintegrators/BndMPoleSymplectic4E2RadPass.c @@ -232,6 +232,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); 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(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -259,7 +266,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; 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, From d462c70a6997c10c25c9346a0125a3556d86591f Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:52:21 +0100 Subject: [PATCH 29/56] check Energy --- atintegrators/BeamLoadingCavityPass.c | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index 5cbaa3c25..ce4394047 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -223,7 +223,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); z_cuts=atGetOptionalDoubleArray(ElemData,"ZCuts"); check_error(); feedback_angle_offset=atGetOptionalDouble(ElemData,"feedback_angle_offset", 0.0); check_error(); - + /* Check energy */ Energy = atEnergy(Param->energy, Energy); if (Energy == 0) { @@ -267,9 +267,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->vbunch_buffer = vbunch_buffer; Elem->feedback_angle_offset = feedback_angle_offset; Elem->fbmode = fbmode; - }else{ - energy = atEnergy(Param->energy, Elem->Energy); }; + energy = atEnergy(Param->energy, Elem->Energy); if(num_particlesnbunch){ atError("Number of particles has to be greater or equal to the number of bunches."); check_error(); From 7eb9905aa9e885d6f3eec7b2a364ece05116a7ec Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:54:38 +0100 Subject: [PATCH 30/56] check Energy --- atintegrators/BndMPoleSymplectic4QuantPass.c | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/atintegrators/BndMPoleSymplectic4QuantPass.c b/atintegrators/BndMPoleSymplectic4QuantPass.c index 0b6922102..d94586046 100644 --- a/atintegrators/BndMPoleSymplectic4QuantPass.c +++ b/atintegrators/BndMPoleSymplectic4QuantPass.c @@ -215,6 +215,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); 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(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -246,7 +253,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; 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, From baf49c1077158aac6afe8e329280b7b23d5e3dfe Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 12:58:06 +0100 Subject: [PATCH 31/56] check Energy --- atintegrators/BndMPoleSymplectic4RadPass.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index 94a71f3fd..7d6042cbb 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -168,6 +168,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); 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(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -199,7 +206,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; 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, @@ -268,7 +274,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); BndMPoleSymplectic4RadPass(r_in, Length, irho, PolynomA, PolynomB, From 04ce66ac6dfc9f8bb9d6c0c2fc82871e129ee321 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:01:21 +0100 Subject: [PATCH 32/56] check Energy --- atintegrators/CrabCavityPass.c | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/atintegrators/CrabCavityPass.c b/atintegrators/CrabCavityPass.c index 6e8474742..0bb3b19ba 100644 --- a/atintegrators/CrabCavityPass.c +++ b/atintegrators/CrabCavityPass.c @@ -101,6 +101,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, SigPhi = atGetOptionalDouble(ElemData,"SigPhi",0.0); check_error(); 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(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->Vx=Voltages[0]; From 48cd9fdfbdc526181e71bf7ba52c6ba7604fb299 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:01:56 +0100 Subject: [PATCH 33/56] remove check_error --- atintegrators/CrabCavityPass.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/CrabCavityPass.c b/atintegrators/CrabCavityPass.c index 0bb3b19ba..d03686752 100644 --- a/atintegrators/CrabCavityPass.c +++ b/atintegrators/CrabCavityPass.c @@ -126,7 +126,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, SigVV=Elem->SigVV; } energy = atEnergy(Param->energy, Elem->Energy); - check_error(); nvx = Elem->Vx; nvy = Elem->Vy; From d2327cce8fee36e83f0632396b09cce2ab51ac52 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:03:34 +0100 Subject: [PATCH 34/56] check Energy --- atintegrators/ExactMultipoleRadPass.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index a39127d20..a1b1a1b37 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -138,6 +138,14 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, atError("NumIntSteps must be positive"); check_error(); } + /* Check energy */ + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + + Elem = (struct elem *)atMalloc(sizeof(struct elem)); Elem->Length = Length; Elem->PolynomA = PolynomA; @@ -158,7 +166,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, Elem->KickAngle = KickAngle; } 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, From c47d255b792c73e28a5f320b091927254509d4b5 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:04:14 +0100 Subject: [PATCH 35/56] remove check_error --- atintegrators/ExactMultipoleRadPass.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index a1b1a1b37..bf4686137 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -216,7 +216,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); multipole_pass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, From 71d94e6b31068f87e60a3e7dda3ed019a9375e44 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:06:08 +0100 Subject: [PATCH 36/56] check Energy --- atintegrators/ExactRectBendRadPass.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index a520584da..723e10047 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -173,6 +173,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, atError("NumIntSteps == 0 not allowed with radiation"); check_error(); } + /* Check energy */ + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -201,7 +208,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); - check_error(); ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, @@ -265,7 +271,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); ExactRectangularBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, From 348216bc66b160696f3bf2b0bb71cda6fafe7d0f Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:08:27 +0100 Subject: [PATCH 37/56] check Energy --- atintegrators/ExactRectangularBendRadPass.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index 8ddb7885e..dff706fb9 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -175,6 +175,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, atError("NumIntSteps must be positive"); check_error(); } + /* Check energy */ + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -203,7 +210,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } 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, @@ -266,7 +272,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); ExactRectangularBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, From 69908c562cbfa33bd4e3e31946c2fd8ab9bceffd Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:09:44 +0100 Subject: [PATCH 38/56] check Energy --- atintegrators/ExactSectorBendRadPass.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index f44e450c7..e194b9438 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -160,6 +160,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, atError("NumIntSteps == 0 not allowed with radiation"); check_error(); } + /* Check energy */ + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -187,7 +194,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } irho = Elem->BendingAngle/Elem->Length; gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); - check_error(); ExactSectorBendRad(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, @@ -251,7 +257,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxDuplicateArray(prhs[1]); irho = BendingAngle/Length; Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); ExactSectorBendRad(r_in, Length, irho, PolynomA, PolynomB, From f975807dce89b30e1c460c469205c4e6fe272fb8 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:13:08 +0100 Subject: [PATCH 39/56] check Energy --- atintegrators/GWigSymplecticPass.c | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 9299cf9d6..9872e52fa 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -145,8 +145,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); /* Check energy */ - gamma = atGamma(Param->energy, Energy, Param->rest_energy); - if (gamma == 0) { + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); check_error(); } @@ -167,9 +167,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->R2=R2; Elem->T1=T1; Elem->T2=T2; - }else{ - gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); }; + gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); GWigSymplecticPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm, From 73408b6015d6a1e61d7e3964874b15c3953feaff Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:14:41 +0100 Subject: [PATCH 40/56] check Energy --- atintegrators/GWigSymplecticRadPass.c | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index b525f6fb0..8f92900a9 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -271,8 +271,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); /* Check energy */ - gamma = atGamma(Param->energy, Energy, Param->rest_energy); - if (gamma == 0) { + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); check_error(); } @@ -293,9 +293,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->R2=R2; Elem->T1=T1; Elem->T2=T2; - }else{ - gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); }; + gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); GWigSymplecticRadPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm, From 88aa4d00622c1b1736c3fe01bed58609c8c62578 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:16:01 +0100 Subject: [PATCH 41/56] remove check_error --- atintegrators/GWigSymplecticPass.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 9872e52fa..2f12c1e2c 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -223,7 +223,6 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); GWigSymplecticPass(r_in, Gamma, Ltot, Lw, Bmax, Nstep, Nmeth, NHharm, From 68984dc21728dbdd45b56170c9086375ac5e1359 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:16:17 +0100 Subject: [PATCH 42/56] remove check_error --- atintegrators/GWigSymplecticRadPass.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index 8f92900a9..cdac5ba27 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -348,7 +348,6 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); GWigSymplecticRadPass(r_in, Gamma, Ltot, Lw, Bmax, Nstep, Nmeth, From 51a498de47818cb98693d81920941de6b438bf1e Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:20:03 +0100 Subject: [PATCH 43/56] check Energy --- atintegrators/RFCavityPass.c | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index 3e63ec3dd..a14656e2e 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -49,8 +49,11 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); /* Check energy */ - energy = atEnergy(energy, Energy); - check_error(); + Energy = atEnergy(Param->energy, Energy); + if (Energy == 0) { + atError("Energy needs to be defined. Check lattice parameters or pass method options.\n"); + check_error(); + } Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; @@ -62,7 +65,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->PhaseLag=PhaseLag; } 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); @@ -101,7 +103,6 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* Check energy */ Energy = atEnergy(Energy, 0); - check_error(); RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); From 1c7fbaabeff0167b62eadac742f4e510f2989a40 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:21:54 +0100 Subject: [PATCH 44/56] check Energy --- atintegrators/StrMPoleSymplectic4QuantPass.c | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/atintegrators/StrMPoleSymplectic4QuantPass.c b/atintegrators/StrMPoleSymplectic4QuantPass.c index 32d5b37b4..00b1eab6c 100644 --- a/atintegrators/StrMPoleSymplectic4QuantPass.c +++ b/atintegrators/StrMPoleSymplectic4QuantPass.c @@ -184,6 +184,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); 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(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -206,7 +213,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } energy = atEnergy(Param->energy, Elem->Energy); - check_error(); StrMPoleSymplectic4QuantPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, From b2eff678e4f50aab46084de9955f82527a6a258a Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 13:23:18 +0100 Subject: [PATCH 45/56] check Energy --- atintegrators/StrMPoleSymplectic4RadPass.c | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index 02c6e69fd..c292829d4 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/atintegrators/StrMPoleSymplectic4RadPass.c @@ -144,6 +144,13 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); 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(); + } + Elem = (struct elem*)atMalloc(sizeof(struct elem)); Elem->Length=Length; Elem->PolynomA=PolynomA; @@ -166,7 +173,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); - check_error(); StrMPoleSymplectic4RadPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, @@ -221,7 +227,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); Gamma = atGamma(Energy, Energy, rest_energy); - check_error(); r_in = mxGetDoubles(plhs[0]); StrMPoleSymplectic4RadPass(r_in, Length, PolynomA, PolynomB, From 3cfdb4e2b2153791e03dd9f0ae10575790fd54ff Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 14:50:10 +0100 Subject: [PATCH 46/56] remove comment --- atintegrators/RFCavityPass.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index a14656e2e..da82fcfa9 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -101,7 +101,6 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); - /* Check energy */ Energy = atEnergy(Energy, 0); RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); From d6139ad7696eb90a9d5259fa5b8afa3330720d22 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 14:50:55 +0100 Subject: [PATCH 47/56] test rf element track --- pyat/test/test_basic_elements.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/pyat/test/test_basic_elements.py b/pyat/test/test_basic_elements.py index be8dae96f..820a4eed0 100644 --- a/pyat/test/test_basic_elements.py +++ b/pyat/test/test_basic_elements.py @@ -511,17 +511,16 @@ def test_quad(rin, func): assert q.PolynomB[1] == 0.1 -@pytest.mark.parametrize("func", [lattice_track, lattice_pass, internal_lpass]) +@pytest.mark.parametrize("func", [element_track, element_pass, internal_epass]) def test_rfcavity(rin: np.ndarray, func: any) -> None: - rf0 = elements.RFCavity("rfcavity", 0.0, 187500, 3.5237e8, 31, 6.0e9) - rflattice = [rf0, rf0, rf0, rf0] + rfc = elements.RFCavity("rfcavity", 0.0, 187500, 3.5237e8, 31, 6.0e9) rin[4, 0] = 1e-6 rin[5, 0] = 1e-6 - if func == lattice_track: - func(rflattice, rin, in_place=True) + if func == element_track: + func(rfc, rin, in_place=True, energy=+6.0e9) else: - func(rflattice, rin) - result = np.array([0.0, 0.0, 0.0, 0.0, 9.990769e-7, 1.0e-6]).reshape(6, 1) + func(rfc, rin, energy=+6.0e9) + result = np.array([0.0, 0.0, 0.0, 0.0, 9.99769215e-07, 1.0e-6]).reshape(6, 1) np.testing.assert_allclose(rin, result, atol=1e-12) From 294d49580ec0113b65a58b6d9f5bb36a4a5c892a Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 15:07:55 +0100 Subject: [PATCH 48/56] flake8 --- pyat/test/test_basic_elements.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/pyat/test/test_basic_elements.py b/pyat/test/test_basic_elements.py index 820a4eed0..30fe460b3 100644 --- a/pyat/test/test_basic_elements.py +++ b/pyat/test/test_basic_elements.py @@ -513,14 +513,22 @@ def test_quad(rin, func): @pytest.mark.parametrize("func", [element_track, element_pass, internal_epass]) def test_rfcavity(rin: np.ndarray, func: any) -> None: + """Test rf cavity tracking. + + Arguments: + func: element track method + rin: np.array of dims (6,1) + """ rfc = elements.RFCavity("rfcavity", 0.0, 187500, 3.5237e8, 31, 6.0e9) rin[4, 0] = 1e-6 rin[5, 0] = 1e-6 - if func == element_track: + if func is element_track: func(rfc, rin, in_place=True, energy=+6.0e9) else: func(rfc, rin, energy=+6.0e9) + result = np.array([0.0, 0.0, 0.0, 0.0, 9.99769215e-07, 1.0e-6]).reshape(6, 1) + np.testing.assert_allclose(rin, result, atol=1e-12) From 55bf4a89255478e58c0e1509d012b41f70bca8a5 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 15:11:06 +0100 Subject: [PATCH 49/56] remove ; --- atintegrators/BeamLoadingCavityPass.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index ce4394047..3f435647a 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -267,7 +267,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->vbunch_buffer = vbunch_buffer; Elem->feedback_angle_offset = feedback_angle_offset; Elem->fbmode = fbmode; - }; + } energy = atEnergy(Param->energy, Elem->Energy); if(num_particlesnbunch){ From 4c8ef7b44853a6aabe96210e7e4daa9df1bc949e Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 15:12:20 +0100 Subject: [PATCH 50/56] remove extra space --- atintegrators/ExactMultipoleRadPass.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index bf4686137..3beab66f4 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -145,7 +145,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, check_error(); } - Elem = (struct elem *)atMalloc(sizeof(struct elem)); Elem->Length = Length; Elem->PolynomA = PolynomA; From f96bf9647d0cbf0ecb7e40374dc78bb92d72b907 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 15:13:35 +0100 Subject: [PATCH 51/56] remove ; --- atintegrators/GWigSymplecticPass.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 2f12c1e2c..dcc40b619 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -167,7 +167,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->R2=R2; Elem->T1=T1; Elem->T2=T2; - }; + } gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); GWigSymplecticPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, From 6f4f9b65a168c8c48137d61980c3db27d2a7ddef Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 15:13:58 +0100 Subject: [PATCH 52/56] remove ; --- atintegrators/GWigSymplecticRadPass.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index cdac5ba27..577389640 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -293,7 +293,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->R2=R2; Elem->T1=T1; Elem->T2=T2; - }; + } gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy); GWigSymplecticRadPass(r_in, gamma, Elem->Length, Elem->Lw, Elem->Bmax, From 9a50c1407b66b1e7e7dc2c789d2713601ca91f08 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 15:17:08 +0100 Subject: [PATCH 53/56] remove atEnergy --- atintegrators/RFCavityPass.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index da82fcfa9..f12e9fe7e 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -101,8 +101,6 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); - Energy = atEnergy(Energy, 0); - RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); } From 31f060def2d3cbe0621636f9edd1ae258479f397 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 17:49:33 +0100 Subject: [PATCH 54/56] hard code pass methods in elempass.m --- atmat/attrack/elempass.m | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m index 8d567a71c..8f13b0bff 100644 --- a/atmat/attrack/elempass.m +++ b/atmat/attrack/elempass.m @@ -26,6 +26,28 @@ props.Particle=particle; end +method_req_energy = { 'BeamLoadingCavityPass' , ... + 'BndMPoleSymplectic4E2RadPass', ... + 'BndMPoleSymplectic4QuantPass', ... + 'BndMPoleSymplectic4RadPass', ... + 'CrabCavityPass', ... + 'ExactMultipoleRadPass', ... + 'ExactRectangularBendRadPass', ... + 'ExactRectBendRadPass', ... + 'ExactSectorBendRadPass', ... + 'GWigSymplecticPass', ... + 'GWigSymplecticRadPass', ... + 'RFCavityPass', ... + 'StrMPoleSymplectic4QuantPass', ... + 'StrMPoleSymplectic4RadPass' ... + }; + +if any(strcmp(method_req_energy, methodname)) + if props.Energy <= 0 + error("Energy parameter must be defined."); + end +end + rout = feval(methodname,elem,rin,props); end \ No newline at end of file From 87e8de7abe1c33a7e40620ce551f7274a1d0d394 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 17:59:00 +0100 Subject: [PATCH 55/56] fix energy --- atmat/atphysics/Radiation/atdiffmat.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atmat/atphysics/Radiation/atdiffmat.m b/atmat/atphysics/Radiation/atdiffmat.m index 1124c542b..8546cb39d 100644 --- a/atmat/atphysics/Radiation/atdiffmat.m +++ b/atmat/atphysics/Radiation/atdiffmat.m @@ -56,7 +56,7 @@ end % Calculate 6-by-6 linear transfer matrix in each element % near the equilibrium orbit - m=findelemm66(elem,passm,orbit); + m=findelemm66(elem,passm,orbit,'Energy',energy); % Cumulative diffusion matrix of the entire ring BCUM = m*BCUM*m' + bdiff; btx=BCUM; From d1778141114f7bcff275f5c1b0d56a6f89b67cd1 Mon Sep 17 00:00:00 2001 From: oscarxblanco Date: Mon, 5 Jan 2026 18:06:29 +0100 Subject: [PATCH 56/56] remove unnecessary check_error --- atintegrators/atelem.c | 1 - 1 file changed, 1 deletion(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 7298699b2..8b597af37 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -83,7 +83,6 @@ double atEnergy(double ringenergy, double elemenergy) return elemenergy; else { atError("Energy not defined."); - check_error(); return 0.0; /* Never reached but makes the compiler happy */ } }