import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# NGC 3198 Data (Representative points from SPARC/Begeman 1989)
# Rad (kpc), Vobs (km/s), errV (km/s), Vgas (km/s), Vdisk (km/s)
# Note: Vdisk is normalized to 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]

# ISL Model
# v_model^2 = (ML_disk * v_disk^2) + v_gas^2 + v_isl^2
# Where v_isl is the Constant Modularity Speed (to be fitted)
def isl_model(r, ml_disk, v_isl):
    # Interpolating v_disk and v_gas for the given r (though we use indices here)
    v_m_sq = (ml_disk * v_disk**2) + (v_gas**2) + (v_isl**2)
    return np.sqrt(v_m_sq)

# Fitting
popt, pcov = curve_fit(isl_model, rad, v_obs, sigma=err_v, p0=[0.8, 150.0])

ml_disk_fit, v_isl_fit = popt
print(f"ISL Fit Results for NGC 3198:")
print(f"  Disk M/L Ratio: {ml_disk_fit:.3f}")
print(f"  ISL Constant Speed (v_isl): {v_isl_fit:.2f} km/s")

# Calculate chi-squared
v_pred = isl_model(rad, *popt)
chi_sq = np.sum(((v_obs - v_pred) / err_v)**2)
red_chi_sq = chi_sq / (len(v_obs) - 2)
print(f"  Reduced Chi-Squared: {red_chi_sq:.3f}")

# Plotting
plt.figure(figsize=(10, 6))
plt.errorbar(rad, v_obs, yerr=err_v, fmt='ko', label='Observed (Begeman 1989)')
plt.plot(rad, v_pred, 'r-', lw=2, label=f'ISL Fit (M/L={ml_disk_fit:.2f}, v_isl={v_isl_fit:.1f})')
plt.plot(rad, np.sqrt(ml_disk_fit * v_disk**2), 'b--', label='Stellar Disk Component')
plt.plot(rad, v_gas, 'g--', label='Gas Component')
plt.axhline(v_isl_fit, color='orange', linestyle=':', label='ISL Modularity Constant')

plt.xlabel('Radius (kpc)')
plt.ylabel('Velocity (km/s)')
plt.title('NGC 3198 Rotation Curve: ISL Quantitative Confrontation')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/NGC3198_ISL_FIT.png')
plt.show()
