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. Capped cylinder (SasView-aligned)¶
A right circular cylinder with spherical end caps whose centres lie inside the cylinder body (cap radius $R \geq$ cylinder radius $r$, protrusion $h = -\sqrt{R^2 - r^2} < 0$). This is the mirror geometry of the barbell.
Reference: https://www.sasview.org/docs/user/models/capped_cylinder.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Contrast Δρ (10⁻⁶ Å⁻²) | contrast |
3 (= sld 4 − sld_solvent 1) |
| Cylinder radius (Å) | radius |
20 |
| Cap radius (Å) | cap_radius |
20 |
| Length (Å) | length |
400 |
| Theta (deg), 2D only | theta |
60 |
| Phi (deg), 2D only | phi |
60 |
Form-factor (SasView capped_cylinder.c)¶
Identical amplitude to barbell but $h = -\sqrt{R^2 - r^2}$ (cap centre inside cylinder):
$$A_{\mathrm{cyl}} = \pi r^2 L\;\mathrm{sinc}\!\left(\tfrac{L}{2}q_\parallel\right)\frac{2J_1(r\,q_\perp)}{r\,q_\perp}$$
$$A_{\mathrm{cap}} = 2\pi R^3\int_{-h/R}^{1}\cos\!\left[q_\parallel\left(Rt+h+\tfrac{L}{2}\right)\right](1-t^2)\frac{2J_1\!\left(R\,q_\perp\sqrt{1-t^2}\right)}{R\,q_\perp\sqrt{1-t^2}}\,dt$$
$$I(q) = \frac{\mathrm{scale}\cdot\Delta\rho^2}{V}\int_0^{\pi/2}(A_{\mathrm{cyl}}+A_{\mathrm{cap}})^2\,\sin\alpha\,d\alpha + \mathrm{background}$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
alpha = esc.var("alpha") # angle between cylinder axis and q
t = esc.var("t") # integration variable over spherical cap surface
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
radius = esc.par("Radius", 20.0, units=esc.angstr) # cylinder radius r
cap_radius = esc.par("Cap Radius", 20.0, units=esc.angstr) # cap radius R
length = esc.par("Length", 400.0, units=esc.angstr)
contrast = esc.par("Contrast", 3.0, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Geometry ───────────────────────────────────────────────────────────────
# h = -sqrt(R^2 - r^2): cap centre lies INSIDE the cylinder (negative for capped cylinder)
h_cap = -esc.sqrt(esc.pow(cap_radius, 2) - esc.pow(radius, 2))
# Total volume: cylinder rod + 2 * (spherical cap volume)
v_cyl = np.pi * esc.pow(radius, 2) * length
v_tot = v_cyl + 2.0 * (
(2.0 * np.pi / 3.0) * esc.pow(cap_radius, 3)
+ np.pi * (esc.pow(cap_radius, 2) * h_cap - esc.pow(h_cap, 3) / 3.0)
)
# ── Oriented amplitude ─────────────────────────────────────────────────────
# q_axial = q * cos(alpha): component along cylinder axis
# q_radial = q * sin(alpha): component perpendicular to axis
q_axial = q * esc.cos(alpha)
q_radial = q * esc.sin(alpha)
# Cylinder amplitude
a_cyl = (np.pi * esc.pow(radius, 2) * length
* esc.sinc(0.5 * length * q_axial)
* 2.0 * esc.j1_over_x(radius * q_radial))
# Cap amplitude: same formula as barbell, but h_cap is negative
cap_arg = cap_radius * q_radial * esc.sqrt(1.0 - esc.pow(t, 2))
a_cap = 2.0 * np.pi * esc.pow(cap_radius, 3) * esc.integral(
esc.cos(q_axial * (cap_radius * t + h_cap + 0.5 * length))
* (1.0 - esc.pow(t, 2))
* 2.0 * esc.j1_over_x(cap_arg),
t, -h_cap / cap_radius, 1.0,
numpoints=61, maxiter=10, epsabs=1e-5)
a_tot = a_cyl + a_cap
# ── Powder average ─────────────────────────────────────────────────────────
i1d = (scale * esc.pow(contrast, 2) / v_tot
* esc.integral(esc.pow(a_tot, 2) * esc.sin(alpha),
alpha, 0.0, np.pi / 2.0,
numpoints=61, maxiter=10, epsabs=1e-5)
+ background)
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 1.0, 300)
i1d.show(coordinates=qs).config(
title="Capped cylinder — powder average (1D)",
xlog=True, ylog=True,
xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")
Out[3]:
2D oriented scattering (qx, qy)¶
For a fixed orientation $(\theta, \phi)$ the amplitude is evaluated directly at detector coordinates. The cylinder axis unit vector is $\hat{\mathbf{u}} = (\sin\theta\cos\phi,\;\sin\theta\sin\phi,\;\cos\theta)$.
$$I_{\mathrm{2D}}(q_x,q_y) = \frac{\mathrm{scale}\cdot\Delta\rho^2}{V}\,(A_{\mathrm{cyl}}+A_{\mathrm{cap}})^2 + \mathrm{background}$$
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
t2 = esc.var("t2") # separate integration variable for 2D (avoids name clash)
theta = esc.par("Theta", 60.0, userlim=[0.0, 180.0], units="deg")
phi = esc.par("Phi", 60.0, userlim=[0.0, 360.0], units="deg")
deg = np.pi / 180.0
sin_t = esc.sin(theta * deg)
ux = sin_t * esc.cos(phi * deg)
uy = sin_t * esc.sin(phi * deg)
q_sq = esc.pow(qx, 2) + esc.pow(qy, 2)
q_par_2d = qx * ux + qy * uy
q_perp_2d = esc.sqrt(q_sq - esc.pow(q_par_2d, 2))
a_cyl_2d = (np.pi * esc.pow(radius, 2) * length
* esc.sinc(0.5 * length * q_par_2d)
* 2.0 * esc.j1_over_x(radius * q_perp_2d))
cap_arg_2d = cap_radius * q_perp_2d * esc.sqrt(1.0 - esc.pow(t2, 2))
a_cap_2d = 2.0 * np.pi * esc.pow(cap_radius, 3) * esc.integral(
esc.cos(q_par_2d * (cap_radius * t2 + h_cap + 0.5 * length))
* (1.0 - esc.pow(t2, 2))
* 2.0 * esc.j1_over_x(cap_arg_2d),
t2, -h_cap / cap_radius, 1.0,
numpoints=61, maxiter=0)
a_tot_2d = a_cyl_2d + a_cap_2d
i2d = scale * esc.pow(contrast, 2) / v_tot * esc.pow(a_tot_2d, 2) + background
i2d.device = "gpu"
xs = np.linspace(-1.0, 1.0, 300); ys = np.linspace(-1.0, 1.0, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Capped cylinder — oriented 2D SAXS (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 |
cylinder radius (Å) |
cap_radius |
radius_cap |
cap sphere radius (Å) |
length |
length |
cylinder length (Å) |
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, 1.0, 300).copy()
kernel = load_model("capped_cylinder")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
sld=4.0, sld_solvent=1.0,
radius=20.0, radius_cap=20.0, length=400.0)
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"Capped cylinder 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_61456\3670665170.py:16: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
SASView GPU : 10 ms ESCAPE GPU : 277 ms ESCAPE CPU : 193 ms (300 q-pts) Max relative diff vs SasView: 5.37%
Out[5]:
In [ ]: