Have added the magnetic_functions.c, micromagnetic_FF_3D.c, and micro…

…magnetic_FF_3D.py from the marketplace. Unsure how to insert the parameter descriptions given in micromagnetic_FF_3D.py to the GUI according to sasview issue #2784

sasmodels/models/magnetic_functions.c

@@ -0,0 +1,193 @@
+static double clipp(double value, double low, double high) //from kernel_iq.c
+{
+  return (value < low ? low : (value > high ? high : value));
+} 
+
+static double length(double x, double y) 
+{
+  return sqrt(x*x + y*y);
+} 
+
+static double fq_core_shell(double q, double sld_core, double radius,
+   double sld_solvent, double fp_n, double sld[], double thickness[])
+{
+  const int n = (int)(fp_n+0.5);
+  double f, r, last_sld;
+  r = radius;
+  last_sld = sld_core;
+  f = 0.;
+  for (int i=0; i<n; i++) {
+    f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sas_3j1x_x(q*r);
+    last_sld = sld[i];
+    r += thickness[i];
+  }
+  f += M_4PI_3 * cube(r) * (sld_solvent - last_sld) * sas_3j1x_x(q*r);
+  return f;
+}
+
+static double langevin(
+    double x) {
+    // cotanh(x)-1\x
+
+    if (x < 0.00001) {
+        // avoid dividing by zero
+        return 1.0/3.0*x;
+    } 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;
+    } 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 = clipp(sqrt(square(in_spin)), 0.0, 1.0);//opencl has ambiguities for abs()
+   out_spin = clipp(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
+ }
+
+
+
+//some basic 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 v0, double v1, double v2)
+{
+ double vec[3];
+ SET_VEC(vec, v0, v1, v2);   
+ return sqrt(SCALAR_VEC(vec,vec));
+}
+
+void ORTH_VEC(double *result_vec, double *vec1, double *vec2)
+{
+ result_vec[0] = vec1[0] - SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2) * vec2[0];
+ result_vec[1] = vec1[1] - SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2) * vec2[1];
+ result_vec[2] = vec1[2] - SCALAR_VEC(vec1,vec2) / SCALAR_VEC(vec2,vec2) * vec2[2];
+}
+
+//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 x, double y, double z,
+  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(x*x + y*y + z*z);
+  SET_VEC(vector, x/q, y/q, z/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
+
+ }

sasmodels/models/micromagnetic_FF_3D.c

@@ -0,0 +1,200 @@
+//Core-shell form factor for anisotropy field (Hkx, Hky and Hkz), Nuc and lonitudinal magnetization Mz
+static double
+form_volume(double radius, double thickness)
+{
+ return M_4PI_3 * cube(radius + thickness);
+}
+
+static double
+radius_effective(int mode, double radius, double thickness)
+{
+ switch (mode) {
+    default:
+      case 1: // outer radius
+        return radius + thickness;
+      case 2: // core radius
+        return radius;
+ }
+}
+
+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)
+{
+  if( Hi > 1.0e-6 ) 
+    return Ms / (Hi + 2.0 * A * 4.0 * M_PI / Ms * q * q * 10.0);//q in 10e10 m-1, A in 10e-12 J/m, mu0 in 1e-7 
+  else 
+    return Ms / (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 x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  const double q=MAG_VEC(x, y, z);  
+  const double f = reduced_field(q, Ms, Hi, A)*(Hkx*(1.0+reduced_field(q, Ms, Hi, A)*y*y/q/q)-Ms*Mz*x*z/q/q*(1.0+reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hky*reduced_field(q, Ms, Hi, A)*x*y/q/q)/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/q/q-square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
+  return f;
+}
+
+static double fqMximag(double x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  const double q=MAG_VEC(x, y, z);      
+  const double f = -reduced_field(q, Ms, Hi, A)*(Ms*Mz*(1.0+reduced_field(q, Ms, Hi, A))*DMI_length(Ms, D,y)+Hky*reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z))/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/q/q -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
+  return f;
+}
+
+static double fqMyreal( double x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  const double q=MAG_VEC(x, y, z);  
+  const double f = reduced_field(q, Ms, Hi, A)*(Hky*(1.0+reduced_field(q, Ms, Hi, A)*x*x/q/q)-Ms*Mz*y*z/q/q*(1.0+reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hkx*reduced_field(q, Ms, Hi, A)*x*y/q/q)/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/q/q -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
+  return f;
+}
+
+static double fqMyimag( double x, double y, double z, double Mz, double Hkx, double Hky, double Hi, double Ms, double A, double D)
+{
+  const double q=MAG_VEC(x, y, z);  
+  const double f = reduced_field(q, Ms, Hi, A)*(Ms*Mz*(1.0+reduced_field(q, Ms, Hi, A))*DMI_length(Ms, D,x)-Hkx*reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z))/(1.0+reduced_field(q, Ms, Hi, A)*(x*x+y*y)/q/q -square(reduced_field(q, Ms, Hi, A)*DMI_length(Ms, D,z)));
+  return f;
+}
+
+//calculate 2D from _fq
+static double
+Iqxy(double qx, double qy, double radius, double thickness,double core_nuc, double shell_nuc, double solvent_nuc, double core_Ms, double shell_Ms, double solvent_Ms, double core_hk,  double Hi, double Ms, double A, double D,  double up_i, double up_f, double alpha, double beta)
+{
+  const double q=MAG_VEC(qx, qy, 0);    
+  if (q > 1.0e-16 ) {
+    const double cos_theta=qx/q;
+    const double sin_theta=qy/q;
+
+    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, core_Ms, shell_Ms, solvent_Ms);
+    double nuc=fq(q, radius, thickness, core_nuc, shell_nuc, solvent_nuc);
+
+    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= fq(q, radius, thickness, core_hk, 0, 0)*sin_gamma;
+      double Hky= fq(q, radius, thickness, core_hk, 0, 0)*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 0.5*1.0e-4*total_F2;
+  }  
+}
+
+static double
+Iq(double q, double radius, double thickness,double core_nuc, double shell_nuc, double solvent_nuc, double core_Ms, double shell_Ms, double solvent_Ms, double core_hk, 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 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, 5=du.imag,  6=ud.real, 7=ud.imag
+    double weights[8];  
+    set_weights(up_i, up_f, weights);
+
+    double mz=fq(q, radius, thickness, core_Ms, shell_Ms, solvent_Ms);
+    double nuc=fq(q, radius, thickness, core_nuc, shell_nuc, solvent_nuc);
+
+    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= fq(q, radius, thickness, core_hk, 0, 0)*sin_gamma;
+      double Hky= fq(q, radius, thickness, core_hk, 0, 0)*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 ;
+    }
+    total_F1D += GAUSS_W[j] * total_F2 ;
+  }
+  //convert from [1e-12 A-1] to [cm-1]
+  return 0.25*1.0e-4*total_F1D;
+}
+
+
+
+
+
+