#!/usr/bin/env python3
"""
Voorspelling #15: Anisotropie op Planck-Schaal (Grid-Roostereffect)
====================================================================
Spectrum Grid Model — Fermi-LAT GRB Analyse

Hypothese:
    Als het universum een discreet grid is op Planck-schaal, dan is ruimte
    niet perfect isotroop. Fotonen van verschillende energieën propageren
    met licht verschillende snelheden door het grid:

        v(E) = c × [1 ± η × (E / E_Planck)^n]

    waar:
        η = orde-1 coëfficiënt (afhankelijk van grid-geometrie)
        E_Planck = √(ℏc⁵/G) ≈ 1.22 × 10¹⁹ GeV
        n = 1 (lineaire LIV) of n = 2 (kwadratische LIV)

    Over kosmologische afstanden (miljarden lichtjaren) accumuleert dit
    verschil tot meetbare aankomsttijdverschillen tussen hoge- en
    lage-energie fotonen uit dezelfde gamma-ray burst (GRB).

Methode:
    1. Simuleer GRB foton-emissie met simultaan vertrek
    2. Pas grid-dispersie toe over kosmologische afstand
    3. Analyseer aankomsttijdverschillen als functie van energie
    4. Fit het Lorentz Invariance Violation (LIV) model
    5. Vergelijk met Fermi-LAT publicaties (Abdo et al., 2009)

Databron:
    - Fermi-LAT GRB catalogus (fermi.gsfc.nasa.gov)
    - Specifiek: GRB 090510 (korte GRB, z=0.903, meest beperkend)
    - GRB 080916C (lange GRB, z=4.35, hoogste energie fotonen)

Auteur: Marald Bes / Spectrum van Alles
Datum: 2026-04-08
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import os

# ─── Natuurconstanten ───────────────────────────────────────
E_PLANCK = 1.22e19  # GeV — Planck-energie
c = 3e8             # m/s
H0 = 67.4           # km/s/Mpc — Hubble constante
Mpc_to_m = 3.086e22 # 1 Mpc in meters

# ─── Kosmologische Afstand ──────────────────────────────────
def comoving_distance(z, H0=67.4, Omega_m=0.315, Omega_L=0.685, n_steps=1000):
    """
    Bereken medebewegende afstand als functie van roodverschuiving z.
    Vlak ΛCDM model.

    D_C = (c/H0) × ∫₀ᶻ dz' / √(Ω_m(1+z')³ + Ω_Λ)
    """
    z_arr = np.linspace(0, z, n_steps)
    dz = z_arr[1] - z_arr[0]

    integrand = 1.0 / np.sqrt(Omega_m * (1 + z_arr)**3 + Omega_L)
    integral = np.trapz(integrand, z_arr)

    D_C = (c / (H0 * 1e3 / Mpc_to_m)) * integral  # meters
    return D_C


# ─── LIV Tijdvertraging ─────────────────────────────────────
def liv_time_delay(E_high, E_low, z, eta=1.0, n=1):
    """
    Bereken Lorentz Invariance Violation tijdvertraging.

    Δt = ± (1+n)/(2×H0) × (E_high^n - E_low^n)/E_QG^n × K(z)

    waar K(z) de kosmologische integraal is die rekening houdt
    met expansie van het universum.

    Parameters:
        E_high, E_low: fotonenergie in GeV
        z: roodverschuiving van bron
        eta: sterkte-parameter (|η| ~ 1 voor Planck-schaal effect)
        n: orde van LIV (1 = lineair, 2 = kwadratisch)

    Returns:
        Δt in seconden
    """
    # Vereenvoudigde K(z) — goede benadering voor z < 5
    K_z = comoving_distance(z) * (1 + z) / c

    # Tijdvertraging
    if abs(eta) < 1e-30:
        return 0.0

    E_QG = E_PLANCK / abs(eta)  # effectieve quantum-gravity schaal

    # Vereenvoudigd: Δt ≈ D/(2c) × (E_h^n - E_l^n)/E_QG^n
    delta_t_simple = (comoving_distance(z) / (2 * c)) * \
                     (E_high**n - E_low**n) / E_QG**n

    return delta_t_simple


# ─── Grid-Geometrie Modellen ────────────────────────────────
def grid_dispersion(E, geometry='cubic', grid_spacing=1.0):
    """
    Bereken dispersierelatie voor verschillende grid-geometrieën.

    Het grid-model voorspelt dat de effectieve η-parameter afhangt
    van de geometrie van het rooster:

    - Kubisch: sterkste anisotropie, η ~ 1
    - Hexagonaal: middelmatig, η ~ 0.3
    - FCC (face-centered cubic): zwak, η ~ 0.1
    - Amorf: willekeurig, η ~ 0 (geen preferente richting)

    Parameters:
        E: fotonenergie in GeV
        geometry: 'cubic', 'hexagonal', 'fcc', 'amorphous'
        grid_spacing: grid-afstand in Planck-lengten

    Returns:
        v/c - 1 (fractionele snelheidsafwijking)
    """
    eta_map = {
        'cubic': 1.0,
        'hexagonal': 0.3,
        'fcc': 0.1,
        'amorphous': 0.01
    }

    eta = eta_map.get(geometry, 1.0)

    # Lineaire dispersie (n=1)
    dv_linear = eta * (E / E_PLANCK) * grid_spacing

    # Kwadratische dispersie (n=2)
    dv_quadratic = eta * (E / E_PLANCK)**2 * grid_spacing

    return dv_linear, dv_quadratic


# ─── GRB Simulatie ───────────────────────────────────────────
def simulate_grb(z=0.903, n_photons=500, E_range=(0.1, 31),
                 eta=1.0, n_liv=1, intrinsic_spread=0.01):
    """
    Simuleer een GRB-observatie met Fermi-LAT.

    Parameters:
        z: roodverschuiving
        n_photons: aantal gedetecteerde fotonen
        E_range: energiebereik in GeV
        eta: LIV parameter
        n_liv: orde van LIV
        intrinsic_spread: intrinsieke tijdspreiding van bron (s)

    Returns:
        energies: fotonenergiën (GeV)
        arrival_times: aankomsttijden (s)
        true_delay: ware LIV-vertraging
    """
    # Genereer energieën (powerlaw-verdeling: dN/dE ∝ E^-2)
    E_min, E_max = E_range
    u = np.random.uniform(0, 1, n_photons)
    energies = E_min * (E_max/E_min)**u  # powerlaw sampling

    # Bereken LIV tijdvertraging per foton
    E_ref = E_min  # referentie-energie (laagste)
    true_delay = np.array([liv_time_delay(E, E_ref, z, eta, n_liv) for E in energies])

    # Intrinsieke tijdspreiding (bron is niet perfect simultaan)
    intrinsic = np.random.normal(0, intrinsic_spread, n_photons)

    # Detectorruis
    detector_noise = np.random.normal(0, 0.001, n_photons)  # 1 ms

    # Totale aankomsttijd
    arrival_times = true_delay + intrinsic + detector_noise

    return energies, arrival_times, true_delay


# ─── Analyse: Fit LIV parameter ─────────────────────────────
def fit_liv_parameter(energies, arrival_times, z, n_liv=1):
    """
    Fit de LIV parameter η uit aankomsttijddata.

    Model: Δt(E) = A × E^n + B (constante offset)

    Returns:
        eta_fit: gefitte η-waarde
        eta_err: onzekerheid op η
        chi2: chi-squared van fit
    """
    def model(E, A, B):
        return A * E**n_liv + B

    try:
        popt, pcov = curve_fit(model, energies, arrival_times,
                               p0=[0, 0], maxfev=10000)
        A_fit, B_fit = popt
        A_err = np.sqrt(pcov[0, 0])

        # Converteer A naar η
        D = comoving_distance(z)
        eta_fit = A_fit * 2 * c * E_PLANCK**n_liv / D
        eta_err = A_err * 2 * c * E_PLANCK**n_liv / D

        # Chi-squared
        residuals = arrival_times - model(energies, *popt)
        chi2 = np.sum(residuals**2) / (len(energies) - 2)

        return eta_fit, eta_err, chi2, popt
    except Exception:
        return 0, np.inf, np.inf, [0, 0]


# ─── Hoofdvisualisatie ───────────────────────────────────────
def create_analysis_figure(output_dir):
    """
    4-paneel analyse:
    A) Dispersierelatie voor verschillende grid-geometrieën
    B) Gesimuleerde GRB met LIV-signaal
    C) Gesimuleerde GRB zonder LIV (null-hypothese)
    D) Bovengrens op η als functie van grid-geometrie
    """
    np.random.seed(42)

    fig, axes = plt.subplots(2, 2, figsize=(16, 12))
    fig.suptitle('Voorspelling #15: Anisotropie op Planck-Schaal (Grid-Roostereffect)\n'
                 'Spectrum Sub-quantum Grid Model — Fermi-LAT GRB Analyse',
                 fontsize=14, fontweight='bold', y=0.98)

    # ─── Panel A: Dispersierelaties ───────────────────────────
    ax = axes[0, 0]
    E_range = np.logspace(-1, 5, 500)  # 0.1 GeV tot 100 TeV

    for geom, color, ls in [('cubic', 'red', '-'), ('hexagonal', 'orange', '--'),
                             ('fcc', 'green', '-.'), ('amorphous', 'blue', ':')]:
        dv_lin, dv_quad = grid_dispersion(E_range, geom)
        ax.loglog(E_range, np.abs(dv_lin), color=color, linestyle=ls,
                  linewidth=2, label=f'{geom} (n=1)')

    ax.axhline(1e-20, color='gray', linestyle=':', alpha=0.5)
    ax.text(1e4, 1.5e-20, 'Fermi-LAT bovengrens', fontsize=8, color='gray')

    ax.set_xlabel('Fotonenergie (GeV)')
    ax.set_ylabel('|Δv/c|')
    ax.set_title('A) Fractionele Snelheidsafwijking per Grid-Geometrie')
    ax.legend(fontsize=9)
    ax.set_ylim(1e-30, 1e-10)
    ax.grid(True, alpha=0.3, which='both')

    # ─── Panel B: GRB met LIV-signaal ────────────────────────
    ax = axes[0, 1]

    # Simuleer GRB 090510-achtig (z=0.903) met klein signaal
    E_sig, t_sig, delay_sig = simulate_grb(
        z=0.903, n_photons=300, eta=0.5, n_liv=1,
        intrinsic_spread=0.005
    )

    scatter = ax.scatter(E_sig, t_sig * 1000, c=E_sig, cmap='hot',
                         s=15, alpha=0.7, edgecolors='none')
    plt.colorbar(scatter, ax=ax, label='Energie (GeV)')

    # Fit lijn
    E_fit = np.linspace(0.1, 31, 100)
    eta_fit, eta_err, chi2, popt = fit_liv_parameter(E_sig, t_sig, z=0.903, n_liv=1)
    t_fit = popt[0] * E_fit + popt[1]
    ax.plot(E_fit, t_fit * 1000, 'r-', linewidth=2, alpha=0.8,
            label=f'Fit: η = {eta_fit:.2e} ± {eta_err:.2e}')

    ax.set_xlabel('Fotonenergie (GeV)')
    ax.set_ylabel('Aankomsttijd (ms)')
    ax.set_title(f'B) GRB z=0.903 MET Grid-Signaal (η=0.5)')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    # ─── Panel C: GRB zonder LIV (null) ──────────────────────
    ax = axes[1, 0]

    E_null, t_null, _ = simulate_grb(
        z=0.903, n_photons=300, eta=0.0, n_liv=1,
        intrinsic_spread=0.005
    )

    scatter2 = ax.scatter(E_null, t_null * 1000, c=E_null, cmap='cool',
                          s=15, alpha=0.7, edgecolors='none')
    plt.colorbar(scatter2, ax=ax, label='Energie (GeV)')

    eta_null, eta_null_err, chi2_null, popt_null = fit_liv_parameter(
        E_null, t_null, z=0.903, n_liv=1
    )
    t_fit_null = popt_null[0] * E_fit + popt_null[1]
    ax.plot(E_fit, t_fit_null * 1000, 'b-', linewidth=2, alpha=0.8,
            label=f'Fit: η = {eta_null:.2e} (consistent met 0)')

    ax.set_xlabel('Fotonenergie (GeV)')
    ax.set_ylabel('Aankomsttijd (ms)')
    ax.set_title('C) GRB z=0.903 ZONDER Grid-Signaal (η=0)')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.3)

    # ─── Panel D: Constraint plot ─────────────────────────────
    ax = axes[1, 1]

    geometries = ['cubic', 'hexagonal', 'fcc', 'amorphous']
    eta_values = [1.0, 0.3, 0.1, 0.01]
    colors = ['red', 'orange', 'green', 'blue']

    # Fermi-LAT bovengrens (Abdo et al. 2009): E_QG > 1.2 × E_Planck voor n=1
    fermi_limit = 1.2  # η < 1/1.2 ≈ 0.83

    bars = ax.barh(geometries, eta_values, color=colors, alpha=0.7, height=0.5)
    ax.axvline(fermi_limit, color='black', linewidth=2, linestyle='--',
               label=f'Fermi-LAT limiet: η < {1/fermi_limit:.2f}')

    for i, (geom, eta) in enumerate(zip(geometries, eta_values)):
        status = '✗ UITGESLOTEN' if eta > 1/fermi_limit else '✓ Consistent'
        ax.text(max(eta + 0.05, 0.15), i, status, va='center', fontsize=10,
                fontweight='bold', color='red' if eta > 1/fermi_limit else 'green')

    ax.set_xlabel('Verwachte η-waarde')
    ax.set_title('D) Grid-Geometrieën vs Fermi-LAT Constraint')
    ax.legend(fontsize=9, loc='lower right')
    ax.set_xlim(0, 1.5)
    ax.grid(True, alpha=0.3, axis='x')

    plt.tight_layout(rect=[0, 0, 1, 0.95])

    filepath = os.path.join(output_dir, 'voorspelling_15_fermi_lat.png')
    plt.savefig(filepath, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f"Saved: {filepath}")

    return filepath


# ─── Main ────────────────────────────────────────────────────
if __name__ == "__main__":
    output_dir = os.path.dirname(os.path.abspath(__file__))

    print("Voorspelling #15: Anisotropie op Planck-Schaal")
    print("=" * 55)

    filepath = create_analysis_figure(output_dir)

    # Samenvatting
    print("\nResultaat-samenvatting:")
    print("-" * 55)
    print("Fermi-LAT GRB 090510 constraint: E_QG > 1.2 × E_Planck")
    print("Dit sluit uit: kubisch grid (η=1.0)")
    print("Nog consistent: hexagonaal (η=0.3), FCC (η=0.1), amorf (η=0.01)")
    print()
    print("Implicatie voor Spectrum Grid Model:")
    print("→ Het grid kan NIET kubisch zijn")
    print("→ Hexagonaal of FCC grid overleeft de huidige data")
    print("→ Kwadratische LIV (n=2) is nauwelijks beperkt — alle geometrieën open")
    print()
    print("Volgende stap: analyseer GRB 190114C (MAGIC telescoop, z=0.42)")
    print("voor TeV-fotonen die sterkere constraints geven")
    print(f"\n✓ Analyse voltooid. Figuur: {filepath}")