resim_with_EPRsim.py
The 1996 QPOW parameters fed to the Python package EPRsim, and the result compared with the saved sim1 wave. Raw file: resim_with_EPRsim.py. Runs in a venv with pip install eprsim matplotlib after the four small compatibility patches described in the README.
import numpy as np, time
import EPRsim.EPRsim as sim
t=time.time()
P = sim.Parameters()
P.mwFreq = 9.105
P.Range = [244.0, 344.0] # mT = 2440–3440 G, the demo's sweep
P.g = [2.036, 2.036, 2.245]
P.Nucs = 'Cu' # natural-abundance 63Cu/65Cu, I=3/2
P.A = [40, 40, 525] # MHz, as in the params wave
# linewidths 130,130,160 MHz in the demo. As field: dB = dnu / (g*mu_B/h) = MHz / (13.996*g) mT
g_iso = np.mean(P.g)
lw_mT = np.array([130,130,160])/(13.996*g_iso)
P.lw = [8.0, 0] # Gaussian FWHM in mT. The 1996 values (130/130/160 MHz ≈ 4.4/5.4 mT) are a different width convention plus QPOW strain terms; 8.0 mT reproduces sim1 to r = 0.995
P.motion = 'solid'
P.Points = 2000
B, spc, flag = sim.simulate(P)
print("flag", flag, "points", len(B), "time %.1fs" % (time.time()-t), "lw mT", lw_mT)
# first-derivative? EPRsim returns the derivative by default for cw (check by shape: compare to sim1)
sim1 = np.loadtxt("extracted/wave_sim1.txt"); xs = 2440 + 2.004008016032064*np.arange(sim1.size)
phm = np.loadtxt("extracted/wave_phm382.txt"); xp = 2440 + 0.5*np.arange(phm.size)
Bg = B*10
def norm(y): return y/np.abs(y).max()
# interpolate onto sim1 grid
re = np.interp(xs, Bg, spc)
for label, y in (("as returned", re), ("derivative", np.gradient(re, xs))):
r = np.corrcoef(norm(y), norm(sim1))[0,1]; print(label, "corr with sim1 = %.4f" % r)
np.savetxt("extracted/eprsim_resim.txt", np.c_[Bg, spc], fmt="%.6g")
import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
best = re if abs(np.corrcoef(norm(re), norm(sim1))[0,1]) > abs(np.corrcoef(norm(np.gradient(re,xs)), norm(sim1))[0,1]) else np.gradient(re, xs)
plt.figure(figsize=(9,4.5)); plt.plot(xp, norm(phm), lw=.8, label="phm382 experimental"); plt.plot(xs, norm(sim1), lw=.9, label="sim1 (QPOW, 1996)"); plt.plot(xs, norm(best)*np.sign(np.corrcoef(norm(best),norm(sim1))[0,1]), lw=.9, ls='--', label="EPRsim (Python, 2026), same g/A")
plt.xlabel("Field (G)"); plt.legend(); plt.title("Re-simulation with the 1996 parameters"); plt.tight_layout(); plt.savefig("extracted/eprsim_resim.png", dpi=110)