#!/usr/bin/env python3
"""
planck_chronotopic_compare.py

Compare standard Planck spectral radiance and chronotopic form.
Also solve for Wien peak constant x_star where:
    5*(1 - exp(-x)) = x  (x = h c / (lambda_max k_B T))

Requires: numpy, scipy
Run: python3 planck_chronotopic_compare.py
"""

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

# constants
h = const.h
hbar = const.hbar
c = const.c
kB = const.k
pi = np.pi

def B_standard_lambda(lmbda, T):
    """Standard Planck spectral radiance per unit wavelength (SI)"""
    # B_lambda = 2 h c^2 / lambda^5 * 1/(exp(hc/(lambda kB T)) - 1)
    x = h * c / (lmbda * kB * T)
    return (2.0 * h * c**2) / (lmbda**5) / (np.exp(x) - 1.0)

def B_chrono_lambda(lmbda, Theta, S_star=hbar, v_sync=c):
    """Chronotopic Planck-type law in wavelength form"""
    y = 2.0 * pi * S_star * v_sync / (lmbda * Theta)
    return (2.0 * S_star * v_sync**2) / (lmbda**5) / (np.exp(y) - 1.0)

# Example evaluation
T = 5778.0           # Sun-like surface temperature (K)
lmbda = 500e-9       # 500 nm
Theta = kB * T

B_std = B_standard_lambda(lmbda, T)
B_chr = B_chrono_lambda(lmbda, Theta)

# Ratio should be ~1 under the mapping S_*->hbar, v_sync->c, Theta->k_B T
ratio = B_chr / B_std

# Solve for Wien x_star root of 5(1-e^{-x}) = x
def wien_func(x):
    return 5.0 * (1.0 - np.exp(-x)) - x

x_star = brentq(wien_func, 0.1, 20.0)
b_constant = (h * c) / (kB * x_star)  # Wien displacement: lambda_max * T = b

# Print results
print("Planck vs Chronotopic comparison")
print("--------------------------------")
print(f"lambda = {lmbda:.3e} m, T = {T} K")
print(f"B_standard(lambda,T) = {B_std:.8e} W m^-2 sr^-1 m^-1")
print(f"B_chrono(lambda,Theta) = {B_chr:.8e} W m^-2 sr^-1 m^-1")
print(f"ratio (chrono / standard) = {ratio:.8e}")
print()
print("Wien displacement root and constant")
print("-----------------------------------")
print(f"x_star (root of 5(1-e^-x)=x) = {x_star:.8f}")
print(f"Wien b constant (lambda_max * T) = {b_constant:.8e} m*K")

# Expected sanity:
# ratio ≈ 1.0 (to machine precision) if S_*=hbar, v_sync=c, Theta=kB*T
# x_star ≈ 4.965114231744276  (typical textbook value ~ 4.965114)
# b ≈ 2.897771955e-3 m*K