Changeset f12357f in sasmodels for sasmodels/models/pearl_necklace.c

Timestamp:

Feb 3, 2016 6:49:57 AM (8 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.

File:

• sasmodels/models/pearl_necklace.c (modified) (5 diffs)

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;
-}