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. Hollow rectangular prism (SasView-aligned)¶
Hollow rectangular parallelepiped with uniform wall thickness $\Delta$. The amplitude is the difference of two solid parallelepipeds (outer minus inner). Matches hollow_rectangular_prism — SasView 6.1.3.
Reference: https://www.sasview.org/docs/user/models/hollow_rectangular_prism.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Contrast Δρ (10⁻⁶ Å⁻²) | contrast |
5.3 (= sld 6.3 − sld_solvent 1) |
| Shortest external side A (Å) | length_a |
35 |
| Ratio b/a | b2a_ratio |
1 |
| Ratio c/a | c2a_ratio |
1 |
| Wall thickness Δ (Å) | thickness |
1 |
| Theta (deg), 2D only | theta |
0 |
| Phi (deg), 2D only | phi |
0 |
| Psi (deg), 2D only | psi |
0 |
Form-factor (Nayuk, 2012)¶
Outer dimensions: $A, B=b/a\cdot A, C=c/a\cdot A$. Inner dimensions: $A'=A-2\Delta$, $B'=B-2\Delta$, $C'=C-2\Delta$.
Shell volume: $V = ABC - A'B'C'$
$$A_{P\Delta}(q) = ABC\,\mathrm{sinc}\!\left(\frac{C}{2}q\cos\theta\right)\mathrm{sinc}\!\left(\frac{A}{2}q\sin\theta\sin\phi\right)\mathrm{sinc}\!\left(\frac{B}{2}q\sin\theta\cos\phi\right)$$ $$\quad - A'B'C'\,\mathrm{sinc}\!\left(\frac{C'}{2}q\cos\theta\right)\mathrm{sinc}\!\left(\frac{A'}{2}q\sin\theta\sin\phi\right)\mathrm{sinc}\!\left(\frac{B'}{2}q\sin\theta\cos\phi\right)$$
$$P(q) = \frac{1}{V^2}\frac{2}{\pi}\int_0^{\pi/2}\int_0^{\pi/2}A_{P\Delta}^2\sin\theta\,d\theta\,d\phi$$
$$I(q) = \mathrm{scale}\cdot V\cdot(\Delta\rho)^2\cdot P(q) + \mathrm{background}$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
theta = esc.var("theta") # polar angle in powder average
phi = esc.var("phi") # azimuthal angle in powder average
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
length_a = esc.par("Length A", 35.0, units=esc.angstr)
b2a_ratio = esc.par("b/a ratio", 1.0, userlim=[0.01, 10.0])
c2a_ratio = esc.par("c/a ratio", 1.0, userlim=[0.01, 10.0])
thickness = esc.par("Thickness", 1.0, units=esc.angstr)
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 ───────────────────────────────────────────────────────────────
A = length_a
B = b2a_ratio * length_a
C = c2a_ratio * length_a
Ap = A - 2.0 * thickness
Bp = B - 2.0 * thickness
Cp = C - 2.0 * thickness
V_shell = A * B * C - Ap * Bp * Cp
# ── Amplitude difference ───────────────────────────────────────────────────
st = esc.sin(theta); ct = esc.cos(theta)
sp = esc.sin(phi); cp = esc.cos(phi)
qC_arg = 0.5 * q * C * ct
qA_arg = 0.5 * q * A * st * sp
qB_arg = 0.5 * q * B * st * cp
qCp_arg = 0.5 * q * Cp * ct
qAp_arg = 0.5 * q * Ap * st * sp
qBp_arg = 0.5 * q * Bp * st * cp
A_outer = A * B * C * esc.sinc(qC_arg) * esc.sinc(qA_arg) * esc.sinc(qB_arg)
A_inner = Ap * Bp * Cp * esc.sinc(qCp_arg) * esc.sinc(qAp_arg) * esc.sinc(qBp_arg)
A_pdelta = A_outer - A_inner
# ── Powder average ─────────────────────────────────────────────────────────
inner_phi = esc.integral(
esc.pow(A_pdelta, 2) * st,
phi, 0.0, np.pi / 2.0,
numpoints=61, maxiter=0,
)
P_q = (1.0 / esc.pow(V_shell, 2)) * (2.0 / np.pi) * esc.integral(
inner_phi, theta, 0.0, np.pi / 2.0,
numpoints=61, maxiter=0,
)
i1d = scale * V_shell * esc.pow(contrast, 2) * P_q + background
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.7, 300)
i1d.show(coordinates=qs).config(
title="Hollow rectangular prism — 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)¶
Same SasView rotation convention: $\theta$ in $xz$-plane, $\phi$ around $z$, $\psi$ around $C$-axis.
$$I_{\mathrm{2D}}(q_x,q_y) = \mathrm{scale}\cdot V\cdot(\Delta\rho)^2\cdot\frac{A_{P\Delta}^2}{V^2} + \mathrm{background}$$
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
theta_2d = esc.par("Theta", 0.0, userlim=[0.0, 180.0], units="deg")
phi_2d = esc.par("Phi", 0.0, userlim=[-180.0, 180.0], units="deg")
psi_2d = esc.par("Psi", 0.0, userlim=[-180.0, 180.0], units="deg")
deg = np.pi / 180.0
th = theta_2d * deg; ph = phi_2d * deg; ps = psi_2d * deg
ct2 = esc.cos(th); st2 = esc.sin(th)
cp2 = esc.cos(ph); sp2 = esc.sin(ph)
ck2 = esc.cos(ps); sk2 = esc.sin(ps)
Cx2, Cy2 = st2 * cp2, st2 * sp2
a0x2, a0y2 = cp2 * ct2, sp2 * ct2
b0x2, b0y2 = -sp2, cp2
Ax2 = ck2 * a0x2 + sk2 * b0x2; Ay2 = ck2 * a0y2 + sk2 * b0y2
Bx2 = -sk2 * a0x2 + ck2 * b0x2; By2 = -sk2 * a0y2 + ck2 * b0y2
qA_2d = qx * Ax2 + qy * Ay2
qB_2d = qx * Bx2 + qy * By2
qC_2d = qx * Cx2 + qy * Cy2
Aout_2d = A * B * C * esc.sinc(0.5*qC_2d*C) * esc.sinc(0.5*qA_2d*A) * esc.sinc(0.5*qB_2d*B)
Ain_2d = Ap * Bp * Cp * esc.sinc(0.5*qC_2d*Cp) * esc.sinc(0.5*qA_2d*Ap) * esc.sinc(0.5*qB_2d*Bp)
Apd_2d = Aout_2d - Ain_2d
i2d = scale * esc.pow(contrast, 2) * esc.pow(Apd_2d, 2) / V_shell + background
i2d.device = "gpu"
xs = np.linspace(-0.7, 0.7, 300); ys = np.linspace(-0.7, 0.7, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Hollow rectangular prism — 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 Å⁻² |
length_a |
length_a |
shortest external side A (Å) |
b2a_ratio |
b2a_ratio |
ratio B/A |
c2a_ratio |
c2a_ratio |
ratio C/A |
thickness |
thickness |
wall thickness Δ (Å) |
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.7, 300).copy()
kernel = load_model("hollow_rectangular_prism")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001, sld=6.3, sld_solvent=1.0,
length_a=35.0, b2a_ratio=1.0, c2a_ratio=1.0, thickness=1.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"Hollow rectangular prism 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_16616\963615288.py:15: UserWarning: Input array does not own its data (e.g. it is a view or slice); data will be copied
SASView GPU : 11 ms ESCAPE GPU : 1 ms ESCAPE CPU : 41 ms (300 q-pts) Max relative diff vs SasView: 0.02%
Out[5]:
In [ ]: