Changeset 6d96b66 in sasmodels for sasmodels/models/fractal_core_shell.c

Timestamp:

Oct 18, 2016 8:32:18 PM (8 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:

f0afad2

Parents:

3e8ea5d

Message:

fractal models: handle limits of q=0 and dim=1

File:

• sasmodels/models/fractal_core_shell.c (modified) (1 diff)

sasmodels/models/fractal_core_shell.c

@@ -27 +27 @@
 
 
-   double intensity = core_shell_kernel(q,
-                              radius,
-                              thickness,
-                              core_sld,
-                              shell_sld,
-                              solvent_sld);
+    const double pq = core_shell_kernel(q, radius, thickness,
+                                        core_sld, shell_sld, solvent_sld);
+
+
     //calculate S(q)
-    double frac_1 = fractal_dim-1.0;
-    double qr = q*radius;
+    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;
+    }
 
-    double t1 = fractal_dim*sas_gamma(frac_1)*sin(frac_1*atan(q*cor_length));
-    double t2 = (1.0 + 1.0/(q*cor_length)/(q*cor_length));
-    double t3 = pow(qr, fractal_dim) * pow(t2, (frac_1/2.0));
-    double sq = t1/t3;
-    sq += 1.0;
-
-    return sq*intensity*volfraction;
+    // Note: core_shell_kernel already performs the 1e-4 unit conversion
+    return volfraction * sq * pq;
 }