Changeset ed246ab in sasmodels for sasmodels/models/lib

Timestamp:

Apr 26, 2016 5:17:30 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

13b99fd

Parents:

639c4e3 (diff), 98cb4d7 (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 branch 'master' into polydisp

Location:

sasmodels/models/lib

Files:

• sas_gamma.c (modified) (2 diffs)
• librefl.c (added)
• sas_JN.c (modified) (1 diff)
• sph_j1c.c (modified) (2 diffs)
• sphere_form.c (modified) (1 diff)
• wrc_cyl.c (modified) (2 diffs)

sasmodels/models/lib/sas_gamma.c

@@ -1 +1 @@
 /*
 The wrapper for gamma function from OpenCL and standard libraries
-The OpenCL gamma function fails misserably on values lower than 1.0
+The OpenCL gamma function fails miserably on values lower than 1.0
 while works fine on larger values.
 We use gamma definition Gamma(t + 1) = t * Gamma(t) to compute
@@ -7 +7 @@
 */
 
+#if defined(NEED_TGAMMA)
+static double cephes_stirf(double x)
+{
+        const double MAXSTIR=143.01608;
+        const double SQTPI=2.50662827463100050242E0;
+        double y, w, v;
+
+        w = 1.0 / x;
+
+        w = ((((
+                7.87311395793093628397E-4*w
+                - 2.29549961613378126380E-4)*w
+                - 2.68132617805781232825E-3)*w
+                + 3.47222221605458667310E-3)*w
+                + 8.33333333333482257126E-2)*w
+                + 1.0;
+        y = exp(x);
+        if (x > MAXSTIR)
+        { /* Avoid overflow in pow() */
+                v = pow(x, 0.5 * x - 0.25);
+                y = v * (v / y);
+        }
+        else
+        {
+                y = pow(x, x - 0.5) / y;
+        }
+        y = SQTPI * y * w;
+        return(y);
+}
+
+static double tgamma(double x) {
+        double p, q, z;
+        int sgngam;
+        int i;
+
+        sgngam = 1;
+        if (isnan(x))
+                return(x);
+        q = fabs(x);
+
+        if (q > 33.0)
+        {
+                if (x < 0.0)
+                {
+                        p = floor(q);
+                        if (p == q)
+                        {
+                                return (NAN);
+                        }
+                        i = p;
+                        if ((i & 1) == 0)
+                                sgngam = -1;
+                        z = q - p;
+                        if (z > 0.5)
+                        {
+                                p += 1.0;
+                                z = q - p;
+                        }
+                        z = q * sin(M_PI * z);
+                        if (z == 0.0)
+                        {
+                                return(NAN);
+                        }
+                        z = fabs(z);
+                        z = M_PI / (z * cephes_stirf(q));
+                }
+                else
+                {
+                        z = cephes_stirf(x);
+                }
+                return(sgngam * z);
+        }
+
+        z = 1.0;
+        while (x >= 3.0)
+        {
+                x -= 1.0;
+                z *= x;
+        }
+
+        while (x < 0.0)
+        {
+                if (x > -1.E-9)
+                        goto small;
+                z /= x;
+                x += 1.0;
+        }
+
+        while (x < 2.0)
+        {
+                if (x < 1.e-9)
+                        goto small;
+                z /= x;
+                x += 1.0;
+        }
+
+        if (x == 2.0)
+                return(z);
+
+        x -= 2.0;
+        p = (((((
+                +1.60119522476751861407E-4*x
+                + 1.19135147006586384913E-3)*x
+                + 1.04213797561761569935E-2)*x
+                + 4.76367800457137231464E-2)*x
+                + 2.07448227648435975150E-1)*x
+                + 4.94214826801497100753E-1)*x
+                + 9.99999999999999996796E-1;
+        q = ((((((
+                -2.31581873324120129819E-5*x
+                + 5.39605580493303397842E-4)*x
+                - 4.45641913851797240494E-3)*x
+                + 1.18139785222060435552E-2)*x
+                + 3.58236398605498653373E-2)*x
+                - 2.34591795718243348568E-1)*x
+                + 7.14304917030273074085E-2)*x
+                + 1.00000000000000000320E0;
+        return(z * p / q);
+
+small:
+        if (x == 0.0)
+        {
+                return (NAN);
+        }
+        else
+                return(z / ((1.0 + 0.5772156649015329 * x) * x));
+}
+#endif
+
 
 inline double sas_gamma( double x) { return tgamma(x+1)/x; }

sasmodels/models/lib/sas_JN.c

@@ -48 +48 @@
 */
 
-static double
-sas_JN( int n, double x ) {
+double sas_JN( int n, double x );
+
+double sas_JN( int n, double x ) {
 
     const double MACHEP = 1.11022302462515654042E-16;

sasmodels/models/lib/sph_j1c.c

@@ -7 +7 @@
 * using double precision that are the source.
 */
+double sph_j1c(double q);
 
 // The choice of the number of terms in the series and the cutoff value for
@@ -43 +44 @@
 #endif
 
-double sph_j1c(double q);
 double sph_j1c(double q)
 {

sasmodels/models/lib/sphere_form.c

@@ -1 @@
-inline double
-sphere_volume(double radius)
+double sphere_volume(double radius);
+double sphere_form(double q, double radius, double sld, double solvent_sld);
+
+double sphere_volume(double radius)
 {
     return M_4PI_3*cube(radius);
 }
 
-inline double
-sphere_form(double q, double radius, double sld, double solvent_sld)
+double sphere_form(double q, double radius, double sld, double solvent_sld)
 {
     const double fq = sphere_volume(radius) * sph_j1c(q*radius);

sasmodels/models/lib/wrc_cyl.c

@@ -2 +2 @@
     Functions for WRC implementation of flexible cylinders
 */
+double Sk_WR(double q, double L, double b);
+
+
 static double
 AlphaSquare(double x)
@@ -363 +366 @@
 }
 
-double Sk_WR(double q, double L, double b);
 double Sk_WR(double q, double L, double b)
 {