Skip to content

verify accuracy of fractal S(q) for small q #722

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
pkienzle wants to merge 1 commit into
base: master
Choose a base branch
from
Open

verify accuracy of fractal S(q) for small q #722

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 43 additions & 0 deletions explore/precision.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
"""
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand Down
2 changes: 1 addition & 1 deletion sasmodels/models/fractal.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------------------------
Expand Down
24 changes: 9 additions & 15 deletions sasmodels/models/lib/fractal_sq.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.));

@pkienzle pkienzle Jun 29, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Regarding the transformation, I multiplied $t_2$ by $(qξ)^{D-1}$ giving intermediate $t_2'' = (ξ/r)^D/(qξ)$, and $t_3$ by $1/(qξ)^{D-1}$ giving the new $t_3'$.

I then rearranged $t_1$ and $t_2''$ forming the new $t_2'$ using $\sin((D-1)\tan^{-1}(qξ))/(D-1)$ from $t_1$ and $1/(qξ)$ from $t_2''$, and forming the new $t_1'$ using $Γ(D+1)$ from $t_1$ and $(ξ/r)^D$ from $t_2''$.

Limit as $q \rightarrow 0$ for $t_2'$ and $t_3'$ are both 1, so only the constant $t_1'$ remains.

Limit as $D \rightarrow 0$ for $t_1'$ is 1. For $t_2'$ it is $\sin (\tan^{-1}(qξ))/(qξ) = 1/\sqrt{1+(qξ)^2}$ and for $t_3'$ it is $\sqrt{1 + (qξ)^2}$ which cancel.

Limit as $D \rightarrow 1$ for $t_1'$ is $ξ/r$, for $t_2'$ is $\tan^{-1}(qξ) / (qξ)$ and for $t_3'$ is 1.

Further check: Limit $q \rightarrow 0$ for $D \rightarrow 1$ equals $ξ/r$ matches limit $D \rightarrow 1$ for $q \rightarrow 0$.

term = t1 * t2 * t3;
}
return 1.0 + term;
return 1. + term;
}
21 changes: 20 additions & 1 deletion sasmodels/special/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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**
Expand Down Expand Up @@ -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:
Expand Down
Loading