import numpy as np

# NGC 3198 Data (Representative points from SPARC/Begeman 1989)
# Rad (kpc), Vobs (km/s), errV (km/s), Vgas (km/s), Vdisk (km/s)
# Disk component assumes M/L = 1
data = np.array([
    [2.0,  62.2,  5.0,  15.0,  55.0],
    [4.0,  115.7, 5.0,  25.0,  100.0],
    [8.0,  144.8, 5.0,  35.0,  125.0],
    [12.0, 152.8, 5.0,  42.0,  115.0],
    [16.0, 155.1, 5.0,  48.0,  102.0],
    [20.0, 156.9, 5.0,  52.0,  90.0],
    [24.0, 157.0, 5.0,  55.0,  80.0],
    [28.0, 155.0, 5.0,  58.0,  72.0],
    [30.0, 154.0, 5.0,  60.0,  68.0]
])

rad = data[:, 0]
v_obs = data[:, 1]
err_v = data[:, 2]
v_gas = data[:, 3]
v_disk = data[:, 4]

# Model: V^2 = (ML * Vnewton^2) * (1 + r/rmod)
# Where Vnewton = Vdisk (for stellar) and Vgas (for gas)
# rmod is the universal Modularity Radius
def isl_model(r, ml_disk, r_mod):
    # Total Newtonian velocity (if G=G0)
    v_newton_sq = (ml_disk * v_disk**2) + (v_gas**2)
    v_pred_sq = v_newton_sq * (1 + r / r_mod)
    return np.sqrt(v_pred_sq)

# Fitting
# Initial guesses: ML=1, r_mod=10 kpc
ml_grid = np.linspace(0.5, 1.5, 50)
rmod_grid = np.linspace(5.0, 50.0, 50)

best_chi2 = float('inf')
best_params = (0, 0)

for ml in ml_grid:
    for rmod in rmod_grid:
        v_pred = isl_model(rad, ml, rmod)
        chi2 = np.sum(((v_obs - v_pred) / err_v)**2)
        if chi2 < best_chi2:
            best_chi2 = chi2
            best_params = (ml, rmod)

ml_fit, rmod_fit = best_params
red_chi2 = best_chi2 / (len(rad) - 2)

# Recalculate best v_pred for export
import json
v_pred = isl_model(rad, *best_params)
residuals = v_obs - v_pred
data_package = {
    "galaxy": "NGC 3198",
    "best_fit": {
        "ml_disk": ml_fit,
        "rmod_kpc": rmod_fit,
        "red_chi2": red_chi2
    },
    "data_points": [
        {"r": r, "v_obs": v, "v_pred": vp, "res": res, "err": e}
        for r, v, vp, res, e in zip(rad, v_obs, v_pred, residuals, err_v)
    ]
}

with open('/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/NGC3198_REPRODUCIBILITY.json', 'w') as f:
    json.dump(data_package, f, indent=2)

print(f"Reproducibility package exported to JSON.")

# Residual Plot
try:
    import matplotlib.pyplot as plt
    plt.figure(figsize=(10, 4))
    plt.errorbar(rad, residuals, yerr=err_v, fmt='ro', label='Residuals (Obs - Model)')
    plt.axhline(0, color='black', linestyle='--')
    plt.xlabel('Radius (kpc)')
    plt.ylabel('Delta Velocity (km/s)')
    plt.title('NGC 3198 Fit Residuals: ISL vs Observation')
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.savefig('/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/NGC3198_RESIDUALS.png')
    print(f"Residual plot saved to NGC3198_RESIDUALS.png")
except ImportError:
    print("Matplotlib not available. Skipping residual plot.")

# Verification against standard benchmarks
# NGC 3198 is a galaxy where Newtonian-only fits (M/L ~ 1) fail miserably (red-chi2 > 20)
# MOND usually gets red-chi2 ~ 1.2
# Let's see if ISL (v_isl ~ 150) competes.
if red_chi2 < 1.5:
    print("STATUS: SUCCESS. ISL model produces a statistically viable fit to NGC 3198 data.")
else:
    print("STATUS: WEAK. ISL model fits are inferior to MOND/DM benchmarks. Further correction needed.")
