import numpy as np
import scipy.constants as const
from scipy.optimize import fsolve

# Physical constants in SI units
c = const.c  # speed of light: 299792458 m/s
e = const.e  # elementary charge: 1.602176634e-19 C
k = 1 / (4 * np.pi * const.epsilon_0)  # Coulomb constant: 8.9875517923e9 N⋅m²/C²
hbar = const.hbar  # reduced Planck constant: 1.054571817e-34 J⋅s
m_e = const.m_e  # electron mass: 9.1093837015e-31 kg
m_p = const.m_p  # proton mass: 1.67262192369e-27 kg
alpha = const.alpha  # fine structure constant: ~1/137

print("=== DIMENSIONAL ANALYSIS ===")
print(f"ℏ = {hbar:.6e} J⋅s = kg⋅m²/s")
print(f"ℏ² = {hbar**2:.6e} J²⋅s² = kg²⋅m⁴/s²")
print(f"Units of ℏ²: [M²L⁴T⁻²]")
print(f"\nInteresting: L⁴ = (L²)² - suggesting 2D area squared!")
print(f"So ℏ² relates mass², area², and time² fundamentally")

print("\n=== OUR EQUATION: c² = ke²γEr/(ℏ²) ===")
print("Rearranged forms:")
print("γ = c²ℏ²/(ke²Er)")
print("E = c²ℏ²/(ke²γr)")
print("r = c²ℏ²/(ke²γE)")

# Case 1: Hydrogen-antihydrogen annihilation
print("\n=== CASE 1: H + anti-H annihilation ===")
E_annihilation = 2 * (m_e + m_p) * c**2  # Total rest mass energy
r_bohr = const.physical_constants['Bohr radius'][0]  # 5.29e-11 m

# Solve for gamma at Bohr radius
gamma_annihilation = (c**2 * hbar**2) / (k * e**2 * E_annihilation * r_bohr)
v_from_gamma = c * np.sqrt(1 - 1/gamma_annihilation**2) if gamma_annihilation > 1 else 0

print(f"Total annihilation energy: {E_annihilation:.6e} J = {E_annihilation/const.eV:.3f} eV")
print(f"At Bohr radius r = {r_bohr:.3e} m:")
print(f"  γ = {gamma_annihilation:.6f}")
print(f"  Implied velocity: {v_from_gamma/c:.6f}c")

# Case 2: Ground state hydrogen (binding energy)
print("\n=== CASE 2: Ground state hydrogen ===")
E_binding = 13.6 * const.eV  # Hydrogen ionization energy
gamma_ground = (c**2 * hbar**2) / (k * e**2 * E_binding * r_bohr)

print(f"Binding energy: {E_binding:.6e} J = 13.6 eV")
print(f"γ = {gamma_ground:.3e}")
print(f"This huge γ suggests we're probing quantum regime!")

# Case 3: What if we set γ = 1 (non-relativistic)?
print("\n=== CASE 3: Non-relativistic case (γ = 1) ===")
gamma_nr = 1
E_nr = (c**2 * hbar**2) / (k * e**2 * gamma_nr * r_bohr)
print(f"Energy when γ = 1: {E_nr:.6e} J = {E_nr/const.eV:.3f} eV")
print(f"This is {E_nr/E_binding:.1f}× the actual binding energy")

# Explore the meaning of the equation parts
print("\n=== EQUATION COMPONENTS ===")
quantum_term = hbar**2
em_term = k * e**2
print(f"Quantum term (ℏ²): {quantum_term:.6e} J²⋅s²")
print(f"EM term (ke²): {em_term:.6e} J⋅m")
print(f"Ratio c²ℏ²/(ke²) = {(c**2 * quantum_term)/em_term:.6e} J⋅m")
print(f"This has units of [Energy × Length]")

# The fine structure connection
print("\n=== FINE STRUCTURE CONSTANT CONNECTION ===")
# α = ke²/(ℏc) in SI units
calculated_alpha = (k * e**2) / (hbar * c)
print(f"α calculated: {calculated_alpha:.10f}")
print(f"α known: {alpha:.10f}")
print(f"Our equation can be rewritten using α:")
print(f"γ = c/(αEr/ℏ)")

# Exploring unit independence
print("\n=== UNIT INDEPENDENCE CHECK ===")
# In natural units where c = ℏ = 1
print("In natural units (c = ℏ = 1):")
print("γ = 1/(αEr)")
print("This shows γ depends only on:")
print("  - The fine structure constant α (dimensionless)")
print("  - Energy × distance (action units)")

# What this means for spacetime
print("\n=== SPACETIME INTERPRETATION ===")
print("Since γ = dt/dτ (time dilation factor):")
print("Our equation suggests coordinate time / proper time = c²ℏ²/(ke²Er)")
print("\nThis means time dilation emerges from the balance between:")
print("  - Quantum action (ℏ²) pushing toward uncertainty")
print("  - EM binding (ke²Er) pulling toward classical behavior")
print("\nAt high energies or small distances, γ becomes huge,")
print("suggesting extreme time dilation in the quantum regime!")

# Function to explore different scenarios
def explore_scenario(E, r, label):
    gamma = (c**2 * hbar**2) / (k * e**2 * E * r)
    print(f"\n{label}:")
    print(f"  E = {E/const.eV:.3e} eV, r = {r:.3e} m")
    print(f"  γ = {gamma:.3e}")
    if gamma > 1 and gamma < 1e6:
        v = c * np.sqrt(1 - 1/gamma**2)
        print(f"  v/c = {v/c:.6f}")
    
# Explore various energy/distance scales
print("\n=== EXPLORING DIFFERENT SCALES ===")
explore_scenario(1*const.eV, 1e-9, "1 eV at 1 nm")
explore_scenario(1*const.eV*1000, 1e-12, "1 keV at 1 pm") 
explore_scenario(511*const.eV*1000, const.physical_constants['classical electron radius'][0], 
                 "Electron rest mass at classical radius")
explore_scenario(const.m_e * c**2, const.physical_constants['Compton wavelength'][0]/(2*np.pi),
                 "Electron at reduced Compton wavelength")