Changeset 0da2a01 in sasmodels for sasmodels/models/lib

Timestamp:

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

7d256c8

Parents:

0420af7 (diff), 790bcc4c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge remote-tracking branch 'origin/master'

Location:

sasmodels/models/lib

Files:

• J1c.c (deleted)
• sph_j1c.c (added)
• Si.c (modified) (3 diffs)

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