diff --git a/scripts/elastic_leash_test.py b/scripts/elastic_leash_test.py new file mode 100644 index 0000000..57af543 --- /dev/null +++ b/scripts/elastic_leash_test.py @@ -0,0 +1,417 @@ +#!/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() \ No newline at end of file diff --git a/scripts/multiscale_high_precision.py b/scripts/multiscale_high_precision.py new file mode 100644 index 0000000..a795177 --- /dev/null +++ b/scripts/multiscale_high_precision.py @@ -0,0 +1,706 @@ +#!/usr/bin/env python3 +""" +High-precision multiscale verification using externally sourced data + +Tests the geometric force principle F = ℏ²/(γmr³) across scales: +- 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_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_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""" + + print("MULTISCALE HIGH-PRECISION VERIFICATION") + print("="*70) + print("External data sources only. No interpretations - 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) + + codata_constants = load_codata_constants() + nuclear_constants = load_pdg_constants() + planetary_data = load_nasa_planetary_data() + galactic_data = load_galactic_data() + + 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)") + + # 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("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.") + +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") + sys.exit(0) + + main() diff --git a/scripts/multiscale_v21.py b/scripts/multiscale_v21.py new file mode 100644 index 0000000..f416b8c --- /dev/null +++ b/scripts/multiscale_v21.py @@ -0,0 +1,829 @@ +#!/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() \ No newline at end of file