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. Core-shell ellipsoid (SasView-aligned)¶
Core-shell spheroid with independent core and shell SLDs. Matches core_shell_ellipsoid — SasView 6.1.3.
Reference: https://www.sasview.org/docs/user/models/core_shell_ellipsoid.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Core equatorial radius (Å) | radius_equat_core |
200 |
| Core axial ratio (Rp_core/Re_core) | x_core |
0.1 |
| Shell thickness at equator (Å) | thick_shell |
30 |
| Polar/equatorial shell-thickness ratio | x_polar_shell |
1.0 |
| Core SLD (10⁻⁶ Å⁻²) | sld_core |
2 |
| Shell SLD (10⁻⁶ Å⁻²) | sld_shell |
1 |
| Solvent SLD (10⁻⁶ Å⁻²) | sld_solvent |
6.3 |
| Theta (deg), 2D only | theta |
0 |
| Phi (deg), 2D only | phi |
0 |
Form-factor¶
For each interface (core–shell and shell–solvent), using $u = \cos\alpha$:
$$r_{core}(u) = \sqrt{R_{e,core}^2(1-u^2) + R_{p,core}^2 u^2}$$ $$r_{shell}(u) = \sqrt{R_{e,shell}^2(1-u^2) + R_{p,shell}^2 u^2}$$
$$f(q,r) = 3\Delta\rho\,V\,\frac{\sin(qr)-qr\cos(qr)}{(qr)^3}$$
Total amplitude: $F = f_{core-shell} + f_{shell-solvent}$
$$I(q) = \frac{\mathrm{scale}}{V_{shell}}\int_0^1 F^2\,du + \mathrm{background}$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
u = esc.var("u") # u = cos(alpha)
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
re_core = esc.par("Core Equatorial Radius", 200.0, units=esc.angstr)
x_core = esc.par("Core Axial Ratio", 0.1, userlim=[0.01, 10.0])
thick_shell = esc.par("Shell Thickness", 30.0, units=esc.angstr)
x_polar_shell = esc.par("Shell Polar Ratio", 1.0, userlim=[0.01, 10.0])
sld_core = esc.par("Core SLD", 2.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_shell = esc.par("Shell SLD", 1.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent = esc.par("Solvent SLD", 6.3, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Derived geometry ───────────────────────────────────────────────────────
rp_core = re_core * x_core
re_shell = re_core + thick_shell
rp_shell = rp_core + thick_shell * x_polar_shell
Vc = 4.0 / 3.0 * np.pi * rp_core * esc.pow(re_core, 2)
Vs = 4.0 / 3.0 * np.pi * rp_shell * esc.pow(re_shell, 2)
# ── Spherical kernels ──────────────────────────────────────────────────────
r_c = esc.sqrt(esc.pow(re_core, 2) * (1.0 - esc.pow(u, 2)) + esc.pow(rp_core, 2) * esc.pow(u, 2))
r_s = esc.sqrt(esc.pow(re_shell, 2) * (1.0 - esc.pow(u, 2)) + esc.pow(rp_shell, 2) * esc.pow(u, 2))
qrc = q * r_c
qrs = q * r_s
kc = esc.conditional(esc.abs(qrc) < 1e-10, 1.0/3.0, (esc.sin(qrc) - qrc*esc.cos(qrc)) / esc.pow(qrc, 3))
ks = esc.conditional(esc.abs(qrs) < 1e-10, 1.0/3.0, (esc.sin(qrs) - qrs*esc.cos(qrs)) / esc.pow(qrs, 3))
# ── Amplitude: core-shell + shell-solvent ──────────────────────────────────
F_cs = 3.0 * Vc * (sld_core - sld_shell) * kc
F_ss = 3.0 * Vs * (sld_shell - sld_solvent) * ks
F = F_cs + F_ss
# ── Powder average ─────────────────────────────────────────────────────────
i1d = (scale / Vs
* esc.integral(esc.pow(F, 2), u, 0.0, 1.0, numpoints=61, maxiter=5, epsabs=1e-5)
+ background)
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.5, 300)
i1d.show(coordinates=qs).config(
title="Core-shell 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)¶
Same angle convention as SasView cylinder/ellipsoid: $\hat{\mathbf{u}} = (\sin\theta\cos\phi,\;\sin\theta\sin\phi,\;\cos\theta)$.
The fraction of $|\mathbf{q}|^2$ along the polar axis selects the effective radius for each interface:
$$I_{\mathrm{2D}}(q_x,q_y) = \frac{\mathrm{scale}}{V_{shell}}\,F^2 + \mathrm{background}$$
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
theta = esc.par("Theta", 0.0, userlim=[0.0, 180.0], units="deg")
phi = esc.par("Phi", 0.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_mag_sq = esc.pow(qx, 2) + esc.pow(qy, 2)
q_parallel = qx * ux + qy * uy
frac = esc.conditional(q_mag_sq < 1e-20, 0.0, esc.pow(q_parallel, 2) / q_mag_sq)
qabs = esc.sqrt(q_mag_sq)
rc_2d = esc.sqrt(esc.pow(re_core, 2) * (1.0 - frac) + esc.pow(rp_core, 2) * frac)
rs_2d = esc.sqrt(esc.pow(re_shell, 2) * (1.0 - frac) + esc.pow(rp_shell, 2) * frac)
qrc2 = qabs * rc_2d
qrs2 = qabs * rs_2d
kc2 = esc.conditional(esc.abs(qrc2) < 1e-10, 1.0/3.0, (esc.sin(qrc2) - qrc2*esc.cos(qrc2)) / esc.pow(qrc2, 3))
ks2 = esc.conditional(esc.abs(qrs2) < 1e-10, 1.0/3.0, (esc.sin(qrs2) - qrs2*esc.cos(qrs2)) / esc.pow(qrs2, 3))
F_2d = 3.0 * Vc * (sld_core - sld_shell) * kc2 + 3.0 * Vs * (sld_shell - sld_solvent) * ks2
i2d = scale / Vs * esc.pow(F_2d, 2) + background
i2d.device = "gpu"
xs = np.linspace(-0.5, 0.5, 300); ys = np.linspace(-0.5, 0.5, 300)
xv, yv = np.meshgrid(xs, ys)
coords_2d = np.column_stack([xv.flatten(), yv.flatten()]).flatten()
i2d.show(coordinates=coords_2d).config(
title="Core-shell 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 |
|---|---|---|
re_core |
radius_equat_core |
equatorial core radius (Å) |
x_core |
x_core |
axial ratio Rp/Re of core |
thick_shell |
thick_shell |
equatorial shell thickness (Å) |
x_polar_shell |
x_polar_shell |
polar/equatorial shell thickness ratio |
sld_core * 1e-6 |
sld_core |
core SLD (Å⁻²) |
sld_shell * 1e-6 |
sld_shell |
shell SLD (Å⁻²) |
sld_solvent * 1e-6 |
sld_solvent |
solvent SLD (Å⁻²) |
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, 300).copy()
kernel = load_model("core_shell_ellipsoid")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
radius_equat_core=200.0, x_core=0.1,
thick_shell=30.0, x_polar_shell=1.0,
sld_core=2.0, sld_shell=1.0, sld_solvent=6.3)
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"Core-shell 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_23776\3515581655.py:17: 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 : 51 ms (300 q-pts) Max relative diff vs SasView: 0.02%
Out[5]:
In [ ]: