import numpy as np

# Constants
phi = (1 + np.sqrt(5)) / 2

PRIMES = [
    2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
    31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
    73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
    127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
    179, 181, 191, 193, 197, 199, 211, 223, 227, 229
]

def fib_real(n):
    """Real-valued Fibonacci extension."""
    from math import cos, pi, sqrt
    phi_inv = 1 / phi
    return (phi**n / sqrt(5)) - (phi_inv**n * cos(pi * n))

def Omega(n, alpha=0.1):
    """Example symbolic Omega function (field tension decay)."""
    return np.exp(-alpha * n)

def D_n(n, r=1.0, k=1.0, alpha=0.1):
    """Symbolic speed of light at recursion depth n."""
    Fn = fib_real(n)
    idx = int(np.floor(n)) % len(PRIMES)
    Pn = PRIMES[idx]
    omega_val = Omega(n, alpha)
    val = phi * Fn * (2 ** n) * Pn * omega_val
    val = max(val, 1e-15)
    return np.sqrt(val) * (r ** k)

def find_reference_n0(c_phys_target=299792.458, n_min=0, n_max=20, steps=100, alpha=0.1):
    """Find recursion depth n0 where symbolic speed D_n best matches physical speed c_phys_target."""
    n_values = np.linspace(n_min, n_max, steps)
    errors = []
    for n in n_values:
        symbolic_c = D_n(n, alpha=alpha)
        # Calculate scale factor S to map symbolic_c to physical c_phys_target
        S = c_phys_target / symbolic_c
        # Check error of mapping (just 0 because we define S exactly)
        error = abs(symbolic_c * S - c_phys_target)
        errors.append(error)
    min_idx = np.argmin(errors)
    best_n0 = n_values[min_idx]
    best_S = c_phys_target / D_n(best_n0, alpha=alpha)
    return best_n0, best_S

# Find reference n0 and scale S
n0, S = find_reference_n0()

print(f"Reference recursion depth n0: {n0:.4f}")
print(f"Emergent physical scale factor S: {S:.4f} km/s per symbolic unit")

def physical_speed_of_light(n, n0=n0, S=S, alpha=0.1):
    """Physical speed of light in km/s at recursion depth n."""
    c_ref = D_n(n0, alpha=alpha)
    c_n = D_n(n, alpha=alpha)
    # Compute relative symbolic speed and scale
    return (c_n / c_ref) * S

# Example: print speeds across recursion depths
print("\nSpeed of light at various recursion depths (km/s):")
for n_test in np.linspace(0, 20, 21):
    c_val = physical_speed_of_light(n_test)
    print(f"  n = {n_test:5.2f} : c = {c_val:10.2f} km/s")
