Changeset b716cc6 in sasmodels for sasmodels/models/surface_fractal.c

Timestamp:

Oct 14, 2016 5:20:17 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:

3a48772

Parents:

6831fa0

Message:

surface_fractal: document usable q range; suppress negative values outside that range

File:

• sasmodels/models/surface_fractal.c (modified) (2 diffs)

sasmodels/models/surface_fractal.c

@@ +1 @@
+// Don't need invalid test since fractal_dim_surf is not polydisperse
+// #define INVALID(v) (v.fractal_dim_surf <= 1.0 || v.fractal_dim_surf >= 3.0)
+
 double form_volume(double radius);
 
@@ -11 +14 @@
     double cutoff_length)
 {
-    double pq, sq, mmo, result;
+    // calculate P(q)
+    const double pq = square(sph_j1c(q*radius));
 
-    //Replaced the original formula with Taylor expansion near zero.
-    //pq = pow((3.0*(sin(q*radius) - q*radius*cos(q*radius))/pow((q*radius),3)),2);
+    // calculate S(q)
+    // Note: lim q->0 S(q) = -gamma(mmo) cutoff_length^mmo (mmo cutoff_length)
+    // however, the surface fractal formula is invalid outside the range
+    const double mmo = 5.0 - fractal_dim_surf;
+    const double sq = sas_gamma(mmo) * pow(cutoff_length, mmo)
+           * pow(1.0 + square(q*cutoff_length), -0.5*mmo)
+           * sin(-mmo * atan(q*cutoff_length)) / q;
 
-    pq = sph_j1c(q*radius);
-    pq = pq*pq;
+    // Empirically determined that the results are valid within this range.
+    // Above 1/r, the form starts to oscillate;  below
+    //const double result = (q > 5./(3-fractal_dim_surf)/cutoff_length) && q < 1./radius
+    //                      ? pq * sq : 0.);
 
-    //calculate S(q)
-    mmo = 5.0 - fractal_dim_surf;
-    sq  = sas_gamma(mmo)*sin(-(mmo)*atan(q*cutoff_length));
-    sq *= pow(cutoff_length, mmo);
-    sq /= pow((1.0 + (q*cutoff_length)*(q*cutoff_length)),(mmo/2.0));
-    sq /= q;
+    double result = pq * sq;
 
-    //combine and return
-    result = pq * sq;
-
-    return result;
+    // exclude negative results
+    return result > 0. ? result : 0.;
 }
-double form_volume(double radius){
-
-    return 1.333333333333333*M_PI*radius*radius*radius;
+double form_volume(double radius)
+{
+    return M_4PI_3*cube(radius);
 }