#!/usr/bin/env python3
"""
Visualize Gravitational Phase Transition
Plots Alpha vs V_flat with shaded phase regions (Laminar vs Turbulent)
"""

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

def sigmoid(x, L, x0, k, b):
    y = L / (1 + np.exp(-k*(x-x0))) + b
    return y

def plot_phase_transition():
    # Load Data
    json_path = "/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/SPARC_FULL_ANALYSIS_RESULTS.json"
    with open(json_path, 'r') as f:
        results = json.load(f)
        
    v_flats = []
    alphas = []
    
    # Extract
    for r in results:
        if r['fit_quality'] == 'failed': continue
        
        # Load raw dat for V_flat
        try:
            dat_path = f"/home/shri/Desktop/MATHTRUTH/sparc_data/{r['name']}_rotmod.dat"
            raw = np.loadtxt(dat_path, comments='#')
            if len(raw) < 3: continue
            v_val = np.mean(raw[:, 1][-3:]) # Last 3 points average
            
            v_flats.append(v_val)
            alphas.append(r['isl_alpha_isl'])
        except:
            continue
            
    v_flats = np.array(v_flats)
    alphas = np.array(alphas)
    
    # Plot Setup
    plt.figure(figsize=(12, 8))
    
    # 1. Define Phases
    dwarf_mask = v_flats < 60
    giant_mask = v_flats > 150
    trans_mask = (v_flats >= 60) & (v_flats <= 150)
    
    # 2. Plot Points by Phase
    plt.scatter(v_flats[dwarf_mask], alphas[dwarf_mask], c='#ef4444', s=80, alpha=0.7, label='Laminar Phase (High Overhead)', edgecolors='k')
    plt.scatter(v_flats[trans_mask], alphas[trans_mask], c='#f59e0b', s=80, alpha=0.7, label='Transition Instability', edgecolors='k')
    plt.scatter(v_flats[giant_mask], alphas[giant_mask], c='#10b981', s=80, alpha=0.7, label='Turbulent Phase (Low Overhead)', edgecolors='k')
    
    # 3. Fit Sigmoid (Phase Transition Curve)
    # We fit to log(V) to make it look like a standard phase diagram
    try:
        # A simple logistic decay function approximating the transition
        # We want High -> Low, so we flip the standard sigmoid
        def phase_curve(v, top, bottom, v_crit, steepness):
            return bottom + (top - bottom) / (1 + (v/v_crit)**steepness)
            
        popt, _ = curve_fit(phase_curve, v_flats, alphas, p0=[0.5, 0.0, 100, 3], maxfev=5000)
        
        v_range = np.logspace(np.log10(min(v_flats)), np.log10(max(v_flats)), 200)
        plt.plot(v_range, phase_curve(v_range, *popt), 'k-', linewidth=3, label=f'Phase Transition ($V_c \\approx {popt[2]:.0f}$ km/s)')
        
        # Mark Critical Velocity
        plt.axvline(x=popt[2], color='k', linestyle='--', alpha=0.5, label='Critical Velocity')
        
    except:
        print("Sigmoid fit failed, plotting Power Law fallback")
        # Fallback to power law visual
        x = np.linspace(min(v_flats), max(v_flats), 100)
        y = 3.75 * (x ** -0.74)
        plt.plot(x, y, 'k--', linewidth=2, label='Power Law Mean')

    # 4. Annotations & Tuning
    plt.xscale('log')
    plt.yscale('linear') # Linear Y shows the "Plateau" better than log-log
    plt.ylim(-0.1, 1.2)
    
    plt.xlabel('Galaxy Velocity $V_{flat}$ (km/s) [Complexity Proxy]', fontsize=13)
    plt.ylabel('Modularity Parameter $\\alpha_{ISL}$ [Overhead Cost]', fontsize=13)
    plt.title('Evidence for Gravitational Phase Transition', fontsize=16, fontweight='bold')
    
    # Shaded Regions
    plt.axvspan(10, 60, color='#ef4444', alpha=0.1)
    plt.axvspan(60, 150, color='#f59e0b', alpha=0.1)
    plt.axvspan(150, 400, color='#10b981', alpha=0.1)
    
    plt.text(20, 1.0, "LAMINAR\n(Viscous Gravity)", fontweight='bold', color='#991b1b')
    plt.text(200, 1.0, "TURBULENT\n(Optimized)", fontweight='bold', color='#065f46')
    
    plt.legend(loc='upper right')
    plt.grid(True, alpha=0.3)
    
    # Save
    out_path = "/home/shri/Desktop/MATHTRUTH/cosmic_synthesis/reports/phase_transition_plot.png"
    plt.savefig(out_path, dpi=300)
    print(f"✅ Phase Transition Plot Saved: {out_path}")

if __name__ == "__main__":
    plot_phase_transition()
