#!/usr/bin/env python3
"""
Elastic Leash Test for Quarks

Tests elastic string tension models where the "leash" force scales with motion:

CONSTANT LEASH (previous):  F = geometric_term + σ_constant
ELASTIC LEASH (new):        F = geometric_term + elastic_tension(v, m, ω)

Key insight: Real strings/leashes have tension proportional to:
- Centripetal force required: T ∝ mv²/r  
- Angular momentum: T ∝ L²/(mr²)
- Rotational energy: T ∝ ½mv²
- Spring-like: T ∝ k(r - r₀)

This creates self-reinforcing dynamics where faster rotation → higher tension
→ stronger binding → stable high-energy configurations.

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

import numpy as np
import sys
import math

# ==============================================================================
# PHYSICAL CONSTANTS
# ==============================================================================

HBAR = 1.054e-34           # J⋅s
C = 2.99792458e8           # m/s
E_CHARGE = 1.602176634e-19 # C

# V21 baseline
V21_PARAMS = {
    'hbar': 1.054e-34,
    'quark_mass': 4e-30,        # kg
    'proton_radius': 1.0e-15,   # m
    'target_force': 8.2e5,      # N
    'quark_spin': 0.5
}

# ==============================================================================
# ELASTIC LEASH MODELS
# ==============================================================================

def constant_leash_model(hbar, s, m, r, sigma_const, gamma=1.0):
    """
    Original model: F = ℏ²s²/(γmr³) + σ_constant
    String tension is independent of motion
    """
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    F_leash = sigma_const
    F_total = F_geometric + F_leash
    
    return {
        'model': 'constant_leash',
        'F_geometric': F_geometric,
        'F_leash': F_leash,
        'F_total': F_total,
        'v': hbar * s / (m * r),  # Velocity from L = mvr = ℏs
        'omega': hbar * s / (m * r**2),  # Angular velocity
        'description': f'Constant σ = {sigma_const:.1e} N'
    }

def centripetal_elastic_model(hbar, s, m, r, k_factor, gamma=1.0):
    """
    Elastic model 1: String tension proportional to centripetal force
    F = ℏ²s²/(γmr³) + k * (mv²/r)
    
    The string provides additional centripetal force proportional to required force
    """
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    
    # Velocity from angular momentum constraint
    v = hbar * s / (m * r)
    
    # Elastic tension proportional to centripetal requirement
    F_centripetal_needed = m * v**2 / r
    F_leash = k_factor * F_centripetal_needed
    
    F_total = F_geometric + F_leash
    
    return {
        'model': 'centripetal_elastic',
        'F_geometric': F_geometric,
        'F_leash': F_leash,
        'F_total': F_total,
        'v': v,
        'v_over_c': v / C,
        'F_centripetal_needed': F_centripetal_needed,
        'description': f'F_leash = {k_factor} × (mv²/r)'
    }

def angular_momentum_elastic_model(hbar, s, m, r, k_factor, gamma=1.0):
    """
    Elastic model 2: String tension proportional to angular momentum squared
    F = ℏ²s²/(γmr³) + k * L²/(mr²)
    
    Higher angular momentum creates higher string tension
    """
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    
    # Angular momentum
    L = hbar * s
    
    # Elastic tension proportional to L²
    F_leash = k_factor * L**2 / (m * r**2)
    
    F_total = F_geometric + F_leash
    
    return {
        'model': 'angular_momentum_elastic',
        'F_geometric': F_geometric, 
        'F_leash': F_leash,
        'F_total': F_total,
        'L': L,
        'description': f'F_leash = {k_factor} × L²/(mr²)'
    }

def rotational_energy_elastic_model(hbar, s, m, r, k_factor, gamma=1.0):
    """
    Elastic model 3: String tension proportional to rotational kinetic energy
    F = ℏ²s²/(γmr³) + k * (½mv²)/r
    
    String tension scales with rotational energy
    """
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    
    # Velocity and rotational energy
    v = hbar * s / (m * r)
    E_rot = 0.5 * m * v**2
    
    # Elastic tension proportional to rotational energy per unit length
    F_leash = k_factor * E_rot / r
    
    F_total = F_geometric + F_leash
    
    return {
        'model': 'rotational_energy_elastic',
        'F_geometric': F_geometric,
        'F_leash': F_leash, 
        'F_total': F_total,
        'E_rot': E_rot,
        'description': f'F_leash = {k_factor} × E_rot/r'
    }

def spring_elastic_model(hbar, s, m, r, k_spring, r0, gamma=1.0):
    """
    Elastic model 4: Hooke's law spring connecting particles
    F = ℏ²s²/(γmr³) + k_spring * (r - r₀)
    
    String acts like spring with natural length r₀
    """
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    
    # Spring force (positive for r > r₀, negative for r < r₀)
    F_leash = k_spring * abs(r - r0)  # Always attractive/restoring
    
    F_total = F_geometric + F_leash
    
    return {
        'model': 'spring_elastic',
        'F_geometric': F_geometric,
        'F_leash': F_leash,
        'F_total': F_total,
        'extension': r - r0,
        'description': f'F_leash = {k_spring:.1e} × |r - {r0*1e15:.1f}fm|'
    }

def velocity_squared_elastic_model(hbar, s, m, r, k_factor, gamma=1.0):
    """
    Elastic model 5: String tension proportional to velocity squared
    F = ℏ²s²/(γmr³) + k * v²
    
    Tension increases with speed (like air resistance, but attractive)
    """
    F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
    
    # Velocity from angular momentum
    v = hbar * s / (m * r)
    
    # Elastic tension proportional to v²
    F_leash = k_factor * v**2
    
    F_total = F_geometric + F_leash
    
    return {
        'model': 'velocity_squared_elastic',
        'F_geometric': F_geometric,
        'F_leash': F_leash,
        'F_total': F_total,
        'v': v,
        'v_squared': v**2,
        'description': f'F_leash = {k_factor:.2e} × v²'
    }

# ==============================================================================
# TEST FUNCTIONS 
# ==============================================================================

def test_elastic_models():
    """Compare different elastic leash models with v21 target"""
    
    print("ELASTIC LEASH MODEL COMPARISON")
    print("="*70)
    print(f"Target force: {V21_PARAMS['target_force']:.1e} N")
    print(f"Goal: Find elastic models that naturally give this force")
    print()
    
    p = V21_PARAMS
    
    # Test parameters for each model
    models = []
    
    # 1. Constant leash (baseline)
    const = constant_leash_model(p['hbar'], p['quark_spin'], p['quark_mass'], 
                                p['proton_radius'], 1.4e5)
    models.append(const)
    
    # 2. Centripetal elastic - try different k factors
    for k in [0.1, 0.2, 0.5, 1.0]:
        cent = centripetal_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'],
                                       p['proton_radius'], k)
        models.append(cent)
    
    # 3. Angular momentum elastic
    for k in [0.5, 1.0, 2.0]:
        ang = angular_momentum_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'],
                                           p['proton_radius'], k)
        models.append(ang)
    
    # 4. Rotational energy elastic
    for k in [1e10, 1e11, 1e12]:
        rot = rotational_energy_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'],
                                            p['proton_radius'], k)
        models.append(rot)
    
    # 5. Spring elastic
    for k in [1e20, 1e21]:
        for r0 in [0.8e-15, 1.2e-15]:
            spring = spring_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'],
                                        p['proton_radius'], k, r0)
            models.append(spring)
    
    # 6. Velocity squared elastic
    for k in [1e-20, 1e-19, 1e-18]:
        vel = velocity_squared_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'],
                                           p['proton_radius'], k)
        models.append(vel)
    
    # Display results
    print(f"{'Model':<25} {'F_total(N)':<12} {'Agree%':<8} {'F_leash(N)':<12} {'Description'}")
    print("-" * 90)
    
    best_agreement = 0
    best_model = None
    
    for model in models:
        agreement = model['F_total'] / p['target_force'] * 100
        
        if abs(agreement - 100) < abs(best_agreement - 100):
            best_agreement = agreement
            best_model = model
        
        print(f"{model['model']:<25} {model['F_total']:<12.2e} {agreement:<8.1f} "
              f"{model['F_leash']:<12.2e} {model['description']}")
    
    print(f"\nBest agreement: {best_model['model']} at {best_agreement:.1f}%")
    return best_model

def test_radius_scaling(best_model_type):
    """Test how elastic models scale with radius"""
    
    print(f"\n" + "="*70)
    print("RADIUS SCALING FOR ELASTIC MODELS") 
    print("="*70)
    
    radii = [0.5e-15, 0.8e-15, 1.0e-15, 1.2e-15, 1.5e-15]
    p = V21_PARAMS
    
    print(f"{'r(fm)':<8} {'Constant':<12} {'Centripetal':<12} {'Angular':<12} {'Energy':<12} {'Spring':<12}")
    print("-" * 75)
    
    for r in radii:
        r_fm = r * 1e15
        
        # Test different models at this radius
        const = constant_leash_model(p['hbar'], p['quark_spin'], p['quark_mass'], r, 1.4e5)
        cent = centripetal_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'], r, 0.2)
        ang = angular_momentum_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'], r, 1.0)
        energy = rotational_energy_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'], r, 1e11)
        spring = spring_elastic_model(p['hbar'], p['quark_spin'], p['quark_mass'], r, 1e21, 1.0e-15)
        
        print(f"{r_fm:<8.1f} {const['F_total']:<12.2e} {cent['F_total']:<12.2e} "
              f"{ang['F_total']:<12.2e} {energy['F_total']:<12.2e} {spring['F_total']:<12.2e}")

def test_mass_dependence():
    """Test how elastic models respond to different quark masses"""
    
    print(f"\n" + "="*70)
    print("MASS DEPENDENCE FOR ELASTIC MODELS")
    print("="*70)
    
    # Different quark masses
    masses = [
        (V21_PARAMS['quark_mass'], "v21 effective"),
        (2.16e-30, "up quark (PDG)"),
        (4.67e-30, "down quark (PDG)"),
        (596e-30, "constituent quark")
    ]
    
    p = V21_PARAMS
    
    print(f"{'Mass Type':<15} {'Mass(MeV)':<10} {'Constant':<12} {'Centripetal':<12} {'Angular':<12}")
    print("-" * 70)
    
    for mass, label in masses:
        mass_mev = mass * C**2 / E_CHARGE / 1e6
        
        const = constant_leash_model(p['hbar'], p['quark_spin'], mass, p['proton_radius'], 1.4e5)
        cent = centripetal_elastic_model(p['hbar'], p['quark_spin'], mass, p['proton_radius'], 0.2)
        ang = angular_momentum_elastic_model(p['hbar'], p['quark_spin'], mass, p['proton_radius'], 1.0)
        
        print(f"{label:<15} {mass_mev:<10.1f} {const['F_total']:<12.2e} "
              f"{cent['F_total']:<12.2e} {ang['F_total']:<12.2e}")

def test_physical_scaling():
    """Test how forces scale with fundamental parameters"""
    
    print(f"\n" + "="*70)
    print("PHYSICAL SCALING ANALYSIS")
    print("="*70)
    
    print("How does each model scale with basic parameters?")
    print()
    
    p = V21_PARAMS
    base_r = p['proton_radius']
    base_m = p['quark_mass']
    
    # Test scaling with radius
    print("RADIUS SCALING (r → 2r):")
    for model_name, model_func, params in [
        ("Constant", constant_leash_model, [p['hbar'], p['quark_spin'], base_m, base_r, 1.4e5]),
        ("Centripetal", centripetal_elastic_model, [p['hbar'], p['quark_spin'], base_m, base_r, 0.2]),
        ("Angular", angular_momentum_elastic_model, [p['hbar'], p['quark_spin'], base_m, base_r, 1.0])
    ]:
        result1 = model_func(*params)
        params_2r = params.copy()
        params_2r[3] = 2 * base_r  # Double radius
        result2 = model_func(*params_2r)
        
        scaling = result2['F_total'] / result1['F_total']
        print(f"  {model_name:<12}: F(2r)/F(r) = {scaling:.3f}")
    
    # Test scaling with mass
    print("\nMASS SCALING (m → 2m):")
    for model_name, model_func, params in [
        ("Constant", constant_leash_model, [p['hbar'], p['quark_spin'], base_m, base_r, 1.4e5]),
        ("Centripetal", centripetal_elastic_model, [p['hbar'], p['quark_spin'], base_m, base_r, 0.2]),
        ("Angular", angular_momentum_elastic_model, [p['hbar'], p['quark_spin'], base_m, base_r, 1.0])
    ]:
        result1 = model_func(*params)
        params_2m = params.copy()
        params_2m[2] = 2 * base_m  # Double mass
        result2 = model_func(*params_2m)
        
        scaling = result2['F_total'] / result1['F_total']
        print(f"  {model_name:<12}: F(2m)/F(m) = {scaling:.3f}")

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

def main():
    """Run elastic leash tests"""
    
    print("ELASTIC LEASH TEST FOR QUARKS")
    print("="*70)
    print("Testing elastic string tension models")
    print("Key insight: String tension scales with rotational dynamics")
    print("Real leashes get tighter when spun faster or loaded heavier")
    print()
    
    best_model = test_elastic_models()
    test_radius_scaling(best_model)
    test_mass_dependence()
    test_physical_scaling()
    
    print(f"\n" + "="*70)
    print("PHYSICAL INSIGHTS")
    print("="*70)
    
    print("Elastic leash models suggest:")
    print("1. String tension is NOT constant but responds to motion")
    print("2. Faster rotation → higher tension → stronger binding")
    print("3. This creates self-stabilizing high-energy configurations")
    print("4. Explains why quarks can't be pulled apart (infinite tension)")
    print("5. The 'leash' is more like a rubber band than a rigid rod")
    
    print(f"\nThis elastic behavior might explain:")
    print(f"- Asymptotic freedom (weak at short distances)")
    print(f"- Confinement (strong at long distances)")
    print(f"- Why QCD string tension is ~1 GeV/fm")
    print(f"- Stable high-velocity quark configurations")

if __name__ == "__main__":
    if len(sys.argv) > 1 and sys.argv[1] in ['-h', '--help']:
        print("Usage: python elastic_leash_test.py")
        print("  Tests elastic string tension models for quark confinement")
        print("  Explores how leash tension scales with rotational dynamics")
        sys.exit(0)
    
    main()