make vector calculations accessible for models

sasmodels/kernel_header.c

@@ -194,6 +194,47 @@ inline double cube(double x) { return x*x*x; }
 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

sasmodels/kernel_iq.c

@@ -70,45 +70,7 @@ typedef union {
 #if defined(MAGNETIC) && NUM_MAGNETIC > 0
 // ===== Helper functions for magnetism =====
 
-// 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];
-}
 
 
 // Compute spin cross sections given in_spin and out_spin

sasmodels/models/lib/magnetic_functions.c

@@ -1,19 +1,6 @@
-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;
-}
+//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) {
@@ -40,6 +27,8 @@ static double langevinoverx(
  }
 }
 
+
+
 //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
 {
@@ -70,48 +59,6 @@ static void set_weights(double in_spin, double out_spin, double weight[8]) //fro
 
 
 
-//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,

sasmodels/models/micromagnetic_FF_3D.c

@@ -48,36 +48,36 @@ static double reduced_field(double q, double Ms, double Hi,
 
 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 q = sqrt(x*x + y*y + z*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 q = sqrt(x*x + y*y + z*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 q = sqrt(x*x + y*y + z*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 q = sqrt(x*x + y*y + z*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;
 }
 
 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)
 {
- const double q=MAG_VEC(qx, qy, 0);  
+ 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;

sasmodels/models/micromagnetic_FF_3D.py

@@ -12,27 +12,24 @@
 saturation the scattering cross-section can be evaluated by means of micromagnetic theory
 
 .. math::
- I(\mathbf{Q}) = I_{nuc} +I_{res}+ I_{mag}(\mathbf{Q},H),
+ I(\mathbf{Q}) = I_{nuc} + I_{mag}(\mathbf{Q},H),
 
-with the field-independent nuclear and magnetic residual SANS cross section (due
-to nanoscale spatial variations of the saturation magnetisation) measured at
-complete magnetic saturation. Commonly, a measurement at a high magnetic field is
-taken as reference and subtracted from the scattering cross-section at lower
-field to obtain the so-called spin-misalignment scattering cross-section 
+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 purely magnetic scattering reflects the
-response of the transversal magnetisation components with to an
+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 alignment of the magnetic moments along
-the magnetic field is disturbed by perturbations in the microstructure. The
-anisotropy-field function $S_K$ depends on the Fourier transform of the magnetic
-anisotropy distribution (strength and orientation) in the material and the
+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