Changeset f728001 in sasmodels

Timestamp:

Oct 23, 2017 12:34:11 AM (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:

ea60e08

Parents:

31eea1f

Message:

paracrystals: document relation to equation numbers in Matsuoka

Files:

• explore/angles.py (modified) (1 diff)
• sasmodels/models/bcc_paracrystal.c (modified) (3 diffs)
• sasmodels/models/fcc_paracrystal.c (modified) (1 diff)
• sasmodels/models/sc_paracrystal.c (modified) (1 diff)

explore/angles.py

@@ -38 +38 @@
 
 # include unicode symbols in output, even if piping to a pager
-sys.stdout = codecs.getwriter('utf8')(sys.stdout)
+if sys.version_info[0] < 3:
+    sys.stdout = codecs.getwriter('utf8')(sys.stdout)
 sp.init_printing(use_unicode=True)
 

sasmodels/models/bcc_paracrystal.c

@@ -2 +2 @@
 bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
-#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;
+    // Equations from Matsuoka 26-27-28, multiplied by |q|
+    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
+    const double a3 = (+qa + qb - qc)/2.0;
 
-#if 1
-    // Numerator: (1 - exp(a)^2)^3
-    //         => (-(exp(2a) - 1))^3
-    //         => -expm1(2a)^3
-    // 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)
+#if 0
+    // Matsuoka 29-30-31
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (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);
@@ -25 +23 @@
           * ((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);
+
+#elif 0
+    // ** Alternate form, which perhaps is more approachable
+    //     Z_k numerator   => -[(exp(2a) - 1) / 2.exp(a)] 2.exp(a)
+    //                     => -[sinh(a)] exp(a)
+    //     Z_k denominator => [(exp(2a) + 1) / 2.exp(a) - cos(d a_k)] 2.exp(a)
+    //                     => [cosh(a) - cos(d a_k)] 2.exp(a)
+    //     => Z_k = -sinh(a) / [cosh(a) - cos(d a_k)]
+    //            = sinh(-a) / [cosh(-a) - cos(d a_k)]
+    //
+    // One more step leads to the form in sasview 3.x for 2d models
+    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)]
+    //
+    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);
@@ -33 +42 @@
                     * sinh_qd/(cosh_qd - cos(dnn*a2))
                     * sinh_qd/(cosh_qd - cos(dnn*a3));
+#else
+    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double tanh_qd = tanh(arg);
+    const double cosh_qd = cosh(arg);
+    const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd);
 #endif
 

sasmodels/models/fcc_paracrystal.c

@@ -2 +2 @@
 fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
-#if 0  // Equations as written in Matsuoka
-    const double a1 = ( qa + qb)/2.0;
-    const double a2 = (-qa + qc)/2.0;
-    const double a3 = (-qa + qb)/2.0;
-#else
+    // Equations from Matsuoka 17-18-19, multiplied by |q|
     const double a1 = ( qa + qb)/2.0;
     const double a2 = ( qa + qc)/2.0;
     const double a3 = ( qb + qc)/2.0;
-#endif
 
-    // Numerator: (1 - exp(a)^2)^3
-    //         => (-(exp(2a) - 1))^3
-    //         => -expm1(2a)^3
-    // 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)
+    // Matsuoka 23-24-25
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (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);

sasmodels/models/sc_paracrystal.c

@@ -2 +2 @@
 sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
+    // Equations from Matsuoka 9-10-11, multiplied by |q|
     const double a1 = qa;
     const double a2 = qb;
     const double a3 = qc;
 
-    // Numerator: (1 - exp(a)^2)^3
-    //         => (-(exp(2a) - 1))^3
-    //         => -expm1(2a)^3
-    // 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)
+    // Matsuoka 13-14-15
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (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);