"""
Spintronisches Spiking-Neuron: SV/MTJ-Macrospin.
Parameter nach Louis et al., IEEE Magnetics Letters 15, 4500705 (2024).

Drei unabhaengige Formulierungen derselben Physik als Kreuzcheck:
  (A) LLGS kartesisch, m in R^3            -> rhs_cart
  (B) LLGS in Kugelwinkeln (theta, phi)    -> rhs_sph   [Gl. 2a/2b]
  (C) Artificial Neuron Equation, Pendel   -> rhs_ane   [Gl. 4]
"""
import numpy as np
from scipy.integrate import solve_ivp

# ------------------------------------------------------------------ Parameter
gamma = 2*np.pi*28e9        # rad/(s T)   gyromagnetisches Verhaeltnis
alpha = 0.1                 #             Gilbert-Daempfung
Bext  = 3e-3                # T           externes Feld, entlang +x
Bd    = 1.0                 # T           Entmagnetisierung, leichte Ebene (x-y)
R0    = 1500.0              # Ohm         mittlerer MTJ-Widerstand
dR    = 500.0               # Ohm         Amplitude der Widerstandsaenderung
Ith   = 0.622e-3            # A           Schwellstrom (Paper, Gl. 5)

wB    = gamma*Bext          # rad/s       Praezession im externen Feld
wM    = gamma*Bd            # rad/s       Praezession im Entmagnetisierungsfeld
sigma = wB/Ith              # rad/(s A)   Spin-Torque-Koeffizient, aus I_th
p1    = np.array([0., 0., 1.])   # Polarisator: senkrecht zur Ebene

def I_in(t, bias=0.6e-3, amp=0.15e-3, t0=10e-9, width=0.35e-9):
    """Bias-Strom plus ein gaussfoermiger Puls."""
    return bias + amp*np.exp(-((t - t0)/width)**2)

# ------------------------------------------------------------------ (A) kartesisch
def rhs_cart(t, m):
    m = m/np.linalg.norm(m)
    Beff = np.array([Bext, 0., -Bd*m[2]])           # externes + Entmag.-Feld
    a = gamma*np.cross(Beff, m) + sigma*I_in(t)*np.cross(m, np.cross(m, p1))
    # implizites m x dm/dt aufloesen:  x = a + alpha m x x  =>  x = (a + alpha m x a)/(1+alpha^2)
    return (a + alpha*np.cross(m, a))/(1 + alpha**2)

# ------------------------------------------------------------------ (B) Kugelwinkel
def rhs_sph(t, y):
    """theta = Winkel aus der Ebene heraus, phi = Azimut. Gl. (2a)/(2b) des Papers,
    nach dtheta/dt und dphi/dt aufgeloest (2x2-System)."""
    th, ph = y
    c, s = np.cos(th), np.sin(th)
    # dth = wB sin(ph) + alpha c dph - sigma I c
    # dph = -wM s - wB tan(th) cos(ph) - alpha dth / c
    A = np.array([[1.0, -alpha*c], [alpha/c, 1.0]])
    b = np.array([wB*np.sin(ph) - sigma*I_in(t)*c,
                  -wM*s - wB*np.tan(th)*np.cos(ph)])
    return np.linalg.solve(A, b)

# ------------------------------------------------------------------ (C) Pendel / ANE
def rhs_ane(t, y):
    """(1/wM) phi'' + alpha phi' + wB sin(phi) = sigma I    -- Gl. (4)"""
    ph, dph = y
    return [dph, wM*(sigma*I_in(t) - alpha*dph - wB*np.sin(ph))]

# ------------------------------------------------------------------ Laeufe
def run(t_end=22e-9, n=6000):
    t = np.linspace(0, t_end, n)
    opt = dict(max_step=2e-12, rtol=1e-9, atol=1e-11, dense_output=True)

    A = solve_ivp(rhs_cart, (0, t_end), [1., 0., 0.], **opt)
    m = A.sol(t); m /= np.linalg.norm(m, axis=0)
    phi_A = np.unwrap(np.arctan2(m[1], m[0]))

    B = solve_ivp(rhs_sph, (0, t_end), [1e-9, 0.], **opt)
    phi_B = np.unwrap(B.sol(t)[1])

    C = solve_ivp(rhs_ane, (0, t_end), [0., 0.], **opt)
    phi_C = np.unwrap(C.sol(t)[0])

    return t, m, phi_A, phi_B, phi_C

def R_of(m, phi_p2_deg):
    """Widerstand fuer einen Analysator, der in der Ebene um phi_p2 gedreht ist."""
    a = np.radians(phi_p2_deg)
    p2 = np.array([np.cos(a), np.sin(a), 0.])
    return R0 + dR*(p2 @ m)

if __name__ == "__main__":
    t, m, pA, pB, pC = run()
    rest = t < 9e-9
    print(f"sigma_j = {sigma:.4e} rad/(s A)     I_th = {Ith*1e3:.3f} mA")
    print(f"{'':22s}{'kartesisch':>12s}{'Kugelwinkel':>14s}{'ANE/Pendel':>13s}   Paper")
    print(f"{'Ruhewinkel phi0 [deg]':22s}"
          f"{np.degrees(pA[rest][-1]):12.2f}{np.degrees(pB[rest][-1]):14.2f}"
          f"{np.degrees(pC[rest][-1]):13.2f}   73")
    print(f"{'Endwinkel phi [deg]':22s}"
          f"{np.degrees(pA).max():12.2f}{np.degrees(pB).max():14.2f}"
          f"{np.degrees(pC).max():13.2f}   433")
    print(f"max. Auslenkung aus der Ebene: {np.degrees(np.arcsin(np.abs(m[2]).max())):.2f} deg")
    print()
    print(f"{'Analysator':>12s}{'R_ruhe':>10s}{'R_min':>9s}{'R_max':>9s}{'Einbruch':>10s}{'Peak':>9s}")
    for name, ang in [("p2 = +x", 0), ("p2 = -y", 270), ("p2 = +y", 90), ("p2 = -x", 180)]:
        R = R_of(m, ang); r0 = R[rest][-1]
        print(f"{name:>12s}{r0:10.0f}{R.min():9.0f}{R.max():9.0f}"
              f"{R.min()-r0:10.0f}{R.max()-r0:9.0f}")