Changeset 217590b in sasmodels

Timestamp:

Oct 20, 2016 2:48:48 PM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

5c94f41

Parents:

2b9e63f

git-author:

Paul Kienzle <pkienzle@…> (10/20/16 14:47:25)

git-committer:

Paul Kienzle <pkienzle@…> (10/20/16 14:48:48)

Message:

fractal, fractal_core_shell: support fractal_dim as low as 0.

Location:

sasmodels/models

Files:

• fractal.c (modified) (3 diffs)
• fractal.py (modified) (2 diffs)
• fractal_core_shell.c (modified) (1 diff)
• fractal_core_shell.py (modified) (1 diff)
• lib/fractal_sq.c (added)
• lib/sas_gamma.c (modified) (1 diff)

sasmodels/models/fractal.c

@@ -1 @@
-#define INVALID(p) (p.fractal_dim <= 0.0)
+#define INVALID(p) (p.fractal_dim < 0.0)
 
 static double
@@ -10 +10 @@
    double sld_solvent)
 {
+    const double sq = fractal_sq(q, radius, fractal_dim, cor_length);
+
     //calculate P(q) for the spherical subunits
-    const double pq = M_4PI_3*cube(radius) * square((sld_block-sld_solvent)*sph_j1c(q*radius));
-    
-    //calculate S(q),  using Teixeira, Eq(15)
-    double sq;
-    if (q > 0. && fractal_dim > 1.) {
-        // q>0, D>0
-        const double D = fractal_dim;
-        const double Dm1 = fractal_dim - 1.0;
-        // Note: for large Dm1, sin(Dm1*atan(q*cor_length) can go negative
-        const double t1 = D*sas_gamma(Dm1)*sin(Dm1*atan(q*cor_length));
-        const double t2 = pow(q*radius, -D);
-        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1);
-        sq = 1.0 + t1 * t2 * t3;
-    } else if (q > 0.) {
-        // q>0, D=1
-        sq = 1.0 + atan(q*cor_length) / (q*radius);
-    } else if (fractal_dim > 1.) {
-        // q=0, D>1
-        const double D = fractal_dim;
-        sq = 1.0 + pow(cor_length/radius, D)*sas_gamma(D+1.0);
-    } else {
-        // q=0, D=1
-        sq = 1.0 + cor_length/radius;
-    }
+    const double V = M_4PI_3*cube(radius);
+    const double pq = V * square((sld_block-sld_solvent)*sph_j1c(q*radius));
 
     // scale to units cm-1 sr-1 (assuming data on absolute scale)
@@ -41 +21 @@
     //    combined: 1e-4 * I(q)
 
-    return 1.e-4 * volfraction * pq * sq;
+    return 1.e-4 * volfraction * sq * pq;
 }
 

sasmodels/models/fractal.py

@@ -84 +84 @@
               ["radius",    "Ang",  5.0, [0.0, inf], "",
                "radius of particles"],
-              ["fractal_dim",      "",  2.0, [1.0, 6.0], "",
+              ["fractal_dim",      "",  2.0, [0.0, 6.0], "",
                "fractal dimension"],
               ["cor_length", "Ang", 100.0, [0.0, inf], "",
@@ -95 +95 @@
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "fractal.c"]
+source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]
 
 demo = dict(volfraction=0.05,

sasmodels/models/fractal_core_shell.c

@@ -1 @@
-double form_volume(double radius, double thickness);
-
-double Iq(double q,
-          double radius,
-          double thickness,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double volfraction,
-          double fractal_dim,
-          double cor_length);
-
-double form_volume(double radius, double thickness)
+static double
+form_volume(double radius, double thickness)
 {
     return M_4PI_3 * cube(radius + thickness);
 }
 
-double Iq(double q,
-          double radius,
-          double thickness,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double volfraction,
-          double fractal_dim,
-          double cor_length) {
-
-
+static double
+Iq(double q,
+   double radius,
+   double thickness,
+   double core_sld,
+   double shell_sld,
+   double solvent_sld,
+   double volfraction,
+   double fractal_dim,
+   double cor_length)
+{
+    const double sq = fractal_sq(q, radius, fractal_dim, cor_length);
     const double pq = core_shell_kernel(q, radius, thickness,
                                         core_sld, shell_sld, solvent_sld);
-
-
-    //calculate S(q)
-    double sq;
-    if (q > 0. && fractal_dim > 1.) {
-        // q>0, D>0
-        const double D = fractal_dim;
-        const double Dm1 = fractal_dim - 1.0;
-        const double t1 = D*sas_gamma(Dm1)*sin((Dm1)*atan(q*cor_length));
-        const double t2 = pow(q*radius, -D);
-        const double t3 = pow(1.0 + 1.0/square(q*cor_length), -0.5*Dm1);
-        sq = 1.0 + t1 * t2 * t3;
-    } else if (q > 0.) {
-        // q>0, D=1
-        sq = 1.0 + atan(q*cor_length) / (q*radius);
-    } else if (fractal_dim > 1.) {
-        // q=0, D>1
-        const double D = fractal_dim;
-        sq = 1.0 + pow(cor_length/radius, D)*sas_gamma(D+1.0);
-    } else {
-        // q=0, D=1
-        sq = 1.0 + cor_length/radius;
-    }
 
     // Note: core_shell_kernel already performs the 1e-4 unit conversion

sasmodels/models/fractal_core_shell.py

@@ -69 +69 @@
     ["sld_solvent", "1e-6/Ang^2", 3.0,  [-inf, inf], "sld",    "Solvent scattering length density"],
     ["volfraction", "",           1.0,  [0.0, inf],  "",       "Volume fraction of building block spheres"],
-    ["fractal_dim",    "",        2.0,  [1.0, 6.0],  "",       "Fractal dimension"],
+    ["fractal_dim",    "",        2.0,  [0.0, 6.0],  "",       "Fractal dimension"],
     ["cor_length",  "Ang",      100.0,  [0.0, inf],  "",       "Correlation length of fractal-like aggregates"],
 ]
 # pylint: enable=bad-whitespace, line-too-long
 
-source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c", "fractal_core_shell.c"]
+source = ["lib/sph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c",
+          "lib/fractal_sq.c", "fractal_core_shell.c"]
 
 demo = dict(scale=0.05,

sasmodels/models/lib/sas_gamma.c

@@ -137 +137 @@
 
 
-inline double sas_gamma( double x) { return tgamma(x+1)/x; }
+inline double sas_gamma(double x) {
+#if 1
+    // Fast reliable gamma in [-n,171], where n is a small integer
+    double norm = 1.;
+    while (x < 1.) {
+        norm *= x;
+        x += 1.0;
+    }
+    return tgamma(x)/norm;
+#else
+    // Fast reliable gamma in [0,170]
+    return tgamma(x+1)/x;
+#endif
+}