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

# ── Variables ──────────────────────────────────────────────────────────────
q = esc.var("Q")

# ── Parameters ─────────────────────────────────────────────────────────────
scale        = esc.par("Scale",       1.0,  scale=1e8, fixed=True)
radius       = esc.par("Core Radius",60.0,  units=esc.angstr)
thickness    = esc.par("Shell Thickness", 10.0, units=esc.angstr)
sld_core     = esc.par("Core SLD",    1.0,  scale=1e-6, units=f"{esc.angstr}^-2")
sld_shell    = esc.par("Shell SLD",   2.0,  scale=1e-6, units=f"{esc.angstr}^-2")
sld_solvent  = esc.par("Solvent SLD", 3.0,  scale=1e-6, units=f"{esc.angstr}^-2")
background   = esc.par("Background",  0.001, userlim=[0.0, 0.03])

# ── Geometry ───────────────────────────────────────────────────────────────
Rs = radius + thickness
Vc = 4.0 / 3.0 * np.pi * esc.pow(radius, 3)
Vs = 4.0 / 3.0 * np.pi * esc.pow(Rs, 3)

QRc = q * radius
QRs = q * Rs

jc = esc.conditional(esc.abs(QRc) < 1e-10, 1.0/3.0,
    (esc.sin(QRc) - QRc * esc.cos(QRc)) / esc.pow(QRc, 3))
js = esc.conditional(esc.abs(QRs) < 1e-10, 1.0/3.0,
    (esc.sin(QRs) - QRs * esc.cos(QRs)) / esc.pow(QRs, 3))

F   = 3.0 * Vc * (sld_core - sld_shell) * jc + 3.0 * Vs * (sld_shell - sld_solvent) * js
i1d = scale / Vs * esc.pow(F, 2) + background

i1d.device = "gpu"

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

qx = esc.var("qx")
qy = esc.var("qy")
q2d = esc.sqrt(esc.pow(qx, 2) + esc.pow(qy, 2))

QRc2 = q2d * radius
QRs2 = q2d * Rs

jc2 = esc.conditional(esc.abs(QRc2) < 1e-10, 1.0/3.0,
    (esc.sin(QRc2) - QRc2 * esc.cos(QRc2)) / esc.pow(QRc2, 3))
js2 = esc.conditional(esc.abs(QRs2) < 1e-10, 1.0/3.0,
    (esc.sin(QRs2) - QRs2 * esc.cos(QRs2)) / esc.pow(QRs2, 3))

F2  = 3.0 * Vc * (sld_core - sld_shell) * jc2 + 3.0 * Vs * (sld_shell - sld_solvent) * js2
i2d = scale / Vs * esc.pow(F2, 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="Core-shell sphere — isotropic 2D (qx, qy)",
    xlabel=f"qx [{esc.angstr}^-1]", ylabel=f"qy [{esc.angstr}^-1]",
    cblog=True, colorscale="jet")

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("core_shell_sphere")
f_sas  = DirectModel(empty_data1D(qs), kernel)
sas_pars = dict(scale=1.0, background=0.001,
                radius=60.0, thickness=10.0,
                sld_core=1.0, sld_shell=2.0, sld_solvent=3.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 sphere I(q) — {len(qs)} pts",
    labels=["SASView", "ESCAPE GPU", "ESCAPE CPU"],
    line_styles=["solid", "dash", "dot"],
    line_widths=[2, 3, 3]
)

SASView GPU : 11 ms
ESCAPE GPU  : 0 ms
ESCAPE CPU  : 4 ms  (300 q-pts)
Max relative diff vs SasView: 0.34%

C:\Users\User\AppData\Local\Temp\ipykernel_60060\3429340380.py:16: UserWarning:

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

Parameter	Variable	Value
Scale	`scale`	1
Background (cm⁻¹)	`background`	0.001
Core radius Rc (Å)	`radius`	60
Shell thickness S (Å)	`thickness`	10
Core SLD (10⁻⁶ Å⁻²)	`sld_core`	1
Shell SLD (10⁻⁶ Å⁻²)	`sld_shell`	2
Solvent SLD (10⁻⁶ Å⁻²)	`sld_solvent`	3

ESCAPE parameter	SasView parameter	Notes
`radius`	`radius`	core radius Rc (Å)
`thickness`	`thickness`	shell thickness S (Å)
`sld_core * 1e-6`	`sld_core`	core SLD (Å⁻²)
`sld_shell * 1e-6`	`sld_shell`	shell SLD (Å⁻²)
`sld_solvent * 1e-6`	`sld_solvent`	solvent SLD (Å⁻²)

SAXS. Form-factors. Core-shell sphere (SasView-aligned)¶

Parameters (SasView defaults)¶

Form-factor¶

2D isotropic scattering (qx, qy)¶

SasView reference model & comparison¶