Changeset 0420af7 in sasmodels

Timestamp:

Jan 21, 2016 4:17:17 PM (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:

0da2a01

Parents:

6458608

Message:

Modified hollow_cylinder to reject invalid inputs, added VR, ER and
tests. Si now has pre-calculated factorials for faster processing.

Location:

sasmodels/models

Files:

• hollow_cylinder.c (modified) (2 diffs)
• hollow_cylinder.py (modified) (4 diffs)
• lib/Si.c (modified) (3 diffs)

sasmodels/models/hollow_cylinder.c

@@ -114 +114 @@
         double norm,volume;     //final calculation variables
         
+        if (core_radius >= radius || radius >= length) {
+                return NAN;
+        }
+        
         delrho = solvent_sld - sld;
         lower = 0.0;
@@ -132 +136 @@
 }
 
-//FIXME: Factor of two difference
 double Iqxy(double qx, double qy, double radius, double core_radius, double length, double sld,
         double solvent_sld, double theta, double phi)

sasmodels/models/hollow_cylinder.py

@@ -64 +64 @@
 """
 
-from numpy import inf
+from numpy import pi, inf
 
 name = "hollow_cylinder"
@@ -92 +92 @@
 source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"]
 
+def ER(radius, core_radius, length):
+    if radius == 0 or length == 0:
+        return 0.0
+    len1 = radius
+    len2 = length/2.0
+    term1 = len1*len1*2.0*len2/2.0
+    term2 = 1.0 + (len2/len1)*(1.0 + 1/len2/2.0)*(1.0 + pi*len1/len2/2.0)
+    ddd = 3.0*term1*term2
+    diam = pow(ddd, (1.0/3.0))
+    return diam
+
+def VR(radius, core_radius, length):
+    vol_core = pi*core_radius*core_radius*length
+    vol_total = pi*radius*radius*length
+    vol_shell = vol_total - vol_core
+    return vol_shell, vol_total
+
 # parameters for demo
 demo = dict(scale=1.0,background=0.0,length=400.0,radius=30.0,core_radius=20.0,
@@ -97 +114 @@
             radius_pd=.2, radius_pd_n=9,
             length_pd=.2, length_pd_n=10,
+            core_radius_pd=.2, core_radius_pd_n=9,
             theta_pd=10, theta_pd_n=5,
             )
@@ -109 +127 @@
 # Parameters for unit tests
 tests = [
-         [{"radius" : 30.0},0.00005,1764.926]
+         [{"radius" : 30.0},0.00005,1764.926],
+         [{},'VR',1.8],
+         [{},0.001,1756.76]
          ]

sasmodels/models/lib/Si.c

@@ -1 @@
-int factorial(int f);
 double Si(double x);
 
-// integral of sin(x)/x: approximated to w/i 1%
+// integral of sin(x)/x Taylor series approximated to w/i 0.1%
 double Si(double x)
 {
         int i;
-        int nmax=6;
+        int nmax=10;
         double out;
         long power;
@@ -16 +15 @@
                 out = pi/2.0;
 
-                for (i=0; i<nmax-2; i+=1) {
-                        out_cos += pow(-1.0, i) * (double)factorial(2*i) / pow(x, 2*i+1);
-                        out_sin += pow(-1.0, i) * (double)factorial(2*i+1) / pow(x, 2*i+2);
-                }
+                // 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(x) * out_cos;
@@ -26 +24 @@
         }
 
-        out = 0.0;
+        // 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;
 
-        for (i=0; i<nmax; i+=1) {
-                if (i==0) {
-                        out += x;
-                        continue;
-                }
-
-                power = pow(x,(2 * i + 1));
-                out += pow(-1.0, i) * power / ((2.0 * (double)i + 1.0) * (double)factorial(2 * i + 1));
-
-                //printf ("Si=%g %g %d\n", x, out, i);
-        }
-
+        //printf ("Si=%g %g\n", x, out);
         return out;
 }
-
-int factorial(int f)
-{
-    if ( f == 0 ) 
-        return 1;
-    return(f * factorial(f - 1));
-}