Changeset 7e0b281 in sasmodels for sasmodels/models/bcc_paracrystal.c

Timestamp:

Apr 17, 2017 6:34:52 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

64ca163

Parents:

cb038a2

Message:

adjust paracrystal equations to better match terminology in Matsuoka paper

File:

• sasmodels/models/bcc_paracrystal.c (modified) (6 diffs)

sasmodels/models/bcc_paracrystal.c

@@ -1 +1 @@
 static double
-_sq_bcc(double qa, double qb, double qc, double dnn, double d_factor)
+bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
-    // Rewriting equations for efficiency, accuracy and readability, and so
-    // code is reusable between 1D and 2D models.
-    const double a1 = +qa - qc + qb;
-    const double a2 = +qa + qc - qb;
-    const double a3 = -qa + qc + qb;
-
-    const double half_dnn = 0.5*dnn;
-    const double arg = 0.5*square(half_dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+#if 0  // Equations as written in Matsuoka
+    const double a1 = (+qa + qb + qc)/2.0;
+    const double a2 = (-qa - qb + qc)/2.0;
+    const double a3 = (-qa + qb - qc)/2.0;
+#else
+    const double a1 = (+qa + qb - qc)/2.0;
+    const double a2 = (+qa - qb + qc)/2.0;
+    const double a3 = (-qa + qb + qc)/2.0;
+#endif
 
 #if 1
@@ -15 +16 @@
     //         => (-(exp(2a) - 1))^3
     //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(xk) exp(a) + exp(a)^2)
-    //         => exp(a)^2 - 2 cos(xk) exp(a) + 1
-    //         => (exp(a) - 2 cos(xk)) * exp(a) + 1
-    const double exp_arg = exp(-arg);
-    const double Sq = -cube(expm1(-2.0*arg))
-        / ( ((exp_arg - 2.0*cos(half_dnn*a1))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(half_dnn*a2))*exp_arg + 1.0)
-          * ((exp_arg - 2.0*cos(half_dnn*a3))*exp_arg + 1.0));
+    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
+    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double exp_arg = exp(arg);
+    const double Zq = -cube(expm1(2.0*arg))
+        / ( ((exp_arg - 2.0*cos(dnn*a1))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
+          * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
 #else
     // Alternate form, which perhaps is more approachable
+    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
     const double sinh_qd = sinh(arg);
     const double cosh_qd = cosh(arg);
-    const double Sq = sinh_qd/(cosh_qd - cos(half_dnn*a1))
-                    * sinh_qd/(cosh_qd - cos(half_dnn*a2))
-                    * sinh_qd/(cosh_qd - cos(half_dnn*a3));
+    const double Zq = sinh_qd/(cosh_qd - cos(dnn*a1))
+                    * sinh_qd/(cosh_qd - cos(dnn*a2))
+                    * sinh_qd/(cosh_qd - cos(dnn*a3));
 #endif
 
-    return Sq;
+    return Zq;
 }
 
@@ -38 +41 @@
 // occupied volume fraction calculated from lattice symmetry and sphere radius
 static double
-_bcc_volume_fraction(double radius, double dnn)
+bcc_volume_fraction(double radius, double dnn)
 {
     return 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
@@ -75 +78 @@
             const double qa = qab*cos_phi;
             const double qb = qab*sin_phi;
-            const double fq = _sq_bcc(qa, qb, qc, dnn, d_factor);
-            inner_sum += Gauss150Wt[j] * fq;
+            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor);
+            inner_sum += Gauss150Wt[j] * form;
         }
         inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx
@@ -82 +85 @@
     }
     outer_sum *= theta_m;
-    const double Sq = outer_sum/(4.0*M_PI);
+    const double Zq = outer_sum/(4.0*M_PI);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-
-    return _bcc_volume_fraction(radius, dnn) * Pq * Sq;
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }
 
@@ -101 +103 @@
 
     q = sqrt(qa*qa + qb*qb + qc*qc);
+    const double Zq = bcc_Zq(qa, qb, qc, dnn, d_factor);
     const double Pq = sphere_form(q, radius, sld, solvent_sld);
-    const double Sq = _sq_bcc(qa, qb, qc, dnn, d_factor);
-    return _bcc_volume_fraction(radius, dnn) * Pq * Sq;
+    return bcc_volume_fraction(radius, dnn) * Pq * Zq;
 }