#!/usr/bin/env python3
"""
test_relativistic_quarks.py - WITH PROPER LORENTZ FACTORS

Critical fix: Include relativistic effects properly!
"""

import numpy as np
import scipy.constants as const

def analyze_with_relativity(name, mass_mev, radius_fm, alpha_s):
    """Properly include relativistic effects"""
    
    # Constants
    hbar = const.hbar
    c = const.c
    e = const.e
    mev_to_kg = e * 1e6 / c**2
    
    # Convert units
    m0 = mass_mev * mev_to_kg  # REST mass
    r = radius_fm * 1e-15
    CF = 4.0/3.0
    
    # This gets tricky - we need to solve self-consistently
    # because v depends on s, but γ depends on v
    
    # Start with non-relativistic estimate
    s_squared_nr = CF * alpha_s * m0 * c * r / hbar
    s_nr = np.sqrt(s_squared_nr)
    
    # Iterate to find self-consistent solution
    s = s_nr
    for i in range(10):
        v = s * hbar / (m0 * r)  # Using rest mass
        beta = v / c
        if beta >= 1.0:
            print(f"ERROR: v > c for {name}!")
            beta = 0.99
        gamma = 1.0 / np.sqrt(1 - beta**2)
        
        # Recalculate s with relativistic correction
        # But how exactly? This is the key question!
        s_new = np.sqrt(CF * alpha_s * m0 * c * r * gamma / hbar)
        
        if abs(s_new - s) < 0.001:
            break
        s = s_new
    
    # Final forces with proper γ
    F_geometric = (hbar**2 * s**2) / (gamma * m0 * r**3)
    sigma = 0.18 * (e * 1e9 / 1e-15)
    F_total = F_geometric + sigma
    
    print(f"\n{name}:")
    print(f"  Rest mass: {mass_mev} MeV/c²")
    print(f"  Velocity: v = {v/c:.3f}c")
    print(f"  Lorentz γ = {gamma:.3f}")
    print(f"  s factor = {s:.3f}")
    print(f"  F_geometric = {F_geometric:.2e} N")
    print(f"  F_total = {F_total:.2e} N")
    
    return s, gamma, F_total

def main():
    print("RELATIVISTIC QUARK ANALYSIS - PROPER LORENTZ FACTORS")
    print("="*60)
    
    systems = [
        ("Pion (qq̄)", 140, 1.0, 0.5),
        ("Light quark", 336, 0.875, 0.4),
        ("J/ψ (cc̄)", 3097, 0.2, 0.3),
    ]
    
    for system in systems:
        analyze_with_relativity(*system)
    
    print("\n" + "="*60)
    print("CRITICAL INSIGHT:")
    print("At v ~ c, relativistic effects DOMINATE!")
    print("This could explain why different systems need different s")
    print("Maybe s encodes relativistic dynamics?")

if __name__ == "__main__":
    main()