Changeset f12357f in sasmodels

Timestamp:

Feb 3, 2016 6:49:57 AM (9 years ago)

Author:

krzywon

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

a17fe40

Parents:

8696a87

Message:

Modified the pearl_necklace model for better double precision accuracy
and added unit tests.

Location:

sasmodels/models

Files:

• lib/Si.c (modified) (2 diffs)
• pearl_necklace.c (modified) (5 diffs)
• pearl_necklace.py (modified) (3 diffs)

sasmodels/models/lib/Si.c

@@ -13 +13 @@
 
                 // Explicitly writing factorial values triples the speed of the calculation
-                out_cos = 1/x - 2/pow(x,3) + 24/pow(x,5) - 720/pow(x,7) + 40320/pow(x,9);
-                out_sin = 1/x - 6/pow(x,4) + 120/pow(x,6) - 5040/pow(x,8) + 362880/pow(x,10);
+                out_cos = 1./x - 2./pow(x,3) + 24./pow(x,5) - 720./pow(x,7);
+                out_sin = 1./pow(x,2) - 6./pow(x,4) + 120./pow(x,6) - 5040./pow(x,8);
 
                 out -= cos(x) * out_cos;
@@ -22 +22 @@
 
         // Explicitly writing factorial values triples the speed of the calculation
-        out = x - pow(x, 3)/18 + pow(x,5)/600 - pow(x,7)/35280 + pow(x,9)/3265920;
+        out = x - pow(x, 3)/18. + pow(x,5)/600. - pow(x,7)/35280. + pow(x,9)/3265920. - pow(x,11)/439084800.;
 
-        //printf ("Si=%g %g\n", x, out);
         return out;
 }

sasmodels/models/pearl_necklace.c

@@ -12 +12 @@
         double string_thickness, double number_of_pearls, double sld, 
         double string_sld, double solvent_sld);
-        
-double ER(double radius, double edge_separation,
-        double string_thickness, double number_of_pearls);
-double VR(double radius, double edge_separation,
-        double string_thickness, double number_of_pearls);
 
 // From Igor library
@@ -22 +17 @@
         double num_pearls, double sld_pearl, double sld_string, double sld_solv)
 {
+        //relative slds
         double contrast_pearl = sld_pearl - sld_solv;
         double contrast_string = sld_string - sld_solv;
+        
+        // number of string segments
+        num_pearls = floor(num_pearls + 0.5); //Force integer number of pearls
+        double num_strings = num_pearls - 1.0;
+        
+        //Pi
+        double pi = 4.0*atan(1.0);
+        
+        // center to center distance between the neighboring pearls
+        double A_s = edge_separation + 2.0 * radius;
+        
+        // Repeated Calculations
+        double sincasq = sinc(q*A_s);
+        double oneminussinc = 1 - sincasq;
+        double q_r = q * radius;
+        double q_edge = q * edge_separation;
+        
+        // each volume
+        double string_vol = edge_separation * pi * thick_string * thick_string / 4.0;
+        double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius;
 
         //total volume
-        double pi = 4.0*atan(1.0);
-        double tot_vol = form_volume(radius, edge_separation, thick_string, num_pearls);
-        double string_vol = edge_separation * pi * pow((thick_string / 2.0), 2);
-        double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3);
-        double num_strings = num_pearls - 1;
-        //each mass
+        double tot_vol;
+        //each masses
         double m_r= contrast_string * string_vol;
         double m_s= contrast_pearl * pearl_vol;
         double psi, gamma, beta;
         //form factors
-        double sss; //pearls
-        double srr; //strings
-        double srs; //cross
-        double A_s, srr_1, srr_2, srr_3;
+        double sss, srr, srs; //cross
+        double srr_1, srr_2, srr_3;
         double form_factor;
+        tot_vol = num_strings * string_vol;
+        tot_vol += num_pearls * pearl_vol;
 
         //sine functions of a pearl
-        psi = sin(q*radius);
-        psi -= q * radius * cos(q * radius);
-        psi /= pow((q * radius), 3);
+        psi = sin(q_r);
+        psi -= q_r * cos(q_r);
         psi *= 3.0;
+        psi /= q_r * q_r * q_r;
 
         // Note take only 20 terms in Si series: 10 terms may be enough though.
-        gamma = Si(q* edge_separation);
-        gamma /= (q* edge_separation);
-        beta = Si(q * (edge_separation + radius));
-        beta -= Si(q * radius);
-        beta /= (q* edge_separation);
-
-        // center to center distance between the neighboring pearls
-        A_s = edge_separation + 2.0 * radius;
+        gamma = Si(q_edge);
+        gamma /= (q_edge);
+        beta = Si(q * (A_s - radius));
+        beta -= Si(q_r);
+        beta /= q_edge;
 
         // form factor for num_pearls
-        sss = 1.0 - pow(sinc(q*A_s), num_pearls );
-        sss /= pow((1.0-sinc(q*A_s)), 2);
-        sss *= -sinc(q*A_s);
-        sss -= num_pearls/2.0;
-        sss += num_pearls/(1.0-sinc(q*A_s));
-        sss *= 2.0 * pow((m_s*psi), 2);
+        sss = 1.0 - pow(sincasq, num_pearls);
+        sss /= oneminussinc * oneminussinc;
+        sss *= -sincasq;
+        sss -= num_pearls / 2.0;
+        sss += num_pearls / oneminussinc;
+        sss *= 2.0 * m_s * psi * m_s * psi;
 
         // form factor for num_strings (like thin rods)
-        srr_1 = -pow(sinc(q*edge_separation/2.0), 2);
+        srr_1 = -sinc(q_edge/2.0) * sinc(q_edge/2.0);
 
         srr_1 += 2.0 * gamma;
         srr_1 *= num_strings;
-        srr_2 = 2.0/(1.0-sinc(q*A_s));
-        srr_2 *= num_strings * pow(beta, 2);
-        srr_3 = 1.0 - pow(sinc(q*A_s), num_strings);
-        srr_3 /= pow((1.0-sinc(q*A_s)), 2);
-        srr_3 *= -2.0 * pow(beta, 2);
+        srr_2 = 2.0/oneminussinc;
+        srr_2 *= num_strings;
+        srr_2 *= beta * beta;
+        srr_3 = 1.0 - pow(sincasq, num_strings);
+        srr_3 /= oneminussinc * oneminussinc;
+        srr_3 *= beta * beta;
+        srr_3 *= -2.0;
 
         // total srr
         srr = srr_1 + srr_2 + srr_3;
-        srr *= pow(m_r, 2);
+        srr *= m_r * m_r;
 
         // form factor for correlations
         srs = 1.0;
-        srs -= pow(sinc(q*A_s), num_strings);
-        srs /= pow((1.0-sinc(q*A_s)), 2);
-        srs *= -sinc(q*A_s);
-        srs += (num_strings/(1.0-sinc(q*A_s)));
-        srs *= 4.0 * (m_r * m_s * beta * psi);
+        srs -= pow(sincasq, num_strings);
+        srs /= oneminussinc * oneminussinc;
+        srs *= -sincasq;
+        srs += num_strings/oneminussinc;
+        srs *= 4.0;
+        srs *= (m_r * m_s * beta * psi);
 
         form_factor = sss + srr + srs;
@@ -103 +115 @@
         double number_of_strings = number_of_pearls - 1.0;
         
-        double string_vol = edge_separation * pi * pow((string_thickness / 2.0), 2);
-        double pearl_vol = 4.0 /3.0 * pi * pow(radius, 3);
+        double string_vol = edge_separation * pi * string_thickness * string_thickness / 4.0;
+        double pearl_vol = 4.0 / 3.0 * pi * radius * radius * radius;
 
         total_vol = number_of_strings * string_vol;
@@ -124 +136 @@
         double string_sld, double solvent_sld)
 {
-        double value = 0.0;
-        double tot_vol = 0.0;
+        double value, tot_vol;
+        
+        if (string_thickness >= radius || number_of_pearls <= 0) {
+                return NAN;
+        }
         
         value = _pearl_necklace_kernel(q, radius, edge_separation, string_thickness,
@@ -142 +157 @@
                 sld, string_sld, solvent_sld));
 }
-
-
-double ER(double radius, double edge_separation,
-        double string_thickness, double number_of_pearls)
-{
-        double tot_vol = form_volume(radius, edge_separation, string_thickness, number_of_pearls);
-        double pi = 4.0*atan(1.0);
-
-    double rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333);
-    
-    return rad_out;
-}
-double VR(double radius, double edge_separation,
-        double string_thickness, double number_of_pearls)
-{
-        return 1.0;
-}

sasmodels/models/pearl_necklace.py

@@ -58 +58 @@
 """
 
-from numpy import inf
+from numpy import inf, pi
 
 name = "pearl_necklace"
@@ -96 +96 @@
 
 source = ["lib/Si.c", "pearl_necklace.c"]
+# new flag to let the compiler know to never use single precision
+single = False
+
+def volume(radius, edge_separation, string_thickness, number_of_pearls):
+    """
+    Calculates the total particle volume of the necklace.
+    Redundant with form_volume.
+    """
+    number_of_strings = number_of_pearls - 1.0
+    string_vol = edge_separation * pi * pow((string_thickness / 2.0), 2.0)
+    pearl_vol = 4.0 /3.0 * pi * pow(radius, 3.0)
+    total_vol = number_of_strings * string_vol
+    total_vol += number_of_pearls * pearl_vol
+    return total_vol
+
+def ER(radius, edge_separation, string_thickness, number_of_pearls):
+    """
+    Calculation for effective radius.
+    """
+    tot_vol = volume(radius, edge_separation, string_thickness, number_of_pearls)
+    rad_out = pow((3.0*tot_vol/4.0/pi), 0.33333)
+    return rad_out
 
 # parameters for demo
@@ -114 +136 @@
                string_thickness='thick_string', sld='sld_pearl',
                string_sld='sld_string', edge_separation='edge_separation')
+
+tests = [[{}, 0.001, 17380.245], [{}, 'ER', 115.39502]]