From e4ed2995cac641cbb6af2b2f6919de9f398bd8f0 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Tue, 13 Aug 2024 16:38:28 +0200 Subject: [PATCH 01/18] Allow more than 2 arguments in mex files --- atintegrators/AperturePass.c | 2 +- atintegrators/BeamLoadingCavityPass.c | 2 +- atintegrators/BeamMomentsPass.c | 2 +- atintegrators/BendLinearPass.c | 2 +- atintegrators/BndMPoleSymplectic4E2Pass.c | 2 +- atintegrators/BndMPoleSymplectic4E2RadPass.c | 2 +- atintegrators/BndMPoleSymplectic4Pass.c | 2 +- atintegrators/BndMPoleSymplectic4QuantPass.c | 2 +- atintegrators/BndMPoleSymplectic4RadPass.c | 2 +- atintegrators/BndStrMPoleSymplectic4Pass.c | 2 +- atintegrators/CavityPass.c | 2 +- atintegrators/ChangePRefPass.c | 2 +- atintegrators/CorrectorPass.c | 2 +- atintegrators/DeltaQPass.c | 2 +- atintegrators/DriftPass.c | 2 +- atintegrators/EAperturePass.c | 2 +- atintegrators/EnergyLossRadPass.c | 2 +- atintegrators/ExactDriftPass.c | 2 +- atintegrators/ExactMultipolePass.c | 2 +- atintegrators/ExactMultipoleRadPass.c | 2 +- atintegrators/ExactRectBendPass.c | 2 +- atintegrators/ExactRectBendRadPass.c | 2 +- atintegrators/ExactRectangularBendPass.c | 2 +- atintegrators/ExactRectangularBendRadPass.c | 2 +- atintegrators/ExactSectorBendPass.c | 2 +- atintegrators/ExactSectorBendRadPass.c | 2 +- atintegrators/GWigSymplecticPass.c | 2 +- atintegrators/GWigSymplecticRadPass.c | 2 +- atintegrators/IdTablePass.c | 2 +- atintegrators/IdentityPass.c | 2 +- atintegrators/ImpedanceTablePass.c | 2 +- atintegrators/Matrix66Pass.c | 2 +- atintegrators/MatrixTijkPass.c | 2 +- atintegrators/QuantDiffPass.c | 2 +- atintegrators/SimpleQuantDiffPass.c | 2 +- atintegrators/SimpleRadiationPass.c | 2 +- atintegrators/SliceMomentsPass.c | 2 +- atintegrators/SolenoidLinearPass.c | 2 +- atintegrators/StrMPoleSymplectic4Pass.c | 2 +- atintegrators/StrMPoleSymplectic4QuantPass.c | 2 +- atintegrators/StrMPoleSymplectic4RadPass.c | 2 +- atintegrators/TestRandomPass.c | 2 +- atintegrators/ThinMPolePass.c | 2 +- atintegrators/VariableThinMPolePass.c | 2 +- atintegrators/WakeFieldPass.c | 2 +- atintegrators/WiggLinearPass.c | 2 +- 46 files changed, 46 insertions(+), 46 deletions(-) diff --git a/atintegrators/AperturePass.c b/atintegrators/AperturePass.c index c8101784d..9d03f4d2b 100644 --- a/atintegrators/AperturePass.c +++ b/atintegrators/AperturePass.c @@ -46,7 +46,7 @@ MODULE_DEF(AperturePass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index 4c95aabae..e549cce43 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -268,7 +268,7 @@ MODULE_DEF(BeamLoadingCavityPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if(nrhs == 2) + if(nrhs >= 2) { double *r_in; diff --git a/atintegrators/BeamMomentsPass.c b/atintegrators/BeamMomentsPass.c index 3af8ab98b..a8974dfe4 100644 --- a/atintegrators/BeamMomentsPass.c +++ b/atintegrators/BeamMomentsPass.c @@ -98,7 +98,7 @@ MODULE_DEF(BeamMomentsPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; diff --git a/atintegrators/BendLinearPass.c b/atintegrators/BendLinearPass.c index 5a96beb3d..bcf4aa44d 100644 --- a/atintegrators/BendLinearPass.c +++ b/atintegrators/BendLinearPass.c @@ -220,7 +220,7 @@ MODULE_DEF(BendLinearPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2 ) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, K, ByError, FringeInt1, FringeInt2, FullGap; double *R1, *R2, *T1, *T2; double *r_in; diff --git a/atintegrators/BndMPoleSymplectic4E2Pass.c b/atintegrators/BndMPoleSymplectic4E2Pass.c index f2a8c82e7..28329151d 100644 --- a/atintegrators/BndMPoleSymplectic4E2Pass.c +++ b/atintegrators/BndMPoleSymplectic4E2Pass.c @@ -232,7 +232,7 @@ MODULE_DEF(BndMPoleSymplectic4E2Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double irho; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2; int MaxOrder, NumIntSteps; diff --git a/atintegrators/BndMPoleSymplectic4E2RadPass.c b/atintegrators/BndMPoleSymplectic4E2RadPass.c index 0756450e8..076a32472 100644 --- a/atintegrators/BndMPoleSymplectic4E2RadPass.c +++ b/atintegrators/BndMPoleSymplectic4E2RadPass.c @@ -275,7 +275,7 @@ MODULE_DEF(BndMPoleSymplectic4E2RadPass) /* Dummy module initialisation * #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double irho; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps; diff --git a/atintegrators/BndMPoleSymplectic4Pass.c b/atintegrators/BndMPoleSymplectic4Pass.c index 91818201e..1aef9d8e2 100644 --- a/atintegrators/BndMPoleSymplectic4Pass.c +++ b/atintegrators/BndMPoleSymplectic4Pass.c @@ -214,7 +214,7 @@ MODULE_DEF(BndMPoleSymplectic4Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/BndMPoleSymplectic4QuantPass.c b/atintegrators/BndMPoleSymplectic4QuantPass.c index b542a70f9..107649e6f 100644 --- a/atintegrators/BndMPoleSymplectic4QuantPass.c +++ b/atintegrators/BndMPoleSymplectic4QuantPass.c @@ -265,7 +265,7 @@ MODULE_DEF(BndMPoleSymplectic4QuantPass) /* Dummy module initialisation * #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index acb47b660..13e43767c 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -212,7 +212,7 @@ MODULE_DEF(BndMPoleSymplectic4RadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/BndStrMPoleSymplectic4Pass.c b/atintegrators/BndStrMPoleSymplectic4Pass.c index 6f79d9e54..b626574b8 100644 --- a/atintegrators/BndStrMPoleSymplectic4Pass.c +++ b/atintegrators/BndStrMPoleSymplectic4Pass.c @@ -271,7 +271,7 @@ MODULE_DEF(BndStrMPoleSymplectic4Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, X0ref, ByError, RefDZ; int MaxOrder, NumIntSteps; diff --git a/atintegrators/CavityPass.c b/atintegrators/CavityPass.c index fe3e202f6..571c93ce7 100644 --- a/atintegrators/CavityPass.c +++ b/atintegrators/CavityPass.c @@ -92,7 +92,7 @@ MODULE_DEF(CavityPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ChangePRefPass.c b/atintegrators/ChangePRefPass.c index 1810a8379..0161f9f66 100644 --- a/atintegrators/ChangePRefPass.c +++ b/atintegrators/ChangePRefPass.c @@ -43,7 +43,7 @@ MODULE_DEF(ChangePRefPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/CorrectorPass.c b/atintegrators/CorrectorPass.c index 5e58cb46f..f5470c986 100644 --- a/atintegrators/CorrectorPass.c +++ b/atintegrators/CorrectorPass.c @@ -132,7 +132,7 @@ MODULE_DEF(CorrectorPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/DeltaQPass.c b/atintegrators/DeltaQPass.c index e92f34082..ef4c53dca 100644 --- a/atintegrators/DeltaQPass.c +++ b/atintegrators/DeltaQPass.c @@ -141,7 +141,7 @@ MODULE_DEF(DeltaQPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/DriftPass.c b/atintegrators/DriftPass.c index bdbeb9fbb..4935dc598 100644 --- a/atintegrators/DriftPass.c +++ b/atintegrators/DriftPass.c @@ -91,7 +91,7 @@ MODULE_DEF(DriftPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/EAperturePass.c b/atintegrators/EAperturePass.c index 3277961cc..66ca6e60f 100644 --- a/atintegrators/EAperturePass.c +++ b/atintegrators/EAperturePass.c @@ -44,7 +44,7 @@ MODULE_DEF(EAperturePass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/EnergyLossRadPass.c b/atintegrators/EnergyLossRadPass.c index 9d639584b..9dcfa840c 100644 --- a/atintegrators/EnergyLossRadPass.c +++ b/atintegrators/EnergyLossRadPass.c @@ -62,7 +62,7 @@ MODULE_DEF(EnergyLossRadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Check if the number of input arguments is correct. */ - if (nrhs == 2) { + if (nrhs >= 2) { /* Get the input arguments. */ double *r_in; const mxArray *ElemData = prhs[0]; diff --git a/atintegrators/ExactDriftPass.c b/atintegrators/ExactDriftPass.c index 2c04e14a5..5f6df8af0 100644 --- a/atintegrators/ExactDriftPass.c +++ b/atintegrators/ExactDriftPass.c @@ -83,7 +83,7 @@ MODULE_DEF(ExactDriftPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactMultipolePass.c b/atintegrators/ExactMultipolePass.c index de1ea5979..f0e3f8030 100644 --- a/atintegrators/ExactMultipolePass.c +++ b/atintegrators/ExactMultipolePass.c @@ -160,7 +160,7 @@ MODULE_DEF(ExactMultipolePass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index 248d4c8e9..0155b32af 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -163,7 +163,7 @@ MODULE_DEF(ExactMultipoleRadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactRectBendPass.c b/atintegrators/ExactRectBendPass.c index d1460dd90..f8ab7a859 100644 --- a/atintegrators/ExactRectBendPass.c +++ b/atintegrators/ExactRectBendPass.c @@ -221,7 +221,7 @@ MODULE_DEF(ExactRectBendPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 26b6d0bb5..0274c4326 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -216,7 +216,7 @@ MODULE_DEF(ExactRectBendRadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactRectangularBendPass.c b/atintegrators/ExactRectangularBendPass.c index 82bc6d6eb..a91fa447b 100644 --- a/atintegrators/ExactRectangularBendPass.c +++ b/atintegrators/ExactRectangularBendPass.c @@ -214,7 +214,7 @@ MODULE_DEF(ExactRectangularBendPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index 99e353c93..c4af47272 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -217,7 +217,7 @@ MODULE_DEF(ExactRectangularBendRadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactSectorBendPass.c b/atintegrators/ExactSectorBendPass.c index 64254cec7..44b910834 100644 --- a/atintegrators/ExactSectorBendPass.c +++ b/atintegrators/ExactSectorBendPass.c @@ -206,7 +206,7 @@ MODULE_DEF(ExactSectorBendPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index df7a1c3f9..02b5eb506 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -200,7 +200,7 @@ MODULE_DEF(ExactSectorBendRadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/GWigSymplecticPass.c b/atintegrators/GWigSymplecticPass.c index 26b5e945c..9d783f56c 100644 --- a/atintegrators/GWigSymplecticPass.c +++ b/atintegrators/GWigSymplecticPass.c @@ -219,7 +219,7 @@ MODULE_DEF(GWigSymplecticPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index 802aab584..ff653c121 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -236,7 +236,7 @@ MODULE_DEF(GWigSymplecticRadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/IdTablePass.c b/atintegrators/IdTablePass.c index 1739e85c9..c697f38a7 100644 --- a/atintegrators/IdTablePass.c +++ b/atintegrators/IdTablePass.c @@ -150,7 +150,7 @@ MODULE_DEF(IdTablePass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/IdentityPass.c b/atintegrators/IdentityPass.c index 1ba8085d4..791da90fe 100644 --- a/atintegrators/IdentityPass.c +++ b/atintegrators/IdentityPass.c @@ -74,7 +74,7 @@ MODULE_DEF(IdentityPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ImpedanceTablePass.c b/atintegrators/ImpedanceTablePass.c index d2bc29b2d..183509c1b 100755 --- a/atintegrators/ImpedanceTablePass.c +++ b/atintegrators/ImpedanceTablePass.c @@ -302,7 +302,7 @@ MODULE_DEF(ImpedanceTablePass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/Matrix66Pass.c b/atintegrators/Matrix66Pass.c index 11e126bfb..8dd7f89ef 100644 --- a/atintegrators/Matrix66Pass.c +++ b/atintegrators/Matrix66Pass.c @@ -60,7 +60,7 @@ MODULE_DEF(Matrix66Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/MatrixTijkPass.c b/atintegrators/MatrixTijkPass.c index 81789fbbb..e7a10aa9a 100644 --- a/atintegrators/MatrixTijkPass.c +++ b/atintegrators/MatrixTijkPass.c @@ -96,7 +96,7 @@ MODULE_DEF(MatrixTijkPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/QuantDiffPass.c b/atintegrators/QuantDiffPass.c index 3920dde9d..0399c3b28 100644 --- a/atintegrators/QuantDiffPass.c +++ b/atintegrators/QuantDiffPass.c @@ -68,7 +68,7 @@ MODULE_DEF(QuantDiffPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double* r_in; const mxArray* ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/SimpleQuantDiffPass.c b/atintegrators/SimpleQuantDiffPass.c index 988fad0be..6887f8036 100644 --- a/atintegrators/SimpleQuantDiffPass.c +++ b/atintegrators/SimpleQuantDiffPass.c @@ -86,7 +86,7 @@ MODULE_DEF(SimpleQuantDiffPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/SimpleRadiationPass.c b/atintegrators/SimpleRadiationPass.c index 2bed3b515..3443aff02 100644 --- a/atintegrators/SimpleRadiationPass.c +++ b/atintegrators/SimpleRadiationPass.c @@ -110,7 +110,7 @@ MODULE_DEF(SimpleRadiationPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/SliceMomentsPass.c b/atintegrators/SliceMomentsPass.c index 8ba718a78..b388ce2de 100644 --- a/atintegrators/SliceMomentsPass.c +++ b/atintegrators/SliceMomentsPass.c @@ -186,7 +186,7 @@ MODULE_DEF(SliceMomentsPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; diff --git a/atintegrators/SolenoidLinearPass.c b/atintegrators/SolenoidLinearPass.c index 5cd1698cb..95fc10547 100644 --- a/atintegrators/SolenoidLinearPass.c +++ b/atintegrators/SolenoidLinearPass.c @@ -100,7 +100,7 @@ MODULE_DEF(SolenoidLinearPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2 ) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/StrMPoleSymplectic4Pass.c b/atintegrators/StrMPoleSymplectic4Pass.c index 6105eee81..9ac88ad4f 100644 --- a/atintegrators/StrMPoleSymplectic4Pass.c +++ b/atintegrators/StrMPoleSymplectic4Pass.c @@ -172,7 +172,7 @@ MODULE_DEF(StrMPoleSymplectic4Pass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/StrMPoleSymplectic4QuantPass.c b/atintegrators/StrMPoleSymplectic4QuantPass.c index 83fa35448..bf8bfcbb0 100644 --- a/atintegrators/StrMPoleSymplectic4QuantPass.c +++ b/atintegrators/StrMPoleSymplectic4QuantPass.c @@ -222,7 +222,7 @@ MODULE_DEF(StrMPoleSymplectic4QuantPass) /* Dummy module initialisation * #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index 88fb00470..25cd330c5 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/atintegrators/StrMPoleSymplectic4RadPass.c @@ -171,7 +171,7 @@ MODULE_DEF(StrMPoleSymplectic4RadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/TestRandomPass.c b/atintegrators/TestRandomPass.c index 858578468..40dea1c2c 100644 --- a/atintegrators/TestRandomPass.c +++ b/atintegrators/TestRandomPass.c @@ -64,7 +64,7 @@ MODULE_DEF(TestRandomPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/ThinMPolePass.c b/atintegrators/ThinMPolePass.c index e83ec304c..af8e7dcfd 100644 --- a/atintegrators/ThinMPolePass.c +++ b/atintegrators/ThinMPolePass.c @@ -129,7 +129,7 @@ MODULE_DEF(ThinMPolePass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/VariableThinMPolePass.c b/atintegrators/VariableThinMPolePass.c index 0c25e9afa..8ba2f88fa 100644 --- a/atintegrators/VariableThinMPolePass.c +++ b/atintegrators/VariableThinMPolePass.c @@ -186,7 +186,7 @@ MODULE_DEF(VariableThinMPolePass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double* r_in; const mxArray* ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/WakeFieldPass.c b/atintegrators/WakeFieldPass.c index 0e705f1af..703f2fac3 100755 --- a/atintegrators/WakeFieldPass.c +++ b/atintegrators/WakeFieldPass.c @@ -165,7 +165,7 @@ MODULE_DEF(WakeFieldPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/WiggLinearPass.c b/atintegrators/WiggLinearPass.c index 20976281d..beeb13659 100644 --- a/atintegrators/WiggLinearPass.c +++ b/atintegrators/WiggLinearPass.c @@ -125,7 +125,7 @@ MODULE_DEF(WiggLinearPass) /* Dummy module initialisation */ #ifdef MATLAB_MEX_FILE void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if (nrhs == 2 ) { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); From 1a10c329bca3782da51d1000defb61ded11bb472 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Tue, 13 Aug 2024 16:39:39 +0200 Subject: [PATCH 02/18] rebqsed on master --- atintegrators/RFCavityPass.c | 10 ++++++---- atintegrators/atelem.c | 11 ++++++++++- atmat/atphysics/LinearOptics/findelemm66.m | 5 ++++- .../atphysics/ParameterSummaryFunctions/atx.m | 3 ++- atmat/atphysics/Radiation/ohmienvelope.m | 18 +++++++----------- 5 files changed, 29 insertions(+), 18 deletions(-) diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index c2ac702f6..bbb39bd50 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -37,11 +37,12 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, { int nturn=Param->nturn; double T0=Param->T0; + double energy; if (!Elem) { double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag; Length=atGetDouble(ElemData,"Length"); check_error(); Voltage=atGetDouble(ElemData,"Voltage"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); @@ -54,7 +55,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->TimeLag=TimeLag; Elem->PhaseLag=PhaseLag; } - RFCavityPass(r_in, Elem->Length, Elem->Voltage/Elem->Energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, + if (Param->energy == 0.0) energy = Elem->Energy; + else energy = Param->energy; + RFCavityPass(r_in, Elem->Length, Elem->Voltage/energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, Elem->PhaseLag, nturn, T0, num_particles); return Elem; } @@ -67,8 +70,7 @@ MODULE_DEF(RFCavityPass) /* Dummy module initialisation */ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { - if(nrhs == 2) - { + if (nrhs >= 2) { double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 6fde04348..7c3587776 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -6,6 +6,7 @@ #define ATELEM_C #include "atcommon.h" +#include "attypes.h" /*----------------------------------------------------*/ /* For the integrator code */ @@ -167,6 +168,14 @@ static double* atGetOptionalDoubleArray(const mxArray *ElemData, const char *fie return atGetOptionalDoubleArraySz(ElemData, fieldname, &msz, &nsz); } +static double atEnergy(double elvalue, struct parameters *prm) +{ + double energy; + if (prm->energy == 0.0) energy = elvalue; + else energy = prm->energy; + return energy; +} + #endif /* MATLAB_MEX_FILE */ /*----------------------------------------------------*/ @@ -180,6 +189,7 @@ typedef PyObject atElem; #define atError(...) return (struct elem *) PyErr_Format(PyExc_ValueError, __VA_ARGS__) #define atWarning(...) if (PyErr_WarnFormat(PyExc_RuntimeWarning, 0, __VA_ARGS__) != 0) return NULL #define atPrintf(...) PySys_WriteStdout(__VA_ARGS__) +#define atEnergy(elvalue,prm) (prm->energy) static int array_imported = 0; @@ -327,7 +337,6 @@ static double *atGetOptionalDoubleArray(const PyObject *element, char *name) #endif /* defined(PYAT) */ #if defined(PYAT) || defined(MATLAB_MEX_FILE) -#include "attypes.h" #ifdef __cplusplus #define C_LINK extern "C" diff --git a/atmat/atphysics/LinearOptics/findelemm66.m b/atmat/atphysics/LinearOptics/findelemm66.m index 7fb7b361f..04d439058 100644 --- a/atmat/atphysics/LinearOptics/findelemm66.m +++ b/atmat/atphysics/LinearOptics/findelemm66.m @@ -14,7 +14,10 @@ [XYStep,varargs]=getoption(varargin,'XYStep'); [R0,varargs]=getoption(varargs,'orbit',zeros(6,1)); +[energy,varargs]=getoption(varargs,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',atparticle('relativistic')); [MethodName,R0]=getargs(varargs,ELEM.PassMethod,R0); +particle=atparticle.loadobj(particle); % Convert class object into struct (for C access) % Build a diagonal matrix of initial conditions %scaling=2*XYStep*[1 0.1 1 0.1 1 1]; @@ -23,6 +26,6 @@ % Add to the orbit_in RIN = R0 + [D6, -D6]; % Propagate through the element -ROUT = feval(MethodName,ELEM,RIN); +ROUT = feval(MethodName,ELEM,RIN,energy,particle); % Calculate numerical derivative M66 = (ROUT(:,1:6)-ROUT(:,7:12))./scaling; diff --git a/atmat/atphysics/ParameterSummaryFunctions/atx.m b/atmat/atphysics/ParameterSummaryFunctions/atx.m index 1278ce643..d460c6f1a 100644 --- a/atmat/atphysics/ParameterSummaryFunctions/atx.m +++ b/atmat/atphysics/ParameterSummaryFunctions/atx.m @@ -192,7 +192,7 @@ if has_cavity try - [envelope,espread,blength,m,T]=ohmienvelope(ron,radindex,refpts); + [envelope,espread,blength,m,T]=ohmienvelope(ron,radindex,refpts,energy); [tns,chi]=atdampingrates(m); fs=abs(tns(3))*cell_revfreq; alpha=chi*cell_revfreq; @@ -201,6 +201,7 @@ if any(radindex) jmt=jmat(3); lindata=cellfun(@process,{envelope.R},reshape(num2cell(T,[1 2]),1,[]),num2cell(lindata)); + disp('step3'); else lindata=arrayfun(@deflt,lindata); end diff --git a/atmat/atphysics/Radiation/ohmienvelope.m b/atmat/atphysics/Radiation/ohmienvelope.m index fb72a4245..af92cc3b0 100644 --- a/atmat/atphysics/Radiation/ohmienvelope.m +++ b/atmat/atphysics/Radiation/ohmienvelope.m @@ -1,4 +1,4 @@ -function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts) +function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts,energy) %OHMIENVELOPE calculates equilibrium beam envelope in a % circular accelerator using Ohmi's beam envelope formalism [1]. % [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994) @@ -33,12 +33,10 @@ check_6d(ring,true,'strict',0); NumElements = length(ring); +if nargin<4, energy=atGetRingProperties(ring, 'Energy'); end if nargin<3, refpts=1; end - -% Erase wigglers from the radiative element list. -% Diffusion matrix to be computed with separate FDW function. -Wig=atgetcells(ring,'Bmax'); -radindex = radindex & ~Wig; +ringrad=ring(radindex); +energy=num2cell(energy(ones(size(ringrad)))); [mring, ms, orbit] = findm66(ring,1:NumElements+1); mt=squeeze(num2cell(ms,[1 2])); @@ -49,9 +47,7 @@ % calculate Radiation-Diffusion matrix B for elements with radiation B(radindex)=cellfun(@findmpoleraddiffmatrix,... - ring(radindex),orb(radindex),'UniformOutput',false); -B(Wig)=cellfun(@FDW,... - ring(Wig),orb(Wig),'UniformOutput',false); + ringrad,orb(radindex),energy,'UniformOutput',false); % Calculate cumulative Radiation-Diffusion matrix for the ring BCUM = zeros(6,6); @@ -63,7 +59,7 @@ % Equation for the moment matrix R is % R = MRING*R*MRING' + BCUM; % We rewrite it in the form of Sylvester-Lyapunov equation -% to use MATLAB's SYLVERTER function: +% to use MATLAB's SYLVESTER function: % AA*R + R*BB = CC % where % AA = inv(MRING) @@ -77,7 +73,7 @@ R = sylvester(AA,BB,CC); % Envelope matrix at the ring entrance rmsdp = sqrt(R(5,5)); % R.M.S. energy spread -rmsbl = sqrt(R(6,6)); % R.M.S. bunch length +rmsbl = sqrt(R(6,6)); % R.M.S. bunch lenght [rr,tt,ss]=cellfun(@propag,mt(refpts),Batbeg(refpts),'UniformOutput',false); envelope=struct('R',rr,'Sigma',ss,'Tilt',tt); From 5c4f56721cd01dec4b6bd14873f50600c81acfa0 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Tue, 13 Aug 2024 19:22:36 +0200 Subject: [PATCH 03/18] energy field no more necessary in Matlab --- atintegrators/BndMPoleSymplectic4RadPass.c | 9 +++-- atintegrators/RFCavityPass.c | 15 ++++++--- atintegrators/StrMPoleSymplectic4RadPass.c | 9 +++-- atintegrators/atelem.c | 33 +++++++++++++++++++ atmat/atphysics/LinearOptics/findelemm66.m | 13 +++++--- .../atphysics/ParameterSummaryFunctions/atx.m | 1 - atmat/atphysics/Radiation/ohmienvelope.m | 12 ++++--- 7 files changed, 73 insertions(+), 19 deletions(-) diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index 13e43767c..f51896364 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -213,8 +213,11 @@ MODULE_DEF(BndMPoleSymplectic4RadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double Energy = 0.0; + double rest_energy = 0.0; + double charge = -1.0; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, - FringeInt1, FringeInt2, Energy; + FringeInt1, FringeInt2; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; @@ -232,8 +235,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Energy); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -252,10 +255,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); irho = BendingAngle/Length; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + BndMPoleSymplectic4RadPass(r_in, Length, irho, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index bbb39bd50..486c7d144 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -37,12 +37,12 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, { int nturn=Param->nturn; double T0=Param->T0; - double energy; + double energy = Param->energy; if (!Elem) { double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag; Length=atGetDouble(ElemData,"Length"); check_error(); Voltage=atGetDouble(ElemData,"Voltage"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); + Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); @@ -55,8 +55,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->TimeLag=TimeLag; Elem->PhaseLag=PhaseLag; } - if (Param->energy == 0.0) energy = Elem->Energy; - else energy = Param->energy; + if (energy == 0.0) energy = Elem->Energy; RFCavityPass(r_in, Elem->Length, Elem->Voltage/energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, Elem->PhaseLag, nturn, T0, num_particles); return Elem; @@ -71,21 +70,27 @@ MODULE_DEF(RFCavityPass) /* Dummy module initialisation */ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double Energy = 0.0; + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); double Length=atGetDouble(ElemData,"Length"); double Voltage=atGetDouble(ElemData,"Voltage"); - double Energy=atGetDouble(ElemData,"Energy"); + Energy=atGetOptionalDouble(ElemData,"Energy",Energy); double Frequency=atGetDouble(ElemData,"Frequency"); double TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); double PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); double T0=1.0/Frequency; /* Does not matter since nturns == 0 */ double HarmNumber=round(Frequency*T0); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + RFCavityPass(r_in, Length, Voltage/Energy, Frequency, HarmNumber, TimeLag, PhaseLag, 0, T0, num_particles); } diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index 25cd330c5..656b61a9e 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/atintegrators/StrMPoleSymplectic4RadPass.c @@ -172,10 +172,13 @@ MODULE_DEF(StrMPoleSymplectic4RadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double Energy = 0.0; + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); - double Length, Energy, Scaling; + double Length, Scaling; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; if (mxGetM(prhs[1]) != 6) mexErrMsgTxt("Second argument must be a 6 x N matrix"); @@ -185,8 +188,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Energy); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -199,10 +202,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + StrMPoleSymplectic4RadPass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit, diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 7c3587776..74d4721a1 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -176,6 +176,39 @@ static double atEnergy(double elvalue, struct parameters *prm) return energy; } +static double atGetOptionalDoubleProp(const mxArray *obj, const char *fieldname, double default_value) +{ + mxArray *field=mxGetProperty(obj, 0, fieldname); + return (field) ? mxGetScalar(field) : default_value; +} + +static void atParticle(const mxArray *opts, double *rest_energy, double *charge) +{ + const mxArray *part = mxGetField(opts, 0, "Particle"); + if (part) { + if (mxIsClass(part, "atparticle")) { /* OK */ + *rest_energy = atGetOptionalDoubleProp(part, "rest_energy", 0.0); + *charge = atGetOptionalDoubleProp(part, "charge", -1.0); + } + else { /* particle is not a Particle object */ + mexErrMsgIdAndTxt("Atpass:WrongParameter","Particle must be an 'atparticle' object"); + } + } +} + +static void atProperties(const mxArray *opts, double *energy, double *rest_energy, double *charge) +{ + mxArray *field; + if (!mxIsStruct(opts)) { + mexErrMsgIdAndTxt("Atpass:WrongParameter","ring properties must be a struct"); + } + field = mxGetField(opts, 0, "Energy"); + if (field) { + *energy = mxGetScalar(field); + atParticle(opts, rest_energy, charge); + } +} + #endif /* MATLAB_MEX_FILE */ /*----------------------------------------------------*/ diff --git a/atmat/atphysics/LinearOptics/findelemm66.m b/atmat/atphysics/LinearOptics/findelemm66.m index 04d439058..2f2dffd1c 100644 --- a/atmat/atphysics/LinearOptics/findelemm66.m +++ b/atmat/atphysics/LinearOptics/findelemm66.m @@ -9,15 +9,20 @@ % M66=FINDELEMM66(...,'orbit',ORBITIN) % ORBITIN - 6-by-1 phase space coordinates at the entrance % (default: zeros(6,1)) +% +% M66=FINDELEMM66(...,'Energy',ENERGY) +% Use ENERGY and ignore the 'Energy' field of elements +% +% M66=FINDELEMM66(...,'Particle',PARTICLE) +% Use PARTICLE (default is relativistic) % % See also FINDELEMM44 [XYStep,varargs]=getoption(varargin,'XYStep'); [R0,varargs]=getoption(varargs,'orbit',zeros(6,1)); -[energy,varargs]=getoption(varargs,'Energy',0.0); -[particle,varargs]=getoption(varargs,'Particle',atparticle('relativistic')); +[props.Energy,varargs]=getoption(varargs,'Energy',0.0); +[props.Particle,varargs]=getoption(varargs,'Particle',atparticle('relativistic')); [MethodName,R0]=getargs(varargs,ELEM.PassMethod,R0); -particle=atparticle.loadobj(particle); % Convert class object into struct (for C access) % Build a diagonal matrix of initial conditions %scaling=2*XYStep*[1 0.1 1 0.1 1 1]; @@ -26,6 +31,6 @@ % Add to the orbit_in RIN = R0 + [D6, -D6]; % Propagate through the element -ROUT = feval(MethodName,ELEM,RIN,energy,particle); +ROUT = feval(MethodName,ELEM,RIN,props); % Calculate numerical derivative M66 = (ROUT(:,1:6)-ROUT(:,7:12))./scaling; diff --git a/atmat/atphysics/ParameterSummaryFunctions/atx.m b/atmat/atphysics/ParameterSummaryFunctions/atx.m index d460c6f1a..d6f2eb968 100644 --- a/atmat/atphysics/ParameterSummaryFunctions/atx.m +++ b/atmat/atphysics/ParameterSummaryFunctions/atx.m @@ -201,7 +201,6 @@ if any(radindex) jmt=jmat(3); lindata=cellfun(@process,{envelope.R},reshape(num2cell(T,[1 2]),1,[]),num2cell(lindata)); - disp('step3'); else lindata=arrayfun(@deflt,lindata); end diff --git a/atmat/atphysics/Radiation/ohmienvelope.m b/atmat/atphysics/Radiation/ohmienvelope.m index af92cc3b0..02bd7e74a 100644 --- a/atmat/atphysics/Radiation/ohmienvelope.m +++ b/atmat/atphysics/Radiation/ohmienvelope.m @@ -35,8 +35,6 @@ NumElements = length(ring); if nargin<4, energy=atGetRingProperties(ring, 'Energy'); end if nargin<3, refpts=1; end -ringrad=ring(radindex); -energy=num2cell(energy(ones(size(ringrad)))); [mring, ms, orbit] = findm66(ring,1:NumElements+1); mt=squeeze(num2cell(ms,[1 2])); @@ -46,8 +44,8 @@ B=zr(ones(NumElements,1)); % B{i} is the diffusion matrix of the i-th element % calculate Radiation-Diffusion matrix B for elements with radiation -B(radindex)=cellfun(@findmpoleraddiffmatrix,... - ringrad,orb(radindex),energy,'UniformOutput',false); +B(radindex)=cellfun(@diffmatrix,... + ring(radindex),orb(radindex),'UniformOutput',false); % Calculate cumulative Radiation-Diffusion matrix for the ring BCUM = zeros(6,6); @@ -84,10 +82,14 @@ if nout>=2, varargout{2}=ms(:,:,refpts); end if nout>=1, varargout{1}=mring; end + function diff=diffmatrix(elem,orbit) + diff=findmpoleraddiffmatrix(elem, orbit, energy); + end + function btx=cumulb(elem,orbit,b) % Calculate 6-by-6 linear transfer matrix in each element % near the equilibrium orbit - m=findelemm66(elem,elem.PassMethod,orbit); + m=findelemm66(elem,elem.PassMethod,'orbit',orbit,'Energy',energy); % Cumulative diffusion matrix of the entire ring BCUM = m*BCUM*m' + b; btx=BCUM; From 561d3228e6a3c3583df4eff437896b2a32d6d71c Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 14 Aug 2024 17:04:37 +0200 Subject: [PATCH 04/18] energy field no more necessary in Matlab --- atintegrators/BndMPoleSymplectic4RadPass.c | 13 +++--- atintegrators/RFCavityPass.c | 1 + atintegrators/StrMPoleSymplectic4RadPass.c | 12 ++--- atintegrators/atelem.c | 45 ++----------------- atintegrators/ringproperties.c | 34 ++++++++++++++ atmat/atphysics/LinearOptics/findelemm66.m | 6 +-- .../Radiation/findmpoleraddiffmatrix.c | 2 +- atmat/attrack/atpass.c | 40 ++--------------- atmat/lattice/at2str.m | 2 +- 9 files changed, 60 insertions(+), 95 deletions(-) create mode 100644 atintegrators/ringproperties.c diff --git a/atintegrators/BndMPoleSymplectic4RadPass.c b/atintegrators/BndMPoleSymplectic4RadPass.c index f51896364..776741648 100644 --- a/atintegrators/BndMPoleSymplectic4RadPass.c +++ b/atintegrators/BndMPoleSymplectic4RadPass.c @@ -128,7 +128,7 @@ void BndMPoleSymplectic4RadPass(double *r, double le, double irho, double *A, do ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { - double irho; + double irho, energy; if (!Elem) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; @@ -143,8 +143,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -193,6 +193,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + BndMPoleSymplectic4RadPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, @@ -201,7 +203,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->fringeIntM0, Elem->fringeIntP0, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -213,11 +215,10 @@ MODULE_DEF(BndMPoleSymplectic4RadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { - double Energy = 0.0; double rest_energy = 0.0; double charge = -1.0; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, - FringeInt1, FringeInt2; + FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; @@ -236,7 +237,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - Energy=atGetOptionalDouble(ElemData,"Energy",Energy); check_error(); + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index 486c7d144..188eb22ec 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -56,6 +56,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->PhaseLag=PhaseLag; } if (energy == 0.0) energy = Elem->Energy; + RFCavityPass(r_in, Elem->Length, Elem->Voltage/energy, Elem->Frequency, Elem->HarmNumber, Elem->TimeLag, Elem->PhaseLag, nturn, T0, num_particles); return Elem; diff --git a/atintegrators/StrMPoleSymplectic4RadPass.c b/atintegrators/StrMPoleSymplectic4RadPass.c index 656b61a9e..6b9971355 100644 --- a/atintegrators/StrMPoleSymplectic4RadPass.c +++ b/atintegrators/StrMPoleSymplectic4RadPass.c @@ -109,6 +109,7 @@ void StrMPoleSymplectic4RadPass(double *r, double le, double *A, double *B, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double Length, Energy, Scaling; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; @@ -118,8 +119,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -154,13 +155,15 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + StrMPoleSymplectic4RadPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->FringeQuadEntrance, Elem->FringeQuadExit, Elem->fringeIntM0, Elem->fringeIntP0, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -172,13 +175,12 @@ MODULE_DEF(StrMPoleSymplectic4RadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { - double Energy = 0.0; double rest_energy = 0.0; double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); - double Length, Scaling; + double Length, Scaling, Energy; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; if (mxGetM(prhs[1]) != 6) mexErrMsgTxt("Second argument must be a 6 x N matrix"); @@ -189,7 +191,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); /*optional fields*/ - Energy=atGetOptionalDouble(ElemData,"Energy",Energy); check_error(); + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 74d4721a1..9ef2de4c6 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -71,6 +71,8 @@ typedef mxArray atElem; #define atError(...) mexErrMsgIdAndTxt("AT:PassError", __VA_ARGS__) #define atWarning(...) mexWarnMsgIdAndTxt("AT:PassWarning", __VA_ARGS__) #define atPrintf(...) mexPrintf(__VA_ARGS__) +#define atEnergy(ringenergy,elemenergy) ((ringenergy)==0.0 ? (elemenergy) : (ringenergy)) +#include "ringproperties.c" static mxArray *get_field(const mxArray *pm, const char *fieldname) { @@ -168,47 +170,6 @@ static double* atGetOptionalDoubleArray(const mxArray *ElemData, const char *fie return atGetOptionalDoubleArraySz(ElemData, fieldname, &msz, &nsz); } -static double atEnergy(double elvalue, struct parameters *prm) -{ - double energy; - if (prm->energy == 0.0) energy = elvalue; - else energy = prm->energy; - return energy; -} - -static double atGetOptionalDoubleProp(const mxArray *obj, const char *fieldname, double default_value) -{ - mxArray *field=mxGetProperty(obj, 0, fieldname); - return (field) ? mxGetScalar(field) : default_value; -} - -static void atParticle(const mxArray *opts, double *rest_energy, double *charge) -{ - const mxArray *part = mxGetField(opts, 0, "Particle"); - if (part) { - if (mxIsClass(part, "atparticle")) { /* OK */ - *rest_energy = atGetOptionalDoubleProp(part, "rest_energy", 0.0); - *charge = atGetOptionalDoubleProp(part, "charge", -1.0); - } - else { /* particle is not a Particle object */ - mexErrMsgIdAndTxt("Atpass:WrongParameter","Particle must be an 'atparticle' object"); - } - } -} - -static void atProperties(const mxArray *opts, double *energy, double *rest_energy, double *charge) -{ - mxArray *field; - if (!mxIsStruct(opts)) { - mexErrMsgIdAndTxt("Atpass:WrongParameter","ring properties must be a struct"); - } - field = mxGetField(opts, 0, "Energy"); - if (field) { - *energy = mxGetScalar(field); - atParticle(opts, rest_energy, charge); - } -} - #endif /* MATLAB_MEX_FILE */ /*----------------------------------------------------*/ @@ -222,7 +183,7 @@ typedef PyObject atElem; #define atError(...) return (struct elem *) PyErr_Format(PyExc_ValueError, __VA_ARGS__) #define atWarning(...) if (PyErr_WarnFormat(PyExc_RuntimeWarning, 0, __VA_ARGS__) != 0) return NULL #define atPrintf(...) PySys_WriteStdout(__VA_ARGS__) -#define atEnergy(elvalue,prm) (prm->energy) +#define atEnergy(ringenergy,elemenergy) (ringenergy) static int array_imported = 0; diff --git a/atintegrators/ringproperties.c b/atintegrators/ringproperties.c new file mode 100644 index 000000000..baabef3db --- /dev/null +++ b/atintegrators/ringproperties.c @@ -0,0 +1,34 @@ + +static double atGetOptionalDoubleProp(const mxArray *obj, const char *fieldname, double default_value) +{ + mxArray *field=mxGetProperty(obj, 0, fieldname); + return (field) ? mxGetScalar(field) : default_value; +} + +static void atParticle(const mxArray *opts, double *rest_energy, double *charge) +{ + const mxArray *part = mxGetField(opts, 0, "Particle"); + if (part) { + if (mxIsClass(part, "atparticle")) { /* OK */ + *rest_energy = atGetOptionalDoubleProp(part, "rest_energy", 0.0); + *charge = atGetOptionalDoubleProp(part, "charge", -1.0); + } + else { /* particle is not a Particle object */ + mexErrMsgIdAndTxt("Atpass:WrongParameter","Particle must be an 'atparticle' object"); + } + } +} + +static void atProperties(const mxArray *opts, double *energy, double *rest_energy, double *charge) +{ + mxArray *field; + if (!mxIsStruct(opts)) { + mexErrMsgIdAndTxt("Atpass:WrongParameter","ring properties must be a struct"); + } + field = mxGetField(opts, 0, "Energy"); + if (field) { + double ener = mxGetScalar(field); + if (ener != 0.0) *energy = ener; + atParticle(opts, rest_energy, charge); + } +} diff --git a/atmat/atphysics/LinearOptics/findelemm66.m b/atmat/atphysics/LinearOptics/findelemm66.m index 2f2dffd1c..4cbf32e58 100644 --- a/atmat/atphysics/LinearOptics/findelemm66.m +++ b/atmat/atphysics/LinearOptics/findelemm66.m @@ -20,8 +20,8 @@ [XYStep,varargs]=getoption(varargin,'XYStep'); [R0,varargs]=getoption(varargs,'orbit',zeros(6,1)); -[props.Energy,varargs]=getoption(varargs,'Energy',0.0); -[props.Particle,varargs]=getoption(varargs,'Particle',atparticle('relativistic')); +[energy,varargs]=getoption(varargs,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',[]); [MethodName,R0]=getargs(varargs,ELEM.PassMethod,R0); % Build a diagonal matrix of initial conditions @@ -31,6 +31,6 @@ % Add to the orbit_in RIN = R0 + [D6, -D6]; % Propagate through the element -ROUT = feval(MethodName,ELEM,RIN,props); +ROUT=elempass(ELEM,RIN,'PassMethod',MethodName,'Energy',energy,'Particle',particle); % Calculate numerical derivative M66 = (ROUT(:,1:6)-ROUT(:,7:12))./scaling; diff --git a/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c b/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c index d3f203d35..92e69de50 100644 --- a/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c +++ b/atmat/atphysics/Radiation/findmpoleraddiffmatrix.c @@ -522,7 +522,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) plhs[0] = mxCreateDoubleMatrix(6,6,mxREAL); bdiff = mxGetDoubles(plhs[0]); for (i=0; i<36; i++) bdiff[i]=0.0; - + diffmatrix(mxElem, orb, energy, bdiff); } #endif /*MATLAB_MEX_FILE*/ diff --git a/atmat/attrack/atpass.c b/atmat/attrack/atpass.c index c46b5362c..6057991a5 100644 --- a/atmat/attrack/atpass.c +++ b/atmat/attrack/atpass.c @@ -11,6 +11,7 @@ #include "attypes.h" #include "elempass.h" #include "atrandom.c" +#include "ringproperties.c" /* Get ready for R2018a C matrix API */ #ifndef mxGetDoubles @@ -235,42 +236,7 @@ static mxDouble *passhook(mxArray *mxPassArg[], mxArray *mxElem, mxArray *func) return tempdoubleptr; } -static double getoptionaldoubleprop(const mxArray *obj, const char *fieldname, double default_value) -{ - mxArray *field=mxGetProperty(obj, 0, fieldname); - return (field) ? mxGetScalar(field) : default_value; -} - -static void getparticle(const mxArray *opts, double *rest_energy, double *charge) -{ - const mxArray *part = mxGetField(opts, 0, "Particle"); - if (part) { - if (mxIsClass(part, "atparticle")) { /* OK */ - *rest_energy = getoptionaldoubleprop(part, "rest_energy", 0.0); - *charge = getoptionaldoubleprop(part, "charge", -1.0); - } - else { /* particle is not a Particle object */ - mexErrMsgIdAndTxt("Atpass:WrongParameter","Particle must be an 'atparticle' object"); - } - } -} - -static void getproperties(const mxArray *opts, double *energy, double *rest_energy, double *charge) -{ - mxArray *field; - if (!mxIsStruct(opts)) { - mexErrMsgIdAndTxt("Atpass:WrongParameter","ring properties must be a struct"); - } - field = mxGetField(opts, 0, "Energy"); - if (field) { - *energy = mxGetScalar(field); - getparticle(opts, rest_energy, charge); - } -} - -/*! - * - +/* @param[in] [0] LATTICE @param[in,out] [1] INITCONDITIONS @param[in] [2] NEWLATTICE @@ -345,7 +311,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) param.nturn = counter; if (nrhs >= 10) { - getproperties(RINGPROPERTIES, ¶m.energy, ¶m.rest_energy, ¶m.charge); + atProperties(RINGPROPERTIES, ¶m.energy, ¶m.rest_energy, ¶m.charge); } if (nlhs >= 2) { diff --git a/atmat/lattice/at2str.m b/atmat/lattice/at2str.m index 631ca57ab..74e0201a2 100644 --- a/atmat/lattice/at2str.m +++ b/atmat/lattice/at2str.m @@ -79,7 +79,7 @@ [options,args]=doptions(elem,create,{'PolynomA','PolynomB'}); case 'RFCavity' create=@atrfcavity; - [options,args]=doptions(elem,create,{'Length','Voltage','Frequency','HarmNumber','Energy'}); + [options,args]=doptions(elem,create,{'Length','Voltage','Frequency','HarmNumber'}); case 'RingParam' create=@atringparam; [options,args]=doptions(elem,create,{'Energy','Periodicity'}); From f53bbc656141ef042945a548100e12a17435d78b Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 14 Aug 2024 18:01:29 +0200 Subject: [PATCH 05/18] rebased on master --- atintegrators/BndMPoleSymplectic4E2RadPass.c | 13 ++++++--- atintegrators/BndMPoleSymplectic4QuantPass.c | 14 ++++++--- atintegrators/ExactMultipoleRadPass.c | 11 +++++-- atintegrators/ExactRectBendRadPass.c | 11 +++++-- atintegrators/ExactRectangularBendRadPass.c | 10 +++++-- atintegrators/ExactSectorBendRadPass.c | 11 +++++-- atintegrators/GWigSymplecticRadPass.c | 30 +++++++++----------- atintegrators/StrMPoleSymplectic4QuantPass.c | 13 +++++++-- 8 files changed, 78 insertions(+), 35 deletions(-) diff --git a/atintegrators/BndMPoleSymplectic4E2RadPass.c b/atintegrators/BndMPoleSymplectic4E2RadPass.c index 076a32472..2c5ae5ec4 100644 --- a/atintegrators/BndMPoleSymplectic4E2RadPass.c +++ b/atintegrators/BndMPoleSymplectic4E2RadPass.c @@ -202,7 +202,7 @@ void BndMPoleSymplectic4E2RadPass(double *r, double le, double irho, double *A, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { - double irho; + double irho, energy; if (!Elem) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, Energy, FringeInt1, FringeInt2; @@ -216,8 +216,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1",0); check_error(); @@ -258,13 +258,15 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + BndMPoleSymplectic4E2RadPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeInt1, Elem->FringeInt2, Elem->FullGap, Elem->h1, Elem->h2, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -276,6 +278,8 @@ MODULE_DEF(BndMPoleSymplectic4E2RadPass) /* Dummy module initialisation * void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double irho; double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; int MaxOrder, NumIntSteps; @@ -293,8 +297,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap", 0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeInt1=atGetOptionalDouble(ElemData,"FringeInt1", 0); check_error(); @@ -309,6 +313,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); irho = BendingAngle/Length; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); diff --git a/atintegrators/BndMPoleSymplectic4QuantPass.c b/atintegrators/BndMPoleSymplectic4QuantPass.c index 107649e6f..b4bb1304f 100644 --- a/atintegrators/BndMPoleSymplectic4QuantPass.c +++ b/atintegrators/BndMPoleSymplectic4QuantPass.c @@ -180,7 +180,7 @@ void BndMPoleSymplectic4QuantPass(double *r, double le, double irho, double *A, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { - double irho; + double irho, energy; if (!Elem) { double Length, BendingAngle, EntranceAngle, ExitAngle, FullGap, Scaling, FringeInt1, FringeInt2, Energy; @@ -195,8 +195,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -245,6 +245,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->KickAngle=KickAngle; } irho = Elem->BendingAngle/Elem->Length; + energy = atEnergy(Param->energy, Elem->Energy); + BndMPoleSymplectic4QuantPass(r_in, Elem->Length, irho, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, @@ -253,7 +255,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->fringeIntM0, Elem->fringeIntP0, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, + Elem->KickAngle, Elem->Scaling, energy, Param->thread_rng, num_particles); return Elem; } @@ -271,6 +273,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) int MaxOrder, NumIntSteps, FringeBendEntrance, FringeBendExit, FringeQuadEntrance, FringeQuadExit; double *PolynomA, *PolynomB, *R1, *R2, *T1, *T2, *EApertures, *RApertures, *fringeIntM0, *fringeIntP0, *KickAngle; + double rest_energy = 0.0; + double charge = -1.0; double irho; double *r_in; const mxArray *ElemData = prhs[0]; @@ -285,8 +289,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) BendingAngle=atGetDouble(ElemData,"BendingAngle"); check_error(); EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); FullGap=atGetOptionalDouble(ElemData,"FullGap",0); check_error(); @@ -305,10 +309,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); irho = BendingAngle/Length; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + BndMPoleSymplectic4QuantPass(r_in, Length, irho, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/ExactMultipoleRadPass.c b/atintegrators/ExactMultipoleRadPass.c index 0155b32af..ce3094005 100644 --- a/atintegrators/ExactMultipoleRadPass.c +++ b/atintegrators/ExactMultipoleRadPass.c @@ -106,6 +106,7 @@ static void multipole_pass( ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double Length = atGetDouble(ElemData, "Length"); check_error(); double *PolynomA = atGetDoubleArray(ElemData, "PolynomA"); check_error(); @@ -148,12 +149,14 @@ ExportMode struct elem *trackFunction(const atElem *ElemData, struct elem *Elem, Elem->RApertures = RApertures; Elem->KickAngle = KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + multipole_pass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->FringeQuadEntrance, Elem->FringeQuadExit, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Energy, Elem->Scaling, num_particles); + Elem->KickAngle, energy, Elem->Scaling, num_particles); return Elem; } @@ -164,6 +167,8 @@ MODULE_DEF(ExactMultipoleRadPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -174,7 +179,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { int MaxOrder = atGetLong(ElemData, "MaxOrder"); check_error(); int NumIntSteps = atGetLong(ElemData, "NumIntSteps"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); check_error(); + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); int FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); int FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -189,6 +194,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (NumIntSteps <= 0) { atError("NumIntSteps must be positive"); check_error(); } + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); diff --git a/atintegrators/ExactRectBendRadPass.c b/atintegrators/ExactRectBendRadPass.c index 0274c4326..3272c1665 100644 --- a/atintegrators/ExactRectBendRadPass.c +++ b/atintegrators/ExactRectBendRadPass.c @@ -139,6 +139,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double Length=atGetDouble(ElemData,"Length"); check_error(); double *PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); @@ -197,6 +198,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, @@ -205,7 +208,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->gK,Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -217,6 +220,8 @@ MODULE_DEF(ExactRectBendRadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -231,7 +236,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); check_error(); + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); int FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); @@ -251,10 +256,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (NumIntSteps == 0) { atError("NumIntSteps == 0 not allowed with radiation"); check_error(); } + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + ExactRectangularBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/ExactRectangularBendRadPass.c b/atintegrators/ExactRectangularBendRadPass.c index c4af47272..f40988877 100644 --- a/atintegrators/ExactRectangularBendRadPass.c +++ b/atintegrators/ExactRectangularBendRadPass.c @@ -141,6 +141,7 @@ static void ExactRectangularBendRad(double *r, double le, double bending_angle, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double Length=atGetDouble(ElemData,"Length"); check_error(); double *PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); @@ -199,6 +200,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + ExactRectangularBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, Elem->FringeBendEntrance,Elem->FringeBendExit, @@ -206,7 +209,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->gK,Elem->x0ref,Elem->refdz, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -218,6 +221,8 @@ MODULE_DEF(ExactRectangularBendRadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -232,7 +237,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); check_error(); + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); int FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); @@ -252,6 +257,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (NumIntSteps <= 0) { atError("NumIntSteps must be positive"); check_error(); } + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); diff --git a/atintegrators/ExactSectorBendRadPass.c b/atintegrators/ExactSectorBendRadPass.c index 02b5eb506..42469e4f6 100644 --- a/atintegrators/ExactSectorBendRadPass.c +++ b/atintegrators/ExactSectorBendRadPass.c @@ -127,6 +127,7 @@ static void ExactSectorBendRad(double *r, double le, double bending_angle, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double Length=atGetDouble(ElemData,"Length"); check_error(); double *PolynomA=atGetDoubleArray(ElemData,"PolynomA"); check_error(); @@ -181,6 +182,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + ExactSectorBendRad(r_in, Elem->Length, Elem->BendingAngle, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->EntranceAngle, Elem->ExitAngle, @@ -189,7 +192,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->gK, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, - Elem->KickAngle, Elem->Scaling, Elem->Energy, num_particles); + Elem->KickAngle, Elem->Scaling, energy, num_particles); return Elem; } @@ -201,6 +204,8 @@ MODULE_DEF(ExactSectorBendRadPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -215,7 +220,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double EntranceAngle=atGetDouble(ElemData,"EntranceAngle"); check_error(); double ExitAngle=atGetDouble(ElemData,"ExitAngle"); check_error(); /*optional fields*/ - double Energy=atGetDouble(ElemData,"Energy"); check_error(); + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); double Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); int FringeBendEntrance=atGetOptionalLong(ElemData,"FringeBendEntrance",1); check_error(); int FringeBendExit=atGetOptionalLong(ElemData,"FringeBendExit",1); check_error(); @@ -233,10 +238,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) if (NumIntSteps == 0) { atError("NumIntSteps == 0 not allowed with radiation"); check_error(); } + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + ExactSectorBendRad(r_in, Length, BendingAngle, PolynomA, PolynomB, MaxOrder, NumIntSteps, EntranceAngle, ExitAngle, FringeBendEntrance, FringeBendExit, diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index ff653c121..6785380bc 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -6,10 +6,6 @@ *--------------------------------------------------------------------------- * Modification Log: * ----------------- - * .03 2024-05-06 J. Arenillas, ALBA, jarenillas@axt.email - * Adding rotations and translations to wiggler. - * Bug fix in wiggler initialisation. - * Energy parameter bug fix. * .02 2003-06-18 J. Li * Cleanup the code * @@ -127,7 +123,7 @@ void GWigInit(struct gwigR *Wig,double design_energy, double Ltot, double Lw, #define second 2 #define fourth 4 -void GWigSymplecticRadPass(double *r, double Energy, double Ltot, double Lw, +void GWigSymplecticRadPass(double *r,double Energy, double Ltot, double Lw, double Bmax, int Nstep, int Nmeth, int NHharm, int NVharm, double *By, double *Bx, double *T1, double *T2, double *R1, double *R2, int num_particles) @@ -145,13 +141,11 @@ void GWigSymplecticRadPass(double *r, double Energy, double Ltot, double Lw, zEndPointV[0] = 0; zEndPointV[1] = Ltot; + GWigInit(&Wig, Energy, Ltot, Lw, Bmax, Nstep, Nmeth, NHharm, NVharm,0, 0, zEndPointH, zEndPointV, By, Bx, T1, T2, R1, R2); + for(c = 0;cenergy); check_error(); Ltot = atGetDouble(ElemData, "Length"); check_error(); Lw = atGetDouble(ElemData, "Lw"); check_error(); Bmax = atGetDouble(ElemData, "Bmax"); check_error(); @@ -196,6 +187,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, By = atGetDoubleArray(ElemData, "By"); check_error(); Bx = atGetDoubleArray(ElemData, "Bx"); check_error(); /* Optional fields */ + Energy = atGetOptionalDouble(ElemData, "Energy",Param->energy); check_error(); R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); @@ -218,7 +210,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T1=T1; Elem->T2=T2; } - GWigSymplecticRadPass(r_in, Elem->Energy, Elem->Length, Elem->Lw, + energy = atEnergy(Param->energy, Elem->Energy); + + GWigSymplecticRadPass(r_in, energy, Elem->Length, Elem->Lw, Elem->Bmax, Elem->Nstep, Elem->Nmeth, Elem->NHharm, Elem->NVharm, Elem->By, Elem->Bx, Elem->T1, Elem->T2, Elem->R1, Elem->R2, num_particles); @@ -237,6 +231,8 @@ MODULE_DEF(GWigSymplecticRadPass) /* Dummy module initialisation */ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -245,8 +241,8 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs double Ltot, Lw, Bmax, Energy; int Nstep, Nmeth; int NHharm, NVharm; + if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); - Energy = atGetDouble(ElemData, "Energy"); check_error(); Ltot = atGetDouble(ElemData, "Length"); check_error(); Lw = atGetDouble(ElemData, "Lw"); check_error(); Bmax = atGetDouble(ElemData, "Bmax"); check_error(); @@ -257,11 +253,13 @@ void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs By = atGetDoubleArray(ElemData, "By"); check_error(); Bx = atGetDoubleArray(ElemData, "Bx"); check_error(); /* Optional fields */ + Energy = atGetOptionalDouble(ElemData, "Energy",0.0); check_error(); R1 = atGetOptionalDoubleArray(ElemData, "R1"); check_error(); R2 = atGetOptionalDoubleArray(ElemData, "R2"); check_error(); T1 = atGetOptionalDoubleArray(ElemData, "T1"); check_error(); T2 = atGetOptionalDoubleArray(ElemData, "T2"); check_error(); - if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); diff --git a/atintegrators/StrMPoleSymplectic4QuantPass.c b/atintegrators/StrMPoleSymplectic4QuantPass.c index bf8bfcbb0..e641f151f 100644 --- a/atintegrators/StrMPoleSymplectic4QuantPass.c +++ b/atintegrators/StrMPoleSymplectic4QuantPass.c @@ -159,6 +159,7 @@ void StrMPoleSymplectic4QuantPass(double *r, double le, double *A, double *B, ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy; if (!Elem) { double Length, Energy, Scaling; int MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit; @@ -168,8 +169,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -204,6 +205,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->RApertures=RApertures; Elem->KickAngle=KickAngle; } + energy = atEnergy(Param->energy, Elem->Energy); + StrMPoleSymplectic4QuantPass(r_in, Elem->Length, Elem->PolynomA, Elem->PolynomB, Elem->MaxOrder, Elem->NumIntSteps, Elem->FringeQuadEntrance, Elem->FringeQuadExit, @@ -211,7 +214,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->T1, Elem->T2, Elem->R1, Elem->R2, Elem->RApertures, Elem->EApertures, Elem->KickAngle, Elem->Scaling, - Elem->Energy, Param->thread_rng, num_particles); + energy, Param->thread_rng, num_particles); return Elem; } @@ -223,6 +226,8 @@ MODULE_DEF(StrMPoleSymplectic4QuantPass) /* Dummy module initialisation * void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -236,8 +241,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) PolynomB=atGetDoubleArray(ElemData,"PolynomB"); check_error(); MaxOrder=atGetLong(ElemData,"MaxOrder"); check_error(); NumIntSteps=atGetLong(ElemData,"NumIntSteps"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); Scaling=atGetOptionalDouble(ElemData,"FieldScaling",1.0); check_error(); FringeQuadEntrance=atGetOptionalLong(ElemData,"FringeQuadEntrance",0); check_error(); FringeQuadExit=atGetOptionalLong(ElemData,"FringeQuadExit",0); check_error(); @@ -250,10 +255,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) EApertures=atGetOptionalDoubleArray(ElemData,"EApertures"); check_error(); RApertures=atGetOptionalDoubleArray(ElemData,"RApertures"); check_error(); KickAngle=atGetOptionalDoubleArray(ElemData,"KickAngle"); check_error(); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + StrMPoleSymplectic4QuantPass(r_in, Length, PolynomA, PolynomB, MaxOrder, NumIntSteps, FringeQuadEntrance, FringeQuadExit, From e48cb8181ae61271a43677ad5b13094740202f84 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Wed, 14 Aug 2024 18:22:39 +0200 Subject: [PATCH 06/18] energy field no more necessary in Matlab --- atintegrators/BeamLoadingCavityPass.c | 25 ++++++++++++++----------- atintegrators/CavityPass.c | 15 ++++++++++++--- atintegrators/RFCavityPass.c | 3 ++- 3 files changed, 28 insertions(+), 15 deletions(-) diff --git a/atintegrators/BeamLoadingCavityPass.c b/atintegrators/BeamLoadingCavityPass.c index e549cce43..ed239a298 100644 --- a/atintegrators/BeamLoadingCavityPass.c +++ b/atintegrators/BeamLoadingCavityPass.c @@ -51,7 +51,8 @@ void write_buffer(double *data, double *buffer, int datasize, int buffersize){ void BeamLoadingCavityPass(double *r_in,int num_particles,int nbunch, double *bunch_spos,double *bunch_currents, - double circumference,int nturn,struct elem *Elem) { + double circumference,int nturn,double energy, + struct elem *Elem) { /* * r_in - 6-by-N matrix of initial conditions reshaped into * 1-d array of 6*N elements @@ -63,7 +64,6 @@ void BeamLoadingCavityPass(double *r_in,int num_particles,int nbunch, long buffersize = Elem->buffersize; double normfact = Elem->normfact; double le = Elem->Length; - double energy = Elem->Energy; double rffreq = Elem->Frequency; double harmn = Elem->HarmNumber; double tlag = Elem->TimeLag; @@ -156,6 +156,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { double rl = Param->RingLength; + double energy; int nturn=Param->nturn; if (!Elem) { long nslice,nturns,blmode,cavitymode, buffersize; @@ -176,7 +177,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, /*attributes for RF cavity*/ Length=atGetDouble(ElemData,"Length"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); /*attributes for resonator*/ @@ -202,7 +202,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, vbeam_buffer=atGetDoubleArray(ElemData,"_vbeam_buffer"); check_error(); vbunch_buffer=atGetDoubleArray(ElemData,"_vbunch_buffer"); check_error(); /*optional attributes*/ - + Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error(); z_cuts=atGetOptionalDoubleArray(ElemData,"ZCuts"); check_error(); int dimsth[] = {Param->nbunch*nslice*nturns, 4}; @@ -239,6 +239,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->vbeam_buffer = vbeam_buffer; Elem->vbunch_buffer = vbunch_buffer; } + energy = atEnergy(Param->energy, Elem->Energy); + if(num_particlesnbunch){ atError("Number of particles has to be greater or equal to the number of bunches."); }else if (num_particles%Param->nbunch!=0){ @@ -257,7 +259,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, } #endif BeamLoadingCavityPass(r_in,num_particles,Param->nbunch,Param->bunch_spos, - Param->bunch_currents,rl,nturn,Elem); + Param->bunch_currents,rl,nturn,energy,Elem); return Elem; } @@ -267,10 +269,10 @@ MODULE_DEF(BeamLoadingCavityPass) /* Dummy module initialisation */ #if defined(MATLAB_MEX_FILE) void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) -{ - if(nrhs >= 2) - { - +{ + if(nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); @@ -293,7 +295,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double *vbunch_buffer; /*attributes for RF cavity*/ Length=atGetDouble(ElemData,"Length"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); /*attributes for resonator*/ @@ -319,6 +320,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) vbeam_buffer=atGetDoubleArray(ElemData,"_vbeam_buffer"); check_error(); vbunch_buffer=atGetDoubleArray(ElemData,"_vbunch_buffer"); check_error(); /*optional attributes*/ + Energy=atGetOptionalDouble(ElemData,"Energy",0.0); check_error(); z_cuts=atGetOptionalDoubleArray(ElemData,"ZCuts"); check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); @@ -348,6 +350,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) Elem->vgen_buffer = vgen_buffer; Elem->vbeam_buffer = vbeam_buffer; Elem->vbunch_buffer = vbunch_buffer; + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); /* ALLOCATE memory for the output array of the same size as the input */ @@ -356,7 +359,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) double bspos = 0.0; double bcurr = 0.0; - BeamLoadingCavityPass(r_in,num_particles,1,&bspos,&bcurr,1,0,Elem); + BeamLoadingCavityPass(r_in,num_particles,1,&bspos,&bcurr,1,0,Energy,Elem); } else if (nrhs == 0) { /* return list of required fields */ diff --git a/atintegrators/CavityPass.c b/atintegrators/CavityPass.c index 571c93ce7..bce190f57 100644 --- a/atintegrators/CavityPass.c +++ b/atintegrators/CavityPass.c @@ -63,12 +63,14 @@ void CavityPass(double *r_in, double le, double nv, double freq, double lag, dou ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double *r_in, int num_particles, struct parameters *Param) { + double energy = Param->energy; if (!Elem) { double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag; Length=atGetDouble(ElemData,"Length"); check_error(); Voltage=atGetDouble(ElemData,"Voltage"); check_error(); - Energy=atGetDouble(ElemData,"Energy"); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); + /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); @@ -79,7 +81,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, Elem->TimeLag=TimeLag; Elem->PhaseLag=PhaseLag; } - CavityPass(r_in,Elem->Length,Elem->Voltage/Elem->Energy,Elem->Frequency, + if (energy == 0.0) energy = Elem->Energy; + + CavityPass(r_in,Elem->Length,Elem->Voltage/energy,Elem->Frequency, Elem->TimeLag,Elem->PhaseLag,num_particles); return Elem; } @@ -93,19 +97,24 @@ MODULE_DEF(CavityPass) /* Dummy module initialisation */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { if (nrhs >= 2) { + double rest_energy = 0.0; + double charge = -1.0; double *r_in; const mxArray *ElemData = prhs[0]; int num_particles = mxGetN(prhs[1]); double Length=atGetDouble(ElemData,"Length"); double Voltage=atGetDouble(ElemData,"Voltage"); - double Energy=atGetDouble(ElemData,"Energy"); double Frequency=atGetDouble(ElemData,"Frequency"); + double Energy=atGetOptionalDouble(ElemData,"Energy",0.0); double TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); double PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); + if (nrhs > 2) atProperties(prhs[2], &Energy, &rest_energy, &charge); + if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix"); /* ALLOCATE memory for the output array of the same size as the input */ plhs[0] = mxDuplicateArray(prhs[1]); r_in = mxGetDoubles(plhs[0]); + CavityPass(r_in,Length,Voltage/Energy,Frequency,TimeLag,PhaseLag,num_particles); } else if (nrhs == 0) { /* return list of required fields */ diff --git a/atintegrators/RFCavityPass.c b/atintegrators/RFCavityPass.c index 188eb22ec..60bf4c72a 100755 --- a/atintegrators/RFCavityPass.c +++ b/atintegrators/RFCavityPass.c @@ -42,8 +42,9 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem, double Length, Voltage, Energy, Frequency, TimeLag, PhaseLag; Length=atGetDouble(ElemData,"Length"); check_error(); Voltage=atGetDouble(ElemData,"Voltage"); check_error(); - Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); Frequency=atGetDouble(ElemData,"Frequency"); check_error(); + /*optional fields*/ + Energy=atGetOptionalDouble(ElemData,"Energy",energy); check_error(); TimeLag=atGetOptionalDouble(ElemData,"TimeLag",0); check_error(); PhaseLag=atGetOptionalDouble(ElemData,"PhaseLag",0); check_error(); Elem = (struct elem*)atMalloc(sizeof(struct elem)); From 31dc40daf91b7caab021f3a13f941dc580e34ad1 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 15 Aug 2024 18:57:23 +0200 Subject: [PATCH 07/18] Added python tests for lattices without energy in elements or no RingParam in .mat files --- atmat/attrack/elempass.m | 31 +++++++++++++++++++++++++++ machine_data/noenergy.mat | Bin 0 -> 9229 bytes pyat/at/load/utils.py | 6 +++++- pyat/machine_data/noenergy.mat | Bin 0 -> 9229 bytes pyat/test/conftest.py | 14 +++++++++++++ pyat/test/test_linopt6.py | 24 +++++++++++++-------- pyat/test/test_physics.py | 37 ++++++++++++++++++++------------- 7 files changed, 88 insertions(+), 24 deletions(-) create mode 100644 atmat/attrack/elempass.m create mode 100644 machine_data/noenergy.mat create mode 100644 pyat/machine_data/noenergy.mat diff --git a/atmat/attrack/elempass.m b/atmat/attrack/elempass.m new file mode 100644 index 000000000..8d567a71c --- /dev/null +++ b/atmat/attrack/elempass.m @@ -0,0 +1,31 @@ +function rout = elempass(elem,rin,varargin) +%ELEMPASS Tracks particles through a single element +% +%ROUT=ELEMPASS(ELEM,RIN) Tracks particles through ELEM +% +% ELEM: lattice element +% RIN 6xN matrix: input coordinates of N particles +% +% ROUT 6xN matrix: output coordinates of N particles +% +% ROUT=ELEMPASS(...,'PassMethod',PASSMETHOD,...) +% Use PASSMETHOD (default: ELEM.PassMethod) +% +% ROUT=ELEMPASS(...,'Energy',ENERGY,...) +% Use ENERGY and ignore the 'Energy' field of elements +% +% ROUT=ELEMPASS(...,'Particle',PARTICLE,...) +% Use PARTICLE (default: relativistic) +% +% See also: RINGPASS, LINEPASS + +[props.Energy,varargs]=getoption(varargin,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',[]); +[methodname,varargs]=getoption(varargs,'PassMethod',elem.PassMethod); %#ok +if ~isempty(particle) + props.Particle=particle; +end + +rout = feval(methodname,elem,rin,props); + +end \ No newline at end of file diff --git a/machine_data/noenergy.mat b/machine_data/noenergy.mat new file mode 100644 index 0000000000000000000000000000000000000000..fa4192612173cff47d3e9187b17bc53ab826247f GIT binary patch literal 9229 zcmb7}1z6L4yT`ZDB_JRmHA=cc1YwNsMoNUyDX9VqV|1x>BLf5mqDZNrbVxTUB1%X~ zr;I#*(C2x5-t(SwowJK^vG{QftU8HFCNOmA~r<6_J37T$2`&lok_(iVBO0L;tTA z&|ep@x;8oRcQ`u;(ydjZnJ5X3Lp?Px_?e47neTstq74{6b*LJ>{ z#s-Y0m8W(m$0t~p$&;$PAk^FlV{u#Km(6ytkCEvXr|V-J)?7`O$_tYNxbpV`n<{di z3i`lse+XcS8}jRh2lyz)*!wXCB2ibwhmsF-GnJ@5H_E3yq@7h7e4!wl(beTuteMG8 zbg|oXU?hQ~(+mXe-de(A&MRFIu;(;r5s3^D>$$hw7k?iD76>a%|vVL4f)^Lp` zFQ(kf%auXrsGIGC#MIYvqm*%9iMd6j6*f)dnj?ehBS5A|#C0h2b!rgA9ux3S)Mq!= zRUCgt38iRUW-1!w7N?SLs?eSv?n**W%oHeQs17G5w|ushiPATotrQIb!XK1_#;8{t{pk*%pY=$nz|;VF;PO~mV$-Zk)T zuB*Xq)~k+{$|-z){jKH1hV!@QZ|vNz%Y|x3vZ-#xbNM6 zc|-57G@YEpE)HLdSx^{oBk+G5+we(7?u31`a@$Gy-L$Vt1GDAG`kFsX7?u;#S))x@ z_mzrF-m*G$g;)RP)*y-EN}u+8LUeNirSr*TFM3>}tV$r+9_2DXObV? z5;Wg21a}hZD6&hXi+Lu#?LOlUp0c|%i_1-??+ptJb+%kFPUtU%6kv)tm4qCmBFR*; z9>Vf4QQw&3S#;UB;+b2ohekgX?R;`LnJ5p;h>HY~>Ty6XH)YOFJry}!yUvt=7 zADUDgKPUVE`z`OO! z-Jklew$ZeQ9lV!b-)DJ6e>hPGa|xWpew-OFF84ZGlF-o&U+?&~T0`aVEQ#st~xWZCRm25coB)%whR)-9DkA=<$%TAv!$h z`VBAb6Z@{0uM-L@YHJDb0<+S75=OK0t*!SgSmtK6;`y)MI2{kro?p&;h(RX3^5B=y z^|+M7Qr`#+kJD2dI~3>&zKU=kK^s227r%{B?=r4wUVsV*=}3oK5iC8B<{@FA7Vb`k zXt}VbB32r{lA6ApmYzkNNUX*$xTh`9)v0T}q28w{)CF?&><1HzHGw-}s9PqZ^U@y+5tf9fShNvblPHx##jw61; zv4IE&xv2Fl)8r@SxI=Qm;nf?o5nHu%`YkN@@iSoHii1BmgNK`Kux%i^J;la5xU5 zpBMu4hAp?$07sL`5YaDx2_j9PE>OEUJ!88_=b>S7K&?7LK68Dpa33O@e%T`9W_(hn zMPh~mXvI}`aKPK`Q*mvUOIs^3Iw!=GpX57P{fU-u0ibQ z<*z-39qOt-PoBSgg((mCtXp={dAaY~ZsF^vb5u0D318fq1km_qSUVZNZZ5~^_I*`2&|Hp$PiVb5zv%Q1M>45hPIs3lXX=>F#*yoQ~F-q^U_~? zxWGT#x#H^ipzi`gQibf&bQ!Cd2V})t?fy6s0i3vRR5M-P!FiCTdty@SFi_IR@z;Tl zAMT#C%PMiHx<;d_&RC-+HFK=eQXOvIQBRB=8kD?#k=g_C9zG$1%aw_wm@MQLd%*#k zruH83X!{zp;&L@fWQ%<01XPRTate0X`Vk&W+V+velo!l@121)kNwfJ;p`*wV7odWmhYc`DziAf}5*=>C2Za*oTOLAmE;ax!|Jw3<(@i|4fSGA;TU&Hn9i zY0n1JCZ&KHoL~K6kFcJyMj#NlGq})un!UT4z1sGmU{u?T&E#m}!y{^Z`~^m6Bp4E4 zusy6P3SB|{t~`n_k`*nD+?=ln)h8(1_W1gCd-;e2)djd6rf{~%Q0Zh*QmAx1mn(e9 z1B5iTVKF2i|7O96$Iw7#+@0vBi=i8#nlS;v?&@S^y#k%;+JaqLSHb=yD;EeL^IrKG z1TzE}a?-T5;_}qCi5@t_7rXEw%pF*b^XPZMD@VU5JDz)@nY&0A!rPDL45=|2`(5iMd<$2*G zQR~!G3?2w2i=&O)C4#?)g0;Sq&i#lb_F3dT@k9B5$)md7#3Q!+6?=t>G$MDQk^?&W zt`XI;K;b||&c>Ufq{l}&n=V`b?&YL2`Opb_~B!f?3G(^z3_TxmjJPYQ?{* zO~EB9K1=6{dmQ$G;@orDx7P{hT!B(UG{yc@z-B`n2)w&>Rq)4CRcppsX2WP@kD@gv zDR|__$mk-FSo*sA$I@kd%S~7taP(@S{3P02?DD;1o~8FOW5-08l8@-bavljyZ{jO_ z=(s)1MjCkpX=`Z!-EAxV7lCu);|&UnEVr#5y?$Aqf>t3_Y&xN?s5o9elFN zr^<0vgUP@!uG%DG4Z05@waSTTi@oFj2eFsBep1zF?1ee0|*BsTW?e71doBXg6YiyPNghR>%9FD%5q|vNFsbK z3aONvRux=}piUVa$UN)2NTN;~$Z*@Dt;$dbX;D|jpT&rO3E!RDiXiIiMW5q~K8ErJ z;OO@Dh7wMCEDveF6#S}EsIH6UeRX`90QjC6v^EJZHt@OH9uXWXD8Cs|q<4G)cgH6m zsHR&nq9YoRsWP_7UzJMzt$gQFcOD8@n$pF84-aP@Y#$|(J+&F(>PKBhLy?~92||O8CCw?dcv2-@oEcqLri>J{ zq`xbOk+%_GONy&OPHZi5vjs`fUu;jl=}Ec9tzB?C$Qaxqfrev7`X=>hX)6+MFpG#% zn3M_UE!u`0W=(8GWpCt1%6Tn4MPj3q)rTsDIp;9+IIG>dHs4FGj(X+?&rv^cNp2J~ z?c|fIrHdAjtDsIbjWNQE)=?w;=hM?J^z4$BCtgq|H5z8+vu-lqu;DSNpsq8GaU>?O zw3!Kg1eSoo187;;F@B0sNbhbdOdWOhoRHL#D>@-KVqDoDjShp`fo*&F zGEYoEw@f8vrdcL}KTc)iP)988>T2@o2!O>xlc;mkQNzVj-9g!y^o_=^k%7sEcQ8pDWUB&!inL0PY} z_G-rRgVuPrz6K{-CG1k&v|J%Id-iVayN0L#X2_o-I9!et+Kxk-+fHC=i`22VM3zw{ z#(UU=2kC3LP!^rd;>KZ%fx5)0#8F5M=QmQirWb^!&kOKVuN2h){*H5=6T?Pw!!8xd zr=r_2HWMBWmIB%f^0hY`ejr1*NaYBBq)Q=JQ_xAfST#D;JrR;UI-A&QN9D=OdgX~% z9w3^bM@B~G@#K6ry0p2{jJ5ePh6AXwbuJPvIx3EMpr~`|zofo$uMt>mf`j4>8$>TI z$z7!mB$SvP=Tj)CA8Vnx#GJgJW)Y}ZnElEQn3OA}+cV3zueyO2&_w*RDODB^+fu7c ziEr>Lzch4VE)JgqdWw&U4#CNmNLQ?+WZZ&!g?evr*_NlO#0sADg}JxSi%9%V$RZae z`(C$Zvzbm(5{xl@;WVb(y_yTt0B&cm-7~(vB^2i}+`1@4Q5vM+x13@&a9?CXTjN3c zw^dE^S}K!zxX>Z@3v!3SIt#Y4V@Zo{PJ?ZcVx@d#f2Dyfy9lG!HXKJ-{H>fj>h-0)U`b3}rHz+_}iU;V z`to;FwUY;ML*;Nt?bwMOCn9e$Zg1wvOu~V4>!Kg$jSfD6h_8NgOsfOd6{|mam|yPf zyvi)Uw+%eTZ(%(xV%|R^=txE%fs{JDB@(po@{wv278Ba$j~F|0hPc!`r7rK}&;vYkmK0Q?`N#!O5efstj2zYS`0)APE@ ze~^)|BNpsWN(yL-(olz5HO%ER7vzIUQ`BIG5}bt(hXhA5`6W)9dvq-_y0!#b@`slL ztQTzor)ah#K$2*}1)h z6XD?LZKK1hB+e}%!YL-i)Y~hh&fnELz%FVrHrRsK=1aWsNK%+n5a~C#6S|+nQ{vB= z!I@IdozHPf;VpMP6X;3CW~GLmr2mG z>zI!k)A;%0sv$lZT-bhsdo8)9`H&^Aj7{`9n z_#E}_mltD)nn3lux?UT_eX=dzEuKYXIa1z8j_%P% z@~6A4PRA(CNvG5MgF7!DsT_WaI+2k0aoc>SNA>WN;;XM3S)0aG5z`Xiq*vxmdLvYV zZ0;b#8;*@C52%`s+|4!K9qG6s=Nigy8!uHD(j9}nyJJb&T9|%SyB%Px&U_HGkFf3k zl)hOgZh}M@D11O7Boh{dXiDA=LK<4&Z3M_NhFU-mmu6%2%;U|#s6`=HJ9%mmZ?Gc#c> z+6>dvX8J>6ow~|}cyKDp6IQM#Ux!W@FMUC>WcUAT_~$ejQ+;&5bO!+ib@Bw!MendM z#K6KXd#}>44qLiS6O)i@9BG!pexy~YVHpIhl5WAkJ!WF~<0)0M$w{+qJ{0TJaih8n zO;|?=iT@lwEqlcrf2snF=FoC-W*UV5CnL~sg*Z6Rm6R^xq(pUU(3M{rG}W->x4*he z;^e)&KEap3TaNLl!KxYi&1#276!{%c1n*LTOR1r(|i~(W3#Rya`VF z%tKW&*dhKqvcAjSl|4x#BhteXupqrE-yg~r1rGr@$c zNgn5yi8Hz`wX~ODG)5jTXIN;V$>hrCy%8I0AWAi|Yv?1UsR8E;<{@VTpgFA_|-i(#^E{!nTjtJ#oM|V?86WQD{sN$wgnU#$@uNp73qGn9;6$iUq zc(W8iXN1!)-UTNzncc8%wqfB*420-4XzJA^;2AaIfRF!Zk($y@UxH3X{)6a83r|=z z?|urja3;Utdu6Ld7UEdLj(637dLL6rgGohYG;fkAA1*m*>v9K2ByTo}zGyZ|KuKM5^Jo$&5% z@95T_k;0|o{OiiPZ;(VvGiFFO!X#OA&*tPS8vna z9jY4|2}vl6F6-F;FmOpdE%Ak1fa{MpPrNMYPKq3`k)?quLWhel;)d@{Ibp*%#N+3Z zx~6#0C!X5#C(B-s*5*bgla3ztB0nUj3LSQR2(mvY`tDQLr`p5C`-4FM9ie{C*D&8k z{8b<7j5G#F50%gU7+Rfu6O;>thq)1DsuA{Ve3f$~_b64jwID$Bvdr8FUf)xDp3k3_ zEdS0i_prD?eh@_Ov#(i7)Lb`*0nY|b5#}cAp5e7|N4s3sfXtJ*q-vczd52s4l|5We zScQAoS7|METfzw5ijO((gPF{Ci+dsIO)97D~07v~r zqB#Di_}VVa>VZ!>Wt^Ommh+#$2YJ-roQk)fcuq{Aoi1+%0=LYKEg?}u6MzF}xDTL5 z0F1NB9WV&x^C9szNYWHd(uQ1&jcs1@gd&*#Id}hoskr`yPs*1Z6-I;r`>i7xqf~YnW=vup`hvW)-fSYSuQhcKu7V<0K@*$KKZP{|;XwnmF@&cHCV>vvoOL}j& ze*tL#L>I~n|7ueN3cRz?oGNHWuQ$RObp|Be2tTaD;c|Zv$+PqV7A1&9+I*L?+da5z zVe&C}yB%_6xaA2!w=SH_cNV33!9O`yRggV_M0M7Hvi2GWhLPl+2xGiA>CJ#M76+xV zwd?pqdm8d6cI@$&r-vk2NC)xbi=^pLZ{9ZzHo6pM*zD-wnt=X?`Em4XLEb%8c?K{# z5z^zWydu<+*XkU@)Q&OS1`GN(Zbt0fj{6rkBQ8(DahOlzNVS_^prjXcbDDn4k{D*Kos**)S8f4s`UpjYO@3LgdNrr_X!Ft03787ydK=vo>Jj zcufNCbFPvHYnX8mm2L;*&3I(nzv@Y97wX z|F|D>$5wgeBy7Zp!7AZjgzWLlDIu0#x(E}j@OVSM`JoHEyKrdlgVUi4cM8CE0Y27y zRNAv|i;V%Zpz(iMDg>d1wpLfKrXwgM_b9O~vG+&@nnWHY@Egj8gCS=dkYYx5BMVUX zxhXFXTl9tIq6yBdY(`?HNurqnJ&@0PL3bm;iBfK+z1B*ikwo2W?~hlU_Cog?u>KjG zSY@TpJz+LNR1=m4yCR~Am^h0=!~roKPCrPw=a9o8cDsXOsZ@cFqrQn{FvLT!+sihIF}Qy%pNYvA*l;S06Z$ z4>zKst34y#;Xz?{51~#M=Gms(h65mk{QZ`euSoo`H0_sGfyh?t}E{6Rr<1%80xh zFjD~6ZUtLEGoE_myxjC}7*6mCfNQH!km`Mqs}k>Z;0B@mjkzd}FpFGxup|z?Cm-bM z_4Cb*#kNMw#qC7Hgm7L_&3Kw5LqWW^bmVHaZQ%~gQRG{Olm^X8?$Si{Ro@rU*E?2* zzRDMWTL^oL$DGsn<6!cpSkf6r8+=_KQYyp7Rkt{;YjIF)g=DigsOHy_$2!|I!QX8q zjrqmY9LC-1bEuS+=$@np!an*e?4wPIHTjjForQhJv#{^+xJg5rLqtzK@_888wK=v> zt=#~iXgNImeqE~sZ};imzGqO_NisZwA+550dYbHOAEb?zcsrT|;6P*$H+x+jPrS}W%=IqkMT61ycy zkocm}QPy|P_K^l@knPAoSOJ=-H?)m~HSi%fuy0U`8dViIvlT!D#rwS+g^aNL!4#${ zA;~NOhMJ;VhIGHyLQR*mec+#QFGO`U{tc+FMEpBTpwMBJJa*SaIzM>qds*dxYn%n@ z3tyFM(mMBMwEC?~`0>CJbP&RV8=sF~Mp*^y#yo9#zEv1Os4uP4Ev2V#2@ke+rc8Lt z))Wvr)~>;uWCX>x;Q{gJ7!8sw>knMDOKvlJ>5Y;`{o_yHh3eI|6Wt{t-FV*UKPnWw z3liv|aY-9F-&m+OZc*_}=FgdEbnCJM5eGe+90WuhL8^fzEY`6Jn;#0ymxbUhD!o>w z!YL&|k%Um@+J{Plj0u{`vorD8-uz>&e(-J)x6Rd?nTCnc-0wRdzdR!~T>7iB4*0W2 z?WAwgj90l`J`?c0iO`Tj^a}6v1ozW^dWk)aEGQkiN1ZaR?;g=w#`J1vw$-M@BaO#h z^+lF2*PH?sQAZGOAr^A+i+=ezb{RW7J%=vw2E5m# zC{wFK4X6xJg>-*xoe&HQK4-v19QxmZOTn>x=RZcdxPZ_gFvN#n2ecG=6wmqs`h=dW Sx2#Af%&d?Gfq2BXKK>VqV_A^^ literal 0 HcmV?d00001 diff --git a/pyat/at/load/utils.py b/pyat/at/load/utils.py index 8d4872ab9..55d1a7b5e 100644 --- a/pyat/at/load/utils.py +++ b/pyat/at/load/utils.py @@ -86,7 +86,11 @@ class RingParam(elt.Element): "Periodicity", ] _conversions = dict( - elt.Element._conversions, Energy=float, Periodicity=int, Particle=_particle + elt.Element._conversions, + Energy=float, + Periodicity=int, + Particle=_particle, + cell_harmnumber=int, ) # noinspection PyPep8Naming diff --git a/pyat/machine_data/noenergy.mat b/pyat/machine_data/noenergy.mat new file mode 100644 index 0000000000000000000000000000000000000000..fa4192612173cff47d3e9187b17bc53ab826247f GIT binary patch literal 9229 zcmb7}1z6L4yT`ZDB_JRmHA=cc1YwNsMoNUyDX9VqV|1x>BLf5mqDZNrbVxTUB1%X~ zr;I#*(C2x5-t(SwowJK^vG{QftU8HFCNOmA~r<6_J37T$2`&lok_(iVBO0L;tTA z&|ep@x;8oRcQ`u;(ydjZnJ5X3Lp?Px_?e47neTstq74{6b*LJ>{ z#s-Y0m8W(m$0t~p$&;$PAk^FlV{u#Km(6ytkCEvXr|V-J)?7`O$_tYNxbpV`n<{di z3i`lse+XcS8}jRh2lyz)*!wXCB2ibwhmsF-GnJ@5H_E3yq@7h7e4!wl(beTuteMG8 zbg|oXU?hQ~(+mXe-de(A&MRFIu;(;r5s3^D>$$hw7k?iD76>a%|vVL4f)^Lp` zFQ(kf%auXrsGIGC#MIYvqm*%9iMd6j6*f)dnj?ehBS5A|#C0h2b!rgA9ux3S)Mq!= zRUCgt38iRUW-1!w7N?SLs?eSv?n**W%oHeQs17G5w|ushiPATotrQIb!XK1_#;8{t{pk*%pY=$nz|;VF;PO~mV$-Zk)T zuB*Xq)~k+{$|-z){jKH1hV!@QZ|vNz%Y|x3vZ-#xbNM6 zc|-57G@YEpE)HLdSx^{oBk+G5+we(7?u31`a@$Gy-L$Vt1GDAG`kFsX7?u;#S))x@ z_mzrF-m*G$g;)RP)*y-EN}u+8LUeNirSr*TFM3>}tV$r+9_2DXObV? z5;Wg21a}hZD6&hXi+Lu#?LOlUp0c|%i_1-??+ptJb+%kFPUtU%6kv)tm4qCmBFR*; z9>Vf4QQw&3S#;UB;+b2ohekgX?R;`LnJ5p;h>HY~>Ty6XH)YOFJry}!yUvt=7 zADUDgKPUVE`z`OO! z-Jklew$ZeQ9lV!b-)DJ6e>hPGa|xWpew-OFF84ZGlF-o&U+?&~T0`aVEQ#st~xWZCRm25coB)%whR)-9DkA=<$%TAv!$h z`VBAb6Z@{0uM-L@YHJDb0<+S75=OK0t*!SgSmtK6;`y)MI2{kro?p&;h(RX3^5B=y z^|+M7Qr`#+kJD2dI~3>&zKU=kK^s227r%{B?=r4wUVsV*=}3oK5iC8B<{@FA7Vb`k zXt}VbB32r{lA6ApmYzkNNUX*$xTh`9)v0T}q28w{)CF?&><1HzHGw-}s9PqZ^U@y+5tf9fShNvblPHx##jw61; zv4IE&xv2Fl)8r@SxI=Qm;nf?o5nHu%`YkN@@iSoHii1BmgNK`Kux%i^J;la5xU5 zpBMu4hAp?$07sL`5YaDx2_j9PE>OEUJ!88_=b>S7K&?7LK68Dpa33O@e%T`9W_(hn zMPh~mXvI}`aKPK`Q*mvUOIs^3Iw!=GpX57P{fU-u0ibQ z<*z-39qOt-PoBSgg((mCtXp={dAaY~ZsF^vb5u0D318fq1km_qSUVZNZZ5~^_I*`2&|Hp$PiVb5zv%Q1M>45hPIs3lXX=>F#*yoQ~F-q^U_~? zxWGT#x#H^ipzi`gQibf&bQ!Cd2V})t?fy6s0i3vRR5M-P!FiCTdty@SFi_IR@z;Tl zAMT#C%PMiHx<;d_&RC-+HFK=eQXOvIQBRB=8kD?#k=g_C9zG$1%aw_wm@MQLd%*#k zruH83X!{zp;&L@fWQ%<01XPRTate0X`Vk&W+V+velo!l@121)kNwfJ;p`*wV7odWmhYc`DziAf}5*=>C2Za*oTOLAmE;ax!|Jw3<(@i|4fSGA;TU&Hn9i zY0n1JCZ&KHoL~K6kFcJyMj#NlGq})un!UT4z1sGmU{u?T&E#m}!y{^Z`~^m6Bp4E4 zusy6P3SB|{t~`n_k`*nD+?=ln)h8(1_W1gCd-;e2)djd6rf{~%Q0Zh*QmAx1mn(e9 z1B5iTVKF2i|7O96$Iw7#+@0vBi=i8#nlS;v?&@S^y#k%;+JaqLSHb=yD;EeL^IrKG z1TzE}a?-T5;_}qCi5@t_7rXEw%pF*b^XPZMD@VU5JDz)@nY&0A!rPDL45=|2`(5iMd<$2*G zQR~!G3?2w2i=&O)C4#?)g0;Sq&i#lb_F3dT@k9B5$)md7#3Q!+6?=t>G$MDQk^?&W zt`XI;K;b||&c>Ufq{l}&n=V`b?&YL2`Opb_~B!f?3G(^z3_TxmjJPYQ?{* zO~EB9K1=6{dmQ$G;@orDx7P{hT!B(UG{yc@z-B`n2)w&>Rq)4CRcppsX2WP@kD@gv zDR|__$mk-FSo*sA$I@kd%S~7taP(@S{3P02?DD;1o~8FOW5-08l8@-bavljyZ{jO_ z=(s)1MjCkpX=`Z!-EAxV7lCu);|&UnEVr#5y?$Aqf>t3_Y&xN?s5o9elFN zr^<0vgUP@!uG%DG4Z05@waSTTi@oFj2eFsBep1zF?1ee0|*BsTW?e71doBXg6YiyPNghR>%9FD%5q|vNFsbK z3aONvRux=}piUVa$UN)2NTN;~$Z*@Dt;$dbX;D|jpT&rO3E!RDiXiIiMW5q~K8ErJ z;OO@Dh7wMCEDveF6#S}EsIH6UeRX`90QjC6v^EJZHt@OH9uXWXD8Cs|q<4G)cgH6m zsHR&nq9YoRsWP_7UzJMzt$gQFcOD8@n$pF84-aP@Y#$|(J+&F(>PKBhLy?~92||O8CCw?dcv2-@oEcqLri>J{ zq`xbOk+%_GONy&OPHZi5vjs`fUu;jl=}Ec9tzB?C$Qaxqfrev7`X=>hX)6+MFpG#% zn3M_UE!u`0W=(8GWpCt1%6Tn4MPj3q)rTsDIp;9+IIG>dHs4FGj(X+?&rv^cNp2J~ z?c|fIrHdAjtDsIbjWNQE)=?w;=hM?J^z4$BCtgq|H5z8+vu-lqu;DSNpsq8GaU>?O zw3!Kg1eSoo187;;F@B0sNbhbdOdWOhoRHL#D>@-KVqDoDjShp`fo*&F zGEYoEw@f8vrdcL}KTc)iP)988>T2@o2!O>xlc;mkQNzVj-9g!y^o_=^k%7sEcQ8pDWUB&!inL0PY} z_G-rRgVuPrz6K{-CG1k&v|J%Id-iVayN0L#X2_o-I9!et+Kxk-+fHC=i`22VM3zw{ z#(UU=2kC3LP!^rd;>KZ%fx5)0#8F5M=QmQirWb^!&kOKVuN2h){*H5=6T?Pw!!8xd zr=r_2HWMBWmIB%f^0hY`ejr1*NaYBBq)Q=JQ_xAfST#D;JrR;UI-A&QN9D=OdgX~% z9w3^bM@B~G@#K6ry0p2{jJ5ePh6AXwbuJPvIx3EMpr~`|zofo$uMt>mf`j4>8$>TI z$z7!mB$SvP=Tj)CA8Vnx#GJgJW)Y}ZnElEQn3OA}+cV3zueyO2&_w*RDODB^+fu7c ziEr>Lzch4VE)JgqdWw&U4#CNmNLQ?+WZZ&!g?evr*_NlO#0sADg}JxSi%9%V$RZae z`(C$Zvzbm(5{xl@;WVb(y_yTt0B&cm-7~(vB^2i}+`1@4Q5vM+x13@&a9?CXTjN3c zw^dE^S}K!zxX>Z@3v!3SIt#Y4V@Zo{PJ?ZcVx@d#f2Dyfy9lG!HXKJ-{H>fj>h-0)U`b3}rHz+_}iU;V z`to;FwUY;ML*;Nt?bwMOCn9e$Zg1wvOu~V4>!Kg$jSfD6h_8NgOsfOd6{|mam|yPf zyvi)Uw+%eTZ(%(xV%|R^=txE%fs{JDB@(po@{wv278Ba$j~F|0hPc!`r7rK}&;vYkmK0Q?`N#!O5efstj2zYS`0)APE@ ze~^)|BNpsWN(yL-(olz5HO%ER7vzIUQ`BIG5}bt(hXhA5`6W)9dvq-_y0!#b@`slL ztQTzor)ah#K$2*}1)h z6XD?LZKK1hB+e}%!YL-i)Y~hh&fnELz%FVrHrRsK=1aWsNK%+n5a~C#6S|+nQ{vB= z!I@IdozHPf;VpMP6X;3CW~GLmr2mG z>zI!k)A;%0sv$lZT-bhsdo8)9`H&^Aj7{`9n z_#E}_mltD)nn3lux?UT_eX=dzEuKYXIa1z8j_%P% z@~6A4PRA(CNvG5MgF7!DsT_WaI+2k0aoc>SNA>WN;;XM3S)0aG5z`Xiq*vxmdLvYV zZ0;b#8;*@C52%`s+|4!K9qG6s=Nigy8!uHD(j9}nyJJb&T9|%SyB%Px&U_HGkFf3k zl)hOgZh}M@D11O7Boh{dXiDA=LK<4&Z3M_NhFU-mmu6%2%;U|#s6`=HJ9%mmZ?Gc#c> z+6>dvX8J>6ow~|}cyKDp6IQM#Ux!W@FMUC>WcUAT_~$ejQ+;&5bO!+ib@Bw!MendM z#K6KXd#}>44qLiS6O)i@9BG!pexy~YVHpIhl5WAkJ!WF~<0)0M$w{+qJ{0TJaih8n zO;|?=iT@lwEqlcrf2snF=FoC-W*UV5CnL~sg*Z6Rm6R^xq(pUU(3M{rG}W->x4*he z;^e)&KEap3TaNLl!KxYi&1#276!{%c1n*LTOR1r(|i~(W3#Rya`VF z%tKW&*dhKqvcAjSl|4x#BhteXupqrE-yg~r1rGr@$c zNgn5yi8Hz`wX~ODG)5jTXIN;V$>hrCy%8I0AWAi|Yv?1UsR8E;<{@VTpgFA_|-i(#^E{!nTjtJ#oM|V?86WQD{sN$wgnU#$@uNp73qGn9;6$iUq zc(W8iXN1!)-UTNzncc8%wqfB*420-4XzJA^;2AaIfRF!Zk($y@UxH3X{)6a83r|=z z?|urja3;Utdu6Ld7UEdLj(637dLL6rgGohYG;fkAA1*m*>v9K2ByTo}zGyZ|KuKM5^Jo$&5% z@95T_k;0|o{OiiPZ;(VvGiFFO!X#OA&*tPS8vna z9jY4|2}vl6F6-F;FmOpdE%Ak1fa{MpPrNMYPKq3`k)?quLWhel;)d@{Ibp*%#N+3Z zx~6#0C!X5#C(B-s*5*bgla3ztB0nUj3LSQR2(mvY`tDQLr`p5C`-4FM9ie{C*D&8k z{8b<7j5G#F50%gU7+Rfu6O;>thq)1DsuA{Ve3f$~_b64jwID$Bvdr8FUf)xDp3k3_ zEdS0i_prD?eh@_Ov#(i7)Lb`*0nY|b5#}cAp5e7|N4s3sfXtJ*q-vczd52s4l|5We zScQAoS7|METfzw5ijO((gPF{Ci+dsIO)97D~07v~r zqB#Di_}VVa>VZ!>Wt^Ommh+#$2YJ-roQk)fcuq{Aoi1+%0=LYKEg?}u6MzF}xDTL5 z0F1NB9WV&x^C9szNYWHd(uQ1&jcs1@gd&*#Id}hoskr`yPs*1Z6-I;r`>i7xqf~YnW=vup`hvW)-fSYSuQhcKu7V<0K@*$KKZP{|;XwnmF@&cHCV>vvoOL}j& ze*tL#L>I~n|7ueN3cRz?oGNHWuQ$RObp|Be2tTaD;c|Zv$+PqV7A1&9+I*L?+da5z zVe&C}yB%_6xaA2!w=SH_cNV33!9O`yRggV_M0M7Hvi2GWhLPl+2xGiA>CJ#M76+xV zwd?pqdm8d6cI@$&r-vk2NC)xbi=^pLZ{9ZzHo6pM*zD-wnt=X?`Em4XLEb%8c?K{# z5z^zWydu<+*XkU@)Q&OS1`GN(Zbt0fj{6rkBQ8(DahOlzNVS_^prjXcbDDn4k{D*Kos**)S8f4s`UpjYO@3LgdNrr_X!Ft03787ydK=vo>Jj zcufNCbFPvHYnX8mm2L;*&3I(nzv@Y97wX z|F|D>$5wgeBy7Zp!7AZjgzWLlDIu0#x(E}j@OVSM`JoHEyKrdlgVUi4cM8CE0Y27y zRNAv|i;V%Zpz(iMDg>d1wpLfKrXwgM_b9O~vG+&@nnWHY@Egj8gCS=dkYYx5BMVUX zxhXFXTl9tIq6yBdY(`?HNurqnJ&@0PL3bm;iBfK+z1B*ikwo2W?~hlU_Cog?u>KjG zSY@TpJz+LNR1=m4yCR~Am^h0=!~roKPCrPw=a9o8cDsXOsZ@cFqrQn{FvLT!+sihIF}Qy%pNYvA*l;S06Z$ z4>zKst34y#;Xz?{51~#M=Gms(h65mk{QZ`euSoo`H0_sGfyh?t}E{6Rr<1%80xh zFjD~6ZUtLEGoE_myxjC}7*6mCfNQH!km`Mqs}k>Z;0B@mjkzd}FpFGxup|z?Cm-bM z_4Cb*#kNMw#qC7Hgm7L_&3Kw5LqWW^bmVHaZQ%~gQRG{Olm^X8?$Si{Ro@rU*E?2* zzRDMWTL^oL$DGsn<6!cpSkf6r8+=_KQYyp7Rkt{;YjIF)g=DigsOHy_$2!|I!QX8q zjrqmY9LC-1bEuS+=$@np!an*e?4wPIHTjjForQhJv#{^+xJg5rLqtzK@_888wK=v> zt=#~iXgNImeqE~sZ};imzGqO_NisZwA+550dYbHOAEb?zcsrT|;6P*$H+x+jPrS}W%=IqkMT61ycy zkocm}QPy|P_K^l@knPAoSOJ=-H?)m~HSi%fuy0U`8dViIvlT!D#rwS+g^aNL!4#${ zA;~NOhMJ;VhIGHyLQR*mec+#QFGO`U{tc+FMEpBTpwMBJJa*SaIzM>qds*dxYn%n@ z3tyFM(mMBMwEC?~`0>CJbP&RV8=sF~Mp*^y#yo9#zEv1Os4uP4Ev2V#2@ke+rc8Lt z))Wvr)~>;uWCX>x;Q{gJ7!8sw>knMDOKvlJ>5Y;`{o_yHh3eI|6Wt{t-FV*UKPnWw z3liv|aY-9F-&m+OZc*_}=FgdEbnCJM5eGe+90WuhL8^fzEY`6Jn;#0ymxbUhD!o>w z!YL&|k%Um@+J{Plj0u{`vorD8-uz>&e(-J)x6Rd?nTCnc-0wRdzdR!~T>7iB4*0W2 z?WAwgj90l`J`?c0iO`Tj^a}6v1ozW^dWk)aEGQkiN1ZaR?;g=w#`J1vw$-M@BaO#h z^+lF2*PH?sQAZGOAr^A+i+=ezb{RW7J%=vw2E5m# zC{wFK4X6xJg>-*xoe&HQK4-v19QxmZOTn>x=RZcdxPZ_gFvN#n2ecG=6wmqs`h=dW Sx2#Af%&d?Gfq2BXKK>VqV_A^^ literal 0 HcmV?d00001 diff --git a/pyat/test/conftest.py b/pyat/test/conftest.py index c8364d24b..fa1afc6a3 100644 --- a/pyat/test/conftest.py +++ b/pyat/test/conftest.py @@ -58,3 +58,17 @@ def hmba_lattice(): with as_file(files(machine_data) / 'hmba.mat') as path: ring = at.load_lattice(path) return ring + + +@pytest.fixture(scope='session') +def noenergy_lattice(): + with as_file(files(machine_data) / 'noenergy.mat') as path: + ring = at.load_lattice(path) + return ring + + +@pytest.fixture(scope='session') +def noringparam_lattice(): + with as_file(files(machine_data) / 'noringparam.mat') as path: + ring = at.load_lattice(path) + return ring diff --git a/pyat/test/test_linopt6.py b/pyat/test/test_linopt6.py index ba3c60862..f22bf1fab 100644 --- a/pyat/test/test_linopt6.py +++ b/pyat/test/test_linopt6.py @@ -4,7 +4,10 @@ from at import linopt2, linopt4, linopt6, get_optics -@pytest.mark.parametrize("lattice", ["dba_lattice", "hmba_lattice"]) +@pytest.mark.parametrize( + "lattice", + ["dba_lattice", "noenergy_lattice", "hmba_lattice", "noringparam_lattice"], +) def test_linopt6_norad(request, lattice): """Compare the results of linopt2 and linopt6 in 4d""" lattice = request.getfixturevalue(lattice) @@ -19,13 +22,16 @@ def test_linopt6_norad(request, lattice): assert_close(ld2.W, ld6.W, atol=1e-6, rtol=0) -@pytest.mark.parametrize("lattice", ["hmba_lattice"]) +@pytest.mark.parametrize( + "lattice", + ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) def test_linopt6_rad(request, lattice): """Compare the results with and without radiation""" lattice = request.getfixturevalue(lattice) refpts = range(len(lattice) + 1) # Turn cavity ON, without radiation - radlattice = lattice.radiation_on(dipole_pass=None, copy=True) + radlattice = lattice.enable_6d(dipole_pass=None, copy=True) ld04, rd4, ld4 = linopt6(lattice, refpts, get_w=True) ld06, rd6, ld6 = linopt6(radlattice, refpts, get_w=True) @@ -47,12 +53,12 @@ def test_linopt6_line(request, lattice, dp, method): refpts = lattice.uint32_refpts(range(len(lattice) + 1)) ld04, rd4, ld4 = get_optics(lattice, refpts, dp=dp, method=method) - twin = dict( - alpha=ld04.alpha, - beta=ld04.beta, - dispersion=ld04.dispersion, - closed_orbit=ld04.closed_orbit, - ) + twin = { + "alpha": ld04.alpha, + "beta": ld04.beta, + "dispersion": ld04.dispersion, + "closed_orbit": ld04.closed_orbit, + } # twin = ld04 tr04, bd4, tr4 = get_optics(lattice, refpts, dp=dp, twiss_in=twin, method=method) for field in ["s_pos", "closed_orbit", "dispersion", "alpha", "beta", "mu"]: diff --git a/pyat/test/test_physics.py b/pyat/test/test_physics.py index cd2519acf..65af3e268 100644 --- a/pyat/test/test_physics.py +++ b/pyat/test/test_physics.py @@ -1,12 +1,12 @@ -import at import numpy +import pytest from numpy.testing import assert_allclose as assert_close from numpy.testing import assert_equal -import pytest + +import at from at import AtWarning, physics -from at import lattice_track from at import lattice_pass, internal_lpass - +from at import lattice_track DP = 1e-5 DP2 = 0.005 @@ -101,10 +101,13 @@ def test_find_m44_returns_same_answer_as_matlab(dba_lattice, refpts): assert mstack.shape == (len(refpts), 4, 4) +@pytest.mark.parametrize( + "lattice", ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) @pytest.mark.parametrize('refpts', ([145], [20], [1, 2, 3])) -def test_find_m66(hmba_lattice, refpts): - hmba_lattice = hmba_lattice.radiation_on(copy=True) - m66, mstack = physics.find_m66(hmba_lattice, refpts=refpts) +def test_find_m66(request, lattice, refpts): + lattice = request.getfixturevalue(lattice).enable_6d(copy=True) + m66, mstack = lattice.find_m66(refpts=refpts) assert_close(m66, M66_MATLAB, rtol=0, atol=1e-8) stack_size = 0 if refpts is None else len(refpts) assert mstack.shape == (stack_size, 6, 6) @@ -139,10 +142,13 @@ def test_find_sync_orbit_finds_zeros(dba_lattice): numpy.testing.assert_equal(sync_orbit, numpy.zeros(6)) -def test_find_orbit6(hmba_lattice): - hmba_lattice = hmba_lattice.radiation_on(copy=True) - refpts = numpy.ones(len(hmba_lattice), dtype=bool) - orbit6, all_points = physics.find_orbit6(hmba_lattice, refpts) +@pytest.mark.parametrize( + "lattice", ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) +def test_find_orbit6(request, lattice): + lattice = request.getfixturevalue(lattice).enable_6d(copy=True) + refpts = numpy.ones(len(lattice), dtype=bool) + orbit6, all_points = lattice.find_orbit6(refpts) assert_close(orbit6, orbit6_MATLAB, rtol=0, atol=1e-12) @@ -350,10 +356,13 @@ def test_simple_ring(): assert_close(ring.get_tune(), [0.1, 0.2], atol=1e-10) +@pytest.mark.parametrize( + "lattice", ["hmba_lattice", "noenergy_lattice", "noringparam_lattice"], +) @pytest.mark.parametrize('refpts', ([121], [0, 40, 121])) -def test_ohmi_envelope(hmba_lattice, refpts): - hmba_lattice = hmba_lattice.radiation_on(copy=True) - emit0, beamdata, emit = hmba_lattice.ohmi_envelope(refpts) +def test_ohmi_envelope(request, lattice, refpts): + lattice = request.getfixturevalue(lattice).enable_6d(copy=True) + emit0, beamdata, emit = lattice.ohmi_envelope(refpts) obs = emit[-1] # All expected values are Matlab results From 6692ed5253f7aab26c85015a53139b172523fa60 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 15 Aug 2024 19:01:16 +0200 Subject: [PATCH 08/18] Added test lattice --- machine_data/noringparam.mat | Bin 0 -> 9008 bytes pyat/machine_data/noringparam.mat | Bin 0 -> 9008 bytes 2 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 machine_data/noringparam.mat create mode 100644 pyat/machine_data/noringparam.mat diff --git a/machine_data/noringparam.mat b/machine_data/noringparam.mat new file mode 100644 index 0000000000000000000000000000000000000000..40185d4700fec78b234b7906c20d65ac0de55aab GIT binary patch literal 9008 zcma)>2T)V{+VzvrYv@GpRn{Fy!f{x`qlfGBzhN6G2r&^@^FbQ%g8UO7tIb z6_fU6=a&5SX{t8M*=O@9%$PPRVMc{erNYmUR*9j2WHUBl2I{x{I8R9n8nN`@xG1$Z zW=d}`O@z&?YQv!+*yD5Hc5inH0f;8o9$CFrXHv(KPlXj__GKfeM=rGa#(Fodg4^)y zLqRdM(r=k-mc4cAxSoo4k<)}SI{DPGsQJ*5Xyj*V_y}5y7naGax(Gy$)m)V0dPm{U zU&3(j@l)ZXh3INS;?#PO=9qMzIFG8rO;B6-(JJvuNPx(p8`t5gX;FMfe!}^R<`ThqDA5}<)EX% z4aRRgTGgK;?R|kRLWDtBk{_#d=jfDK%6o;p?C}u!-$N zWZ0oQxV7)vkDAzyn$(z!Jg`~vtl6c^G2aao^E?R&W%2uz8+1oGf7eNWq%AmEv{CH- z8;Q8R^Y46XXKlJ-$jObLO1mAy1?!17Y>8Xy0(s9{{QBOmPn@f9>r`;!A*PQjs z!^NoaRNIoHpZK#sOwHYiv;-a;#e42cu2`BE&1)pZ_*QR7r-XxpxDibiYo=;Y*Ym6K zsTXg_T!MD{>3buDSfA>c!;_C1*V)VtXv7a}UOhMWwMyD_o~3?u?gw@MH#mW8Msm;Y z?nj#E8c(a~K8PIpQ`TaYpNE|6c&$$R+_ADKudcTSFyQq|N^3TZ`sz}^)Q=}~6_RU=G5NJaCiSUm^ z_=zmlY0MP|TQRXOlYQ86@jG`+6>p(_TwMPAczI-gZZ5K_p|?{ZVq5Fhk3O7_SbuMM z)$6+njW36;`%Wv2&e}j3byyNBFQ5Ctwih z6B!0d+epLAoln80p0mg@Pf5+TKxetprI=atG?AoLBPElm68)337Q9s3HZB^R91MSd za!IQZ4E|A(UmF|UEFdEmZ4vm~nbi~J{^JWkw+ZTPTDBaN>-M>CNe7N% z#vf!Dm=3pV6*LTNk!d@_y5V1%O%Mtc8&=Y5h@Dg(~)v=s;u8R%4e;AnG>X?+Y$oy%+h$OfH}%=<#EOPIoCpW3UVE&Tkb9 zE)G%^#H4g_*LcTe9(E_#V}U15$KYX2JhFjAD#NowXr_z#rT*pt%`$iod9+T2uDZ&m zN^C&~jzM(Hk(`@q==qVvp^ zY_BwbFv}kGIPdeah=Ye-5?vAujO^@M)T3<0`8YBJ_jsP9n&(laZ#74Rhu{bRTRb_j zwYRpmCbHEaCZ!46x|tw3q!T91L1PY*I8Wn)bRjM~M_q={;5o`M+`N1|b?~N_eLQl9 zS)6l6inhP6tEr`l3a1G?bK7TDw(f#jg0RZkF}W8!+{)3-a2c_;JOodh6-f9g2Uf5d zS95Ma9EsuZg(6k8NMm{v$?d`BoHPZ7L0WhzQULWCne1< zbB`i28z356gSpfM&wOLk#!}F*9)#*1U?*~KCg>dK9wUA|U%b)j zaKz!TGKfgNk;9r#WHP$oFtKcwnm}t7N}#m~B5>XWf6T<1WMQx?Q9SdH%x109&^MyW z0pWE#mmZaR&Qr_t0^}b5Q+iB7<>AHg0>9PUp?WIbgn=Q6mwq-j$Op_!F-SRj%@hO_ zvfs*T_=Ua|Dk`u=v@a}EEBOirm#O2C`b3zc6Q}KD;k*#O2UBKDXPSSfbH#h8vvo}( z)T1oth2}kyf}mcqI~0_Do0j}olcrT93g6j*oY*3rLeP@8B(fj37JI~{VPl9XFE<`i zTD}co8%@766I2=2x=Z&o_hHz%b$%kacYOSXoel>c@7XvDK`!Rl`82y`cMSsL2R4Zt zGs6jVMYc(+duaEkaZ^~&arP;zu^p+xWsG$bA)Rcoqc*D99_PbHJh!DFm4I-pJ?Rwb zi_yVbcZk<_#SZK;j*MSJ6|Mp1x5Pn{SL+BC{o7X*I*6hsL{xC#`fqF3!lAiXs5-Ax=m8byN zkEd!I+&vi|cp3Nf(rd+vWhd1ku`q@_lw7FhO^en5&?O79x=6APTl=OYsOqj%xPR;N zNjA?80tgL?Q!dZ-G8O!s8LARK z!2!F6T*Ef2-l(OxjJHVQ|96EzQ04vtOT8T4?3t@M=#T2F%Lg$@&qjT2L%>RH7>#saz(GrW9=F&4uw1g-ZrCA<^l zqL$#Lhj+Sko39_EQyX?!S)EC-4ekc0OOOOj@3r)g@uI#~dC?D_#gAaOk9}i9v?)LB zLrADLmG5m_4aI^VxWnD@~!ZfLUA#Z}KOIxR!Y zex=dbHXv%qx+dtH2v(yTPJ5!0*J~lKf!>-NzV&8f-~P9$YcG(U9`Da?q`>uX4RC|Fk(a zz~)}YTt0dAeL6o7KP>w7v43G{^x#&s+mRxEq>o7(W4jj!+)eG4OIeR@rWEC2ed|cy zTFBRjkIE&*jGELeYc*3?){k9W+R9T@rR5c9uU3#0^F_E4h)Xr{RPsnECwJ^MUv0QZ zwDAu zQf}ntm8~fuEyoj{34#uA3BV>W`F@uvV7(-*)W{y<)(gm<77+NUA%dI6N7aSqw}G6_ zdBGgaN0>D8U4kREP|pS7{ZAE3PHKaFWQ^0SRi^i*+S7MoM`D?@`&pz*PDHN_S-KiY zTixU9LgMSnDYZcJiqpgBq(g?<6d&rAcEg#}Y>Xv3oQomlzdVFb<@%?B)tR1j2&-JB zBxfcdKHO|+3A-RBOd0O(_)D)+`S0Oz92Um&6g$^lJssRtP_EwEj#Xsx=`^Vy{^jzZ z?6_}Gh2{CdCoa4u*&Y!?+7R%OoAW+ZN=x;rOuZY-KNI^TE@Ocy;Eb*b%NoB#( zdoY$@`_Rhh)BWA~VYz%ip3%>) zn#VXJ{ea#1uUMMT0S+IlEh+EOn^1>?B=$26cO(WZIO^70uZZ<`BXm|GS*67=akoyW zFy8UQL6m|IosHNBs+F=hN5nUqvY_@D_)TLmDfmK9v5S9%E|vfs+<_y5ORw6@c$TlO zaeN5vnbc~)`CqmUZr z4Tj`lgYz2gLgg1NT_jfD3vD=uQYYSSG}!bjE`2@{_2Q9cyO7N+IYsh-H`6_ZTOEMC zeA-CBD~8@yf_zYB(<_Qrgr;H7tWMawGvU7!bPK z%k*Je!Q0!Cx?72HRT#21OTOB-MPKdgc?>w}cAIw?#kqDQXi2YKwI=lj<_o!+k-37& zd_<0hsHc0lYHNg$NBU4pkF&7+Y#l7tVVQf@4#@jZnrUAC&mAmJro+QRP|`rM(lv6t zp$(e6GXI9ao|@?>YLFXlTJM?aqI7SED=_7lOt>ie4bW-Y;^@}J!x)am!w15*+a%%E$PPdYuJWd+yL?s%1c9M9xU)EX zLtrv&*4?)px;|+g}puIS8+INV6X|G^VHA{H%&{(umB0ob2hA1+PSSr0+r#9-% zD-QqR-zrkoonBJ^h!=eSyR^tP&t)H2KD5qaHY(;4CBugol!}J#HrS9Lu{|Kn6Xrqyl-Y+%`BWfEt~WSQPc0t&GnAR{Gd;uN$zu#oAVA_ zWAWP3gaPkg%HW=jlH!`~jn;*0>l71?>FkQp101l@`VJPz?E3u5EX%y4^p-?{S2$wc6O>*N^?WEt#~#`lQM=H1ph4CFyfIde5F z3zK)wEIV^9`yt zS^+38isRUKwi6o@kvnjj{=S%!X^?=@rQpB$o=yobH6cOyv5K7cB*E8{)b^99)HzhC z>UI-oT*O{_Mn{wI9`jza0&}Z!(1U}YHZinQ;Z}Y%9#MIt@WIl5l*RvdM7V*gWq>xZ zPAyj~^yN->S6Ft4N)TY0+)`yBB#?cD)7~Mfp5?*|_Pq9RgE~dQY={1QpGv)nlCkn_ z&?5H#6xM$@`~P7A64s;&Tx58@t9#7FMNFJz`1FN*tp5mQsHsts6QFNkaL9++QQ zuAR5PpXf0*$9r5rjL1nrAGPmkop&N6_+}t4XP=lMd4Ak{t)o%tiFRL}#Swjr`SYHD zkBPJHwq2YwFbxj+1O#993ioJk?N=foTe$>_M=s##GE_-*jbtpkc#XFu%G4gK`>XYB z&1+Wof1mujyHnHvWw4?(s&qkhXp0bo0KSEkZhDmS*t@1#<4+UHIVbC-PIYSTA*Drc zbrO`zah7I#riO3dCYbK!GwQd_U#8_bAk06Udz{wO>U%Ui*@LDtBrCw$T0^Az!o&QP zPD@grK)F9}7H@h5pD^#NCQ7kWu9E+}sp*4{(8BKoJ#1;iP{!Ns>1uf8l$**ODB;=* zrQ%LCA+VOI%YTRD4`U0Dz7B-Q1O!~5JaTm%9<5nQPTu`Xv;F$R_UpZ@F~{;Nh1+v0 zw`GyPHmq>#K{X&%>&XXeo4EW*uD-dp`5`7^~03ZJGp2qBQ-ZcJD^hG^Sr zJBF~Qke%u0lx$V}?j1kLR)Ketp%<8c`b@();o;DZyRdvLiN12&Mix5z8v431D*(7C z=}5Bd!tANfGx@MQ?6HZc+Kco2FY0p?f)6s=MNX*#&=r!nsh{Xd#W=rg>;J71;qPGt zlR&?2?g~4{0a%A&{X<>?ujP39(;G2P%eqPz=1Klqk+c|&D%<-NVb?l$NjlgJ51cAe znATei5zO2=weLX2d%#<%33MYUEX>5!QpeHqwVh=Ebq{j19{sqre%_QvVl<(L&Yl-` ze%nWId&?&`rDsd1_j`{JnW(l1zlI4euBc0{vZaxq_XVkT$0IO7EhG45`srx{@u9Aw z=~vKljgN=A;nEV5AL`Z%s5O#+kLFHuacMn+RX&NyA7!nnJZW<0#ox32C%2E}RtL7F zWY1V2(FH&u6F<{B5%*2{ze3@S`=7X$jjO)Ldgvd$l=?vNZhBbwq<@8kmg_qZ9|BhF ziVUZn@lia|DY(?m83F9Er~bn6ir4*h+|YFDVNy@X2v9%{I!8hs_DjC_)c0z3!v(*d z7s1A8GsZ=k>U*m0Ox?KNb}PSVt_(4?qr|O8I*mDORYiq794B| zS9Qak@v3+6X+Jwv&soVy%P_cXZ`m>yW;3(M`>PZdbk(>1cjU~3`Trzm@e_8(pA8jm zC9hMf;+WvWBH}Z?+xLds)Sr?#Uw>x!CWjARzcg!bBLfSF;gnigTa#)Di|+oMMWD`B zu{JfG5vNEn@m|;%B<@pKCPCcw{JGP1lGaPg$PCb}h`IYz1!rVQ4IoQz7#vE{37c-J zHz*<(tOxpH)Rno>A*eFHj<4TWYM(u=CRw&>it6uELn;gPM=KeIKxp3ZF6{FzltkcC zAnXeBMiB_uuWsAAtN-+yn2xrHD&0N8$|D`0OWQgTy2ca<&BszO8iX+a?7J44Rof-V z;Gzl=D}SRlVy~H3y#e{Ei1^w9(|i#h#-8s{m|f!J8PgLG0i7ddgv&KF>@mApwBDWk zlfFDOsr$|t7(tyME5`Xl9DH5RFqCGnJ?U{W>w2j18h!X#zsNn!f0X>oX$!bW+IRh8 z9IoVCGkF4IDis_6x<4S$N)f-YmSo zS`5Yc8?|A3C~(6jF;i6kkI55gL8svutHqpl4ACS$jxiqsO z>o-o&P^`e4W$Fyo>z6VYYN;8@%RG#}z^-SAp|caKLO7Pm6Gk0zl}_3L@aKCY8Sqa* z)7(#avG?`Q%}#4WA$OrSn#|pm^?+7gddp*t6!$40Z~L84ua}($R0#!~Vg94eu(Zp+*&F8_y$*di{j2uFw*}rGPr2?_e}L;+ zHzod1#pggANxZvAG>FGv8hEB_&POp{I*iwyHWX0S2Z|h?RVk=usxOgVb9}hJ`x?GM@$^E;d9o28r9TSmNF}A0)W0{RI}lXiQuyT>GoQG$dP! z@@PwH`?=_sWSmgU;$97IX>Vs~Px|}T8=}*1GxqS!k9&_U@r+?1KSgnha-(c%tR1Xv zX2Q4V!^T7;M#UuL7=@VuA6R{K1N{)&(lnDpFpvDcSF})!u}6-9@e1kw_&OBqJWGV_mNMait{98!}fMdNVD$hD@p{&qDa z2mkLWXAsEZDmI>1>k4Bv7@KXvG+_x1B3ygJAkAFgxJo|n>$!ctgWFe8FlcT>-% zSywNpLL&b?N!-Sq-O@P^(g6b$ekBbPQGLh@556Qg={fdg@J^ujS6V|Q^sVP7SDtbPUh`u;VXa(ME~9l6#|5sEU8>F(+o-gM>F1> zE)!t&j|41s;mKPHvk1krJDIM|d$U!k!5h;Ec%~<|Gs`^CQL_AE)2WLAx@t;x(laXHpMYGO;?6E{>#p@Mg z?f>kcY)Sw34$Ai)4Y?SQhhENOUM~8Teo9?~rYPp`{;&0eb2t>DWkfLg->x72KjbSZ z??Gj>5CwsEZIw9TYnHl*9e&l;6eqkQ@h?z$7i5aL(WO`KRe_pJ>-Vx z`DFQ_`uda!J~00qr1_8g2=@19`-oqVMndBkq*;sawCbVkz$dc0iE2D5HcJ?ST(M_X+b2+M;2h-_d{H<_8GBRFwXid`&3IbWWb(lcIU7`0G>9p7_VWV( zoW}+zb^!|m4EC=zC@W8n1^EYNBA%Cp*0b#4Qj?WGPL75K6HGXZO-~)t2|%nRAYF-n zF;-8T%M{EsXBoUdtafkq&z)7v<2dv#3FnXM!cWk-gM~$hjCtmPntJ4DdxzG<78$0Y Z06l$b6(K(Y+HAR;-+u)rvY@8?e*pPz%@Y6s literal 0 HcmV?d00001 diff --git a/pyat/machine_data/noringparam.mat b/pyat/machine_data/noringparam.mat new file mode 100644 index 0000000000000000000000000000000000000000..40185d4700fec78b234b7906c20d65ac0de55aab GIT binary patch literal 9008 zcma)>2T)V{+VzvrYv@GpRn{Fy!f{x`qlfGBzhN6G2r&^@^FbQ%g8UO7tIb z6_fU6=a&5SX{t8M*=O@9%$PPRVMc{erNYmUR*9j2WHUBl2I{x{I8R9n8nN`@xG1$Z zW=d}`O@z&?YQv!+*yD5Hc5inH0f;8o9$CFrXHv(KPlXj__GKfeM=rGa#(Fodg4^)y zLqRdM(r=k-mc4cAxSoo4k<)}SI{DPGsQJ*5Xyj*V_y}5y7naGax(Gy$)m)V0dPm{U zU&3(j@l)ZXh3INS;?#PO=9qMzIFG8rO;B6-(JJvuNPx(p8`t5gX;FMfe!}^R<`ThqDA5}<)EX% z4aRRgTGgK;?R|kRLWDtBk{_#d=jfDK%6o;p?C}u!-$N zWZ0oQxV7)vkDAzyn$(z!Jg`~vtl6c^G2aao^E?R&W%2uz8+1oGf7eNWq%AmEv{CH- z8;Q8R^Y46XXKlJ-$jObLO1mAy1?!17Y>8Xy0(s9{{QBOmPn@f9>r`;!A*PQjs z!^NoaRNIoHpZK#sOwHYiv;-a;#e42cu2`BE&1)pZ_*QR7r-XxpxDibiYo=;Y*Ym6K zsTXg_T!MD{>3buDSfA>c!;_C1*V)VtXv7a}UOhMWwMyD_o~3?u?gw@MH#mW8Msm;Y z?nj#E8c(a~K8PIpQ`TaYpNE|6c&$$R+_ADKudcTSFyQq|N^3TZ`sz}^)Q=}~6_RU=G5NJaCiSUm^ z_=zmlY0MP|TQRXOlYQ86@jG`+6>p(_TwMPAczI-gZZ5K_p|?{ZVq5Fhk3O7_SbuMM z)$6+njW36;`%Wv2&e}j3byyNBFQ5Ctwih z6B!0d+epLAoln80p0mg@Pf5+TKxetprI=atG?AoLBPElm68)337Q9s3HZB^R91MSd za!IQZ4E|A(UmF|UEFdEmZ4vm~nbi~J{^JWkw+ZTPTDBaN>-M>CNe7N% z#vf!Dm=3pV6*LTNk!d@_y5V1%O%Mtc8&=Y5h@Dg(~)v=s;u8R%4e;AnG>X?+Y$oy%+h$OfH}%=<#EOPIoCpW3UVE&Tkb9 zE)G%^#H4g_*LcTe9(E_#V}U15$KYX2JhFjAD#NowXr_z#rT*pt%`$iod9+T2uDZ&m zN^C&~jzM(Hk(`@q==qVvp^ zY_BwbFv}kGIPdeah=Ye-5?vAujO^@M)T3<0`8YBJ_jsP9n&(laZ#74Rhu{bRTRb_j zwYRpmCbHEaCZ!46x|tw3q!T91L1PY*I8Wn)bRjM~M_q={;5o`M+`N1|b?~N_eLQl9 zS)6l6inhP6tEr`l3a1G?bK7TDw(f#jg0RZkF}W8!+{)3-a2c_;JOodh6-f9g2Uf5d zS95Ma9EsuZg(6k8NMm{v$?d`BoHPZ7L0WhzQULWCne1< zbB`i28z356gSpfM&wOLk#!}F*9)#*1U?*~KCg>dK9wUA|U%b)j zaKz!TGKfgNk;9r#WHP$oFtKcwnm}t7N}#m~B5>XWf6T<1WMQx?Q9SdH%x109&^MyW z0pWE#mmZaR&Qr_t0^}b5Q+iB7<>AHg0>9PUp?WIbgn=Q6mwq-j$Op_!F-SRj%@hO_ zvfs*T_=Ua|Dk`u=v@a}EEBOirm#O2C`b3zc6Q}KD;k*#O2UBKDXPSSfbH#h8vvo}( z)T1oth2}kyf}mcqI~0_Do0j}olcrT93g6j*oY*3rLeP@8B(fj37JI~{VPl9XFE<`i zTD}co8%@766I2=2x=Z&o_hHz%b$%kacYOSXoel>c@7XvDK`!Rl`82y`cMSsL2R4Zt zGs6jVMYc(+duaEkaZ^~&arP;zu^p+xWsG$bA)Rcoqc*D99_PbHJh!DFm4I-pJ?Rwb zi_yVbcZk<_#SZK;j*MSJ6|Mp1x5Pn{SL+BC{o7X*I*6hsL{xC#`fqF3!lAiXs5-Ax=m8byN zkEd!I+&vi|cp3Nf(rd+vWhd1ku`q@_lw7FhO^en5&?O79x=6APTl=OYsOqj%xPR;N zNjA?80tgL?Q!dZ-G8O!s8LARK z!2!F6T*Ef2-l(OxjJHVQ|96EzQ04vtOT8T4?3t@M=#T2F%Lg$@&qjT2L%>RH7>#saz(GrW9=F&4uw1g-ZrCA<^l zqL$#Lhj+Sko39_EQyX?!S)EC-4ekc0OOOOj@3r)g@uI#~dC?D_#gAaOk9}i9v?)LB zLrADLmG5m_4aI^VxWnD@~!ZfLUA#Z}KOIxR!Y zex=dbHXv%qx+dtH2v(yTPJ5!0*J~lKf!>-NzV&8f-~P9$YcG(U9`Da?q`>uX4RC|Fk(a zz~)}YTt0dAeL6o7KP>w7v43G{^x#&s+mRxEq>o7(W4jj!+)eG4OIeR@rWEC2ed|cy zTFBRjkIE&*jGELeYc*3?){k9W+R9T@rR5c9uU3#0^F_E4h)Xr{RPsnECwJ^MUv0QZ zwDAu zQf}ntm8~fuEyoj{34#uA3BV>W`F@uvV7(-*)W{y<)(gm<77+NUA%dI6N7aSqw}G6_ zdBGgaN0>D8U4kREP|pS7{ZAE3PHKaFWQ^0SRi^i*+S7MoM`D?@`&pz*PDHN_S-KiY zTixU9LgMSnDYZcJiqpgBq(g?<6d&rAcEg#}Y>Xv3oQomlzdVFb<@%?B)tR1j2&-JB zBxfcdKHO|+3A-RBOd0O(_)D)+`S0Oz92Um&6g$^lJssRtP_EwEj#Xsx=`^Vy{^jzZ z?6_}Gh2{CdCoa4u*&Y!?+7R%OoAW+ZN=x;rOuZY-KNI^TE@Ocy;Eb*b%NoB#( zdoY$@`_Rhh)BWA~VYz%ip3%>) zn#VXJ{ea#1uUMMT0S+IlEh+EOn^1>?B=$26cO(WZIO^70uZZ<`BXm|GS*67=akoyW zFy8UQL6m|IosHNBs+F=hN5nUqvY_@D_)TLmDfmK9v5S9%E|vfs+<_y5ORw6@c$TlO zaeN5vnbc~)`CqmUZr z4Tj`lgYz2gLgg1NT_jfD3vD=uQYYSSG}!bjE`2@{_2Q9cyO7N+IYsh-H`6_ZTOEMC zeA-CBD~8@yf_zYB(<_Qrgr;H7tWMawGvU7!bPK z%k*Je!Q0!Cx?72HRT#21OTOB-MPKdgc?>w}cAIw?#kqDQXi2YKwI=lj<_o!+k-37& zd_<0hsHc0lYHNg$NBU4pkF&7+Y#l7tVVQf@4#@jZnrUAC&mAmJro+QRP|`rM(lv6t zp$(e6GXI9ao|@?>YLFXlTJM?aqI7SED=_7lOt>ie4bW-Y;^@}J!x)am!w15*+a%%E$PPdYuJWd+yL?s%1c9M9xU)EX zLtrv&*4?)px;|+g}puIS8+INV6X|G^VHA{H%&{(umB0ob2hA1+PSSr0+r#9-% zD-QqR-zrkoonBJ^h!=eSyR^tP&t)H2KD5qaHY(;4CBugol!}J#HrS9Lu{|Kn6Xrqyl-Y+%`BWfEt~WSQPc0t&GnAR{Gd;uN$zu#oAVA_ zWAWP3gaPkg%HW=jlH!`~jn;*0>l71?>FkQp101l@`VJPz?E3u5EX%y4^p-?{S2$wc6O>*N^?WEt#~#`lQM=H1ph4CFyfIde5F z3zK)wEIV^9`yt zS^+38isRUKwi6o@kvnjj{=S%!X^?=@rQpB$o=yobH6cOyv5K7cB*E8{)b^99)HzhC z>UI-oT*O{_Mn{wI9`jza0&}Z!(1U}YHZinQ;Z}Y%9#MIt@WIl5l*RvdM7V*gWq>xZ zPAyj~^yN->S6Ft4N)TY0+)`yBB#?cD)7~Mfp5?*|_Pq9RgE~dQY={1QpGv)nlCkn_ z&?5H#6xM$@`~P7A64s;&Tx58@t9#7FMNFJz`1FN*tp5mQsHsts6QFNkaL9++QQ zuAR5PpXf0*$9r5rjL1nrAGPmkop&N6_+}t4XP=lMd4Ak{t)o%tiFRL}#Swjr`SYHD zkBPJHwq2YwFbxj+1O#993ioJk?N=foTe$>_M=s##GE_-*jbtpkc#XFu%G4gK`>XYB z&1+Wof1mujyHnHvWw4?(s&qkhXp0bo0KSEkZhDmS*t@1#<4+UHIVbC-PIYSTA*Drc zbrO`zah7I#riO3dCYbK!GwQd_U#8_bAk06Udz{wO>U%Ui*@LDtBrCw$T0^Az!o&QP zPD@grK)F9}7H@h5pD^#NCQ7kWu9E+}sp*4{(8BKoJ#1;iP{!Ns>1uf8l$**ODB;=* zrQ%LCA+VOI%YTRD4`U0Dz7B-Q1O!~5JaTm%9<5nQPTu`Xv;F$R_UpZ@F~{;Nh1+v0 zw`GyPHmq>#K{X&%>&XXeo4EW*uD-dp`5`7^~03ZJGp2qBQ-ZcJD^hG^Sr zJBF~Qke%u0lx$V}?j1kLR)Ketp%<8c`b@();o;DZyRdvLiN12&Mix5z8v431D*(7C z=}5Bd!tANfGx@MQ?6HZc+Kco2FY0p?f)6s=MNX*#&=r!nsh{Xd#W=rg>;J71;qPGt zlR&?2?g~4{0a%A&{X<>?ujP39(;G2P%eqPz=1Klqk+c|&D%<-NVb?l$NjlgJ51cAe znATei5zO2=weLX2d%#<%33MYUEX>5!QpeHqwVh=Ebq{j19{sqre%_QvVl<(L&Yl-` ze%nWId&?&`rDsd1_j`{JnW(l1zlI4euBc0{vZaxq_XVkT$0IO7EhG45`srx{@u9Aw z=~vKljgN=A;nEV5AL`Z%s5O#+kLFHuacMn+RX&NyA7!nnJZW<0#ox32C%2E}RtL7F zWY1V2(FH&u6F<{B5%*2{ze3@S`=7X$jjO)Ldgvd$l=?vNZhBbwq<@8kmg_qZ9|BhF ziVUZn@lia|DY(?m83F9Er~bn6ir4*h+|YFDVNy@X2v9%{I!8hs_DjC_)c0z3!v(*d z7s1A8GsZ=k>U*m0Ox?KNb}PSVt_(4?qr|O8I*mDORYiq794B| zS9Qak@v3+6X+Jwv&soVy%P_cXZ`m>yW;3(M`>PZdbk(>1cjU~3`Trzm@e_8(pA8jm zC9hMf;+WvWBH}Z?+xLds)Sr?#Uw>x!CWjARzcg!bBLfSF;gnigTa#)Di|+oMMWD`B zu{JfG5vNEn@m|;%B<@pKCPCcw{JGP1lGaPg$PCb}h`IYz1!rVQ4IoQz7#vE{37c-J zHz*<(tOxpH)Rno>A*eFHj<4TWYM(u=CRw&>it6uELn;gPM=KeIKxp3ZF6{FzltkcC zAnXeBMiB_uuWsAAtN-+yn2xrHD&0N8$|D`0OWQgTy2ca<&BszO8iX+a?7J44Rof-V z;Gzl=D}SRlVy~H3y#e{Ei1^w9(|i#h#-8s{m|f!J8PgLG0i7ddgv&KF>@mApwBDWk zlfFDOsr$|t7(tyME5`Xl9DH5RFqCGnJ?U{W>w2j18h!X#zsNn!f0X>oX$!bW+IRh8 z9IoVCGkF4IDis_6x<4S$N)f-YmSo zS`5Yc8?|A3C~(6jF;i6kkI55gL8svutHqpl4ACS$jxiqsO z>o-o&P^`e4W$Fyo>z6VYYN;8@%RG#}z^-SAp|caKLO7Pm6Gk0zl}_3L@aKCY8Sqa* z)7(#avG?`Q%}#4WA$OrSn#|pm^?+7gddp*t6!$40Z~L84ua}($R0#!~Vg94eu(Zp+*&F8_y$*di{j2uFw*}rGPr2?_e}L;+ zHzod1#pggANxZvAG>FGv8hEB_&POp{I*iwyHWX0S2Z|h?RVk=usxOgVb9}hJ`x?GM@$^E;d9o28r9TSmNF}A0)W0{RI}lXiQuyT>GoQG$dP! z@@PwH`?=_sWSmgU;$97IX>Vs~Px|}T8=}*1GxqS!k9&_U@r+?1KSgnha-(c%tR1Xv zX2Q4V!^T7;M#UuL7=@VuA6R{K1N{)&(lnDpFpvDcSF})!u}6-9@e1kw_&OBqJWGV_mNMait{98!}fMdNVD$hD@p{&qDa z2mkLWXAsEZDmI>1>k4Bv7@KXvG+_x1B3ygJAkAFgxJo|n>$!ctgWFe8FlcT>-% zSywNpLL&b?N!-Sq-O@P^(g6b$ekBbPQGLh@556Qg={fdg@J^ujS6V|Q^sVP7SDtbPUh`u;VXa(ME~9l6#|5sEU8>F(+o-gM>F1> zE)!t&j|41s;mKPHvk1krJDIM|d$U!k!5h;Ec%~<|Gs`^CQL_AE)2WLAx@t;x(laXHpMYGO;?6E{>#p@Mg z?f>kcY)Sw34$Ai)4Y?SQhhENOUM~8Teo9?~rYPp`{;&0eb2t>DWkfLg->x72KjbSZ z??Gj>5CwsEZIw9TYnHl*9e&l;6eqkQ@h?z$7i5aL(WO`KRe_pJ>-Vx z`DFQ_`uda!J~00qr1_8g2=@19`-oqVMndBkq*;sawCbVkz$dc0iE2D5HcJ?ST(M_X+b2+M;2h-_d{H<_8GBRFwXid`&3IbWWb(lcIU7`0G>9p7_VWV( zoW}+zb^!|m4EC=zC@W8n1^EYNBA%Cp*0b#4Qj?WGPL75K6HGXZO-~)t2|%nRAYF-n zF;-8T%M{EsXBoUdtafkq&z)7v<2dv#3FnXM!cWk-gM~$hjCtmPntJ4DdxzG<78$0Y Z06l$b6(K(Y+HAR;-+u)rvY@8?e*pPz%@Y6s literal 0 HcmV?d00001 From 233340047c200766d8e38e7130951d43907a0703 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Thu, 15 Aug 2024 21:05:11 +0200 Subject: [PATCH 09/18] fix NaN in cell_harmnumber --- pyat/at/load/matfile.py | 4 +++- pyat/at/load/utils.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/pyat/at/load/matfile.py b/pyat/at/load/matfile.py index 66acea2e3..721c3af5d 100644 --- a/pyat/at/load/matfile.py +++ b/pyat/at/load/matfile.py @@ -11,6 +11,7 @@ from os.path import abspath, basename, splitext from typing import Any from collections.abc import Sequence, Generator +from math import isfinite from warnings import warn import numpy as np @@ -147,7 +148,8 @@ def ringparam_filter( for k, v in elem.items(): k2 = _m2p.get(k, k) if k2 is not None: - params.setdefault(k2, v) + if k2 != "cell_harmnumber" or isfinite(v): + params.setdefault(k2, v) if keep_all: pars = vars(elem).copy() name = pars.pop("FamName") diff --git a/pyat/at/load/utils.py b/pyat/at/load/utils.py index 55d1a7b5e..7aa6eab42 100644 --- a/pyat/at/load/utils.py +++ b/pyat/at/load/utils.py @@ -90,7 +90,7 @@ class RingParam(elt.Element): Energy=float, Periodicity=int, Particle=_particle, - cell_harmnumber=int, + cell_harmnumber=float, ) # noinspection PyPep8Naming From 4989b8778387228236fbb462198c7b82fe8926c2 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Fri, 16 Aug 2024 20:35:36 +0200 Subject: [PATCH 10/18] findelemm44 --- atmat/atphysics/LinearOptics/findelemm44.m | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/atmat/atphysics/LinearOptics/findelemm44.m b/atmat/atphysics/LinearOptics/findelemm44.m index 3e5976d03..d67278714 100644 --- a/atmat/atphysics/LinearOptics/findelemm44.m +++ b/atmat/atphysics/LinearOptics/findelemm44.m @@ -12,10 +12,19 @@ % The transverse matrix is momentum-dependent, % the 5-th component of ORBITIN is used as the DP value % +% M66=FINDELEMM44(...,'Energy',ENERGY) +% Use ENERGY and ignore the 'Energy' field of elements +% +% M66=FINDELEMM44(...,'Particle',PARTICLE) +% Use PARTICLE (default is relativistic) +% +% % See also FINDELEMM66 [XYStep,varargs]=getoption(varargin,'XYStep'); [R0,varargs]=getoption(varargs,'orbit',zeros(6,1)); +[energy,varargs]=getoption(varargs,'Energy',0.0); +[particle,varargs]=getoption(varargs,'Particle',[]); [MethodName,R0]=getargs(varargs,ELEM.PassMethod,R0); % Build a diagonal matrix of initial conditions @@ -25,6 +34,6 @@ % Add to the orbit_in RIN = R0 + [D4 -D4]; % Propagate through the element -ROUT = feval(MethodName,ELEM,RIN); +ROUT=elempass(ELEM,RIN,'PassMethod',MethodName,'Energy',energy,'Particle',particle); % Calculate numerical derivative M44 = ((ROUT(1:4,1:4)-ROUT(1:4,5:8))./scaling); From b8835388773ef74977134f00ab3e330543c53fcc Mon Sep 17 00:00:00 2001 From: Nicola Carmignani Date: Thu, 29 Aug 2024 15:43:35 +0200 Subject: [PATCH 11/18] atEnergy for matlab checks that energy is not 0 --- atintegrators/atelem.c | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 9ef2de4c6..84253fcf2 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -71,9 +71,19 @@ typedef mxArray atElem; #define atError(...) mexErrMsgIdAndTxt("AT:PassError", __VA_ARGS__) #define atWarning(...) mexWarnMsgIdAndTxt("AT:PassWarning", __VA_ARGS__) #define atPrintf(...) mexPrintf(__VA_ARGS__) -#define atEnergy(ringenergy,elemenergy) ((ringenergy)==0.0 ? (elemenergy) : (ringenergy)) #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."); +} + static mxArray *get_field(const mxArray *pm, const char *fieldname) { mxArray *field; From 24eaef4db76b982a98408565be01d5b22264a110 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Fri, 30 Aug 2024 16:53:44 +0200 Subject: [PATCH 12/18] update GWigSymplecticRadPass.c --- atintegrators/GWigSymplecticRadPass.c | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/atintegrators/GWigSymplecticRadPass.c b/atintegrators/GWigSymplecticRadPass.c index 6785380bc..af306484a 100644 --- a/atintegrators/GWigSymplecticRadPass.c +++ b/atintegrators/GWigSymplecticRadPass.c @@ -6,6 +6,10 @@ *--------------------------------------------------------------------------- * Modification Log: * ----------------- + * .03 2024-05-06 J. Arenillas, ALBA, jarenillas@axt.email + * Adding rotations and translations to wiggler. + * Bug fix in wiggler initialisation. + * Energy parameter bug fix. * .02 2003-06-18 J. Li * Cleanup the code * @@ -123,7 +127,7 @@ void GWigInit(struct gwigR *Wig,double design_energy, double Ltot, double Lw, #define second 2 #define fourth 4 -void GWigSymplecticRadPass(double *r,double Energy, double Ltot, double Lw, +void GWigSymplecticRadPass(double *r, double Energy, double Ltot, double Lw, double Bmax, int Nstep, int Nmeth, int NHharm, int NVharm, double *By, double *Bx, double *T1, double *T2, double *R1, double *R2, int num_particles) @@ -141,11 +145,13 @@ void GWigSymplecticRadPass(double *r,double Energy, double Ltot, double Lw, zEndPointV[0] = 0; zEndPointV[1] = Ltot; - GWigInit(&Wig, Energy, Ltot, Lw, Bmax, Nstep, Nmeth, NHharm, NVharm,0, 0, zEndPointH, zEndPointV, By, Bx, T1, T2, R1, R2); - for(c = 0;c Date: Fri, 30 Aug 2024 17:13:39 +0200 Subject: [PATCH 13/18] update ohmienvelope for wigglers --- atmat/atphysics/Radiation/ohmienvelope.m | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/atmat/atphysics/Radiation/ohmienvelope.m b/atmat/atphysics/Radiation/ohmienvelope.m index 02bd7e74a..854828e31 100644 --- a/atmat/atphysics/Radiation/ohmienvelope.m +++ b/atmat/atphysics/Radiation/ohmienvelope.m @@ -83,7 +83,11 @@ if nout>=1, varargout{1}=mring; end function diff=diffmatrix(elem,orbit) - diff=findmpoleraddiffmatrix(elem, orbit, energy); + if isfield(elem,'Bmax') + diff=FDW(elem,orbit,energy); % For wigglers + else + diff=findmpoleraddiffmatrix(elem, orbit, energy); % For other elements + end end function btx=cumulb(elem,orbit,b) From ca508d874da34a7406b6aafe3402ca3b9b82d5d4 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Sun, 1 Sep 2024 10:43:39 +0200 Subject: [PATCH 14/18] compile FDW and integrate wiggler diffusion matrix --- atintegrators/atelem.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index 84253fcf2..da14ef3d9 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -80,8 +80,10 @@ double atEnergy(double ringenergy, double elemenergy) else if (elemenergy!=0.0) return elemenergy; - else + else { atError("Energy not defined."); + return 0.0; /* Never reached but makes the compiler happy */ + } } static mxArray *get_field(const mxArray *pm, const char *fieldname) From 80c34be7f24c413f50595c7b1f1751111f9b7224 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Mon, 2 Sep 2024 12:02:33 +0200 Subject: [PATCH 15/18] factorise the computation of the diffusion matrix used in ohmienvelope.m and quantumDiff.m --- atmat/atphysics/Radiation/atdiffmat.m | 51 ++++++++++++++++++++ atmat/atphysics/Radiation/ohmienvelope.m | 39 ++------------- atmat/atphysics/Radiation/quantumDiff.m | 61 ++++++------------------ 3 files changed, 70 insertions(+), 81 deletions(-) create mode 100644 atmat/atphysics/Radiation/atdiffmat.m diff --git a/atmat/atphysics/Radiation/atdiffmat.m b/atmat/atphysics/Radiation/atdiffmat.m new file mode 100644 index 000000000..0dc61611c --- /dev/null +++ b/atmat/atphysics/Radiation/atdiffmat.m @@ -0,0 +1,51 @@ +function [BCUM,Batbeg] = atdiffmat(ring,energy,varargin) +%quantumDiff Compute the radiation-diffusion matrix +% +%[BCUM,BS]=ATDIFFMAT(RING,ENERGY) +% RING: Closed ring AT structure, containing radiative elements and +% RF cavity. Radiative elements are identified by a +% PassMethod ending with 'RadPass'. +% ENERGY: Lattice energy +% +% BCUM: Cumulated diffusion matrix +% BS: Cumulative diffusion matrix at the beginning of each element +% +%[BCUM,BS]=ATDIFFMAT(...,'orbit',ORBITIN) +% ORBITIN: Initial 6-D closed orbit. +% In this mode, RING may be a section of a ring. + +NumElements=length(ring); + +[orb0,varargs]=getoption(varargin, 'orbit', []); %#ok +if isempty(orb0) + orbit=num2cell(findorbit6(ring, 1:NumElements),1)'; +else + orbit=num2cell(linepass(ring, orb0, 1:NumElements),1)'; +end + +zr=zeros(6,6); +% Calculate cumulative Radiation-Diffusion matrix for the ring +BCUM = zeros(6,6); +% Batbeg{i} is the cumulative diffusion matrix from +% 0 to the beginning of the i-th element +Batbeg=[zr;cellfun(@cumulb,ring,orbit,'UniformOutput',false)]; + + function btx=cumulb(elem,orbit) + if endsWith(elem.PassMethod,'RadPass') + if isfield(elem,'Bmax') + b=FDW(elem,orbit,energy); % For wigglers + else + b=findmpoleraddiffmatrix(elem, orbit, energy); % For other elements + end + else + b=zr; + end + % Calculate 6-by-6 linear transfer matrix in each element + % near the equilibrium orbit + m=findelemm66(elem,elem.PassMethod,orbit); + % Cumulative diffusion matrix of the entire ring + BCUM = m*BCUM*m' + b; + btx=BCUM; + end +end + diff --git a/atmat/atphysics/Radiation/ohmienvelope.m b/atmat/atphysics/Radiation/ohmienvelope.m index 854828e31..e1d5f4e18 100644 --- a/atmat/atphysics/Radiation/ohmienvelope.m +++ b/atmat/atphysics/Radiation/ohmienvelope.m @@ -1,4 +1,4 @@ -function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts,energy) +function [envelope, rmsdp, rmsbl, varargout] = ohmienvelope(ring,radindex,refpts,energy) %#ok %OHMIENVELOPE calculates equilibrium beam envelope in a % circular accelerator using Ohmi's beam envelope formalism [1]. % [1] K.Ohmi et al. Phys.Rev.E. Vol.49. (1994) @@ -36,22 +36,9 @@ if nargin<4, energy=atGetRingProperties(ring, 'Energy'); end if nargin<3, refpts=1; end -[mring, ms, orbit] = findm66(ring,1:NumElements+1); -mt=squeeze(num2cell(ms,[1 2])); -orb=num2cell(orbit,1)'; - -zr={zeros(6,6)}; -B=zr(ones(NumElements,1)); % B{i} is the diffusion matrix of the i-th element - -% calculate Radiation-Diffusion matrix B for elements with radiation -B(radindex)=cellfun(@diffmatrix,... - ring(radindex),orb(radindex),'UniformOutput',false); - -% Calculate cumulative Radiation-Diffusion matrix for the ring -BCUM = zeros(6,6); -% Batbeg{i} is the cumulative diffusion matrix from -% 0 to the beginning of the i-th element -Batbeg=[zr;cellfun(@cumulb,ring,orb(1:end-1),B,'UniformOutput',false)]; +orb0=findorbit(ring); +[BCUM,Batbeg]=atdiffmat(ring,energy,'orbit',orb0); +[mring, ms, orbit] = findm66(ring,1:NumElements+1,'orbit',orb0); % ------------------------------------------------------------------------ % Equation for the moment matrix R is @@ -73,6 +60,7 @@ rmsdp = sqrt(R(5,5)); % R.M.S. energy spread rmsbl = sqrt(R(6,6)); % R.M.S. bunch lenght +mt=squeeze(num2cell(ms,[1 2])); [rr,tt,ss]=cellfun(@propag,mt(refpts),Batbeg(refpts),'UniformOutput',false); envelope=struct('R',rr,'Sigma',ss,'Tilt',tt); @@ -82,23 +70,6 @@ if nout>=2, varargout{2}=ms(:,:,refpts); end if nout>=1, varargout{1}=mring; end - function diff=diffmatrix(elem,orbit) - if isfield(elem,'Bmax') - diff=FDW(elem,orbit,energy); % For wigglers - else - diff=findmpoleraddiffmatrix(elem, orbit, energy); % For other elements - end - end - - function btx=cumulb(elem,orbit,b) - % Calculate 6-by-6 linear transfer matrix in each element - % near the equilibrium orbit - m=findelemm66(elem,elem.PassMethod,'orbit',orbit,'Energy',energy); - % Cumulative diffusion matrix of the entire ring - BCUM = m*BCUM*m' + b; - btx=BCUM; - end - function [r,tilt,sigma]=propag(m,cumb) r=m*R*m'+cumb; [u,dr] = eig(r([1 3],[1 3])); diff --git a/atmat/atphysics/Radiation/quantumDiff.m b/atmat/atphysics/Radiation/quantumDiff.m index 23aa0233e..3b37f4ebc 100644 --- a/atmat/atphysics/Radiation/quantumDiff.m +++ b/atmat/atphysics/Radiation/quantumDiff.m @@ -1,65 +1,32 @@ function DiffMat = quantumDiff(elems, varargin) -%quantumDiff Compute the radiation-diffusion matrix +%QUANTUMDIFF Compute the radiation-diffusion matrix % %DIFFMAT=QUANTUMDIFF(RING) % RING: Closed ring AT structure, containing radiative elements and % RF cavity. Radiative elements are identified by a % PassMethod ending with 'RadPass'. % -%DIFFMAT=QUANTUMDIFF(RING,RADINDEX) -% RADINDEX: Indices of elements where diffusion occurs, typically dipoles -% and quadrupoles. -% %DIFFMAT=QUANTUMDIFF(LINE,RADINDEX,ORBITIN) (Deprecated syntax) %DIFFMAT=QUANTUMDIFF(...,'orbit',ORBITIN) +% RADINDEX: Ignored % ORBITIN: Initial 6-D closed orbit. % In this mode, LINE may be a section of a ring. +% +%DIFFMAT=QUANTUMDIFF(...,'Energy',ENERGY) +% ENERGY: Lattice energy -NumElements=length(elems); - -[orb0,varargs]=getoption(varargin, 'orbit', []); -if length(varargs) >= 2 % QUANTUMDIFF(RING,RADINDEX,ORBITIN) - orb0 = varargs{2}; -end -if length(varargs) >= 1 % QUANTUMDIFF(RING,RADINDEX,...) - radindex=varargs{1}; -else % QUANTUMDIFF(RING) - radindex=atgetcells(elems,'PassMethod',@(elem, pass) endsWith(pass,'RadPass')); -end +[orb0,varargs]=getoption(varargin, 'orbit',[]); +[energy,varargs]=getoption(varargs,'Energy',[]); +[radindex,orb0,varargs]=getargs(varargs,[],orb0); %#ok -%[mring, ms, orbit] = findm66(ring,1:NumElements+1); if isempty(orb0) - orb=num2cell(findorbit6(elems, 1:NumElements),1)'; -else - orb=num2cell(linepass(elems, orb0, 1:NumElements),1)'; + orb0=findorbit6(elems); +end +if isempty(energy) + energy=atGetRingProperties(elems,'Energy'); end -zr={zeros(6,6)}; -B=zr(ones(NumElements,1)); % B{i} is the diffusion matrix of the i-th element - -% calculate Radiation-Diffusion matrix B for elements with radiation -B(radindex)=cellfun(@findmpoleraddiffmatrix,... - elems(radindex),orb(radindex),'UniformOutput',false); - -% Calculate cumulative Radiation-Diffusion matrix for the ring -BCUM = zeros(6,6); -% Batbeg{i} is the cumulative diffusion matrix from -% 0 to the beginning of the i-th element -Batbeg=[zr;cellfun(@cumulb,elems,orb,B,'UniformOutput',false)]; %#ok - -DiffCum = BCUM; - -DiffMat=(DiffCum+DiffCum')/2; - -%Lmat=chol((DiffCum+DiffCum')/2); +BCUM=atdiffmat(elems,energy,'orbit',orb0); - function btx=cumulb(elem,orbit,b) - % Calculate 6-by-6 linear transfer matrix in each element - % near the equilibrium orbit - m=findelemm66(elem,elem.PassMethod,orbit); - % Cumulative diffusion matrix of the entire ring - BCUM = m*BCUM*m' + b; - btx=BCUM; - end +DiffMat=(BCUM+BCUM')/2; end - From 29f75cb472cd71444d76904ad5ab1c1b7fc48157 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Mon, 2 Sep 2024 15:38:17 +0200 Subject: [PATCH 16/18] merged from mster --- atmat/attests/pytests.m | 61 +++++++++++++++++------------------------ 1 file changed, 25 insertions(+), 36 deletions(-) diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m index ed8c9d547..28bb07665 100644 --- a/atmat/attests/pytests.m +++ b/atmat/attests/pytests.m @@ -1,6 +1,4 @@ classdef pytests < matlab.unittest.TestCase - % ">> run(pytests)" to run all tests - % ">> run(pytests, 'orbit4')" to run a the single orbit4 method properties(Constant) mlist=[... @@ -166,25 +164,23 @@ function offmomdf(testCase, dp) function tunechrom4(testCase,lat,dp) % test on and off-momentum tunes of 4D lattices lattice=testCase.ring4.(lat); - periodicity = atGetRingProperties(lattice.m,'Periodicity'); [mtune,mchrom]=tunechrom(lattice.m,'get_chrom',dp=dp); ptune=double(lattice.p.get_tune(pyargs(dp=dp))); pchrom=double(lattice.p.get_chrom(pyargs(dp=dp))); - testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=2.e-9); - testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=3.e-4); + testCase.verifyEqual(mtune,ptune,AbsTol=2.e-9); + testCase.verifyEqual(mchrom,pchrom,AbsTol=3.e-4); end function tunechrom6(testCase,lat2,dp) % test on and off-momentum tunes of 6D lattices lattice=testCase.ring6.(lat2); - periodicity = atGetRingProperties(lattice.m,'Periodicity'); mlat=atsetcavity(lattice.m,frequency='nominal',dp=dp); plat=lattice.p.set_rf_frequency(pyargs(dp=dp,copy=true)); [mtune,mchrom]=tunechrom(mlat,'get_chrom'); ptune=double(plat.get_tune()); pchrom=double(plat.get_chrom()); - testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=1.e-9); - testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=2.e-4); + testCase.verifyEqual(mtune,ptune,AbsTol=1.e-9); + testCase.verifyEqual(mchrom,pchrom,AbsTol=2.e-4); end function linopt1(testCase,dp) @@ -209,6 +205,7 @@ function linopt1(testCase,dp) end function avlin1(testCase, lat, dp) + % Check average beta, dispersion, phase advance lattice=testCase.ring4.(lat); mrefs=true(1,length(lattice.m)); prefs=py.numpy.array(mrefs); @@ -243,34 +240,26 @@ function ringparameters(testCase,lat2) testCase.verifyEqual(mmcf,lattice.p.mcf,RelTol=1.E-8); end - function fastring(testCase,lat2) - lattice=testCase.ring4.(lat2); - [rm, rmrad]=atfastring(lattice.m); - rpy=cell(py.at.fast_ring(lattice.p)); - rp=cell(rpy{1}); - rprad=cell(rpy{2}); - checkattr(rm, rp); - checkattr(rmrad,rprad); - - function checkattr(rm, rp) - testCase.verifyEqual(rm{2}.Frequency, rp{1}.Frequency, RelTol=1.0e-20); - testCase.verifyEqual(rm{2}.Voltage, rp{1}.Voltage, RelTol=1.0e-20); - testCase.verifyEqual(rm{3}.I2, rp{2}.I2, RelTol=1.0e-15); - testCase.verifyEqual(rm{3}.Length, rp{2}.Length, RelTol=1.0e-20); - testCase.verifyEqual(rm{3}.M66, double(rp{2}.M66), AbsTol=1.0e-7); - testCase.verifyEqual(rm{end}.A1, rp{3}.A1, RelTol=0.01); - testCase.verifyEqual(rm{end}.A2, rp{3}.A2, RelTol=0.02); - testCase.verifyEqual(rm{end}.A3, rp{3}.A3, RelTol=0.01); - testCase.verifyEqual(rm{end}.Alphax, rp{3}.Alphax, AbsTol=1.e-10); - testCase.verifyEqual(rm{end}.Alphay, rp{3}.Alphay, AbsTol=1.e-10); - testCase.verifyEqual(rm{end}.Betax, rp{3}.Betax, RelTol=1.e-10); - testCase.verifyEqual(rm{end}.Betay, rp{3}.Betay, RelTol=1.e-10); - testCase.verifyEqual(rm{end}.Qpx, rp{3}.Qpx, RelTol=1.e-8); - testCase.verifyEqual(rm{end}.Qpy, rp{3}.Qpy, RelTol=1.e-8); - if length(rm) >= 5 - testCase.verifyEqual(rm{end-1}.Lmatp, double(rp{4}.Lmatp), AbsTol=2.e-7); - end - end + function emittances(testCase, lat2) + % Check emittances, tunes and damping rates + lattice=testCase.ring6.(lat2); + % python + pdata=cell(lattice.p.ohmi_envelope()); + pemit=double(py.getattr(pdata{2},'mode_emittances')); + ptunes=double(py.getattr(pdata{2},'tunes')); + pdamprate=double(py.getattr(pdata{2},'damping_rates')); + %matlab + [mdata,~,~,m]=ohmienvelope(lattice.m); + jmt=jmat(3); + aa=amat(m); + nn=-aa'*jmt*mdata.R*jmt*aa; + memit=0.5*[nn(1,1)+nn(2,2) nn(3,3)+nn(4,4) nn(5,5)+nn(6,6)]; + [mtunes,mdamprate]=atdampingrates(m); + %check + testCase.verifyEqual(memit,pemit,AbsTol=1.e-30,RelTol=1.e-6); + testCase.verifyEqual(mtunes,ptunes,AbsTol=1.e-10); + testCase.verifyEqual(mdamprate,pdamprate,AbsTol=1.e-10); end + end end From e235a7864c8c667871517ce68ad9513216afd13b Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Tue, 29 Oct 2024 16:05:05 +0100 Subject: [PATCH 17/18] merged pytests.m from master --- atmat/attests/pytests.m | 41 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 37 insertions(+), 4 deletions(-) diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m index 28bb07665..c4818e0a8 100644 --- a/atmat/attests/pytests.m +++ b/atmat/attests/pytests.m @@ -1,4 +1,6 @@ classdef pytests < matlab.unittest.TestCase + % ">> run(pytests)" to run all tests + % ">> run(pytests, 'orbit4')" to run a the single orbit4 method properties(Constant) mlist=[... @@ -164,23 +166,25 @@ function offmomdf(testCase, dp) function tunechrom4(testCase,lat,dp) % test on and off-momentum tunes of 4D lattices lattice=testCase.ring4.(lat); + periodicity = atGetRingProperties(lattice.m,'Periodicity'); [mtune,mchrom]=tunechrom(lattice.m,'get_chrom',dp=dp); ptune=double(lattice.p.get_tune(pyargs(dp=dp))); pchrom=double(lattice.p.get_chrom(pyargs(dp=dp))); - testCase.verifyEqual(mtune,ptune,AbsTol=2.e-9); - testCase.verifyEqual(mchrom,pchrom,AbsTol=3.e-4); + testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=2.e-9); + testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=3.e-4); end function tunechrom6(testCase,lat2,dp) % test on and off-momentum tunes of 6D lattices lattice=testCase.ring6.(lat2); + periodicity = atGetRingProperties(lattice.m,'Periodicity'); mlat=atsetcavity(lattice.m,frequency='nominal',dp=dp); plat=lattice.p.set_rf_frequency(pyargs(dp=dp,copy=true)); [mtune,mchrom]=tunechrom(mlat,'get_chrom'); ptune=double(plat.get_tune()); pchrom=double(plat.get_chrom()); - testCase.verifyEqual(mtune,ptune,AbsTol=1.e-9); - testCase.verifyEqual(mchrom,pchrom,AbsTol=2.e-4); + testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=1.e-9); + testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=2.e-4); end function linopt1(testCase,dp) @@ -240,6 +244,35 @@ function ringparameters(testCase,lat2) testCase.verifyEqual(mmcf,lattice.p.mcf,RelTol=1.E-8); end + function fastring(testCase,lat2) + lattice=testCase.ring4.(lat2); + [rm, rmrad]=atfastring(lattice.m); + rpy=cell(py.at.fast_ring(lattice.p)); + rp=cell(rpy{1}); + rprad=cell(rpy{2}); + checkattr(rm, rp); + checkattr(rmrad,rprad); + + function checkattr(rm, rp) + testCase.verifyEqual(rm{2}.Frequency, rp{1}.Frequency, RelTol=1.0e-20); + testCase.verifyEqual(rm{2}.Voltage, rp{1}.Voltage, RelTol=1.0e-20); + testCase.verifyEqual(rm{3}.I2, rp{2}.I2, RelTol=1.0e-15); + testCase.verifyEqual(rm{3}.Length, rp{2}.Length, RelTol=1.0e-20); + testCase.verifyEqual(rm{3}.M66, double(rp{2}.M66), AbsTol=1.0e-7); + testCase.verifyEqual(rm{end}.A1, rp{3}.A1, RelTol=0.01); + testCase.verifyEqual(rm{end}.A2, rp{3}.A2, RelTol=0.02); + testCase.verifyEqual(rm{end}.A3, rp{3}.A3, RelTol=0.01); + testCase.verifyEqual(rm{end}.Alphax, rp{3}.Alphax, AbsTol=1.e-10); + testCase.verifyEqual(rm{end}.Alphay, rp{3}.Alphay, AbsTol=1.e-10); + testCase.verifyEqual(rm{end}.Betax, rp{3}.Betax, RelTol=1.e-10); + testCase.verifyEqual(rm{end}.Betay, rp{3}.Betay, RelTol=1.e-10); + testCase.verifyEqual(rm{end}.Qpx, rp{3}.Qpx, RelTol=1.e-8); + testCase.verifyEqual(rm{end}.Qpy, rp{3}.Qpy, RelTol=1.e-8); + if length(rm) >= 5 + testCase.verifyEqual(rm{end-1}.Lmatp, double(rp{4}.Lmatp), AbsTol=2.e-7); + end + end + function emittances(testCase, lat2) % Check emittances, tunes and damping rates lattice=testCase.ring6.(lat2); From c790d9b9c9eb1776f8ecff5438a773f486339cc3 Mon Sep 17 00:00:00 2001 From: Laurent Farvacque Date: Tue, 29 Oct 2024 16:21:47 +0100 Subject: [PATCH 18/18] merged pytests.m from master --- atintegrators/atelem.c | 11 +++++++++++ atmat/attests/pytests.m | 3 ++- atmat/lattice/at2str.m | 2 +- 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/atintegrators/atelem.c b/atintegrators/atelem.c index da14ef3d9..04af7ccca 100755 --- a/atintegrators/atelem.c +++ b/atintegrators/atelem.c @@ -7,6 +7,7 @@ #include "atcommon.h" #include "attypes.h" +#include "atconstants.h" /*----------------------------------------------------*/ /* For the integrator code */ @@ -86,6 +87,15 @@ double atEnergy(double ringenergy, double elemenergy) } } +double atGamma(double ringenergy, double elemenergy, double rest_energy) +{ + double energy = atEnergy(ringenergy, elemenergy); + if (rest_energy == 0.0) + return 1.0E-9 * energy / __E0; + else + return energy / rest_energy; +} + static mxArray *get_field(const mxArray *pm, const char *fieldname) { mxArray *field; @@ -196,6 +206,7 @@ typedef PyObject atElem; #define atWarning(...) if (PyErr_WarnFormat(PyExc_RuntimeWarning, 0, __VA_ARGS__) != 0) return NULL #define atPrintf(...) PySys_WriteStdout(__VA_ARGS__) #define atEnergy(ringenergy,elemenergy) (ringenergy) +#define atGamma(ringenergy,elemenergy,rest_energy) ((rest_energy) == 0.0 ? 1.0E-9*(ringenergy)/__E0 : (ringenergy)/(rest_energy)) static int array_imported = 0; diff --git a/atmat/attests/pytests.m b/atmat/attests/pytests.m index c4818e0a8..567aa3471 100644 --- a/atmat/attests/pytests.m +++ b/atmat/attests/pytests.m @@ -184,7 +184,7 @@ function tunechrom6(testCase,lat2,dp) ptune=double(plat.get_tune()); pchrom=double(plat.get_chrom()); testCase.verifyEqual(mod(mtune*periodicity,1),ptune,AbsTol=1.e-9); - testCase.verifyEqual(mchrom*periodicity,pchrom,AbsTol=2.e-4); + testCase.verifyEqual(mchrom*periodicity,pchrom,RelTol=1.e-4,AbsTol=3.e-4); end function linopt1(testCase,dp) @@ -272,6 +272,7 @@ function checkattr(rm, rp) testCase.verifyEqual(rm{end-1}.Lmatp, double(rp{4}.Lmatp), AbsTol=2.e-7); end end + end function emittances(testCase, lat2) % Check emittances, tunes and damping rates diff --git a/atmat/lattice/at2str.m b/atmat/lattice/at2str.m index 74e0201a2..631ca57ab 100644 --- a/atmat/lattice/at2str.m +++ b/atmat/lattice/at2str.m @@ -79,7 +79,7 @@ [options,args]=doptions(elem,create,{'PolynomA','PolynomB'}); case 'RFCavity' create=@atrfcavity; - [options,args]=doptions(elem,create,{'Length','Voltage','Frequency','HarmNumber'}); + [options,args]=doptions(elem,create,{'Length','Voltage','Frequency','HarmNumber','Energy'}); case 'RingParam' create=@atringparam; [options,args]=doptions(elem,create,{'Energy','Periodicity'});