In [1]:
import os
os.environ["SAS_OPENCL"] = "cuda" # use CUDA GPU backend for sasmodels
import escape as esc
import numpy as np
esc.require("0.9.8")
Loading material database from C:\dev\escape-core\python\src\escape\scattering\..\data\mdb\materials.db
SAXS. Form-factors. Pringle (SasView-aligned)¶
A saddle-shaped (hyperbolic paraboloid) disc whose curvature is controlled by two parameters $\alpha$ and $\beta$. The model is isotropic.
Reference: https://www.sasview.org/docs/user/models/pringle.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Radius (Å) | radius |
60 |
| Thickness (Å) | thickness |
10 |
| Curvature α | alpha_c |
0.001 |
| Curvature β | beta_c |
0.02 |
| Contrast Δρ (10⁻⁶ Å⁻²) | contrast |
−5.3 (= sld 1 − sld_solvent 6.3) |
Form-factor (SasView pringle.c, Rautu derivation)¶
$$I(q) = \text{scale}\cdot\Delta\rho^2 V\int_0^{\pi/2}4\sin\psi\;\mathrm{sinc}^2\!\left(\frac{qd\cos\psi}{2}\right)\left[(S_0^2+C_0^2)+2\sum_{n=1}^{3}(S_n^2+C_n^2)\right]d\psi + \text{background}$$
where the Bessel integrals are:
$$C_n = \frac{1}{R^2}\int_0^R r\cos(\alpha q r^2\cos\psi)\,J_n(\beta q r^2\cos\psi)\,J_{2n}(qr\sin\psi)\,dr$$
$$S_n = \frac{1}{R^2}\int_0^R r\sin(\alpha q r^2\cos\psi)\,J_n(\beta q r^2\cos\psi)\,J_{2n}(qr\sin\psi)\,dr$$
The 2D scattering pattern is the same as the 1D curve evaluated at $q = \sqrt{q_x^2+q_y^2}$ (isotropic model).
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
psi = esc.var("psi") # orientation angle for the disc normal vs q
r = esc.var("r") # radial integration variable over the disc
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True) # 1e8: Å⁻¹ → cm⁻¹ unit conversion
radius = esc.par("Radius", 60.0, units=esc.angstr)
thickness = esc.par("Thickness", 10.0, units=esc.angstr)
alpha_c = esc.par("Alpha", 0.001) # saddle curvature parameter
beta_c = esc.par("Beta", 0.02) # saddle curvature parameter
contrast = esc.par("Contrast", -5.3, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Geometry ───────────────────────────────────────────────────────────────
volume = np.pi * esc.pow(radius, 2) * thickness
# ── Orientation-dependent q components ─────────────────────────────────────
q_sin_psi = q * esc.sin(psi) # in-plane component
q_cos_psi = q * esc.cos(psi) # normal component
# ── Inner radial integrals ──────────────────────────────────────────────────
# qrrc = r^2 * q * cos(psi) (argument of curvature Bessel functions)
# qrs = r * q * sin(psi) (argument of in-plane Bessel function)
qrrc = esc.pow(r, 2) * q_cos_psi
qrs = r * q_sin_psi
def cn_sn(n):
"""Build ESCAPE functors for Cn and Sn at Bessel order n."""
jn_beta = esc.cyl_bessel_j(n, beta_c * qrrc)
j2n_qrs = esc.cyl_bessel_j(2*n, qrs)
kernel = r * jn_beta * j2n_qrs
Cn = (esc.integral(kernel * esc.cos(alpha_c * qrrc),
r, 0.0, radius, numpoints=61, maxiter=0)
/ esc.pow(radius, 2))
Sn = (esc.integral(kernel * esc.sin(alpha_c * qrrc),
r, 0.0, radius, numpoints=61, maxiter=2)
/ esc.pow(radius, 2))
return Cn, Sn
# Build the Bessel sum: S0^2+C0^2 + 2*(S1^2+C1^2 + S2^2+C2^2 + S3^2+C3^2)
C0, S0 = cn_sn(0)
C1, S1 = cn_sn(1)
C2, S2 = cn_sn(2)
C3, S3 = cn_sn(3)
bessel_sum = (esc.pow(S0, 2) + esc.pow(C0, 2)
+ 2.0 * (esc.pow(S1, 2) + esc.pow(C1, 2))
+ 2.0 * (esc.pow(S2, 2) + esc.pow(C2, 2))
+ 2.0 * (esc.pow(S3, 2) + esc.pow(C3, 2)))
# sinc^2(qd*cos(psi)/2) — axial thickness factor
sinc2 = esc.pow(esc.sinc(0.5 * thickness * q_cos_psi), 2)
# ── Outer integral over psi in [0, pi/2] ───────────────────────────────────
psi_integrand = 4.0 * esc.sin(psi) * sinc2 * bessel_sum
psi_integral = esc.integral(psi_integrand, psi, 0.0, np.pi / 2.0,
numpoints=61, maxiter=2, epsabs=1e-5)
i1d = scale * esc.pow(contrast, 2) * volume * psi_integral + background
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.5, 128)
i1d.show(coordinates=qs).config(
title="Pringle — isotropic 1D",
xlog=True, ylog=True,
xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")
Out[3]:
2D scattering (isotropic)¶
The pringle model is isotropic: the 2D scattering pattern is the same as the 1D curve evaluated at $q = \sqrt{q_x^2+q_y^2}$. No orientation parameters are defined.
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
psi2 = esc.var("psi2") # separate variable name to avoid clash with 1D psi
r2 = esc.var("r2") # separate variable name to avoid clash with 1D r
q_mag = esc.sqrt(esc.pow(qx, 2) + esc.pow(qy, 2))
q_sin_psi2 = q_mag * esc.sin(psi2)
q_cos_psi2 = q_mag * esc.cos(psi2)
qrrc2 = esc.pow(r2, 2) * q_cos_psi2
qrs2 = r2 * q_sin_psi2
def cn_sn_2d(n):
jn_beta = esc.cyl_bessel_j(n, beta_c * qrrc2)
j2n_qrs = esc.cyl_bessel_j(2*n, qrs2)
kernel = r2 * jn_beta * j2n_qrs
Cn = (esc.integral(kernel * esc.cos(alpha_c * qrrc2),
r2, 0.0, radius, numpoints=61, maxiter=0)
/ esc.pow(radius, 2))
Sn = (esc.integral(kernel * esc.sin(alpha_c * qrrc2),
r2, 0.0, radius, numpoints=61, maxiter=0)
/ esc.pow(radius, 2))
return Cn, Sn
C0_2d, S0_2d = cn_sn_2d(0)
C1_2d, S1_2d = cn_sn_2d(1)
C2_2d, S2_2d = cn_sn_2d(2)
C3_2d, S3_2d = cn_sn_2d(3)
bessel_sum_2d = (esc.pow(S0_2d, 2) + esc.pow(C0_2d, 2)
+ 2.0 * (esc.pow(S1_2d, 2) + esc.pow(C1_2d, 2))
+ 2.0 * (esc.pow(S2_2d, 2) + esc.pow(C2_2d, 2))
+ 2.0 * (esc.pow(S3_2d, 2) + esc.pow(C3_2d, 2)))
sinc2_2d = esc.pow(esc.sinc(0.5 * thickness * q_cos_psi2), 2)
psi_integrand_2d = 4.0 * esc.sin(psi2) * sinc2_2d * bessel_sum_2d
psi_integral_2d = esc.integral(psi_integrand_2d, psi2, 0.0, np.pi / 2.0,
numpoints=61, maxiter=0)
i2d = scale * esc.pow(contrast, 2) * volume * psi_integral_2d + background
i2d.device = "gpu"
xs = np.linspace(-0.5, 0.5, 200); ys = np.linspace(-0.5, 0.5, 200)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Pringle — isotropic 2D (qx, qy)",
xlabel=f"qx [{esc.angstr}^-1]", ylabel=f"qy [{esc.angstr}^-1]",
cblog=True, colorscale="jet")
Out[4]:
SasView reference model & comparison¶
| ESCAPE parameter | SasView parameter | Notes |
|---|---|---|
contrast * 1e-6 |
sld - sld_solvent |
contrast in Å⁻² |
radius |
radius |
disc radius (Å) |
thickness |
thickness |
disc thickness (Å) |
alpha_c |
alpha |
saddle curvature α |
beta_c |
beta |
saddle curvature β |
In [5]:
import time
import matplotlib.pyplot as plt
from sasmodels.core import load_model
from sasmodels.data import empty_data1D
from sasmodels.direct_model import DirectModel
qs = np.linspace(0.001, 0.5, 128)
kernel = load_model("pringle")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
sld=1.0, sld_solvent=6.3,
radius=60.0, thickness=10.0, alpha=0.001, beta=0.02)
f_sas(**sas_pars)
i1d.device = "gpu"; i1d(qs[:5])
def timeit(fn, n=5):
t0 = time.perf_counter()
for _ in range(n): result = fn()
return (time.perf_counter() - t0) / n * 1e3, result
t_sas, Iq_sas = timeit(lambda: f_sas(**sas_pars))
i1d.device = "gpu"
t_gpu, Iq_gpu = timeit(lambda: i1d(qs), n=3)
i1d.device = "cpu"
t_cpu, Iq_cpu = timeit(lambda: i1d(qs))
i1d.device = "gpu"
print(f"SASView GPU : {t_sas:.0f} ms")
print(f"ESCAPE GPU : {t_gpu:.0f} ms")
print(f"ESCAPE CPU : {t_cpu:.0f} ms ({len(qs)} q-pts)")
rel = np.max(np.abs((Iq_gpu - Iq_sas) / Iq_sas)) * 100
print(f"Max relative diff vs SasView: {rel:.2f}%")
esc.overlay(Iq_sas, Iq_gpu, Iq_cpu, coordinates=qs).config(
xlabel="Q (1/A)", ylabel="I(q) (1/cm)",
xlog=True, ylog=True,
fig_title=f"Pringle I(q) — {len(qs)} pts",
labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
line_styles=["solid", "dash", "dot"],
line_widths=[2, 3, 3]
)
C:\Users\User\AppData\Local\Temp\ipykernel_50852\140823826.py:16: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied C:\Users\User\AppData\Local\Temp\ipykernel_50852\140823826.py:26: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied C:\Users\User\AppData\Local\Temp\ipykernel_50852\140823826.py:28: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
SASView GPU : 31 ms ESCAPE GPU : 106 ms ESCAPE CPU : 282 ms (128 q-pts) Max relative diff vs SasView: 0.04%
Out[5]:
In [ ]: