refactor model

sasmodels/models/micromagnetic_FF_3D.c

@@ -1,11 +1,11 @@
 //Core-shell form factor for anisotropy field (Hkx, Hky and Hkz), Nuc and lonitudinal magnetization Mz
 static double
-form_volume(double nuc_radius, double nuc_thickness, double mag_radius, double mag_thickness, double hk_radius)
+form_volume(double radius, double thickness)
 {
-  if( nuc_radius+nuc_thickness>mag_radius+ mag_thickness)
-    return M_4PI_3 * cube(nuc_radius + nuc_thickness);
+  if( radius+thickness>radius+ thickness)
+    return M_4PI_3 * cube(radius + thickness);
   else
-    return M_4PI_3 * cube(mag_radius + mag_thickness);
+    return M_4PI_3 * cube(radius + thickness);
 
 }
 
@@ -27,10 +27,9 @@ static double fq(double q, double radius,
 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); 
+    // 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)
@@ -51,7 +50,9 @@ static double fqMxreal( double qx, double qy, double qz, double Mz, double Hkx,
   const double qsq = qx*qx + qy*qy + qz*qz;
   const double q = sqrt(qsq); 
   const double Hr = reduced_field(q, Ms, Hi, A);
-  const double f = Hr*(Hkx*(qsq+Hr*qy*qy)-Ms*Mz*qx*qz*(1.0+Hr*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hky*Hr*qx*qy)/(qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI_length(Ms, D,qz)*q));
+  const double DMI = DMI_length(Ms, D,q);
+  const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
+  const double f = (Hkx*(qsq+Hr*qy*qy)-Ms*Mz*qx*qz*(1.0+Hr*square(DMI))-Hky*Hr*qx*qy)/denominator;
   return f;
 }
 
@@ -60,7 +61,9 @@ static double fqMximag(double qx, double qy, double qz, double Mz, double Hkx, d
   const double qsq = qx*qx + qy*qy + qz*qz;
   const double q = sqrt(qsq);   
   double Hr = reduced_field(q, Ms, Hi, A);
-  const double f = -Hr*qsq*(Ms*Mz*(1.0+Hr)*DMI_length(Ms, D,qy)+Hky*Hr*DMI_length(Ms, D,qz))/(qsq+Hr*(qx*qx+qy*qy) -square(Hr*DMI_length(Ms, D,qz)*q));
+  const double DMI = DMI_length(Ms, D,q);
+  const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
+  const double f = -qsq*(Ms*Mz*(1.0+Hr)*DMI+Hky*Hr*DMI)/denominator;
   return f;
 }
 
@@ -69,7 +72,9 @@ static double fqMyreal( double qx, double qy, double qz, double Mz, double Hkx,
   const double qsq = qx*qx + qy*qy + qz*qz;
   const double q = sqrt(qsq); 
   const double Hr = reduced_field(q, Ms, Hi, A);
-  const double f = Hr*(Hky*(qsq+Hr*qx*qx)-Ms*Mz*qy*qz*(1.0+Hr*DMI_length(Ms, D,q)*DMI_length(Ms, D,q))-Hkx*Hr*qx*qy)/(qsq+Hr*(qx*qx+qy*qy) -square(Hr*DMI_length(Ms, D,qz)*q));
+  const double DMI = DMI_length(Ms, D,q);
+  const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
+  const double f = (Hky*(qsq+Hr*qx*qx)-Ms*Mz*qy*qz*(1.0+Hr*square(DMI))-Hkx*Hr*qx*qy)/denominator;
   return f;
 }
 
@@ -78,27 +83,24 @@ static double fqMyimag( double qx, double qy, double qz, double Mz, double Hkx,
   const double qsq = qx*qx + qy*qy + qz*qz;
   const double q = sqrt(qsq);  
   const double Hr = reduced_field(q, Ms, Hi, A);
-  const double f = Hr*qsq*(Ms*Mz*(1.0+Hr)*DMI_length(Ms, D,qx)-Hkx*Hr*DMI_length(Ms, D,qz))/(qsq+Hr*(qx*qx+qy*qy) -square(Hr*DMI_length(Ms, D,qz)*q));
+  const double DMI = DMI_length(Ms, D,q);
+  const double denominator = (qsq+Hr*(qx*qx+qy*qy)-square(Hr*DMI*q))/Hr;
+  const double f = qsq*(Ms*Mz*(1.0+Hr)*DMI-Hkx*Hr*DMI)/denominator;
   return f;
 }
 
 static double
-Iqxy(double qx, double qy, double nuc_radius, double nuc_thickness, double mag_radius, double mag_thickness, double hk_radius, 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)
+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)
 {
-  const double q = sqrt(qx*qx + qy*qy);  
-  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, mag_radius, mag_thickness, mag_sld_core, mag_sld_shell, mag_sld_solvent);
-    double nuc=fq(q, nuc_radius, nuc_thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent);
-
+    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];
@@ -111,9 +113,9 @@ Iqxy(double qx, double qy, double nuc_radius, double nuc_thickness, double mag_r
       //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, hk_radius, 0, hk_sld_core, 0, 0)*sin_gamma;
-      double Hky= fq(q, hk_radius, 0, hk_sld_core, 0, 0)*cos_gamma;
-	  
+      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);
@@ -131,14 +133,32 @@ Iqxy(double qx, double qy, double nuc_radius, double nuc_thickness, double mag_r
         }
       }
       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)
+{
+  const double q = sqrt(qx*qx + qy*qy);  
+  if (q > 1.0e-16 ) {
+    const double cos_theta=qx/q;
+    const double sin_theta=qy/q;
+  } else {
+    const double cos_theta=0.0;
+    const double 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 nuc_radius, double nuc_thickness, double mag_radius, double mag_thickness, double hk_radius,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)
+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; 
@@ -148,49 +168,7 @@ Iq(double q, double nuc_radius, double nuc_thickness, double mag_radius, double
     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, mag_radius, mag_thickness, mag_sld_core, mag_sld_shell, mag_sld_solvent);
-    double nuc=fq(q, nuc_radius, nuc_thickness, nuc_sld_core, nuc_sld_shell, nuc_sld_solvent);
-
-
-    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, hk_radius, 0, hk_sld_core, 0, 0)*sin_gamma;
-      double Hky= fq(q, hk_radius, 0, hk_sld_core, 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 ;
-    }
+    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]

sasmodels/models/micromagnetic_FF_3D.py

@@ -136,11 +136,8 @@
 
 # pylint: disable=bad-whitespace, line-too-long
 #             ["name", "units", default, [lower, upper], "type","description"],
-parameters = [["nuc_radius", "Ang", 50., [0, inf], "volume", "Structural radius of the core"],
-              ["nuc_thickness", "Ang", 40., [0, inf], "volume", "Structural thickness of shell"],  
-              ["mag_radius", "Ang", 50., [0, inf], "volume", "Magnetic radius of the core"],
-              ["mag_thickness", "Ang", 40., [0, inf], "volume", "Magnetic thickness of shell"], 
-              ["hk_radius", "Ang", 50., [0, inf], "volume", "Anisotropy radius of the core"],            
+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"],
@@ -177,11 +174,8 @@ def random():
     radius = np.random.beta(0.5, 0.5)*(outer_radius-2) + 1
     thickness = outer_radius - radius
     pars = dict(
-        nuc_radius=radius,
-        nuc_thickness=thickness,
-        mag_radius=radius,
-        mag_thickness=thickness,
-        hk_radius=radius
+        radius=radius,
+        thickness=thickness,
 
     )
     return pars