#!/usr/bin/env python3
"""
Subatomic Scale Leash Test

Tests whether quarks in protons need s² multiplier or if simple geometric binding works better:
1. F = ℏ²s²/(γmr³) + σ (v21 formula with spin)
2. F = ℏ²/(γmr³) + σ (simple geometric + leash, like atoms)

Uses v21 parameters but tests both approaches systematically.
Goal: Determine if the "leash" metaphor works better with or without s².

Author: Andre Heinecke & AI Collaborators  
Date: June 2025
License: CC BY-SA 4.0
"""

import numpy as np
import sys

# Try to import scipy for constants
try:
    import scipy.constants as const
    SCIPY_AVAILABLE = True
except ImportError:
    SCIPY_AVAILABLE = False
    print("WARNING: scipy not available, using hardcoded constants")

# ==============================================================================
# PHYSICAL CONSTANTS AND PARAMETERS
# ==============================================================================

# Basic constants
if SCIPY_AVAILABLE:
    HBAR = const.hbar
    C = const.c
    E_CHARGE = const.e
else:
    HBAR = 1.054571817e-34  # J⋅s
    C = 2.99792458e8        # m/s  
    E_CHARGE = 1.602176634e-19  # C

# V21 parameters (exact from paper)
V21_PARAMS = {
    'hbar': 1.054e-34,           # J⋅s
    'quark_spin': 0.5,           # s = 1/2
    'quark_mass': 4e-30,         # kg (few MeV/c²)
    'proton_radius': 1.0e-15,    # m
    'string_tension': 1.4e5,     # N
    'expected_force': 8.2e5      # N (target from v21)
}

# Alternative parameter sets for testing
ALT_PARAMS = {
    'pdg_current_quarks': {
        'up_mass': 2.16e-30,     # kg (2.16 MeV/c²) 
        'down_mass': 4.67e-30,   # kg (4.67 MeV/c²)
        'avg_mass': 3.415e-30,   # kg (average)
        'source': 'PDG 2022 current masses'
    },
    'constituent_quarks': {
        'light_mass': 596e-30,   # kg (~336 MeV/c²)
        'source': 'Constituent quark masses'
    },
    'radii_range': [0.5e-15, 0.8e-15, 1.0e-15, 1.2e-15, 1.5e-15],  # m
    'sigma_range': [0.5e5, 1.0e5, 1.4e5, 2.0e5, 3.0e5]             # N
}

# ==============================================================================
# CALCULATION FUNCTIONS
# ==============================================================================

def calculate_with_spin(hbar, s, m, r, sigma, gamma=1.0):
    """Calculate force using v21 formula: F = ℏ²s²/(γmr³) + σ"""
    
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    F_confinement = sigma
    F_total = F_geometric + F_confinement
    
    return {
        'F_geometric': F_geometric,
        'F_confinement': F_confinement, 
        'F_total': F_total,
        'geometric_fraction': F_geometric / F_total,
        'confinement_fraction': F_confinement / F_total
    }

def calculate_without_spin(hbar, m, r, sigma, gamma=1.0):
    """Calculate force using simple formula: F = ℏ²/(γmr³) + σ"""
    
    F_geometric = hbar**2 / (gamma * m * r**3)
    F_confinement = sigma
    F_total = F_geometric + F_confinement
    
    return {
        'F_geometric': F_geometric,
        'F_confinement': F_confinement,
        'F_total': F_total, 
        'geometric_fraction': F_geometric / F_total,
        'confinement_fraction': F_confinement / F_total
    }

def test_v21_baseline():
    """Test exact v21 parameters with both formulas"""
    
    print("V21 BASELINE TEST")
    print("="*50)
    print("Parameters from paper v21:")
    
    p = V21_PARAMS
    print(f"  ℏ = {p['hbar']:.3e} J⋅s")
    print(f"  s = {p['quark_spin']}")
    print(f"  m = {p['quark_mass']:.3e} kg")
    print(f"  r = {p['proton_radius']*1e15:.1f} fm") 
    print(f"  σ = {p['string_tension']:.1e} N")
    print(f"  Target = {p['expected_force']:.1e} N")
    
    # Test with spin (s²)
    result_spin = calculate_with_spin(
        p['hbar'], p['quark_spin'], p['quark_mass'], 
        p['proton_radius'], p['string_tension']
    )
    
    # Test without spin 
    result_no_spin = calculate_without_spin(
        p['hbar'], p['quark_mass'], p['proton_radius'], p['string_tension']
    )
    
    print(f"\nFORMULA 1: F = ℏ²s²/(γmr³) + σ")
    print(f"  F_geometric = {result_spin['F_geometric']:.2e} N")
    print(f"  F_confinement = {result_spin['F_confinement']:.1e} N")
    print(f"  F_total = {result_spin['F_total']:.2e} N")
    print(f"  Agreement = {result_spin['F_total']/p['expected_force']*100:.1f}%")
    print(f"  Geometric fraction = {result_spin['geometric_fraction']*100:.1f}%")
    
    print(f"\nFORMULA 2: F = ℏ²/(γmr³) + σ (no s²)")
    print(f"  F_geometric = {result_no_spin['F_geometric']:.2e} N")
    print(f"  F_confinement = {result_no_spin['F_confinement']:.1e} N")
    print(f"  F_total = {result_no_spin['F_total']:.2e} N")
    print(f"  Agreement = {result_no_spin['F_total']/p['expected_force']*100:.1f}%")
    print(f"  Geometric fraction = {result_no_spin['geometric_fraction']*100:.1f}%")
    
    return result_spin, result_no_spin

def test_mass_dependence():
    """Test how results depend on quark mass choice"""
    
    print(f"\n" + "="*50)
    print("MASS DEPENDENCE TEST")
    print("="*50)
    
    masses = [
        (V21_PARAMS['quark_mass'], "v21 effective"),
        (ALT_PARAMS['pdg_current_quarks']['avg_mass'], "PDG average"), 
        (ALT_PARAMS['constituent_quarks']['light_mass'], "constituent")
    ]
    
    print(f"{'Mass Type':<15} {'Mass(MeV)':<10} {'With s²':<12} {'Without s²':<12} {'Ratio':<8}")
    print("-" * 60)
    
    for mass, label in masses:
        mass_mev = mass * C**2 / E_CHARGE / 1e6  # Convert to MeV/c²
        
        # Test both formulas
        with_spin = calculate_with_spin(
            V21_PARAMS['hbar'], V21_PARAMS['quark_spin'], mass,
            V21_PARAMS['proton_radius'], V21_PARAMS['string_tension']
        )
        
        without_spin = calculate_without_spin(
            V21_PARAMS['hbar'], mass, V21_PARAMS['proton_radius'], 
            V21_PARAMS['string_tension']
        )
        
        # Agreement with target
        agree_with = with_spin['F_total'] / V21_PARAMS['expected_force'] * 100
        agree_without = without_spin['F_total'] / V21_PARAMS['expected_force'] * 100
        ratio = without_spin['F_total'] / with_spin['F_total']
        
        print(f"{label:<15} {mass_mev:<10.1f} {agree_with:<12.1f} {agree_without:<12.1f} {ratio:<8.2f}")

def test_radius_scaling():
    """Test how forces scale with proton radius"""
    
    print(f"\n" + "="*50)
    print("RADIUS SCALING TEST")
    print("="*50)
    
    print(f"{'r(fm)':<8} {'F_with_s²':<12} {'F_no_s²':<12} {'Ratio':<8} {'Better?':<10}")
    print("-" * 55)
    
    for r in ALT_PARAMS['radii_range']:
        r_fm = r * 1e15
        
        with_spin = calculate_with_spin(
            V21_PARAMS['hbar'], V21_PARAMS['quark_spin'], V21_PARAMS['quark_mass'],
            r, V21_PARAMS['string_tension']
        )
        
        without_spin = calculate_without_spin(
            V21_PARAMS['hbar'], V21_PARAMS['quark_mass'], r, 
            V21_PARAMS['string_tension']
        )
        
        # Which is closer to target?
        error_with = abs(with_spin['F_total'] - V21_PARAMS['expected_force'])
        error_without = abs(without_spin['F_total'] - V21_PARAMS['expected_force'])
        better = "no s²" if error_without < error_with else "with s²"
        
        ratio = without_spin['F_total'] / with_spin['F_total']
        
        print(f"{r_fm:<8.1f} {with_spin['F_total']:<12.2e} {without_spin['F_total']:<12.2e} "
              f"{ratio:<8.2f} {better:<10}")

def test_sigma_optimization():
    """Find optimal σ for both formulas"""
    
    print(f"\n" + "="*50)
    print("SIGMA OPTIMIZATION TEST")
    print("="*50)
    
    target = V21_PARAMS['expected_force']
    
    print(f"Finding σ that gives F_total = {target:.1e} N")
    print(f"{'σ(N)':<10} {'With s²':<15} {'Error%':<10} {'Without s²':<15} {'Error%':<10}")
    print("-" * 65)
    
    best_with_spin = (float('inf'), 0)
    best_without_spin = (float('inf'), 0)
    
    for sigma in ALT_PARAMS['sigma_range']:
        
        with_spin = calculate_with_spin(
            V21_PARAMS['hbar'], V21_PARAMS['quark_spin'], V21_PARAMS['quark_mass'],
            V21_PARAMS['proton_radius'], sigma
        )
        
        without_spin = calculate_without_spin(
            V21_PARAMS['hbar'], V21_PARAMS['quark_mass'], 
            V21_PARAMS['proton_radius'], sigma
        )
        
        error_with = abs(with_spin['F_total'] - target) / target * 100
        error_without = abs(without_spin['F_total'] - target) / target * 100
        
        if error_with < best_with_spin[0]:
            best_with_spin = (error_with, sigma)
        if error_without < best_without_spin[0]:
            best_without_spin = (error_without, sigma)
        
        print(f"{sigma:<10.1e} {with_spin['F_total']:<15.2e} {error_with:<10.2f} "
              f"{without_spin['F_total']:<15.2e} {error_without:<10.2f}")
    
    print(f"\nBest results:")
    print(f"  With s²: σ = {best_with_spin[1]:.1e} N, error = {best_with_spin[0]:.2f}%")
    print(f"  Without s²: σ = {best_without_spin[1]:.1e} N, error = {best_without_spin[0]:.2f}%")

def test_crossover_analysis():
    """Analyze where geometric vs confinement dominates"""
    
    print(f"\n" + "="*50)
    print("CROSSOVER ANALYSIS")
    print("="*50)
    
    radii = np.logspace(-1, 1, 15)  # 0.1 to 10 fm
    
    print("Where does geometric binding = confinement force?")
    print(f"{'r(fm)':<8} {'F_geom(s²)':<12} {'F_geom(no s²)':<15} {'F_conf':<12} {'Dominant':<10}")
    print("-" * 65)
    
    for r_fm in radii:
        r = r_fm * 1e-15
        
        # Geometric forces
        F_geom_spin = (V21_PARAMS['hbar']**2 * V21_PARAMS['quark_spin']**2) / (V21_PARAMS['quark_mass'] * r**3)
        F_geom_no_spin = V21_PARAMS['hbar']**2 / (V21_PARAMS['quark_mass'] * r**3)
        F_conf = V21_PARAMS['string_tension']
        
        # Which dominates?
        if F_geom_no_spin > F_conf:
            dominant = "geometric"
        else:
            dominant = "confinement"
        
        print(f"{r_fm:<8.2f} {F_geom_spin:<12.2e} {F_geom_no_spin:<15.2e} "
              f"{F_conf:<12.2e} {dominant:<10}")

def test_spin_values():
    """Test different spin values to see if s=1/2 is optimal"""
    
    print(f"\n" + "="*50)
    print("SPIN VALUE TEST")
    print("="*50)
    
    spin_values = [0.25, 0.5, 0.75, 1.0, 1.5, 2.0]
    target = V21_PARAMS['expected_force']
    
    print(f"{'s':<6} {'F_total':<12} {'Agreement%':<12} {'Error%':<10}")
    print("-" * 45)
    
    for s in spin_values:
        result = calculate_with_spin(
            V21_PARAMS['hbar'], s, V21_PARAMS['quark_mass'],
            V21_PARAMS['proton_radius'], V21_PARAMS['string_tension']
        )
        
        agreement = result['F_total'] / target * 100
        error = abs(result['F_total'] - target) / target * 100
        
        marker = " ← v21" if s == 0.5 else ""
        print(f"{s:<6.2f} {result['F_total']:<12.2e} {agreement:<12.1f} {error:<10.2f}{marker}")

# ==============================================================================
# MAIN TEST ROUTINE
# ==============================================================================

def main():
    """Run all subatomic leash tests"""
    
    print("SUBATOMIC SCALE LEASH TEST")
    print("="*70)
    print("Testing whether quarks need s² multiplier or if simple leash works better")
    print("Based on v21 parameters but testing formula variations")
    print()
    
    # Run all tests
    test_v21_baseline()
    test_mass_dependence()
    test_radius_scaling()
    test_sigma_optimization()
    test_crossover_analysis()
    test_spin_values()
    
    # Summary
    print(f"\n" + "="*70)
    print("SUMMARY")
    print("="*70)
    
    print("Key questions answered:")
    print("1. Does F = ℏ²/(γmr³) + σ work as well as F = ℏ²s²/(γmr³) + σ?")
    print("2. Which formula is more robust to parameter changes?")
    print("3. What does this tell us about the quark 'leash' mechanism?")
    print("4. Should we follow the atomic scale pattern (no s²)?")
    
    print(f"\nThe data speaks for itself - observer should evaluate which approach")
    print(f"better captures the physics of quark confinement in protons.")

if __name__ == "__main__":
    if len(sys.argv) > 1 and sys.argv[1] in ['-h', '--help']:
        print("Usage: python subatomic_leash_test.py")
        print("  Tests s² dependence in quark confinement formula")
        print("  Compares F = ℏ²s²/(γmr³) + σ vs F = ℏ²/(γmr³) + σ")
        sys.exit(0)
    
    main()