#!/usr/bin/env python3
"""
Generate α_ISL vs V_flat Correlation Plot
"""

import json
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

def generate_plot():
    # Load data
    json_path = "/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/SPARC_CORRELATION_RESULTS.json"
    with open(json_path, 'r') as f:
        data = json.load(f)
    
    # Reconstruct data points from the report logic (since JSON summary doesn't have raw points, 
    # I need to re-read the full analysis results)
    
    full_results_path = "/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/SPARC_FULL_ANALYSIS_RESULTS.json"
    with open(full_results_path, 'r') as f:
        full_results = json.load(f)
    
    alphas = []
    v_flats = []
    names = []
    qualities = []
    
    # We need to re-calculate v_flat effectively or grab it if we saved it?
    # The correlation script didn't save the per-galaxy v_flat in the JSON output, only the summary.
    # So I will re-calculate V_flat quickly here.
    
    print("Re-calculating V_flat for plot...")
    
    for r in full_results:
        if r['fit_quality'] == 'failed':
            continue
            
        # Load raw data for V_flat
        dat_path = f"/home/shri/Desktop/MATHTRUTH/sparc_data/{r['name']}_rotmod.dat"
        try:
            raw_data = np.loadtxt(dat_path, comments='#')
            if len(raw_data) == 0: continue
            
            v_obs = raw_data[:, 1]
            if len(v_obs) >= 3:
                v_flat = np.mean(v_obs[-3:])
            else:
                v_flat = np.max(v_obs)
            
            alphas.append(r['isl_alpha_isl'])
            v_flats.append(v_flat)
            names.append(r['name'])
            qualities.append(r['fit_quality'])
            
        except Exception:
            continue
            
    # Plotting
    plt.figure(figsize=(12, 8))
    
    # Scatter points
    plt.scatter(v_flats, alphas, c='#10b981', alpha=0.6, edgecolors='k', s=60, label='SPARC Galaxies')
    
    # Power Law Fit
    # From correlation report: A = 3.7501, beta = -0.7432
    A = 3.7501
    beta = -0.7432
    
    x_line = np.linspace(min(v_flats), max(v_flats), 100)
    y_line = A * (x_line ** beta)
    
    plt.plot(x_line, y_line, 'r--', linewidth=3, label=f'Power Law: $α \propto V^{{{beta:.2f}}}$')
    
    # Decorations
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('Galaxy Flat Velocity $V_{flat}$ (km/s)', fontsize=14)
    plt.ylabel('ISL Modularity Parameter $α_{ISL}$', fontsize=14)
    plt.title('Discovery of the Galaxy Scaling Law: Information Overhead vs Complexity', fontsize=16, fontweight='bold')
    plt.grid(True, which="both", ls="-", alpha=0.2)
    plt.legend(fontsize=12)
    
    # Annotations
    plt.text(20, 0.8, "High Overhead\n(Dark Matter Dominated)", fontsize=11, color='#666')
    plt.text(200, 0.05, "Low Overhead\n(Newtonian)", fontsize=11, color='#666')
    
    # Save
    output_path = "/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/alpha_vs_vflat_scan.png"
    plt.savefig(output_path, dpi=300, bbox_inches='tight')
    print(f"Plot saved to {output_path}")

if __name__ == "__main__":
    generate_plot()
