Fix wiggler radiation

atintegrators/gwigR.c

@@ -124,7 +124,6 @@ static void GWigPass_4th(struct gwigR *pWig, double *X)
   double dl, dl1, dl0;
   double B[3];
   double ax, ay, axpy, aypx;
-  double factorvect;
 
   Nstep = pWig->PN*(pWig->Nw);
   dl = pWig->Lw/(pWig->PN);
@@ -450,7 +449,7 @@ static void GWigB(struct gwigR *pWig, double *Xvec, double *B)
 {
   int    i;
   double x, y, z;
-  double kx, ky, kz, tz, kw;
+  double kx, ky, kz, tz;
   double cx, sx, chx, shx;
   double cy, sy, chy, shy;
   double cz, sz;
@@ -461,8 +460,6 @@ static void GWigB(struct gwigR *pWig, double *Xvec, double *B)
   y = Xvec[2];
   z = pWig->Zw;
 
-  kw   = 2e0*PI/(pWig->Lw);
-
   B[0] = 0;
   B[1] = 0;
   B[2] = 0;
@@ -596,9 +593,9 @@ static void GWigRadiationKicks(struct gwigR *pWig, double *X, double *Bxy, doubl
 static void Hessian(struct gwigR *pWig, double *Xvec, double *H2)
 {
 	int i,j,k;
-	double x,px,y,py,D;
-	double ax,axx,axxx,axy,axyy,axxy;
-	double ay, ayx, ayxx, ayy, ayyy, ayxy;
+	double px, py,D;
+	double ax,axx,axxx,axy,axxy;
+	double ay, ayx, ayxx, ayy, ayxy;
 	double Pax[6];
 	double Pay[6];
 	double H[6][6];
@@ -610,19 +607,15 @@ static void Hessian(struct gwigR *pWig, double *Xvec, double *H2)
 	axx =Pax[1];
 	axxx=Pax[2];
 	axy =Pax[3];
-	axyy=Pax[4];
 	axxy=Pax[5];
 
 	ay  =Pay[0];
 	ayx =Pay[1];
 	ayxx=Pay[2];
 	ayy =Pay[3];
-	ayyy=Pay[4];
 	ayxy=Pay[5];
 
-	x =Xvec[0];
 	px=Xvec[1];
-	y =Xvec[2];
 	py=Xvec[3];
 	D =Xvec[4];
 

atmat/atmexall.m

@@ -101,6 +101,7 @@ function atmexall(varargin)
 % Diffusion matrices
 cdir=fullfile(atroot,'atphysics','Radiation');
 compile([alloptions, {passinclude}], fullfile(cdir,'findmpoleraddiffmatrix.c'));
+compile([alloptions, {passinclude}], fullfile(cdir,'FDW.c'));
 
 % RDTs
 cdir=fullfile(atroot,'atphysics','NonLinearDynamics');

atmat/atphysics/Radiation/FDW.c

@@ -23,12 +23,10 @@
  * Journal of Instrumentation, May 2021.
  */
 
-#include "atlalib.c"
 #include "atelem.c"
+#include "atlalib.c"
 #include "gwigR.c"
 #include <math.h>
-#include "mex.h"
-#include "matrix.h"
 #include "gwig.h"
 
 /* Fourth order-symplectic integrator constants */
@@ -121,13 +119,11 @@ static void wigglerB(struct gwigR *pWig, double* orbit_in, double L, double *B66
 	double Bxyz[3];
 	double p_norm = 1/(1+orbit_in[4]);
 	double p_norm2 = SQR(p_norm);
-	double x,px,y,py,D;
+	double px,py,D;
 	double DLDS;
   int i;
 
-	x = orbit_in[0];
 	px= orbit_in[1];
-	y = orbit_in[2];
 	py= orbit_in[3];
 	D = orbit_in[4];
 	GWigAx(pWig, orbit_in, &ax, &axpy);
@@ -204,7 +200,6 @@ static void FindElemB(double *orbit_in, double le, double Lw, double Bmax,
 	double ax,ay,axpy,aypx;
 	double B[3];
 	double E0G;
-	double factorvect;
   struct gwigR pWig;
   int flag;
 
@@ -327,7 +322,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 	const mxArray *mxElem, *mxOrbit;
 	double *bdiff;
 	double orb[6];
-    double energy;
+    double energy = mxGetNaN(); /* Default energy value */
     int i, m, n;
 
 	if (nrhs < 2)
@@ -349,11 +344,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
 	    if (!(mxIsDouble(mxEnergy) && mxIsScalar(mxEnergy)))
 	        mexErrMsgIdAndTxt("AT:WrongArg", "Energy must be a double scalar");
 	    energy=mxGetScalar(mxEnergy);
-		}
+	}
 
-    else {
-        energy=mxGetNaN();
-    }
 	/* ALLOCATE memory for the output array */
 	plhs[0] = mxCreateDoubleMatrix(6,6,mxREAL);
 	bdiff = mxGetDoubles(plhs[0]);
@@ -374,7 +366,7 @@ static PyObject *compute_wiggdiffmatrix(PyObject *self, PyObject *args) {
     PyArrayObject *pyOrbit;
     double *orb0, *bdiff, *retval;
     double orb[6];
-    double energy;
+    double energy = NAN; /* Default energy value */
     npy_intp outdims[2] = {6, 6};
     int i;
 

pyat/at/lattice/elements.py

@@ -1301,17 +1301,16 @@ class Wiggler(Radiative, LongElement):
     See atwiggler.m
     """
 
-    _BUILD_ATTRIBUTES = LongElement._BUILD_ATTRIBUTES + ['Lw', 'Bmax',
-                                                         'Energy']
-    _conversions = dict(Element._conversions, Lw=float, Bmax=float,
-                        Energy=float,
+    _BUILD_ATTRIBUTES = LongElement._BUILD_ATTRIBUTES + ['Lw', 'Bmax']
+    _conversions = dict(Element._conversions, Lw=_float, Bmax=_float,
+                        Energy=_float,
                         Bx=lambda v: _array(v, (6, -1)),
                         By=lambda v: _array(v, (6, -1)),
                         Nstep=int, Nmeth=int, NHharm=int, NVharm=int)
 
     # noinspection PyPep8Naming
     def __init__(self, family_name: str, length: float, wiggle_period: float,
-                 b_max: float, energy: float, Nstep: Optional[int] = 5,
+                 b_max: float, energy: float = 0.0, *, Nstep: Optional[int] = 5,
                  Nmeth: Optional[int] = 4,
                  By=(1, 1, 0, 1, 1, 0), Bx=(), **kwargs):
         """

pyat/at/physics/__init__.py

@@ -12,6 +12,7 @@
 from .matrix import *
 from .linear import *
 from .diffmatrix import find_mpole_raddiff_matrix
+from .wiggdiffmatrix import FDW
 from .radiation import *
 from .ring_parameters import *
 from .nonlinear import *