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 bicelle elliptical (SasView-aligned)¶
A core-shell bicelle with an elliptical cross-section. The in-plane anisotropy is averaged over the azimuthal angle $\psi$ for the 1D powder average.
Reference: https://www.sasview.org/docs/user/models/core_shell_bicelle_elliptical.html
Parameters (SasView defaults)¶
| Parameter | Variable | Value |
|---|---|---|
| Scale | scale |
1 |
| Background (cm⁻¹) | background |
0.001 |
| Minor radius (Å) | radius |
30 |
| Axis ratio ν = major/minor | x_core |
3 |
| Rim thickness (Å) | thick_rim |
8 |
| Face thickness (Å) | thick_face |
14 |
| Core length (Å) | length |
40 |
| Core SLD (10⁻⁶ Å⁻²) | sld_core |
4 |
| Face SLD (10⁻⁶ Å⁻²) | sld_face |
7 |
| Rim SLD (10⁻⁶ Å⁻²) | sld_rim |
1 |
| Solvent SLD (10⁻⁶ Å⁻²) | sld_solvent |
1 |
| Theta (deg), 2D only | theta |
90 |
| Phi (deg), 2D only | phi |
0 |
Form-factor¶
Same three-term bicelle amplitude as core_shell_bicelle, but the radial Bessel argument uses the $\psi$-dependent effective radius:
$$r'(\psi) = \frac{r_{\min}}{\sqrt{2}}\sqrt{(1+\nu^2)+(1-\nu^2)\cos\psi}$$
$$I(q) = \frac{\mathrm{scale}}{V_t}\cdot\frac{2}{\pi}\int_0^{\pi/2}d\psi\int_0^{\pi/2}F^2\,\sin\alpha\,d\alpha + \mathrm{background}$$
In [2]:
# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
alpha = esc.var("alpha") # tilt angle between disc normal and q
psi = esc.var("psi") # in-plane azimuthal angle (averages ellipse orientation)
# ── Parameters ─────────────────────────────────────────────────────────────
scale = esc.par("Scale", 1.0, scale=1e8, fixed=True)
radius = esc.par("Radius", 30.0, units=esc.angstr) # minor semi-axis
x_core = esc.par("Axis ratio", 3.0, userlim=[1.0, 10.0])
thick_rim = esc.par("Thick rim", 8.0, units=esc.angstr)
thick_face = esc.par("Thick face", 14.0, units=esc.angstr)
length = esc.par("Length", 40.0, units=esc.angstr)
sld_core = esc.par("SLD core", 4.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_face = esc.par("SLD face", 7.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_rim = esc.par("SLD rim", 1.0, scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent = esc.par("SLD solvent", 1.0, scale=1e-6, units=f"{esc.angstr}^-2")
background = esc.par("Background", 0.001, userlim=[0.0, 0.03])
# ── Geometry ───────────────────────────────────────────────────────────────
# Geometric-mean radius for volumes: r_eff = r_minor * sqrt(x_core)
r_eff = radius * esc.sqrt(x_core) # geometric-mean radius for volumes
h_core = 0.5 * length
h_face = h_core + thick_face
v_core = np.pi * x_core * esc.pow(radius, 2) * length
v_cf = np.pi * x_core * esc.pow(radius, 2) * (length + 2.0 * thick_face)
# Outer ellipse semi-axes: (R + t_r) and (x_core*R + t_r)
v_tot = np.pi * (radius + thick_rim) * (x_core * radius + thick_rim) * (length + 2.0 * thick_face)
# ψ-dependent effective radius for Bessel arguments
r_psi = (radius / esc.sqrt(2.0)) * esc.sqrt(
(1.0 + esc.pow(x_core, 2)) + (1.0 - esc.pow(x_core, 2)) * esc.cos(psi))
r_psi_rim = r_psi + thick_rim
# ── Oriented amplitude ─────────────────────────────────────────────────────
q_axial = q * esc.cos(alpha)
q_radial = q * esc.sin(alpha)
bj_core_psi = 2.0 * esc.j1_over_x(r_psi * q_radial)
bj_rim_psi = 2.0 * esc.j1_over_x(r_psi_rim * q_radial)
f_core = (sld_core - sld_face) * v_core * bj_core_psi * esc.sinc(h_core * q_axial)
f_face = (sld_face - sld_rim) * v_cf * bj_core_psi * esc.sinc(h_face * q_axial)
f_rim = (sld_rim - sld_solvent) * v_tot * bj_rim_psi * esc.sinc(h_face * q_axial)
f_tot = f_core + f_face + f_rim
# ── Powder average: integrate over alpha then psi ──────────────────────────
alpha_integral = esc.integral(esc.pow(f_tot, 2) * esc.sin(alpha),
alpha, 0.0, np.pi / 2.0,
numpoints=61, maxiter=0, epsabs=1e-5)
i1d = (scale / v_tot
* esc.integral(alpha_integral, psi, 0.0, np.pi,
numpoints=61, maxiter=0, epsabs=1e-5)
* (1.0 / np.pi) # normalise psi average over [0, pi] (full ellipse)
+ background)
In [3]:
i1d.device = "gpu"
qs = np.linspace(0.001, 0.8, 300)
i1d.show(coordinates=qs).config(
title="Core-shell bicelle elliptical — 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 the oriented 2D case the geometric-mean radius $r_{\mathrm{eff}} = \sqrt{ab}$ is used as an approximation (the ψ average is replaced by a fixed effective radius).
$$I_{\mathrm{2D}}(q_x,q_y) = \frac{\mathrm{scale}}{V_t}\,F^2(q_\parallel, q_\perp) + \mathrm{background}$$
In [4]:
qx = esc.var("qx")
qy = esc.var("qy")
theta = esc.par("Theta", 90.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_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))
r_rim_2d = r_eff + thick_rim # geometric-mean outer radius for oriented 2D
bj_core_2d = 2.0 * esc.j1_over_x(r_eff * q_perp_2d)
bj_rim_2d = 2.0 * esc.j1_over_x(r_rim_2d * q_perp_2d)
f_core_2d = (sld_core - sld_face) * v_core * bj_core_2d * esc.sinc(h_core * q_par_2d)
f_face_2d = (sld_face - sld_rim) * v_cf * bj_core_2d * esc.sinc(h_face * q_par_2d)
f_rim_2d = (sld_rim - sld_solvent) * v_tot * bj_rim_2d * esc.sinc(h_face * q_par_2d)
f_tot_2d = f_core_2d + f_face_2d + f_rim_2d
i2d = scale / v_tot * esc.pow(f_tot_2d, 2) + background
i2d.device = "gpu"
xs = np.linspace(-0.8, 0.8, 300); ys = np.linspace(-0.8, 0.8, 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 bicelle elliptical — 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 |
|---|---|---|
sld_core * 1e-6 |
sld_core |
core SLD (Å⁻²) |
sld_face * 1e-6 |
sld_face |
face SLD (Å⁻²) |
sld_rim * 1e-6 |
sld_rim |
rim SLD (Å⁻²) |
sld_solvent * 1e-6 |
sld_solvent |
solvent SLD (Å⁻²) |
radius |
radius |
minor semi-axis (Å) |
x_core |
x_core |
major/minor ratio |
thick_rim |
thick_rim |
rim thickness (Å) |
thick_face |
thick_face |
face thickness (Å) |
length |
length |
core 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, 0.8, 300).copy()
kernel = load_model("core_shell_bicelle_elliptical")
f_sas = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
sld_core=4.0, sld_face=7.0, sld_rim=1.0, sld_solvent=1.0,
radius=30.0, x_core=3.0, thick_rim=8.0, thick_face=14.0, length=40.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"Core shell bicelle elliptical 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_58120\1352006094.py:16: 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 : 3 ms ESCAPE CPU : 82 ms (300 q-pts) Max relative diff vs SasView: 0.00%
Out[5]:
In [ ]: