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. Triaxial ellipsoid (SasView-aligned)¶
Ellipsoid of uniform scattering length density with three independent axes $R_a$, $R_b$, $R_c$. Matches triaxial_ellipsoid — SasView 6.1.3.
Reference: https://www.sasview.org/docs/user/models/triaxial_ellipsoid.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Contrast Δρ (10⁻⁶ Å⁻²) | contrast |
3 (= sld 4 − sld_solvent 1) |
| Minor equatorial radius Ra (Å) | radius_equat_minor |
20 |
| Major equatorial radius Rb (Å) | radius_equat_major |
400 |
| Polar radius Rc (Å) | radius_polar |
10 |
| Theta (deg), 2D only | theta |
60 |
| Phi (deg), 2D only | phi |
60 |
| Psi (deg), 2D only | psi |
60 |
Form-factor (SasView triaxial_ellipsoid.c)¶
Using the substitution $u = \sin\gamma$ (so $du = \cos\gamma\,d\gamma$), the double integral over orientations reduces to:
$$\langle\Phi^2\rangle = 8\int_0^{\pi/2}\int_0^1 \Phi^2(qr)\,du\,d\phi$$
where
$$r^2 = R_b^2\bigl(p_a\sin^2\phi\,(1-u^2) + 1 + p_c\,u^2\bigr), \quad p_a = \frac{R_a^2}{R_b^2}-1, \quad p_c = \frac{R_c^2}{R_b^2}-1$$
$$\Phi(qr) = \frac{3(\sin(qr)-qr\cos(qr))}{(qr)^3} = 3\,j_1(qr)/(qr), \quad V = \frac{4\pi}{3}R_a R_b R_c$$
$$I(q) = \mathrm{scale}\cdot(\Delta\rho)^2\cdot\frac{V}{4\pi}\cdot\langle\Phi^2\rangle + \mathrm{background}$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
phi = esc.var("phi") # equatorial angle in powder average
u = esc.var("u") # u = sin(gamma), polar substitution
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
radius_equat_minor = esc.par("Minor Equatorial Ra", 20.0, units=esc.angstr)
radius_equat_major = esc.par("Major Equatorial Rb", 400.0, units=esc.angstr)
radius_polar = esc.par("Polar Radius Rc", 10.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 ───────────────────────────────────────────────────────────────
Ra = radius_equat_minor
Rb = radius_equat_major
Rc = radius_polar
volume = 4.0 / 3.0 * np.pi * Ra * Rb * Rc
# Dimensionless ratios (SasView notation)
pa = esc.pow(Ra, 2) / esc.pow(Rb, 2) - 1.0
pc = esc.pow(Rc, 2) / esc.pow(Rb, 2) - 1.0
# ── Effective radius ───────────────────────────────────────────────────────
r_sq = esc.pow(Rb, 2) * (pa * esc.pow(esc.sin(phi), 2) * (1.0 - esc.pow(u, 2)) + 1.0 + pc * esc.pow(u, 2))
r = esc.sqrt(r_sq)
qr = q * r
# ── Phi(qr) = 3*j1(qr)/(qr) = 3*(sin(qr)-qr*cos(qr))/(qr)^3 ─────────────
Phi = 3.0 * esc.conditional(
esc.abs(qr) < 1e-10,
1.0 / 3.0,
(esc.sin(qr) - qr * esc.cos(qr)) / esc.pow(qr, 3),
)
# ── Double powder average: <Phi^2> = 8 * integral (factor 8 from symmetry) ─
inner_u = esc.integral(esc.pow(Phi, 2), u, 0.0, 1.0, numpoints=61, maxiter=5, epsabs=1e-5)
outer_phi = esc.integral(inner_u, phi, 0.0, np.pi / 2.0, numpoints=61, maxiter=5, epsabs=1e-5)
# I = scale * (Δρ)^2 * V / (4π) * <Phi^2> + background
i1d = (scale * esc.pow(contrast, 2) * volume / (4.0 * np.pi)
* 8.0 * outer_phi
+ background)
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.7, 300)
i1d.show(coordinates=qs).config(
title="Triaxial ellipsoid — 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_{ori}, \psi)$ the amplitude is evaluated directly at detector coordinates $(q_x, q_y)$. The rotation sequence follows SasView: first $\theta$ in the $xz$-plane, then $\phi_{ori}$ around $z$, then $\psi$ around the $C$-axis.
The effective radius for the oriented particle is:
$$r^2 = \left(\frac{q_A}{q}\right)^2 R_a^2 + \left(\frac{q_B}{q}\right)^2 R_b^2 + \left(\frac{q_C}{q}\right)^2 R_c^2$$
$$I_{\mathrm{2D}}(q_x,q_y) = \mathrm{scale}\cdot(\Delta\rho)^2\cdot\frac{V}{4\pi}\cdot\Phi^2(qr) + \mathrm{background}$$
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
theta_ori = esc.par("Theta", 60.0, userlim=[0.0, 180.0], units="deg")
phi_ori = esc.par("Phi", 60.0, userlim=[0.0, 360.0], units="deg")
psi_ori = esc.par("Psi", 60.0, userlim=[-180.0, 180.0], units="deg")
deg = np.pi / 180.0
th = theta_ori * deg
ph = phi_ori * deg
ps = psi_ori * deg
ct = esc.cos(th); st = esc.sin(th)
cp = esc.cos(ph); sp = esc.sin(ph)
ck = esc.cos(ps); sk = esc.sin(ps)
# C-axis (polar)
Cx, Cy = st * cp, st * sp
# Transverse axes before psi roll
a0x, a0y = cp * ct, sp * ct
b0x, b0y = -sp, cp
# Roll around C by psi
Ax = ck * a0x + sk * b0x; Ay = ck * a0y + sk * b0y
Bx = -sk * a0x + ck * b0x; By = -sk * a0y + ck * b0y
# Projections of q onto particle axes (qz=0 for detector plane)
qA = qx * Ax + qy * Ay
qB = qx * Bx + qy * By
qC = qx * Cx + qy * Cy
# Effective radius from direction cosines
q_mag_sq = esc.pow(qx, 2) + esc.pow(qy, 2)
r_2d_sq = esc.conditional(
q_mag_sq < 1e-20,
esc.pow(Rb, 2),
(esc.pow(qA, 2) * esc.pow(Ra, 2) + esc.pow(qB, 2) * esc.pow(Rb, 2) + esc.pow(qC, 2) * esc.pow(Rc, 2)) / q_mag_sq,
)
qr_2d = esc.sqrt(q_mag_sq) * esc.sqrt(r_2d_sq)
Phi_2d = 3.0 * esc.conditional(
esc.abs(qr_2d) < 1e-10,
1.0 / 3.0,
(esc.sin(qr_2d) - qr_2d * esc.cos(qr_2d)) / esc.pow(qr_2d, 3),
)
i2d = (scale * esc.pow(contrast, 2) * volume / (4.0 * np.pi)
* esc.pow(Phi_2d, 2)
+ 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="Triaxial ellipsoid — 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_equat_minor |
radius_equat_minor |
Ra — minor equatorial radius (Å) |
radius_equat_major |
radius_equat_major |
Rb — major equatorial radius (Å) |
radius_polar |
radius_polar |
Rc — polar radius (Å) |
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("triaxial_ellipsoid")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001, sld=4.0, sld_solvent=1.0,
radius_equat_minor=20.0, radius_equat_major=400.0, radius_polar=10.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"Triaxial ellipsoid 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_52868\2604801801.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 : 344 ms ESCAPE CPU : 1807 ms (300 q-pts) Max relative diff vs SasView: 0.81%
Out[5]:
In [ ]: