verify accuracy of fractal S(q) for small q

explore/precision.py

@@ -183,6 +183,7 @@ def compare(self, x, precision, target, linear=False, diff="relative"):
             pylab.semilogx(x, target, '-', label="true value")
         if linear:
             pylab.xscale('linear')
+        #pylab.yscale('linear')
 
 def plotdiff(x, target, actual, label, diff):
     """
@@ -192,10 +193,20 @@ def plotdiff(x, target, actual, label, diff):
     """
     if diff == "relative":
         err = np.array([(abs((t-a)/t) if t != 0 else a) for t, a in zip(target, actual)], 'd')
+        if (err == 0).any():
+            if (err > 0).any():
+                err[err == 0] = err[err > 0].min()/2
+            else:
+                err[err == 0] = 1e-20
         #err = np.clip(err, 0, 1)
         pylab.loglog(x, err, '-', label=label, alpha=0.7)
     elif diff == "absolute":
         err = np.array([abs(t-a) for t, a in zip(target, actual)], 'd')
+        if (err == 0).any():
+            if (err > 0).any():
+                err[err == 0] = err[err > 0].min()/2
+            else:
+                err[err == 0] = 1e-20
         pylab.loglog(x, err, '-', label=label, alpha=0.7)
     else:
         limits = np.min(target), np.max(target)
@@ -608,6 +619,37 @@ def np_star_polymer(x, arms=3):
     xaxis="$(Q R_g)^2$ (unitless)",
 )
 
+from sasmodels.special import fractal_sq
+
+
+def mp_fractal_sq(q, radius, fractal_dim, cor_length):
+    # Eq 16 from Texiera (1988), DOI:10.1107/S0021889888000263
+    D = fractal_dim
+    if D == 0:
+        term = 1
+    elif D == 1:
+        term = 1 if q == 0 else mp.atan(q*cor_length)/(q*radius)
+    elif q == 0:
+        term = mp.power(cor_length/radius, D)*mp.gamma(D+1)
+    else:
+        t1 = mp.power(q*radius, -D)
+        t2 = D*mp.gamma(D-1) / mp.power(1 + 1/(q*cor_length)**2, (D-1)/2)
+        t3 = mp.sin((D-1)*mp.atan(q*cor_length))
+        term = t1*t2*t3
+
+    return 1 + term
+
+FRACTAL_RADIUS=30
+FRACTAL_DIM=1.6
+FRACTAL_CORR=421
+add_function(
+    name="fractal",
+    mp_function=lambda x: mp_fractal_sq(x, FRACTAL_RADIUS, FRACTAL_DIM, FRACTAL_CORR),
+    np_function=lambda x: fractal_sq(x, FRACTAL_RADIUS, FRACTAL_DIM, FRACTAL_CORR),
+    ocl_function=make_ocl(f"return fractal_sq(q, {FRACTAL_RADIUS}, {FRACTAL_DIM}, {FRACTAL_CORR});",
+                          "fractal", source=["lib/sas_gamma.c", "lib/fractal_sq.c"]),
+)
+
 
 RADIUS=3000
 LENGTH=30
@@ -686,6 +728,7 @@ def np_cyl(x):
     ocl_function=make_ocl(lanczos_gamma, "lgamma"),
 )
 
+
 replacement_expm1 = """\
       double x = (double)q;  // go back to float for single precision kernels
       // Adapted from the cephes math library.

sasmodels/models/fractal.py

@@ -43,7 +43,7 @@
 References
 ----------
 
-#.  J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785
+#. J Teixeira, *J. Appl. Cryst.*, 21 (1988) 781-785 DOI:10.1107/S0021889888000263
 
 Authorship and Verification
 ----------------------------

sasmodels/models/lib/fractal_sq.c

@@ -1,28 +1,22 @@
 static double
 fractal_sq(double q, double radius, double fractal_dim, double cor_length)
 {
-    //calculate S(q),  using Teixeira, Eq(15)
-    // mathematica query to check limiting conditions:
-    //    lim x->0 of [ x gamma(x-1) sin(arctan(q c (x-1))) (q r)^(-x) (1 + 1/(q c)^2)^((1-x)/2) ]
-    // Note: gamma(x) may be unreliable for x<0, so the gamma(D-1) is risky.
-    // We instead transform D*gamma(D-1) into gamma(D+1)/(D-1).
+    // Calculate S(q) using Eq(16) from Teixeira (1988), DOI:10.1107/S0021889888000263
+    // with terms rearranged so that limiting cases are easier to calculate.
+    // Values checked against 500 bit floating point using the original equations.
     double term;
     if (q == 0.) {
         const double D = fractal_dim;
-        term = pow(cor_length/radius, D)*sas_gamma(D+1.);
+        term = pow(cor_length/radius, fractal_dim)*sas_gamma(fractal_dim+1.);
     } else if (fractal_dim == 0.) {
-        term = 1.0;
+        term = 1.;
     } else if (fractal_dim == 1.) {
         term = atan(q*cor_length)/(q*radius);
     } else {
-        // 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 = sas_gamma(D+1.)/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);
+        const double t1 = pow(cor_length/radius, fractal_dim) * sas_gamma(fractal_dim+1.);
+        const double t2 = sin((fractal_dim-1.)*atan(q*cor_length))/((fractal_dim-1)*q*cor_length);
+        const double t3 = pow(1. + square(q*cor_length), -0.5*(fractal_dim-1.));
         term = t1 * t2 * t3;
     }
-    return 1.0 + term;
+    return 1. + term;
 }

sasmodels/special/__init__.py

@@ -14,7 +14,8 @@
 
 # Functions to add to our standard set
 # C99 standard math library functions
-from numpy import cos, inf, nan, sin
+from numpy import acos, asin, atan, cos, inf, nan, pow, sin, tan
+from scipy.special import gamma as sas_gamma
 from scipy.special import j1 as sas_J1
 
 # erf, erfc, tgamma, lgamma  **do not use**
@@ -98,6 +99,24 @@ def sas_2J1x_x(x):
         retvalue[x == 0] = 1.
     return retvalue
 
+def fractal_sq(q, radius, fractal_dim, cor_length):
+    # Calculate S(q) using Eq(16) from Teixeira (1988), DOI:10.1107/S0021889888000263
+    # with terms rearranged so that limiting cases are easier to calculate.
+    # Values checked against 500 bit floating point using the original equations.
+    with np.errstate(all='ignore'):
+        if fractal_dim == 0:
+            term = np.ones_like(q)
+        elif fractal_dim == 1:
+            term = atan(q*cor_length)/(q*radius)
+            term[q == 0] = 1
+        else:
+            t1 = pow(cor_length/radius, fractal_dim) * sas_gamma(fractal_dim+1)
+            t2 = sin((fractal_dim-1)*atan(q*cor_length))/((fractal_dim-1)*q*cor_length)
+            t2[q == 0] = 1
+            t3 = pow(1 + square(q*cor_length), -(fractal_dim-1)/2)
+            term = t1*t2*t3
+
+    return 1 + term
 
 # Gaussians
 class Gauss: