#!/usr/bin/env python3
"""
High-precision multiscale verification using externally sourced data
Version 21 - Includes quark confinement test from paper v21

Tests the geometric force principle F = ℏ²s²/(γmr³) + σ across scales:
- Subatomic: Quark confinement in proton (v21 parameters)
- Atomic: F = ℏ²/(γmr³) = ke²/r²
- Nuclear: F = ℏ²/(γm_q r³) + σr (geometric + confinement)
- Planetary: F = s²ℏ²/(γmr³) = mv²/r where s = mvr/ℏ  
- Galactic: Same as planetary but expected breakdown

Uses only external data sources. No conclusions drawn - data speaks for itself.

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

import numpy as np
import sys
import json
import urllib.request
import urllib.error
from decimal import Decimal, getcontext
import time
import warnings

# Set high precision for calculations
getcontext().prec = 100

# Try to import required libraries
try:
    import scipy.constants as const
    from scipy.constants import physical_constants
    SCIPY_AVAILABLE = True
except ImportError:
    SCIPY_AVAILABLE = False
    print("WARNING: scipy.constants not available")

try:
    import requests
    REQUESTS_AVAILABLE = True
except ImportError:
    REQUESTS_AVAILABLE = False
    print("WARNING: requests not available")

# ==============================================================================
# EXTERNAL DATA LOADERS
# ==============================================================================

def load_codata_constants():
    """Load CODATA constants via scipy.constants with error handling"""
    
    if not SCIPY_AVAILABLE:
        return None
    
    constants = {}
    
    try:
        # Direct access to exact constants
        constants['c'] = {
            'value': const.c,
            'uncertainty': 0.0,
            'unit': 'm/s',
            'source': 'SI definition (exact)',
            'relative_uncertainty': 0.0
        }
        
        constants['h'] = {
            'value': const.h,
            'uncertainty': 0.0, 
            'unit': 'J⋅s',
            'source': 'SI definition (exact)',
            'relative_uncertainty': 0.0
        }
        
        constants['hbar'] = {
            'value': const.hbar,
            'uncertainty': 0.0,
            'unit': 'J⋅s', 
            'source': 'SI definition (exact)',
            'relative_uncertainty': 0.0
        }
        
        constants['e'] = {
            'value': const.e,
            'uncertainty': 0.0,
            'unit': 'C',
            'source': 'SI definition (exact)',
            'relative_uncertainty': 0.0
        }
        
        # Measured constants with uncertainties
        for name, scipy_name in [
            ('me', 'electron mass'),
            ('alpha', 'fine-structure constant'),
            ('a0', 'Bohr radius'),
            ('epsilon_0', 'electric constant'),
            ('G', 'Newtonian constant of gravitation')
        ]:
            if scipy_name in physical_constants:
                value, unit, uncertainty = physical_constants[scipy_name]
                constants[name] = {
                    'value': value,
                    'uncertainty': uncertainty,
                    'unit': unit,
                    'source': f'CODATA 2022 via scipy.constants',
                    'relative_uncertainty': uncertainty/value if value != 0 else 0.0
                }
        
        # Calculate Coulomb constant (avoid const.k confusion)
        if 'epsilon_0' in constants:
            eps0 = constants['epsilon_0']['value']
            eps0_unc = constants['epsilon_0']['uncertainty']
            k_val = 1.0 / (4 * np.pi * eps0)
            k_unc = k_val * (eps0_unc / eps0) if eps0 != 0 else 0.0
            
            constants['k'] = {
                'value': k_val,
                'uncertainty': k_unc,
                'unit': 'N⋅m²⋅C⁻²',
                'source': 'Calculated from electric constant',
                'relative_uncertainty': k_unc / k_val if k_val != 0 else 0.0
            }
        
        return constants
        
    except Exception as e:
        print(f"ERROR fetching CODATA constants: {e}")
        return None

def load_v21_parameters():
    """Load exact parameters from paper version 21"""
    
    # Exact values from examples_v21.tex
    v21_params = {
        'subatomic_proton': {
            'hbar': 1.054e-34,                # J·s (Planck's constant)
            'quark_spin': 0.5,                # s = 1/2 for quarks
            'effective_quark_mass': 4e-30,    # kg (few MeV/c²)
            'effective_proton_radius': 1.0e-15,  # m (proton radius)
            'string_tension': 1.4e5,          # N (QCD flux-tube tension)
            'source': 'Paper v21 examples_v21.tex subatomic scale'
        },
        'nuclear_general': {
            'sigma_coefficient': 1.4e5,       # N (from v21)
            'typical_nuclear_radius': 1e-15,  # m
            'source': 'Paper v21 nuclear scale parameters'
        }
    }
    
    return v21_params

def load_pdg_constants():
    """Load particle physics constants from PDG or use authoritative values"""
    
    # PDG 2022 values for nuclear physics
    # Source: https://pdg.lbl.gov/2022/reviews/rpp2022-rev-phys-constants.pdf
    
    constants = {
        'proton_mass': {
            'value': 1.67262192595e-27,  # kg
            'uncertainty': 5.2e-37,
            'unit': 'kg',
            'source': 'PDG 2022'
        },
        'neutron_mass': {
            'value': 1.67492750056e-27,  # kg  
            'uncertainty': 8.5e-37,
            'unit': 'kg',
            'source': 'PDG 2022'
        },
        'up_quark_mass': {
            'value': 2.16e-30,  # kg (2.16 MeV/c²)
            'uncertainty': 0.49e-30,
            'unit': 'kg',
            'source': 'PDG 2022 (current mass)',
            'note': 'Current quark mass, not constituent'
        },
        'down_quark_mass': {
            'value': 4.67e-30,  # kg (4.67 MeV/c²)
            'uncertainty': 0.48e-30,
            'unit': 'kg',
            'source': 'PDG 2022 (current mass)',
            'note': 'Current quark mass, not constituent'
        },
        'qcd_string_tension': {
            'value': 0.18e9 * const.e / 1e-15 if SCIPY_AVAILABLE else 2.88e-11,  # N
            'uncertainty': 0.02e9 * const.e / 1e-15 if SCIPY_AVAILABLE else 3.2e-12,
            'unit': 'N',
            'source': 'Lattice QCD (FLAG 2021)',
            'note': 'σ = 0.18 ± 0.02 GeV²/fm'
        },
        'alpha_s': {
            'value': 0.1179,  # at M_Z scale
            'uncertainty': 0.0010,
            'unit': 'dimensionless',
            'source': 'PDG 2022 (at M_Z = 91.2 GeV)',
            'note': 'Strong coupling constant'
        }
    }
    
    return constants

def load_nasa_planetary_data():
    """Load planetary data from NASA JPL or use authoritative cached values"""
    
    # High-precision values from NASA JPL Planetary Fact Sheet
    # https://nssdc.gsfc.nasa.gov/planetary/factsheet/
    # and IAU 2015 nominal values
    
    planets = {
        'Mercury': {
            'mass': 3.30104e23,           # kg ± 1e18 (NASA JPL 2021)
            'semimajor_axis': 5.7909050e10,  # m (IAU 2015 exact)
            'eccentricity': 0.20563593,      # (NASA JPL)
            'orbital_period': 87.9691,       # days
            'mean_orbital_velocity': 47362,  # m/s
            'mass_uncertainty': 1e18,
            'source': 'NASA JPL Planetary Fact Sheet 2021, IAU 2015'
        },
        'Venus': {
            'mass': 4.86732e24,
            'semimajor_axis': 1.0820893e11,
            'eccentricity': 0.00677672,
            'orbital_period': 224.701,
            'mean_orbital_velocity': 35020,
            'mass_uncertainty': 1e19,
            'source': 'NASA JPL Planetary Fact Sheet 2021'
        },
        'Earth': {
            'mass': 5.97219e24,          # kg (CODATA 2018)
            'semimajor_axis': 1.495978707e11,  # m (IAU 2012 exact)
            'eccentricity': 0.01671123,
            'orbital_period': 365.256363004,
            'mean_orbital_velocity': 29784.77,
            'mass_uncertainty': 6e19,    # CODATA uncertainty
            'source': 'CODATA 2018, IAU 2012'
        },
        'Mars': {
            'mass': 6.41693e23,
            'semimajor_axis': 2.2793664e11,
            'eccentricity': 0.0933941,
            'orbital_period': 686.980,
            'mean_orbital_velocity': 24077,
            'mass_uncertainty': 1e18,
            'source': 'NASA JPL Planetary Fact Sheet 2021'
        },
        'Jupiter': {
            'mass': 1.89813e27,
            'semimajor_axis': 7.7857e11,
            'eccentricity': 0.0489,
            'orbital_period': 4332.589,
            'mean_orbital_velocity': 13056,
            'mass_uncertainty': 1.9e23,
            'source': 'NASA JPL Planetary Fact Sheet 2021'
        }
    }
    
    # Add solar mass from IAU 2015
    planets['_solar_constants'] = {
        'M_sun': 1.9884754153381438e30,  # kg (IAU 2015 nominal)
        'M_sun_uncertainty': 0,  # Defined value
        'AU': 1.495978707e11,    # m (IAU 2012 exact)
        'source': 'IAU 2015 nominal values'
    }
    
    return planets

def load_galactic_data():
    """Load galactic rotation data from literature or use authoritative compilation"""
    
    # Compilation from multiple authoritative sources:
    # - Reid et al. (2019) ApJ 885, 131 - Solar position
    # - Sofue (2020) PASJ 72, 4 - Galaxy rotation curve
    # - Gaia Collaboration (2021) A&A 649, A1 - Gaia EDR3
    # - McMillan (2017) MNRAS 465, 76 - Mass model
    
    galactic_data = {
        'solar_position': {
            'R0': 8178,  # pc ± 13 (Reid et al. 2019)
            'R0_uncertainty': 13,
            'V0': 220,   # km/s ± 3 (Reid et al. 2019) 
            'V0_uncertainty': 3,
            'source': 'Reid et al. (2019) ApJ 885, 131'
        },
        'rotation_curve': [
            # R (kpc), V_circular (km/s), V_error (km/s), Source
            {'R_kpc': 1.0, 'V_kms': 180, 'V_error': 15, 'source': 'Inner Galaxy (Sofue 2020)'},
            {'R_kpc': 2.0, 'V_kms': 220, 'V_error': 10, 'source': 'Gaia DR3 + APOGEE'},
            {'R_kpc': 4.0, 'V_kms': 235, 'V_error': 8, 'source': 'Gaia DR3'},
            {'R_kpc': 6.0, 'V_kms': 240, 'V_error': 10, 'source': 'Gaia DR3'},
            {'R_kpc': 8.178, 'V_kms': 220, 'V_error': 3, 'source': 'Solar position (Reid 2019)'},
            {'R_kpc': 10.0, 'V_kms': 225, 'V_error': 12, 'source': 'Outer disk tracers'},
            {'R_kpc': 15.0, 'V_kms': 220, 'V_error': 20, 'source': 'Globular clusters'},
            {'R_kpc': 20.0, 'V_kms': 210, 'V_error': 25, 'source': 'Satellite galaxies'},
            {'R_kpc': 25.0, 'V_kms': 200, 'V_error': 30, 'source': 'Extended halo'}
        ],
        'mass_model': {
            'bulge_mass': 1.5e10,        # Solar masses (McMillan 2017)
            'bulge_uncertainty': 0.3e10,
            'disk_mass': 6.43e10,        # Solar masses
            'disk_uncertainty': 0.5e10,
            'dark_halo_mass': 1.3e12,    # Solar masses (within 200 kpc)
            'dark_halo_uncertainty': 0.3e12,
            'source': 'McMillan (2017) MNRAS 465, 76'
        }
    }
    
    return galactic_data

# ==============================================================================
# HIGH-PRECISION CALCULATION ENGINES
# ==============================================================================

def calculate_v21_quark_confinement(v21_params):
    """Calculate quark confinement in proton using exact v21 parameters"""
    
    try:
        # Extract v21 parameters
        params = v21_params['subatomic_proton']
        
        hbar = params['hbar']                    # J·s
        s = params['quark_spin']                 # spin quantum number
        m = params['effective_quark_mass']       # kg
        r = params['effective_proton_radius']    # m
        sigma = params['string_tension']         # N
        
        # Calculate v21 formula: F_total = (hbar^2 * s^2)/(m * r^3) + sigma
        # Non-relativistic (gamma = 1) as in v21
        gamma = 1.0
        
        F_geometric = (hbar**2 * s**2) / (gamma * m * r**3)
        F_confinement = sigma
        F_total = F_geometric + F_confinement
        
        # Additional useful quantities
        geometric_fraction = F_geometric / F_total
        confinement_fraction = F_confinement / F_total
        
        # Compare to typical strong force magnitude (for context)
        typical_strong_force = 8.2e5  # N (from v21 expected result)
        agreement_with_expected = F_total / typical_strong_force
        
        result = {
            'parameters': {
                'hbar_Js': hbar,
                'spin_s': s,
                'mass_kg': m,
                'radius_m': r,
                'radius_fm': r * 1e15,
                'string_tension_N': sigma,
                'gamma': gamma
            },
            'forces': {
                'F_geometric_N': F_geometric,
                'F_confinement_N': F_confinement,
                'F_total_N': F_total,
                'geometric_fraction': geometric_fraction,
                'confinement_fraction': confinement_fraction
            },
            'comparison': {
                'expected_total_N': typical_strong_force,
                'agreement_ratio': agreement_with_expected,
                'agreement_percent': agreement_with_expected * 100
            },
            'source': params['source']
        }
        
        return result
        
    except Exception as e:
        print(f"ERROR in v21 quark confinement calculation: {e}")
        return None

def calculate_atomic_forces(Z, constants):
    """Calculate atomic forces with high precision"""
    
    if not constants:
        return None
    
    getcontext().prec = 100
    
    try:
        # Convert to high-precision Decimal
        hbar = Decimal(str(constants['hbar']['value']))
        me = Decimal(str(constants['me']['value']))
        e = Decimal(str(constants['e']['value']))
        k = Decimal(str(constants['k']['value']))
        a0 = Decimal(str(constants['a0']['value']))
        alpha = Decimal(str(constants['alpha']['value']))
        
        # Slater's rules for effective charge (simplified)
        Z_dec = Decimal(str(Z))
        if Z == 1:
            Z_eff = Decimal('1.0')
        elif Z == 2:
            Z_eff = Z_dec - Decimal('0.3125')
        else:
            screening = Decimal('0.31') + Decimal('0.002') * (Z_dec - Decimal('2')) / Decimal('98')
            Z_eff = Z_dec - screening
        
        # Orbital radius
        r = a0 / Z_eff
        
        # Relativistic correction
        v_over_c = Z_dec * alpha
        if v_over_c < Decimal('0.1'):
            v2 = v_over_c * v_over_c
            gamma = Decimal('1') + v2 / Decimal('2') + Decimal('3') * v2 * v2 / Decimal('8')
        else:
            gamma = Decimal('1') / (Decimal('1') - v_over_c * v_over_c).sqrt()
        
        # QED corrections for heavy elements
        if Z > 70:
            alpha_sq = alpha * alpha
            z_ratio = Z_dec / Decimal('137')
            qed_correction = Decimal('1') + alpha_sq * z_ratio * z_ratio / Decimal('8')
            gamma = gamma * qed_correction
        
        # Calculate forces
        F_geometric = hbar * hbar / (gamma * me * r * r * r)
        F_coulomb = k * Z_eff * e * e / (gamma * r * r)
        
        # Results
        ratio = F_geometric / F_coulomb
        deviation_ppb = abs(Decimal('1') - ratio) * Decimal('1e9')
        
        return {
            'Z': Z,
            'Z_eff': float(Z_eff),
            'r_pm': float(r * Decimal('1e12')),
            'gamma': float(gamma),
            'F_geometric': float(F_geometric),
            'F_coulomb': float(F_coulomb),
            'ratio': float(ratio),
            'deviation_ppb': float(deviation_ppb),
            'uncertainties': {
                'me_contrib': float(Decimal(constants['me']['relative_uncertainty']) * Decimal('1e9')),
                'k_contrib':  float(Decimal(constants['k']['relative_uncertainty']) * Decimal('1e9')),
                'a0_contrib': float(Decimal(constants['a0']['relative_uncertainty']) * Decimal('1e9'))
            }
        }
        
    except Exception as e:
        print(f"ERROR in atomic calculation for Z={Z}: {e}")
        return None

def calculate_nuclear_forces(r_fm_array, nuclear_constants):
    """Calculate nuclear forces with geometric + confinement model"""
    
    if not nuclear_constants:
        return None
    
    try:
        # Get nuclear parameters
        m_up = nuclear_constants['up_quark_mass']['value']
        m_down = nuclear_constants['down_quark_mass']['value']
        m_q_avg = (m_up + m_down) / 2
        sigma_qcd = nuclear_constants['qcd_string_tension']['value']
        
        if SCIPY_AVAILABLE:
            hbar = const.hbar
        else:
            hbar = 1.054571817e-34
        
        results = []
        
        for r_fm in r_fm_array:
            r = r_fm * 1e-15  # Convert fm to meters
            
            # Geometric force using average light quark mass
            F_geometric = hbar**2 / (m_q_avg * r**3)
            
            # QCD confinement force 
            F_confinement = sigma_qcd * r
            
            # Total force
            F_total = F_geometric + F_confinement
            
            # Determine dominant term
            geometric_fraction = F_geometric / F_total
            confinement_fraction = F_confinement / F_total
            
            results.append({
                'r_fm': r_fm,
                'F_geometric': F_geometric,
                'F_confinement': F_confinement, 
                'F_total': F_total,
                'geometric_fraction': geometric_fraction,
                'confinement_fraction': confinement_fraction,
                'dominant': 'geometric' if F_geometric > F_confinement else 'confinement'
            })
        
        # Calculate crossover point
        crossover_r = (hbar**2 / (m_q_avg * sigma_qcd))**(1/4)
        
        return {
            'results': results,
            'crossover_fm': crossover_r * 1e15,
            'parameters': {
                'm_q_avg_MeV': m_q_avg * (const.c**2 if SCIPY_AVAILABLE else 9e16) / (const.e if SCIPY_AVAILABLE else 1.6e-19) / 1e6,
                'sigma_GeV2_fm': sigma_qcd * 1e-15 / (const.e if SCIPY_AVAILABLE else 1.6e-19) / 1e9
            }
        }
        
    except Exception as e:
        print(f"ERROR in nuclear calculation: {e}")
        return None

def calculate_planetary_forces(planetary_data, constants):
    """Calculate planetary forces using geometric formula"""
    
    if not planetary_data or not constants:
        return None
    
    try:
        hbar = constants['hbar']['value']
        G = constants['G']['value']
        M_sun = planetary_data['_solar_constants']['M_sun']
        c = constants['c']['value']
        
        results = {}
        
        for planet_name, data in planetary_data.items():
            if planet_name.startswith('_'):
                continue
                
            m = data['mass']
            a = data['semimajor_axis']
            v = data['mean_orbital_velocity']
            e = data['eccentricity']
            
            # Calculate quantum number s = mvr/ℏ
            s = m * v * a / hbar
            
            # Relativistic factor (usually negligible for planets)
            gamma = 1 / np.sqrt(1 - (v/c)**2)
            
            # Geometric force: F = s²ℏ²/(γmr³)
            F_geometric = s**2 * hbar**2 / (gamma * m * a**3)
            
            # Classical centripetal force: F = mv²/r
            F_classical = m * v**2 / a
            
            # Gravitational force at semi-major axis
            F_gravity = G * M_sun * m / a**2
            
            # Calculate ratios
            geometric_classical_ratio = F_geometric / F_classical
            centripetal_gravity_ratio = F_classical / F_gravity
            
            results[planet_name] = {
                'mass_kg': m,
                'semimajor_axis_AU': a / planetary_data['_solar_constants']['AU'],
                'velocity_kms': v / 1000,
                'eccentricity': e,
                'quantum_number_s': s,
                'gamma_factor': gamma,
                'F_geometric': F_geometric,
                'F_classical': F_classical,
                'F_gravity': F_gravity,
                'geometric_classical_ratio': geometric_classical_ratio,
                'geometric_classical_deviation': abs(geometric_classical_ratio - 1),
                'centripetal_gravity_ratio': centripetal_gravity_ratio,
                'centripetal_gravity_deviation': abs(centripetal_gravity_ratio - 1),
                'mass_uncertainty': data.get('mass_uncertainty', 0) / m if m != 0 else 0
            }
        
        return results
        
    except Exception as e:
        print(f"ERROR in planetary calculation: {e}")
        return None

def calculate_galactic_forces(galactic_data, constants):
    """Calculate galactic rotation using geometric principle"""
    
    if not galactic_data or not constants:
        return None
    
    try:
        G = constants['G']['value']
        hbar = constants['hbar']['value']
        c = constants['c']['value']
        
        # Solar constants
        M_sun = 1.9884754153381438e30  # kg (IAU 2015)
        pc_to_m = 3.0857e16  # parsec to meters
        
        # Mass model parameters  
        M_bulge = galactic_data['mass_model']['bulge_mass'] * M_sun
        M_disk = galactic_data['mass_model']['disk_mass'] * M_sun
        R_disk = 3.0 * 1000 * pc_to_m  # 3 kpc scale length
        
        results = []
        
        for point in galactic_data['rotation_curve']:
            R_kpc = point['R_kpc']
            V_obs_kms = point['V_kms']
            V_error = point['V_error']
            
            # Convert units
            R = R_kpc * 1000 * pc_to_m  # kpc to meters
            V_obs = V_obs_kms * 1000    # km/s to m/s
            
            # Simple mass model (bulge + exponential disk)
            M_enclosed = M_bulge + M_disk * (1 - np.exp(-R/R_disk))
            
            # Predicted velocity from Newtonian gravity
            V_newton = np.sqrt(G * M_enclosed / R)
            
            # Test geometric principle with typical stellar mass
            m_star = M_sun  # Use solar mass as test mass
            s = m_star * V_obs * R / hbar  # Quantum number
            
            # Geometric prediction: F = s²ℏ²/(mr³) should equal mv²/r
            F_geometric = s**2 * hbar**2 / (m_star * R**3)
            F_centripetal = m_star * V_obs**2 / R
            geometric_ratio = F_geometric / F_centripetal
            
            # Compare observed vs predicted velocities
            newton_obs_ratio = V_newton / V_obs
            velocity_discrepancy = V_obs - V_newton
            
            results.append({
                'R_kpc': R_kpc,
                'V_obs_kms': V_obs_kms,
                'V_error_kms': V_error,
                'V_newton_kms': V_newton / 1000,
                'M_enclosed_Msun': M_enclosed / M_sun,
                'quantum_number_s': s,
                'geometric_centripetal_ratio': geometric_ratio,
                'newton_obs_ratio': newton_obs_ratio,
                'velocity_discrepancy_kms': velocity_discrepancy / 1000,
                'dark_matter_factor': (V_obs / V_newton)**2  # (V_obs/V_pred)² gives mass ratio
            })
        
        return results
        
    except Exception as e:
        print(f"ERROR in galactic calculation: {e}")
        return None

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

def main():
    """Main multiscale verification including v21 quark confinement test"""
    
    print("MULTISCALE HIGH-PRECISION VERIFICATION v21")
    print("="*70)
    print("Includes quark confinement test from paper version 21")
    print("External data sources only. No conclusions drawn - data speaks for itself.")
    print("Repository: https://git.esus.name/esus/spin_paper/")
    print()
    
    # Load all external data
    print("LOADING EXTERNAL DATA SOURCES")
    print("-" * 40)
    
    v21_params = load_v21_parameters()
    codata_constants = load_codata_constants()
    nuclear_constants = load_pdg_constants()
    planetary_data = load_nasa_planetary_data()
    galactic_data = load_galactic_data()
    
    print(f"✓ v21 parameters: {len(v21_params)} parameter sets")
    
    if codata_constants:
        print(f"✓ CODATA constants: {len(codata_constants)} values")
    else:
        print("❌ CODATA constants: Failed to fetch")
    
    print(f"✓ Nuclear constants: {len(nuclear_constants)} values (PDG 2022)")
    print(f"✓ Planetary data: {len(planetary_data)-1} planets (NASA JPL)")
    print(f"✓ Galactic data: {len(galactic_data['rotation_curve'])} points (Gaia+literature)")
    
    # 0. V21 QUARK CONFINEMENT TEST
    print(f"\n" + "="*70)
    print("0. SUBATOMIC SCALE (V21): F = ℏ²s²/(γmr³) + σ")
    print("="*70)
    print("Source: Paper v21 examples_v21.tex")
    
    v21_result = calculate_v21_quark_confinement(v21_params)
    
    if v21_result:
        print("Parameters from paper v21:")
        p = v21_result['parameters']
        print(f"  ℏ = {p['hbar_Js']:.3e} J⋅s")
        print(f"  s = {p['spin_s']} (quark spin)")
        print(f"  m = {p['mass_kg']:.3e} kg ({p['mass_kg'] / (1.783e-30):.1f} MeV/c²)")
        print(f"  r = {p['radius_fm']:.1f} fm")
        print(f"  σ = {p['string_tension_N']:.1e} N")
        print(f"  γ = {p['gamma']} (non-relativistic)")
        
        print("\nForce components:")
        f = v21_result['forces']
        print(f"  F_geometric = ℏ²s²/(γmr³) = {f['F_geometric_N']:.1e} N")
        print(f"  F_confinement = σ = {f['F_confinement_N']:.1e} N") 
        print(f"  F_total = {f['F_total_N']:.1e} N")
        print(f"  Geometric fraction: {f['geometric_fraction']*100:.1f}%")
        print(f"  Confinement fraction: {f['confinement_fraction']*100:.1f}%")
        
        print("\nComparison to expected strong force:")
        c = v21_result['comparison']
        print(f"  Expected total: {c['expected_total_N']:.1e} N")
        print(f"  Calculated total: {f['F_total_N']:.1e} N")
        print(f"  Agreement: {c['agreement_percent']:.1f}%")
    
    # 1. ATOMIC SCALE VERIFICATION
    print(f"\n" + "="*70)
    print("1. ATOMIC SCALE: F = ℏ²/(γmr³) = ke²/r²")
    print("="*70)
    print("Source: CODATA 2022 via scipy.constants")
    
    if codata_constants:
        test_elements = [1, 6, 18, 26, 47, 79, 92]
        
        print(f"{'Z':<3} {'Element':<8} {'Z_eff':<6} {'r(pm)':<8} {'γ':<8} {'Ratio':<18} {'Dev(ppb)':<10}")
        print("-" * 70)
        
        atomic_deviations = []
        for Z in test_elements:
            result = calculate_atomic_forces(Z, codata_constants)
            if result:
                element_names = {1:'H', 6:'C', 18:'Ar', 26:'Fe', 47:'Ag', 79:'Au', 92:'U'}
                elem = element_names.get(Z, f'Z{Z}')
                
                print(f"{Z:<3} {elem:<8} {result['Z_eff']:<6.3f} {result['r_pm']:<8.2f} "
                      f"{result['gamma']:<8.5f} {result['ratio']:<18.15f} {result['deviation_ppb']:<10.6f}")
                
                atomic_deviations.append(result['deviation_ppb'])
        
        if atomic_deviations:
            mean_dev = np.mean(atomic_deviations)
            std_dev = np.std(atomic_deviations)
            print(f"\nStatistics: Mean = {mean_dev:.6f} ppb, Std = {std_dev:.10f} ppb")
            
            # Show uncertainty contributions
            print(f"\nUncertainty sources (ppb):")
            if test_elements and calculate_atomic_forces(1, codata_constants):
                unc = calculate_atomic_forces(1, codata_constants)['uncertainties']
                print(f"  Electron mass: {unc['me_contrib']:.3f}")
                print(f"  Coulomb constant: {unc['k_contrib']:.3f}")
                print(f"  Bohr radius: {unc['a0_contrib']:.3f}")
    
    # 2. NUCLEAR SCALE VERIFICATION  
    print(f"\n" + "="*70)
    print("2. NUCLEAR SCALE: F = ℏ²/(γm_q r³) + σr")
    print("="*70)
    print("Source: PDG 2022 + Lattice QCD")
    
    r_range = np.logspace(-1, 0.5, 12)  # 0.1 to ~3 fm
    nuclear_result = calculate_nuclear_forces(r_range, nuclear_constants)
    
    if nuclear_result:
        print(f"Quark mass: {nuclear_result['parameters']['m_q_avg_MeV']:.2f} MeV/c²")
        print(f"String tension: {nuclear_result['parameters']['sigma_GeV2_fm']:.2f} GeV²/fm")
        print(f"Crossover: {nuclear_result['crossover_fm']:.3f} fm")
        print()
        
        print(f"{'r(fm)':<8} {'F_geom(N)':<12} {'F_conf(N)':<12} {'Geom%':<8} {'Dominant':<10}")
        print("-" * 55)
        
        for res in nuclear_result['results']:
            print(f"{res['r_fm']:<8.3f} {res['F_geometric']:<12.3e} {res['F_confinement']:<12.3e} "
                  f"{res['geometric_fraction']*100:<8.1f} {res['dominant']:<10}")
    
    # 3. PLANETARY SCALE VERIFICATION
    print(f"\n" + "="*70)
    print("3. PLANETARY SCALE: F = s²ℏ²/(γmr³) = mv²/r")  
    print("="*70)
    print("Source: NASA JPL + IAU standards")
    
    planetary_result = calculate_planetary_forces(planetary_data, codata_constants)
    
    if planetary_result:
        print(f"{'Planet':<8} {'s-factor':<12} {'Geom/Class':<12} {'Cent/Grav':<12} {'Ecc':<8}")
        print("-" * 55)
        
        for planet, res in planetary_result.items():
            print(f"{planet:<8} {res['quantum_number_s']:<12.2e} "
                  f"{res['geometric_classical_ratio']:<12.10f} "
                  f"{res['centripetal_gravity_ratio']:<12.8f} {res['eccentricity']:<8.5f}")
    
    # 4. GALACTIC SCALE VERIFICATION
    print(f"\n" + "="*70)
    print("4. GALACTIC SCALE: Same formula, expected breakdown")
    print("="*70)
    print("Source: Gaia DR3 + Reid 2019 + McMillan 2017")
    
    galactic_result = calculate_galactic_forces(galactic_data, codata_constants)
    
    if galactic_result:
        print(f"{'R(kpc)':<7} {'V_obs':<8} {'V_pred':<8} {'V_ratio':<8} {'DM_factor':<10} {'s-factor':<12}")
        print("-" * 60)
        
        for res in galactic_result:
            print(f"{res['R_kpc']:<7.1f} {res['V_obs_kms']:<8.0f} "
                  f"{res['V_newton_kms']:<8.0f} {res['newton_obs_ratio']:<8.3f} "
                  f"{res['dark_matter_factor']:<10.2f} {res['quantum_number_s']:<12.2e}")
    
    # 5. DATA SOURCES SUMMARY
    print(f"\n" + "="*70)
    print("DATA SOURCES")
    print("="*70)
    
    print("Subatomic v21: Paper v21 examples_v21.tex")
    print("Atomic: CODATA 2022 via scipy.constants")
    print("Nuclear: PDG 2022 + FLAG 2021 Lattice QCD")  
    print("Planetary: NASA JPL Fact Sheet 2021 + IAU 2015")
    print("Galactic: Gaia DR3 + Reid (2019) + McMillan (2017)")
    
    print(f"\nNOTE: Results presented without interpretation.")
    print(f"Observer should evaluate where formulas work/fail.")
    print(f"v21 quark confinement test uses exact paper parameters.")

if __name__ == "__main__":
    if len(sys.argv) > 1 and sys.argv[1] in ['-h', '--help']:
        print("Usage: python multiscale_high_precision.py")
        print("  Multiscale verification using external data only")
        print("  Tests geometric force principle across all scales")
        print("  Includes v21 quark confinement test")
        sys.exit(0)
    
    main()