# Tiny Appendix Notebook — Spin→Projection → Observable (GR baselines + overlay)
# Requirements: numpy, matplotlib
# Run: python spin_projection_lensing_demo.py
# Or copy into a Colab cell (set backend inline).

import numpy as np
import matplotlib.pyplot as plt

# ------------------------------
# Constants (SI unless noted)
# ------------------------------
G   = 6.67430e-11              # m^3 kg^-1 s^-2
c   = 299_792_458.0            # m s^-1
Msun= 1.98847e30               # kg
Rsun= 6.957e8                  # m
AU  = 1.495978707e11           # m
arcsec = np.pi/(180.*3600.)    # rad per arcsec

# ------------------------------
# 1) GR Baselines (compact)
# ------------------------------
def alpha_GR_point_mass(M, b):
    """
    GR deflection for null geodesic near point mass (leading PN):
    alpha(b) = 4GM/(c^2 b) [radians]
    """
    return 4*G*M/(c**2 * b)

def shapiro_delay_approx(M, rE, rR, b):
    """
    Shapiro time delay (schematic, leading log approx):
    Δt ≈ (2GM/c^3) * ln( 4 r_E r_R / b^2 )
    Inputs in meters; output in seconds.
    """
    return (2*G*M/c**3)*np.log(4*rE*rR/(b**2))

def thetaE_SIS(sigma_v, Dds_over_Ds):
    """
    SIS Einstein angle: θ_E = 4π (σ_v^2 / c^2) * (D_ds / D_s)
    Returns radians.
    """
    return 4*np.pi*(sigma_v**2/c**2)*Dds_over_Ds

# ------------------------------
# 2) Spin→Projection placeholder
#    (replace this with your actual mapping)
# ------------------------------
def projection_kernel(b, params):
    """
    YOUR MODEL lives here.
    Idea: map spin-phase/flow → an effective deflection alpha_proj(b).
    'b' is impact parameter array [m].
    'params' can hold your spin frequency, elliptic constraint knobs, etc.

    This default is a toy: GR baseline plus a small phase-y modulation
    that vanishes as b grows (so far field matches GR).
    """
    M      = params.get('M', 1.0*Msun)
    eps    = params.get('eps', 0.0)      # 0.0 by default: no deviation
    kphase = params.get('kphase', 1e-9)  # set your scale
    base   = alpha_GR_point_mass(M, b)
    mod    = 1.0 + eps*np.sin(kphase*b)/(1.0 + (kphase*b)**2)
    return base*mod

# ------------------------------
# 3) One-figure traction demo
# ------------------------------
def demo_overlay_point_mass():
    # Scenario: grazing the Sun & out to 5 R_sun
    M = 1.0*Msun
    b = np.linspace(1.0*Rsun, 5.0*Rsun, 400)  # impact parameter [m]

    # GR baseline curve
    alpha_gr = alpha_GR_point_mass(M, b)     # radians
    alpha_gr_arcsec = alpha_gr/arcsec

    # Your projection curve (edit params)
    params = dict(M=M, eps=0.0, kphase=1e-9)  # set eps>0 to show deviation
    alpha_proj = projection_kernel(b, params)/arcsec

    # Plot
    fig, ax = plt.subplots(figsize=(7.5,4.5), dpi=140)
    ax.plot(b/Rsun, alpha_gr_arcsec,  lw=2.5, label='GR (point mass)', color='#4444aa')
    ax.plot(b/Rsun, alpha_proj,       lw=2.0, ls='--', label='Projection (spin→deflection)', color='#d9534f')
    ax.set_xlabel('impact parameter  b / R$_\\odot$')
    ax.set_ylabel('deflection  α  (arcsec)')
    ax.set_title('Light Bending: projection kernel overlaid on GR baseline')
    ax.grid(alpha=0.2)
    ax.legend()
    plt.tight_layout()
    plt.show()

def demo_shapiro():
    # Earth–Mars schematic at superior conjunction (toy inputs)
    M  = Msun
    rE = 1.0*AU
    rR = 1.52*AU
    b  = 1.0*Rsun
    dt = shapiro_delay_approx(M, rE, rR, b)
    print(f"Shapiro delay (approx, Earth–Mars, b=1 R_sun): {dt:.3e} s")

def demo_SIS():
    # SIS Einstein angle for cluster-like numbers
    sigma_v = 200_000.0         # m/s (200 km/s)
    Dds_over_Ds = 0.5           # unitless ratio
    thetaE = thetaE_SIS(sigma_v, Dds_over_Ds)
    print(f"SIS Einstein angle θ_E ≈ {thetaE/arcsec:.3f} arcsec")

if __name__ == "__main__":
    print("== GR checks ==")
    demo_shapiro()
    demo_SIS()
    print("\n== Overlay demo ==")
    demo_overlay_point_mass()
    print("\nEdit projection_kernel(...) to map your spin/flow parameters to α(b).")