Repository: SasView/sasmodels (language: python)

Issue:
Insert marketplacemodel bulk ferromagnets
I have added the marketplace model "SANS in bulk ferromagnets". It consists of the files "magnetic_functions.c", "micromagnetic_FF_3D.c", and "micromagnetic_FF_3D.py". In the issue  "Insert model "SANS of bulk ferromagnets" from the model marketplace #2784", I have mentioned that it would be good to (i) add the parameter descriptions for the choosable parameters in the GUI (e.g. radius = "Radius of the core", thickness = "Thickness of shell"), and (ii) to add the description of the whole model (given in the .py file) to the sasview webpage. Unfortunately, for both I do not know how to do this - could you have a look on that?


Discussion / hints:

Fix patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 4f6285832..49eae7462 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,10 +1,12 @@
 Release notes
 =============
 
-v1.0.8 2024-06-??
+v1.0.8 2024-09-26
 -----------------
+* New model: Bulk ferromagnets model from marketplace
 * Doc update: Archive built docs on Github
 * Doc update: Display math correctly
+* Fix error in FCC paracrystaline models
 * Fix parameter name checking in kernel call
 
 v1.0.7 2023-03-23
diff --git a/doc/guide/plugin.rst b/doc/guide/plugin.rst
index e7c0ba4bc..589d7872a 100644
--- a/doc/guide/plugin.rst
+++ b/doc/guide/plugin.rst
@@ -851,6 +851,8 @@ Some non-standard constants and functions are also provided:
         $x^2$
     cube(x):
         $x^3$
+    clip(a, a_min, a_max):
+        $\min(\max(a, a_\text{min}), a_\text{max})$, or NaN if $a$ is NaN.
     sas_sinx_x(x):
         $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
     powr(x, y):
diff --git a/explore/precision.py b/explore/precision.py
index 47eb71a85..cb19ce729 100755
--- a/explore/precision.py
+++ b/explore/precision.py
@@ -190,13 +190,13 @@ def plotdiff(x, target, actual, label, diff):
     if diff == "relative":
         err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
         #err = np.clip(err, 0, 1)
-        pylab.loglog(x, err, '-', label=label)
+        pylab.loglog(x, err, '-', label=label, alpha=0.7)
     elif diff == "absolute":
         err = np.array([abs((t-a)) for t, a in zip(target, actual)], 'd')
-        pylab.loglog(x, err, '-', label=label)
+        pylab.loglog(x, err, '-', label=label, alpha=0.7)
     else:
         limits = np.min(target), np.max(target)
-        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label)
+        pylab.semilogx(x, np.clip(actual, *limits), '-', label=label, alpha=0.7)
 
 def make_ocl(function, name, source=[]):
     class Kernel(object):
@@ -412,6 +412,54 @@ def add_function(name, mp_function, np_function, ocl_function,
     np_function=lambda x: np.fmod(x, 2*np.pi),
     ocl_function=make_ocl("return fmod(q, 2*M_PI);", "sas_fmod"),
 )
+
+def sas_langevin(x):
+    scalar = np.isscalar(x)
+    if scalar:
+        x = np.array([x]) # should inherit dtype for single if given single
+    f = np.empty_like(x)
+    cutoff = 0.1 if f.dtype == np.float64 else 1.0
+    #cutoff *= 10
+    index = x < cutoff
+    xp = x[index]
+    xpsq = xp*xp
+    f[index] = xp / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9.))))
+    # 4 terms gets to 1e-7 single, 1e-14 double. Can get to 1e-15 double by adding
+    # another 4 terms and setting cutoff at 1.0. Not worthwhile. Instead we would
+    # need an expansion about x somewhere between 1 and 10 for the interval [0.1, 100.]
+    #f[index] = xp / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9. + xpsq/(11.0 + xpsq/(13. + xpsq/(15. + xpsq/17.)))))))
+    xp = x[~index]
+    f[~index] = 1/np.tanh(xp) - 1/xp
+    return f[0] if scalar else f
+
+def sas_langevin_x(x):
+    scalar = np.isscalar(x)
+    if scalar:
+        x = np.array([x]) # should inherit dtype for single if given single
+    f = np.empty_like(x)
+    cutoff = 0.1 if f.dtype == np.float64 else 1.0
+    index = x < cutoff
+    xp = x[index]
+    xpsq = xp*xp
+    f[index] = 1. / (3. + xpsq / (5. + xpsq/(7. + xpsq/(9.))))
+    xp = x[~index]
+    f[~index] = (1/np.tanh(xp) - 1/xp)/xp
+    return f[0] if scalar else f
+
+add_function(
+    name="langevin(x)",
+    mp_function=lambda x: (1/mp.tanh(x) - 1/x),
+    np_function=sas_langevin,
+    #ocl_function=make_ocl("return q < 0.7 ? q*(1./3. + q*q*(-1./45. + q*q*(2./945. + q*q*(-1./4725.) + q*q*(2./93555.)))) : 1/tanh(q) - 1/q;", "sas_langevin"),
+    ocl_function=make_ocl("return q < 1e-5 ? q/3. : 1/tanh(q) - 1/q;", "sas_langevin"),
+)
+add_function(
+    name="langevin(x)/x",
+    mp_function=lambda x: (1/mp.tanh(x) - 1/x)/x,
+    #np_function=lambda x: sas_langevin(x)/x, # Note: need to test for x=0
+    np_function=sas_langevin_x,
+    ocl_function=make_ocl("return q < 1e-5 ? 1./3. : (1/tanh(q) - 1/q)/q;", "sas_langevin_x"),
+)
 add_function(
     name="gauss_coil",
     mp_function=lambda x: 2*(mp.exp(-x**2) + x**2 - 1)/x**4,
diff --git a/sasmodels/kernel_header.c b/sasmodels/kernel_header.c
index f27584d58..a65b1d7d0 100644
--- a/sasmodels/kernel_header.c
+++ b/sasmodels/kernel_header.c
@@ -187,8 +187,54 @@
 #endif
 inline double square(double x) { return x*x; }
 inline double cube(double x) { return x*x*x; }
+// clip() follows numpy.clip() semantics, returning (x < low ? low : x > high ? high : x)
+// OpenCL/CUDA clamp() returns fmin(fmax(x, low), high)
+// C++(17) clamp() matches numpy.clip()
+// If x is NaN numpy.clip() returns NaN but OpenCL clamp() returns low.
+inline double clip(double x, double low, double high) { return x < low ? low : x > high ? high : x; }
 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
+// vector algebra
+static void SET_VEC(double *vector, double v0, double v1, double v2)
+{
+    vector[0] = v0;
+    vector[1] = v1;
+    vector[2] = v2;
+}
+
+static void SCALE_VEC(double *vector, double a)
+{
+    vector[0] = a*vector[0];
+    vector[1] = a*vector[1];
+    vector[2] = a*vector[2];
+}
+
+static void ADD_VEC(double *result_vec, double *vec1, double *vec2)
+{
+    result_vec[0] = vec1[0] + vec2[0];
+    result_vec[1] = vec1[1] + vec2[1];
+    result_vec[2] = vec1[2] + vec2[2];
+}
+
+static double SCALAR_VEC( double *vec1, double *vec2)
+{
+    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
+}
+
+static double MAG_VEC( double *vec)
+{
+    return sqrt(SCALAR_VEC(vec,vec));
+}
+
+
+static void ORTH_VEC(double *result_vec, double *vec1, double *vec2)
+{
+    double scale =  SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2);
+    result_vec[0] = vec1[0] - scale * vec2[0];
+    result_vec[1] = vec1[1] - scale * vec2[1];
+    result_vec[2] = vec1[2] - scale * vec2[2];
+}
+
 // CRUFT: support old style models with orientation received qx, qy and angles
 
 // To rotate from the canonical position to theta, phi, psi, first rotate by
diff --git a/sasmodels/kernel_iq.c b/sasmodels/kernel_iq.c
index f7e73add0..b8e3394be 100644
--- a/sasmodels/kernel_iq.c
+++ b/sasmodels/kernel_iq.c
@@ -70,52 +70,8 @@ typedef union {
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
 // ===== Helper functions for magnetism =====
 
-// vector algebra
-void SET_VEC(double *vector, double v0, double v1, double v2)
-{
-    vector[0] = v0;
-    vector[1] = v1;
-    vector[2] = v2;
-}
-
-void SCALE_VEC(double *vector, double a)
-{
-    vector[0] = a*vector[0];
-    vector[1] = a*vector[1];
-    vector[2] = a*vector[2];
-}
 
-void ADD_VEC(double *result_vec, double *vec1, double *vec2)
-{
-    result_vec[0] = vec1[0] + vec2[0];
-    result_vec[1] = vec1[1] + vec2[1];
-    result_vec[2] = vec1[2] + vec2[2];
-}
-
-static double SCALAR_VEC( double *vec1, double *vec2)
-{
-    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2];
-}
 
-static double MAG_VEC( double *vec)
-{
-    return sqrt(SCALAR_VEC(vec,vec));
-}
-
-void ORTH_VEC(double *result_vec, double *vec1, double *vec2)
-{
-    double scale =  SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2);
-    result_vec[0] = vec1[0] - scale * vec2[0];
-    result_vec[1] = vec1[1] - scale * vec2[1];
-    result_vec[2] = vec1[2] - scale * vec2[2];
-}
-
-
-// Return value restricted between low and high
-static double clip(double value, double low, double high)
-{
-  return (value < low ? low : (value > high ? high : value));
-}
 
 // Compute spin cross sections given in_spin and out_spin
 // To convert spin cross sections to sld b:
@@ -178,11 +134,10 @@ static double mag_sld(
 )
 {
   double Mvector[3];
-  double Pvector[3]; 
+  double Pvector[3];
   double qvector[3];
-  double rhom[3]; 
   double perpy[3];
-  double perpz[3];  
+  double perpz[3];
   double Mperp[3];
 
   const double qsq = sqrt(qx*qx + qy*qy); 
diff --git a/sasmodels/models/fcc_paracrystal.c b/sasmodels/models/fcc_paracrystal.c
index 76ef61beb..e3f9788e8 100644
--- a/sasmodels/models/fcc_paracrystal.c
+++ b/sasmodels/models/fcc_paracrystal.c
@@ -5,7 +5,7 @@ fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
     const double a1 = ( qa + qb)/2.0;
     const double a2 = ( qa + qc)/2.0;
     const double a3 = ( qb + qc)/2.0;
-    const double d_a = dnn/sqrt(2);
+    const double d_a = dnn*sqrt(2);
 
     // Matsuoka 23-24-25
     //     Z_k numerator: 1 - exp(a)^2
diff --git a/sasmodels/models/fcc_paracrystal.py b/sasmodels/models/fcc_paracrystal.py
index 35b19b16f..b8711cfcf 100644
--- a/sasmodels/models/fcc_paracrystal.py
+++ b/sasmodels/models/fcc_paracrystal.py
@@ -213,9 +213,10 @@ def random():
 #
 # October 26, 2022, PDB fixed unit tests to conform to new maths
 # TODO: fix the 2d tests
-q = 4.*pi/220.
+# Test position of first Bragg Peak at q = 2.0 * pi * np.sqrt(3) / (np.sqrt(2) * 220)
+q=0.035
 tests = [
-    [{}, [0.001, q, 0.25], [1.9839734995338474, 0.3352457353010224, 0.005804136688760973]],
+    [{}, [0.01, q, 0.25], [0.0213498915484313, 69.03586722100925, 0.005713137547083953]],
     #[{}, (-0.047, -0.007), 238.103096286],
     #[{}, (0.053, 0.063), 0.863609587796],
 ]
diff --git a/sasmodels/models/img/micromagnetic_FF.png b/sasmodels/models/img/micromagnetic_FF.png
new file mode 100644
index 000000000..b316ccb56
Binary files /dev/null and b/sasmodels/models/img/micromagnetic_FF.png differ
diff --git a/sasmodels/models/lib/magnetic_functions.c b/sasmodels/models/lib/magnetic_functions.c
new file mode 100644
index 000000000..c37b3f70c
--- /dev/null
+++ b/sasmodels/models/lib/magnetic_functions.c
@@ -0,0 +1,130 @@
+//These functions are required for magnetic analysis models. They are copies 
+//from sasmodels/kernel_iq.c, to enables magnetic parameters for 1D and 2D models.
+
+
+static double langevin(
+    double x) {
+    // cotanh(x)-1\x
+
+    if (x < 0.00001) {
+        // avoid dividing by zero
+        return 1.0/3.0*x-1.0/45.0 * pow(x, 3) + 2.0/945.0 * pow(x, 5) - 1.0/4725.0 * pow(x, 7);
+    } else {
+        return 1.0/tanh(x)-1/x;
+    }
+}
+
+static double langevinoverx(
+    double x) 
+{
+    // cotanh(x)-1\x
+
+    if (x < 0.00001) {
+        // avoid dividing by zero
+        return 1.0/3.0-1.0/45.0 * pow(x, 2) + 2.0/945.0 * pow(x, 4) - 1.0/4725.0 * pow(x, 6);
+    } else {
+        return langevin(x)/x;
+    }
+}
+
+
+
+//weighting of spin resolved cross sections to reconstruct partially polarised beam with imperfect optics using up_i/up_f.
+static void set_weights(double in_spin, double out_spin, double weight[8]) //from kernel_iq.c
+{
+ double norm=out_spin;
+
+
+   in_spin = clip(sqrt(square(in_spin)), 0.0, 1.0);//opencl has ambiguities for abs()
+   out_spin = clip(sqrt(square(out_spin)), 0.0, 1.0);
+
+   if (out_spin < 0.5){norm=1-out_spin;}
+
+// The norm is needed to make sure that the scattering cross sections are
+//correctly weighted, such that the sum of spin-resolved measurements adds up to
+// the unpolarised or half-polarised scattering cross section. No intensity weighting
+// needed on the incoming polariser side assuming that a user has normalised
+// to the incoming flux with polariser in for SANSPOl and unpolarised beam, respectively.
+
+
+   weight[0] = (1.0-in_spin) * (1.0-out_spin) / norm; // dd.real
+   weight[1] = weight[0]; // dd.imag
+   weight[2] = in_spin * out_spin / norm;             // uu.real
+   weight[3] = weight[2];             // uu.imag
+   weight[4] = (1.0-in_spin) * out_spin / norm;       // du.real
+   weight[5] = weight[4];       // du.imag 
+   weight[6] = in_spin * (1.0-out_spin) / norm;       // ud.real
+   weight[7] = weight[6];       // ud.imag
+ }
+
+
+
+//transforms scattering vector q in polarisation/magnetisation coordinate system
+ static void set_scatvec(double *qrot, double q,
+  double cos_theta, double sin_theta,
+  double alpha, double beta)
+{
+  double cos_alpha, sin_alpha;
+  double cos_beta, sin_beta;
+
+  SINCOS(alpha*M_PI/180, sin_alpha, cos_alpha);
+  SINCOS(beta*M_PI/180, sin_beta, cos_beta);
+  //field is defined along (0,0,1), orientation of detector
+  //is precessing in a cone around B with an inclination of theta
+
+  qrot[0] = q*(cos_alpha * cos_theta);
+  qrot[1] = q*(cos_theta * sin_alpha*sin_beta + 
+  cos_beta * sin_theta);
+  qrot[2] = q*(-cos_beta * cos_theta* sin_alpha + 
+  sin_beta * sin_theta);
+}
+
+
+
+//Evaluating the magnetic scattering vector (Halpern Johnson vector) for general orientation of q and collecting terms for the spin-resolved (POLARIS) cross sections. Mz is along the applied magnetic field direction, which is also the polarisation direction.
+static void mag_sld(
+   // 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag,  4=du.real, 5=du.imag,  6=ud.real, 7=ud.imag
+ double qx, double qy, double qz,
+  double mxreal, double mximag, double myreal,  double myimag, double mzreal,double mzimag, double nuc, double sld[8])
+{
+  double vector[3];
+  //The (transversal) magnetisation and hence the magnetic scattering sector is here a complex quantity. The spin-flip (magnetic) scattering amplitude is given with
+  //  MperpPperpQ \pm i MperpP  (Moon-Riste-Koehler Phys Rev 181, 920, 1969) with Mperp and MperpPperpQ the
+  //magnetisation scattering vector components perpendicular to the polarisation/field direction. Collecting terms in SF that are real (MperpPperpQreal + SCALAR_VEC(MperpPimag,qvector) )  and imaginary (MperpPperpQimag \pm SCALAR_VEC(MperpPreal,qvector) )
+  double Mvectorreal[3];
+  double Mvectorimag[3];
+  double Pvector[3];
+  double perpy[3];
+  double perpx[3]; 
+
+  double Mperpreal[3];
+  double Mperpimag[3];
+  
+  const double q = sqrt(qx*qx + qy*qy + qz*qz);
+  SET_VEC(vector, qx/q, qy/q, qz/q); 
+
+  //Moon-Riste-Koehler notation choose z as pointing along field/polarisation axis
+  //totally different to what is used in SASview (historical reasons)
+  SET_VEC(Mvectorreal, mxreal, myreal, mzreal);
+  SET_VEC(Mvectorimag, mximag, myimag, mzimag);
+
+  SET_VEC(Pvector, 0, 0, 1);
+  SET_VEC(perpy, 0, 1, 0);  
+  SET_VEC(perpx, 1, 0, 0);   
+   //Magnetic scattering vector Mperp could be simplified like in Moon-Riste-Koehler
+   //leave the generic computation just to check
+  ORTH_VEC(Mperpreal, Mvectorreal, vector);
+  ORTH_VEC(Mperpimag, Mvectorimag, vector);
+  
+  
+   sld[0] = nuc - SCALAR_VEC(Pvector,Mperpreal); // dd.real => sld - D Pvector \cdot Mperp 
+   sld[1] = +SCALAR_VEC(Pvector,Mperpimag); //dd.imag  = nuc_img - SCALAR_VEC(Pvector,Mperpimg); nuc_img only exist for noncentrosymmetric nuclear structures; 
+   sld[2] = nuc + SCALAR_VEC(Pvector,Mperpreal);              // uu => sld + D Pvector \cdot Mperp
+   sld[3] = -SCALAR_VEC(Pvector,Mperpimag); //uu.imag
+
+   sld[4] = SCALAR_VEC(perpy,Mperpreal)+SCALAR_VEC(perpx,Mperpimag);       // du.real => real part along y + imaginary part along x
+   sld[5] = SCALAR_VEC(perpy,Mperpimag)-SCALAR_VEC(perpx,Mperpreal);       // du.imag => imaginary component along y - i *real part along x
+   sld[6] = SCALAR_VEC(perpy,Mperpreal)-SCALAR_VEC(perpx,Mperpimag);      // ud.real =>  real part along y - imaginary part along x
+   sld[7] = SCALAR_VEC(perpy,Mperpimag)+SCALAR_VEC(perpx,Mperpreal);         // ud.imag => imaginary component along y + i * real part along x
+
+ }
\ No newline at end of file
diff --git a/sasmodels/models/micromagnetic_FF_3D.c b/sasmodels/models/micromagnetic_FF_3D.c
new file mode 100644
index 000000000..e64154f00
--- /dev/null
+++ b/sasmodels/models/micromagnetic_FF_3D.c
@@ -0,0 +1,190 @@
+//Core-shell form factor for anisotropy field (Hkx, Hky and Hkz), Nuc and lonitudinal magnetization Mz
+static double
+form_volume(double radius, double thickness)
+{
+  if( radius + thickness > radius + thickness)
+    return M_4PI_3 * cube(radius + thickness);
+  else
+    return M_4PI_3 * cube(radius + thickness);
+
+}
+
+
+
+static double fq(double q, double radius,
+ double thickness, double core_sld, double shell_sld, double solvent_sld)
+{
+ const double form = core_shell_fq(q,
+  radius,
+  thickness,
+  core_sld,
+  shell_sld,
+  solvent_sld);
+ return form;
+}
+
+
+static double reduced_field(double q, double Ms, double Hi,
+ double A)
+{
+    // q in 10e10 m-1, A in 10e-12 J/m, mu0 in 1e-7 
+    return Ms / (fmax(Hi, 1.0e-6) + 2.0 * A * 4.0 * M_PI / Ms * q * q * 10.0);
+
+ }
+
+ static double DMI_length(double Ms, double D, double qval)
+ {
+   return 2.0 * D * 4.0 * M_PI / Ms / Ms * qval ; //q in 10e10 m-1, A in 10e-3 J/m^2, mu0 in 4 M_PI 1e-7 
+ }
+
+//Mz is defined as the longitudinal magnetisation component along the magnetic field.
+//In the approach to saturation this component is (almost) constant with magnetic 
+//field and simplfy reflects the nanoscale variations in the saturation magnetisation
+//value in the sample. The misalignment of the magnetisation due to perturbing
+//magnetic anisotropy or dipolar magnetic fields in the sample enter Mx and My,
+//the two transversal magnetisation components, reacting to a magnetic field.
+//The micromagnetic solution for the magnetisation are from Michels et al. PRB 94, 054424 (2016).
+
+static double fqMxreal( double qx, double qy, double qz, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  double qsq = qx * qx + qy * qy + qz * qz;
+  double q = sqrt(qsq); 
+  double Hr = reduced_field(q, Ms, Hi, A);
+  double DMI = DMI_length(Ms, D,q);
+  double DMIz = DMI_length(Ms, D,qz);
+  double denominator = (qsq + Hr * (qx * qx + qy * qy) - square(Hr * DMIz * q)) / Hr;
+  double f = (Hkx * (qsq + Hr * qy * qy) - Ms * Mz * qx * qz * (1.0 + Hr * square(DMI)) - Hky * Hr * qx * qy) / denominator;
+  return f;
+}
+
+static double fqMximag(double qx, double qy, double qz, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  double qsq = qx * qx + qy * qy + qz * qz;
+  double q = sqrt(qsq);   
+  double Hr = reduced_field(q, Ms, Hi, A);
+  double DMI = DMI_length(Ms, D,q);
+  double DMIz = DMI_length(Ms, D,qz);
+  double DMIy = DMI_length(Ms, D,qy);
+  double denominator = (qsq + Hr * (qx * qx + qy * qy) - square(Hr * DMIz * q)) / Hr;
+  double f = - qsq * (Ms * Mz * (1.0 + Hr) * DMIy + Hky * Hr * DMIz) / denominator;
+  return f;
+}
+
+static double fqMyreal( double qx, double qy, double qz, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  double qsq = qx * qx + qy * qy + qz * qz;
+  double q = sqrt(qsq); 
+  double Hr = reduced_field(q, Ms, Hi, A);
+  double DMI = DMI_length(Ms, D,q);
+  double DMIz = DMI_length(Ms, D,qz);
+  double denominator = (qsq + Hr * (qx * qx + qy * qy) - square(Hr * DMIz * q))/Hr;
+  double f = (Hky * (qsq + Hr * qx * qx) - Ms * Mz * qy * qz * (1.0 + Hr * square(DMI)) - Hkx * Hr * qx * qy)/denominator;
+  return f;
+}
+
+static double fqMyimag( double qx, double qy, double qz, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  double qsq = qx * qx + qy * qy + qz * qz;
+  double q = sqrt(qsq);  
+  double Hr = reduced_field(q, Ms, Hi, A);
+  double DMI = DMI_length(Ms, D,q);
+  double DMIx = DMI_length(Ms, D,qx);
+  double DMIz = DMI_length(Ms, D,qz);
+  double denominator = (qsq + Hr * (qx * qx + qy * qy) - square(Hr * DMIz * q))/Hr;
+  double f = qsq * (Ms * Mz * (1.0 + Hr) * DMIx - Hkx * Hr * DMIz)/denominator;
+  return f;
+}
+
+static double
+Calculate_Scattering(double q, double cos_theta, double sin_theta, double radius, double thickness, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D,  double up_i, double up_f, double alpha, double beta)
+{
+    double qrot[3];
+    set_scatvec(qrot,q,cos_theta, sin_theta, alpha, beta);
+    // 0=dd.real, 1=dd.imag, 2=uu.real, 3=uu.imag,  4=du.real, 6=du.imag,  7=ud.real, 5=ud.imag
+    double weights[8];
+    set_weights(up_i, up_f, weights);
+   
+    double mz = fq(q, radius, thickness, mag_sld_core, mag_sld_shell, mag_sld_solvent);
+    double nuc = fq(q, radius, thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent);
+	double hk = fq(q, radius, 0, hk_sld_core, 0, 0);
+
+    double cos_gamma, sin_gamma;
+    double sld[8];
+    //loop over random anisotropy axis with isotropic orientation gamma for Hkx and Hky
+    //To be modified for textured material see also Weissmueller et al. PRB 63, 214414 (2001)
+    double total_F2 = 0.0;
+    for (int i = 0; i<GAUSS_N ;i++) {
+      const double gamma = M_PI * (GAUSS_Z[i] + 1.0); // 0 .. 2 pi
+      SINCOS(gamma, sin_gamma, cos_gamma);	
+      //Only the core of the defect/particle in the matrix has an effective
+      //anisotropy (for simplicity), for the effect of different, more complex
+      // spatial profile of the anisotropy see Michels PRB 82, 024433 (2010)
+      double Hkx = hk * sin_gamma;
+      double Hky = hk * cos_gamma;
+
+      double mxreal = fqMxreal(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
+      double mximag = fqMximag(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
+      double myreal = fqMyreal(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
+      double myimag = fqMyimag(qrot[0], qrot[1], qrot[2], mz, Hkx, Hky, Hi, Ms, A, D);
+
+      mag_sld(qrot[0], qrot[1], qrot[2], mxreal, mximag, myreal, myimag, mz, 0, nuc, sld);
+      double form = 0.0;
+      for (unsigned int xs = 0; xs<8; xs++) {
+        if (weights[xs] > 1.0e-8) {
+          // Since the cross section weight is significant, set the slds
+          // to the effective slds for this cross section, call the
+          // kernel, and add according to weight.
+          // loop over uu, ud real, du real, dd, ud imag, du imag 
+          form += weights[xs] * sld[xs] * sld[xs];
+        }
+      }
+      total_F2 += GAUSS_W[i] * form ;      
+    }
+    return total_F2;
+}
+
+
+static double
+Iqxy(double qx, double qy, double radius, double thickness, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D,  double up_i, double up_f, double alpha, double beta)
+{
+  double sin_theta, cos_theta;
+  const double q = sqrt(qx * qx + qy * qy);  
+  if (q > 1.0e-16 ) {
+    cos_theta = qx/q;
+    sin_theta = qy/q;
+  } else {
+    cos_theta = 0.0;
+    sin_theta = 0.0;
+  }  
+
+  double total_F2 = Calculate_Scattering(q, cos_theta, sin_theta, radius, thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent, mag_sld_core, mag_sld_shell, mag_sld_solvent, hk_sld_core, Hi, Ms, A, D,  up_i, up_f, alpha, beta);
+	
+  return 0.5 * 1.0e-4 * total_F2;
+    
+}
+
+
+
+static double
+Iq(double q, double radius, double thickness, double nuc_sld_core, double nuc_sld_shell, double nuc_sld_solvent, double mag_sld_core, double mag_sld_shell, double mag_sld_solvent, double hk_sld_core, double Hi, double Ms, double A, double D,  double up_i, double up_f, double alpha, double beta)
+{
+  // slots to hold sincos function output of the orientation on the detector plane
+  double sin_theta, cos_theta; 
+  double total_F1D = 0.0;
+  for (int j = 0; j<GAUSS_N ;j++) {
+
+    const double theta = M_PI * (GAUSS_Z[j] + 1.0); // 0 .. 2 pi
+    SINCOS(theta, sin_theta, cos_theta);
+
+    double total_F2 = Calculate_Scattering(q, cos_theta, sin_theta, radius, thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent, mag_sld_core, mag_sld_shell, mag_sld_solvent, hk_sld_core, Hi, Ms, A, D,  up_i, up_f, alpha, beta);	
+    total_F1D += GAUSS_W[j] * total_F2 ;
+  }
+  //convert from [1e-12 A-1] to [cm-1]
+  return 0.25 * 1.0e-4 * total_F1D;
+}
+
+
+
+
+
+
diff --git a/sasmodels/models/micromagnetic_FF_3D.py b/sasmodels/models/micromagnetic_FF_3D.py
new file mode 100644
index 000000000..8e0973e7b
--- /dev/null
+++ b/sasmodels/models/micromagnetic_FF_3D.py
@@ -0,0 +1,181 @@
+r"""
+Definition
+----------
+This model is a micromagnetic approach to analyse the SANS that arises from
+nanoscale variations in the magnitude and orientation of the magnetization in
+bulk ferromagnets in the approach to magnetic saturation (single domain state).
+Typical materials are cold-worked elemental magnets, hard and soft magnetic 
+nanocomposites, amorphous alloys and precipitates in magnetic steel [#Michels2014]_.
+The magnetic SANS depends on the magnetic interactions, the magnetic microstructure
+(defect/particle size, magnetocrystalline anisotropy, saturation magnetisation)
+and on the applied magnetic field. As shown in [#Michels2016]_ near magnetic 
+saturation the scattering cross-section can be evaluated by means of micromagnetic theory
+
+.. math::
+    I(\mathbf{Q}) = I_{nuc} + I_{mag}(\mathbf{Q},H),
+
+with the field-independent nuclear and magnetic SANS cross section (due
+to nanoscale spatial variations of the magnetisation).
+
+.. math::
+    I_{mag}(\mathbf{Q},H)= S_K(Q) R_K(\mathbf{Q}, H_i) + S_M(Q) R_M(\mathbf{Q}, H_i),
+
+with $H_i$ the internal field, i.e. the external magnetic field corrected for
+demagnetizing effects and the influence of the magnetodipolar field and of the
+magnetic anisotropy [#Bick2013]_. This magnetic field dependence of the scattering 
+reflects the increasing magnetisation misalignment with decreasing
+externally applied magnetic field with a contribution $S_K \times R_K$ due to
+perturbations around magnetic anisotropy fields and a term $S_M \times R_M$
+related to magnetostatic fields. The magnetic moments decorate perturbations in the
+microstructure (precipitates, grain boundaries etc).
+The anisotropy-field function $S_K$ depends on the Fourier transform of the magnetic
+anisotropy distribution (strength and orientation) in the material, and the
+scattering function of the longitudinal magnetisation $S_M$ reflects the
+variations of the saturation magnetisation, e.g. jumps at the particle-matrix
+interface. $R_K$ and $R_M$ denote the micromagnetic response functions that
+describe the magnetisation distribution around a perturbation in magnetic
+anisotropy and flucutations in the saturation magnetisation value.
+
+.. figure:: img/micromagnetic_FF.png
+
+    Magnetisation distribution around (left) a particle with magnetic easy axis
+    in the vertical direction and (right) a precipitation with a magnetisation
+    that is higher than the matrix phase.  
+
+The micromagnetic response functions depend on magnetic material parameters $M_S$:
+average saturation magnetisation of the material, $H_i$: the internal magnetic
+field, $A$ the average exchange-stiffness constant. In the vicinity of lattice 
+imperfection in ferromagnetic materials, antisymmetric Dzyaloshinskii–Moriya
+interaction (DMI) can occur due to the local structural inversion symmetry
+breaking [#Arrott1963]_. DMI with strength $D$ can give rise to nonuniform spin
+textures resulting in a polarization-dependent asymmetric scattering term for 
+polycrystalline ferromagnetic with a centrosymmetric crystal structure [#Michels2016]_.
+We assume (for simplicity) an isotropic microstructure (for $S_M$) and random 
+orientation of magnetic easy axes (additionally for $S_K$) such that the 
+contributions of the magnetic microstructure only depend on the magnitude of $q$.
+Considerations for a microstructure with a prefered orientation (texture) can be
+found in [#Weissmueller2001]_. In the code the averaging procedure over the random
+anisotropy is explicitely performed. A specific orientation distribution can be
+implemented by rewriting the model. 
+
+The magnetic field is oriented with an inclination of $\alpha$ to the neutron beam
+and rotated by $\beta$. The model for the nuclear scattering amplitude, saturation
+magnetisation is based on spherical particles with a core shell structure. For
+simplicity, only the core has an effective anisotropy, that is varying randomly
+in direction from particle to particle. The effect of different, more complex 
+spatial profiles of the anisotropy can be seen in [#Michels2010]_.
+The magnetic scattering length density (SLD) is defined as 
+$\rho_{\mathrm{mag}}=b_H M_S$, where $b_H= 2.91*10^{8}A^{-1}m^{-1}$ and $M_S$
+is the saturation magnetisation (in $A/m$).
+
+The fraction of "upward" neutrons before ('up_frac_i') and after the sample 
+('up_frac_f') must range between 0 to 1, with 0.5 denoting an unpolarised beam.
+Note that a fit may result in a negative magnetic SLD, and hence magnetisation, 
+when the polarisation state is inverted, i.e. if you have analysed for a $I_{00}$
+state wheras your data are $I_{11}$. The model allows to construct the 4 
+spin-resolved cross sections (non-spin-flip $I_{00}$, $I_{11}$ and spin-flip, here
+$I_{01}=I_{10}$), half-polarised SANS (SANSpol, incoming polarised beam $I_0$ and
+$I_1$, no analysis after sample 'up_frac_f'$=0.5$), and unpolarised beam 
+('up_frac_i'$=$'up_frac_f'$=0.5$). Differences and other combinations between
+polarised scattering cross section, e.g. to obtain the nuclear-magnetic 
+interference scattering, or subtraction of the residual scattering of the high
+field reference state can be constructed with a custom model (Fitting>
+Add/Multiply Model) and using approbriate scales. For dense systems, special
+care has to be taken as the nculear structure factor (arrangement of particles)
+does not need to be identical with the magnetic microstructure e.g. local 
+textures and correlations between easy axes (see [#Honecker2020]_ for further
+details). The use of structure model is therefore strongly discouraged. Better
+$I_{nuc}$, $S_K$ and $S_M$ are fit independent from each other in a model-free way. 
+
+
+
+References
+----------
+
+.. [#Arrott1963] A. Arrott, J. Appl. Phys. 34, 1108 (1963).
+.. [#Weissmueller2001] J. Weissmueller et al., *Phys. Rev. B* 63, 214414 (2001).
+.. [#Bick2013] J.-P. Bick et al., *Appl. Phys. Lett.* 102, 022415 (2013).
+.. [#Michels2010] A. Michels et al., *Phys. Rev. B* 82, 024433 (2010).
+.. [#Michels2014] A. Michels, *J. Phys.: Condens. Matter* 26, 383201 (2014).
+.. [#Michels2016] A. Michels et al., *Phys. Rev. B* 94, 054424 (2016).
+.. [#Honecker2020] D. Honecker, L. Fernandez Barguin, and P. Bender, *Phys. Rev. B* 101, 134401 (2020).
+
+
+
+Authorship and Verification
+----------------------------
+
+* **Author:** Dirk Honecker **Date:** January 14, 2021
+* **Last Modified by:** Dirk Honecker **Date:** September 23, 2024
+* **Last Reviewed by:**
+
+"""
+
+import numpy as np
+from numpy import pi, inf
+
+name = "micromagnetic_FF_3D"
+title = "Field-dependent magnetic microstructure around imperfections in bulk ferromagnets"
+description = """
+    I(q) = A (F_N^2(q)+ C F_N F_M + D F_M^2) +B(H) I_mag(q,H) 
+            A: weighting function =1 for unpolarised beam and non-neutron-spin-flip scattering, zero for spin-flip scattering. The terms in the bracket are the residual scattering at perfectly saturating magnetic field.
+            B(H): weighting function for purely magnetic scattering I_mag(q,H)  due to misaligned magnetic moments, different for the various possible spin-resolved scattering cross sections
+            C: weighting function for nuclear-magnetic interference scattering
+            F_N: nuclear form factor
+            F_M: magnetic form factor
+The underlying defect can have a core-shell structure.
+"""
+category = "shape:sphere"
+
+# pylint: disable=bad-whitespace, line-too-long
+#             ["name", "units", default, [lower, upper], "type","description"],
+parameters = [["radius", "Ang", 50., [0, inf], "volume", "Structural radius of the core"],
+              ["thickness", "Ang", 40., [0, inf], "volume", "Structural thickness of shell"],            
+              ["nuc_sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", "Core scattering length density"],
+              ["nuc_sld_shell", "1e-6/Ang^2", 1.7, [-inf, inf], "", "Scattering length density of shell"],
+              ["nuc_sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", "Solvent scattering length density"],
+              ["mag_sld_core",    "1e-6/Ang^2", 1.0,  [-inf, inf], "",    "Magnetic scattering length density of core"],
+              ["mag_sld_shell", "1e-6/Ang^2", 1.7, [-inf, inf], "", "Magnetic scattering length density of shell"],
+              ["mag_sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "",  "Magnetic scattering length density of solvent"],
+              ["hk_sld_core", "1e-6/Ang^2", 1.0,  [0, inf], "",    "Anisotropy field of defect"],
+              ["Hi", "T", 2.0,  [0, inf], "",    "Effective field inside the material"],
+              ["Ms", "T", 1.0,  [0, inf], "",    "Volume averaged saturation magnetisation"], 
+              ["A", "pJ/m", 10.0,  [0, inf], "",    "Average exchange stiffness constant"],  
+              ["D", "mJ/m^2", 0.0,  [0, inf], "",    "Average DMI constant"],                          
+              ["up_i",          "None",  0.5, [0, 1], "", "Polarisation incoming beam"],
+              ["up_f",          "None",  0.5, [0, 1], "", "Polarisation outgoing beam"],
+              ["alpha",          "None",  90, [0, 180], "", "Inclination of field to neutron beam"],
+              ["beta",          "None",  0, [0, 360], "", "Rotation of field around neutron beam"],             
+             ]
+# pylint: enable=bad-whitespace, line-too-long
+
+
+
+
+source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "lib/gauss76.c", "lib/magnetic_functions.c", "micromagnetic_FF_3D.c"]
+structure_factor = False
+have_Fq = False
+single=False
+opencl = False
+
+
+def random():
+    """Return a random parameter set for the model."""
+    outer_radius = 10**np.random.uniform(1.3, 4.3)
+    # Use a distribution with a preference for thin shell or thin core
+    # Avoid core,shell radii < 1
+    radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
+    thickness = outer_radius - radius
+    pars = dict(
+        radius=radius,
+        thickness=thickness,
+        
+    )
+    return pars
+
+
+
+tests = [
+     [{},1.002266990452620e-03, 7.461046163627724e+03],
+     [{},(0.0688124,  -0.0261013),  22.024],
+]
diff --git a/sasmodels/models/sphere.py b/sasmodels/models/sphere.py
index 1088e3e5a..f64070778 100644
--- a/sasmodels/models/sphere.py
+++ b/sasmodels/models/sphere.py
@@ -75,7 +75,6 @@
 have_Fq = True
 radius_effective_modes = ["radius"]
 #single = False
-
 def random():
     """Return a random parameter set for the model."""
     radius = 10**np.random.uniform(1.3, 4)
@@ -106,6 +105,21 @@ def random():
      0.1, 482.93824329, 29763977.79867414, 120.0, 8087664.122641933, 1.0],
     [{"radius": 120., "radius_pd": 0.2, "radius_pd_n": 45},
      0.2, 1.23330406, 1850806.1197361, 120.0, 8087664.122641933, 1.0],
+
+    # For 2-D data use (qx, qy) pairs. Since sphere is radial, just need the
+    # correct |q| value for the test, so use 3-4-5 triangle. The test code
+    # looks for tuples to detect 2-D data, so can't use simple numpy cheats.
+    [{"scale": 1., "background": 0., "sld": 6., "sld_solvent": 1.,
+      "radius": 120.},
+     [(0.006, 0.008), (0.06,0.08), (0.12, 0.16)],
+     [1.34836265e+04, 6.20114062e+00, 1.04733914e-01]],
+
+    # TODO: magnetism smoke test. Values not validated.
+    [dict(radius=120, sld_M0=4, sld_mphi=20, sld_mtheta=60,
+      up_frac_i=0.05, up_frac_f=0.1, up_theta=-15, up_phi=10),
+     [(0.0, 0.01), (0.0, 0.1), (0.0, 0.2)],
+     [20247.206006297125, 9.312720770235483, 0.15826993186001856]],
+
     #  But note P(Q) = F2/volume
     #  F and F^2 are "unscaled", with for  n <F F*>S(q) or for beta approx
     #          I(q) = n [<F F*> + <F><F*> (S(q) - 1)]
diff --git a/sasmodels/special.py b/sasmodels/special.py
index 1b8b34833..3b9632706 100644
--- a/sasmodels/special.py
+++ b/sasmodels/special.py
@@ -58,6 +58,9 @@
     cube(x):
         $x^3$
 
+    clip(a, a_min, a_max):
+        $\min(\max(a, a_\text{min}), a_\text{max})$, or NaN if $a$ is NaN.
+
     sas_sinx_x(x):
         $\sin(x)/x$, with limit $\sin(0)/0 = 1$.
 
@@ -215,7 +218,7 @@
 from numpy import sin, cos, tan, arcsin as asin, arccos as acos, arctan as atan
 from numpy import sinh, cosh, tanh, arcsinh as asinh, arccosh as acosh, arctanh as atanh
 from numpy import arctan2 as atan2
-from numpy import fabs, fmin, fmax, trunc, rint
+from numpy import fabs, fmin, fmax, clip, trunc, rint
 from numpy import pi, nan, inf
 from scipy.special import gamma as sas_gamma
 from scipy.special import gammaln as sas_gammaln


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index a07707134..4953c9e45 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -29,11 +29,7 @@ jobs:
     - name: setup apt dependencies for Linux
       if: ${{ matrix.os == 'ubuntu-latest' }}
       run: |
-        wget https://repo.radeon.com/amdgpu-install/6.1.1/ubuntu/jammy/amdgpu-install_6.1.60101-1_all.deb
-        sudo apt install ./amdgpu-install_6.1.60101-1_all.deb
         sudo apt update
-        sudo apt install amdgpu-dkms
-        sudo apt install rocm
 
     - name: Install Python dependencies
       run: |
@@ -56,9 +52,11 @@ jobs:
         choco install opencl-intel-cpu-runtime
         python -m pip install --only-binary=pyopencl --find-links http://www.silx.org/pub/wheelhouse/ --trusted-host www.silx.org pyopencl
 
-    - name: Test with pytest
+    - name: Test with pytest (only on Windows for now since PoCL is failing on Ubuntu)
+
       env:
         PYOPENCL_COMPILER_OUTPUT: 1
+        SAS_OPENCL: none
       run: |
         # other CI uses the following, but `setup.py test` is a deprecated way
         # of running tests
@@ -68,6 +66,8 @@ jobs:
 
     - name: check that the docs build (linux only)
       if: ${{ matrix.os == 'ubuntu-latest' }}
+      env:
+        SAS_OPENCL: none
       run: |
         make -j 4 -C doc SPHINXOPTS="-W --keep-going -n" html