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. Ellipsoid (SasView-aligned)

Uniform ellipsoid of revolution with polar radius Rp and equatorial radius Re. Matches ellipsoid — SasView 6.1.3.

Reference: https://www.sasview.org/docs/user/models/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)
Polar radius (Å)	radius_polar	20
Equatorial radius (Å)	radius_equatorial	400
Theta (deg), 2D only	theta	60
Phi (deg), 2D only	phi	60

Form-factor

The oriented amplitude at angle $\alpha$ between the polar axis and $\mathbf{q}$ (substitution $u = \cos\alpha$):

$$r(u) = \sqrt{R_e^2(1-u^2) + R_p^2 u^2}$$

$$F(q,u) = \Delta\rho\,V\,\frac{3(\sin(qr) - qr\cos(qr))}{(qr)^3}, \quad V = \frac{4\pi}{3}R_p R_e^2$$

Powder average:

$$I(q) = \frac{\mathrm{scale}}{V}\int_0^1 F^2(q,u)\,du + \mathrm{background}$$

# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")
u = esc.var("u")   # u = cos(alpha), integration variable

# ── Parameters ─────────────────────────────────────────────────────────────
scale            = esc.par("Scale",              1.0,   scale=1e8, fixed=True)
radius_polar     = esc.par("Polar Radius",      20.0,   units=esc.angstr)
radius_equatorial= esc.par("Equatorial Radius", 400.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 ───────────────────────────────────────────────────────────────
volume = 4.0 / 3.0 * np.pi * radius_polar * esc.pow(radius_equatorial, 2)

# ── Effective radius and spherical kernel ──────────────────────────────────
r_eff = esc.sqrt(esc.pow(radius_equatorial, 2) * (1.0 - esc.pow(u, 2))
                 + esc.pow(radius_polar, 2) * esc.pow(u, 2))
qr    = q * r_eff

kern = esc.conditional(
    esc.abs(qr) < 1e-10,
    1.0 / 3.0,
    (esc.sin(qr) - qr * esc.cos(qr)) / esc.pow(qr, 3),
)
F = 3.0 * volume * contrast * kern

# ── Powder average ─────────────────────────────────────────────────────────
i1d = (scale / volume
       * esc.integral(esc.pow(F, 2), u, 0.0, 1.0, numpoints=61, maxiter=5, epsabs=1e-5)
       + background)

i1d.device = "gpu"

qs = np.linspace(0.001, 0.7, 300)
i1d.show(coordinates=qs).config(
    title="Ellipsoid — powder average (1D)",
    xlog=True, ylog=True,
    xlabel=f"Q [{esc.angstr}^-1]", ylabel="I(q) [cm^-1]")

2D oriented scattering (qx, qy)

For a fixed ellipsoid orientation $(\theta, \phi)$ the amplitude is evaluated directly at detector coordinates $(q_x, q_y)$. The polar axis unit vector is $\hat{\mathbf{u}} = (\sin\theta\cos\phi,\;\sin\theta\sin\phi,\;\cos\theta)$.

The fraction of $|\mathbf{q}|^2$ along the polar axis determines the effective radius:

$$\mathrm{frac} = \frac{(\mathbf{q}\cdot\hat{\mathbf{u}})^2}{|\mathbf{q}|^2}, \qquad r = \sqrt{R_e^2(1-\mathrm{frac}) + R_p^2\,\mathrm{frac}}$$

$$I_{\mathrm{2D}}(q_x,q_y) = \frac{\mathrm{scale}}{V}\,F^2 + \mathrm{background}$$

qx = esc.var("qx")
qy = esc.var("qy")

theta = esc.par("Theta", 60.0, userlim=[0.0, 180.0], units="deg")
phi   = esc.par("Phi",   60.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)

r_2d  = esc.sqrt(esc.pow(radius_equatorial, 2) * (1.0 - frac) + esc.pow(radius_polar, 2) * frac)
qr_2d = esc.sqrt(q_mag_sq) * r_2d

kern_2d = 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),
)
F_2d = 3.0 * volume * contrast * kern_2d
i2d  = scale / volume * esc.pow(F_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="Ellipsoid — oriented 2D SAXS (qx, qy)",
    xlabel=f"qx [{esc.angstr}^-1]", ylabel=f"qy [{esc.angstr}^-1]",
    cblog=True, colorscale="jet")

SasView reference model & comparison

ESCAPE parameter	SasView parameter	Notes
contrast * 1e-6	sld - sld_solvent	contrast in Å⁻²
radius_polar	radius_polar	polar radius (Å)
radius_equatorial	radius_equatorial	equatorial radius (Å)

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()

# ── SasView model ───────────────────────────────────────────────────────────
kernel = load_model("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_polar=20.0, radius_equatorial=400.0)

# ── Warm-up ─────────────────────────────────────────────────────────────────
f_sas(**sas_pars)
i1d.device = "gpu"; i1d(qs[:5])

# ── Timing helper ────────────────────────────────────────────────────────────
def timeit(fn, n=5):
    t0 = time.perf_counter()
    for _ in range(n): result = fn()
    return (time.perf_counter() - t0) / n * 1e3, result

# ── Evaluate ────────────────────────────────────────────────────────────────
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"Ellipsoid I(q) — {len(qs)} pts",
    labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
    line_styles=["solid", "dash", "dot"],
    line_widths=[2, 3, 3]
)

SASView GPU : 10 ms
ESCAPE GPU  : 1 ms
ESCAPE CPU  : 20 ms  (300 q-pts)
Max relative diff vs SasView: 131.39%

C:\Users\User\AppData\Local\Temp\ipykernel_29872\555725375.py:17: UserWarning:

Input array does not own its data (e.g. it is a view or slice); data will be copied

The discrepency appears due to the way how a numerical integration is performed. SASView uses a sum over 76-gaussian quadratures, in ESCAPE adaptive integration is used.