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 — thin walls (SasView-aligned)¶
Hollow rectangular parallelepiped with infinitely thin walls. The scattering amplitude is the sum of face and edge contributions $A_L + A_T$. This model computes 1D only (no 2D per SasView). Matches hollow_rectangular_prism_thin_walls — SasView 6.1.3.
Reference: https://www.sasview.org/docs/user/models/hollow_rectangular_prism_thin_walls.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 side A (Å) | length_a |
35 |
| Ratio b/a | b2a_ratio |
1 |
| Ratio c/a | c2a_ratio |
1 |
Form-factor (Nayuk & Huber, 2012)¶
Surface area: $V = 2AB + 2AC + 2BC$ (used as normalization).
$$A_L = 8\frac{\sin(\tfrac{1}{2}qA\sin\phi\sin\theta)\,\sin(\tfrac{1}{2}qB\cos\phi\sin\theta)\,\cos(\tfrac{1}{2}qC\cos\theta)}{q^2\sin^2\theta\,\sin\phi\cos\phi}$$
$$A_F = 4\frac{\cos(\tfrac{1}{2}qA\sin\phi\sin\theta)\,\sin(\tfrac{1}{2}qB\cos\phi\sin\theta)}{q\cos\phi\sin\theta} + 4\frac{\sin(\tfrac{1}{2}qA\sin\phi\sin\theta)\,\cos(\tfrac{1}{2}qB\cos\phi\sin\theta)}{q\sin\phi\sin\theta}$$
$$A_T = A_F\cdot\frac{2\sin(\tfrac{1}{2}qC\cos\theta)}{q\cos\theta}$$
$$P(q) = \frac{1}{V^2}\frac{2}{\pi}\int_0^{\pi/2}\int_0^{\pi/2}(A_L+A_T)^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
phi = esc.var("phi") # azimuthal angle
# ── 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])
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
# Surface area (thin-wall normalization)
V_surf = 2.0 * (A * B + A * C + B * C)
# ── Trig factors ──────────────────────────────────────────────────────────
st = esc.sin(theta); ct = esc.cos(theta)
sp = esc.sin(phi); cp = esc.cos(phi)
half_qA_st_sp = 0.5 * q * A * st * sp
half_qB_st_cp = 0.5 * q * B * st * cp
half_qC_ct = 0.5 * q * C * ct
# ── A_L: face contribution ─────────────────────────────────────────────────
# Numerically stable: avoid 0/0 at theta=0 or phi=0,pi/2
denom_L = q * q * esc.pow(st, 2) * sp * cp
A_L = esc.conditional(
esc.abs(denom_L) < 1e-20,
0.0,
8.0 * esc.sin(half_qA_st_sp) * esc.sin(half_qB_st_cp) * esc.cos(half_qC_ct) / denom_L,
)
# ── A_F: edge amplitude helper ─────────────────────────────────────────────
denom_cp = q * cp * st
denom_sp = q * sp * st
term1 = esc.conditional(
esc.abs(denom_cp) < 1e-20,
0.0,
4.0 * esc.cos(half_qA_st_sp) * esc.sin(half_qB_st_cp) / denom_cp,
)
term2 = esc.conditional(
esc.abs(denom_sp) < 1e-20,
0.0,
4.0 * esc.sin(half_qA_st_sp) * esc.cos(half_qB_st_cp) / denom_sp,
)
A_F = term1 + term2
# ── A_T: top/bottom contribution ──────────────────────────────────────────
denom_ct = q * ct
A_T = esc.conditional(
esc.abs(denom_ct) < 1e-20,
0.0,
A_F * 2.0 * esc.sin(half_qC_ct) / denom_ct,
)
A_total = A_L + A_T
# ── Powder average ─────────────────────────────────────────────────────────
inner_phi = esc.integral(
esc.pow(A_total, 2) * st,
phi, 0.0, np.pi / 2.0,
numpoints=61, maxiter=0,
)
P_q = (1.0 / esc.pow(V_surf, 2)) * (2.0 / np.pi) * esc.integral(
inner_phi, theta, 0.0, np.pi / 2.0,
numpoints=61, maxiter=0,
)
i1d = scale * V_surf * 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 (thin walls) — powder average (1D)",
xlog=True, ylog=True,
xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")
Out[3]:
2D scattering¶
This model is isotropic in 2D (same as 1D with $q = \sqrt{q_x^2 + q_y^2}$). SasView does not provide a 2D oriented version.
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
q2d = esc.sqrt(esc.pow(qx, 2) + esc.pow(qy, 2))
# Re-use i1d expression substituting q -> q2d
# Build 2D isotropic version directly
st2 = esc.sin(theta); ct2 = esc.cos(theta)
sp2 = esc.sin(phi); cp2 = esc.cos(phi)
hqA2 = 0.5 * q2d * A * st2 * sp2
hqB2 = 0.5 * q2d * B * st2 * cp2
hqC2 = 0.5 * q2d * C * ct2
dL2 = q2d * q2d * esc.pow(st2, 2) * sp2 * cp2
AL2 = esc.conditional(esc.abs(dL2) < 1e-20, 0.0,
8.0 * esc.sin(hqA2) * esc.sin(hqB2) * esc.cos(hqC2) / dL2)
dc2 = q2d * cp2 * st2; ds2 = q2d * sp2 * st2
t12 = esc.conditional(esc.abs(dc2) < 1e-20, 0.0, 4.0 * esc.cos(hqA2) * esc.sin(hqB2) / dc2)
t22 = esc.conditional(esc.abs(ds2) < 1e-20, 0.0, 4.0 * esc.sin(hqA2) * esc.cos(hqB2) / ds2)
AF2 = t12 + t22
dct2 = q2d * ct2
AT2 = esc.conditional(esc.abs(dct2) < 1e-20, 0.0, AF2 * 2.0 * esc.sin(hqC2) / dct2)
inner2 = esc.integral(esc.pow(AL2 + AT2, 2) * st2, phi, 0.0, np.pi/2.0, numpoints=61, maxiter=0)
P2d = (1.0 / esc.pow(V_surf, 2)) * (2.0 / np.pi) * esc.integral(inner2, theta, 0.0, np.pi/2.0, numpoints=61, maxiter=0)
i2d = scale * V_surf * esc.pow(contrast, 2) * P2d + background
i2d.device = "gpu"
xs = np.linspace(-0.7, 0.7, 200); ys = np.linspace(-0.7, 0.7, 200)
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 (thin walls) — isotropic 2D",
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 side A (Å) |
b2a_ratio |
b2a_ratio |
ratio B/A |
c2a_ratio |
c2a_ratio |
ratio C/A |
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_thin_walls")
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)
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 (thin walls) 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_20208\2857610645.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 : 5 ms ESCAPE CPU : 228 ms (300 q-pts) Max relative diff vs SasView: 0.02%
Out[5]:
In [ ]: