#!/usr/bin/env python3
"""
RADIAL ACCELERATION RELATION (RAR) TEST
========================================

The RAR (McGaugh+ 2016) is one of the tightest correlations in astrophysics:
  g_obs = g_bar / μ(g_bar/a₀)

Where μ is the MOND interpolating function.

If LFM is correct, it should ALSO predict the RAR, but with a different
physical interpretation: g_obs = g_bar × (c/c_eff)

Let's test LFM against the RAR.
"""

import numpy as np
import json
from pathlib import Path
from datetime import datetime
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# Physical constants
G = 6.67430e-11
c = 2.998e8
a0 = 1.2e-10

def rar_mond_simple(g_bar):
    """Simple MOND RAR formula (McGaugh+ 2016)"""
    return g_bar / (1 - np.exp(-np.sqrt(g_bar / a0)))

def rar_mond_standard(g_bar):
    """Standard MOND interpolating function"""
    x = g_bar / a0
    mu = x / np.sqrt(1 + x**2)
    return g_bar / mu

def rar_lfm(g_bar):
    """LFM prediction for RAR"""
    # c_eff = c × sqrt(g_bar / (g_bar + a0))
    # Enhancement = (c/c_eff)² = (g_bar + a0) / g_bar
    # So g_obs = g_bar × sqrt(enhancement) = g_bar × sqrt((g_bar + a0)/g_bar)
    # = sqrt(g_bar × (g_bar + a0))
    
    # Wait - need to think about this more carefully.
    # In rotation curves, v² = g_obs × r
    # If g_obs = g_bar × (c/c_eff), and c_eff = c × sqrt(a/(a+a0))
    # Then g_obs = g_bar × sqrt((a+a0)/a)
    
    # But a = g_bar in the circular orbit approximation
    enhancement = np.sqrt((g_bar + a0) / g_bar)
    return g_bar * enhancement

def rar_lfm_v2(g_bar):
    """LFM v2: velocity enhancement, not acceleration"""
    # v_obs = v_bar × (c/c_eff)^0.5
    # v² = g × r, so v = sqrt(g × r)
    # v_obs² = v_bar² × (c/c_eff)
    # g_obs × r = g_bar × r × (c/c_eff)
    # g_obs = g_bar × (c/c_eff)
    # (c/c_eff)² = (g_bar + a0) / g_bar
    # (c/c_eff) = sqrt((g_bar + a0) / g_bar)
    
    enhancement = np.sqrt((g_bar + a0) / g_bar)
    return g_bar * enhancement


def generate_rar_data():
    """
    Generate synthetic RAR data based on SPARC results.
    
    Real data from McGaugh+ 2016 shows incredibly tight correlation.
    """
    
    # Generate g_bar spanning many orders of magnitude
    log_g_bar = np.linspace(-12.5, -8.5, 100)
    g_bar = 10**log_g_bar
    
    # "True" g_obs follows the MOND-like RAR (empirical fit)
    g_obs_true = rar_mond_simple(g_bar)
    
    # Add realistic scatter (very small - about 0.1 dex)
    scatter = 0.1  # dex
    g_obs_data = g_obs_true * 10**(np.random.normal(0, scatter, len(g_bar)))
    
    # Also generate individual galaxy data points
    n_galaxies = 27
    data_points = []
    
    for i in range(n_galaxies):
        n_points = np.random.randint(5, 15)
        log_g_range = np.random.uniform(-12, -9), np.random.uniform(-10, -8.5)
        log_g = np.linspace(min(log_g_range), max(log_g_range), n_points)
        g_b = 10**log_g
        g_o = rar_mond_simple(g_b) * 10**(np.random.normal(0, scatter, n_points))
        
        for j in range(n_points):
            data_points.append({
                'g_bar': g_b[j],
                'g_obs': g_o[j],
                'galaxy': f'Galaxy_{i}'
            })
    
    return data_points


def analyze_rar():
    """Test LFM vs MOND on RAR data."""
    
    print("="*80)
    print("RADIAL ACCELERATION RELATION (RAR) TEST")
    print("="*80)
    print(f"\nDate: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
    
    print("""
The RAR is one of astrophysics' tightest correlations:
  - Observed radial acceleration vs baryonic acceleration
  - Scatter of only ~0.1 dex across 175 SPARC galaxies
  - MOND predicts this exactly
  - Dark matter models struggle to explain the tightness

QUESTION: Does LFM also predict the RAR?
""")
    
    # Generate data
    data = generate_rar_data()
    g_bar_data = np.array([d['g_bar'] for d in data])
    g_obs_data = np.array([d['g_obs'] for d in data])
    
    # Model predictions
    g_obs_mond = rar_mond_simple(g_bar_data)
    g_obs_lfm = rar_lfm(g_bar_data)
    
    # Calculate residuals
    resid_mond = np.log10(g_obs_data) - np.log10(g_obs_mond)
    resid_lfm = np.log10(g_obs_data) - np.log10(g_obs_lfm)
    
    rms_mond = np.sqrt(np.mean(resid_mond**2))
    rms_lfm = np.sqrt(np.mean(resid_lfm**2))
    
    print(f"\nRAR Analysis ({len(data)} data points):")
    print(f"  MOND RMS residual: {rms_mond:.3f} dex")
    print(f"  LFM RMS residual:  {rms_lfm:.3f} dex")
    print(f"  Ratio: {rms_lfm/rms_mond:.2f}x")
    
    # Compare formulas analytically
    print("\n" + "="*80)
    print("ANALYTICAL COMPARISON")
    print("="*80)
    
    print("""
MOND formula (simple):
  g_obs = g_bar / (1 - exp(-sqrt(g_bar/a₀)))
  
LFM formula:
  g_obs = g_bar × sqrt((g_bar + a₀) / g_bar) = sqrt(g_bar × (g_bar + a₀))

Let's compare in different regimes:
""")
    
    test_g = np.array([1e-12, 1e-11, 1e-10, 1e-9, 1e-8])
    
    print(f"{'g_bar (m/s²)':<15} {'MOND g_obs':<15} {'LFM g_obs':<15} {'Ratio':<10}")
    print("-" * 60)
    
    for g in test_g:
        g_m = rar_mond_simple(g)
        g_l = rar_lfm(g)
        print(f"{g:.2e}       {g_m:.2e}       {g_l:.2e}       {g_l/g_m:.2f}")
    
    # Asymptotic behavior
    print("\n" + "="*80)
    print("ASYMPTOTIC BEHAVIOR")
    print("="*80)
    
    print("""
HIGH ACCELERATION (g >> a₀):
  MOND: g_obs → g_bar (no modification)
  LFM:  g_obs → g_bar × sqrt(1 + a₀/g) → g_bar (same!)
  
LOW ACCELERATION (g << a₀):
  MOND: g_obs → sqrt(g_bar × a₀)
  LFM:  g_obs → sqrt(g_bar × a₀)   (SAME!)
  
THE ASYMPTOTES ARE IDENTICAL!
""")
    
    # Transition region
    print("TRANSITION REGION (g ~ a₀):")
    print("This is where MOND and LFM differ most.")
    
    g_trans = np.logspace(-11, -9.5, 20)
    
    diffs = []
    for g in g_trans:
        g_m = rar_mond_simple(g)
        g_l = rar_lfm(g)
        diff = (g_l - g_m) / g_m * 100
        diffs.append(diff)
        if abs(g/a0 - 1) < 0.3:
            print(f"  g/a₀ = {g/a0:.2f}: LFM is {diff:+.1f}% vs MOND")
    
    max_diff = max(abs(d) for d in diffs)
    print(f"\nMaximum difference in transition: {max_diff:.1f}%")
    
    # Create plot
    print("\n" + "="*80)
    print("GENERATING RAR PLOT")
    print("="*80)
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # Plot 1: RAR
    ax1 = axes[0]
    g_range = np.logspace(-13, -8, 200)
    
    ax1.scatter(g_bar_data, g_obs_data, c='gray', alpha=0.3, s=5, label='SPARC-like data')
    ax1.plot(g_range, rar_mond_simple(g_range), 'b-', lw=2, label='MOND')
    ax1.plot(g_range, rar_lfm(g_range), 'r--', lw=2, label='LFM')
    ax1.plot(g_range, g_range, 'k:', lw=1, label='1:1 (no DM)')
    ax1.axvline(a0, color='green', ls='-.', alpha=0.5, label=f'a₀ = {a0:.1e}')
    
    ax1.set_xscale('log')
    ax1.set_yscale('log')
    ax1.set_xlabel(r'$g_{bar}$ (m/s²)', fontsize=12)
    ax1.set_ylabel(r'$g_{obs}$ (m/s²)', fontsize=12)
    ax1.set_title('Radial Acceleration Relation', fontsize=14)
    ax1.legend(loc='lower right')
    ax1.grid(True, alpha=0.3)
    
    # Plot 2: Residuals
    ax2 = axes[1]
    ax2.scatter(g_bar_data, resid_mond, c='blue', alpha=0.5, s=10, label='MOND residuals')
    ax2.scatter(g_bar_data, resid_lfm, c='red', alpha=0.5, s=10, label='LFM residuals')
    ax2.axhline(0, color='k', ls='-', lw=1)
    ax2.axvline(a0, color='green', ls='-.', alpha=0.5)
    
    ax2.set_xscale('log')
    ax2.set_xlabel(r'$g_{bar}$ (m/s²)', fontsize=12)
    ax2.set_ylabel('log(observed/model) [dex]', fontsize=12)
    ax2.set_title(f'Residuals: MOND={rms_mond:.3f}, LFM={rms_lfm:.3f} dex RMS', fontsize=14)
    ax2.legend(loc='upper left')
    ax2.grid(True, alpha=0.3)
    ax2.set_ylim(-0.5, 0.5)
    
    plt.tight_layout()
    
    output_dir = Path(__file__).parent / "results"
    output_dir.mkdir(exist_ok=True)
    fig.savefig(output_dir / "rar_comparison.png", dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"Plot saved to: {output_dir / 'rar_comparison.png'}")
    
    # Conclusion
    print("\n" + "="*80)
    print("CONCLUSIONS")
    print("="*80)
    
    print(f"""
KEY FINDING: LFM and MOND give NEARLY IDENTICAL RAR predictions!

Both have:
  • Same high-a asymptote: g_obs → g_bar
  • Same low-a asymptote: g_obs → √(g_bar × a₀)
  • Maximum difference: {max_diff:.1f}% in transition region

This means:
  1. LFM naturally explains the RAR - just like MOND
  2. The two theories are nearly degenerate for rotation curves
  3. Need HIGH-PRECISION data in transition region to distinguish
  4. Time dilation tests may be the ONLY way to tell them apart!

IMPLICATIONS:
  • LFM inherits all of MOND's rotation curve successes
  • But LFM provides a PHYSICAL MECHANISM (c_eff from χ field)
  • And LFM makes ADDITIONAL predictions (time dilation, GW speed)
""")
    
    # Save results
    summary = {
        'experiment': 'Radial Acceleration Relation Test',
        'date': datetime.now().isoformat(),
        'n_data_points': len(data),
        'rms_mond': float(rms_mond),
        'rms_lfm': float(rms_lfm),
        'max_difference_percent': float(max_diff),
        'conclusion': 'LFM and MOND nearly identical for RAR'
    }
    
    with open(output_dir / "rar_analysis.json", 'w') as f:
        json.dump(summary, f, indent=2)
    
    return summary


if __name__ == "__main__":
    analyze_rar()
