Noether’s Theorem Reinterpreted in QFunity

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace, convolve
import matplotlib.animation as animation

# Paramètres de la simulation
NX, NY = 64, 64  # Taille de la grille
DX = 1.0          # Pas d'espace
DT = 0.01         # Pas de temps
STEPS = 500       # Nombre d'itérations
EPS = 0.01        # Paramètre de régularisation (principe "zéro n'existe pas")

# Initialisation du champ Psi (complexe) avec des fluctuations aléatoires
np.random.seed(42)
Psi = 0.1 * (np.random.randn(NX, NY) + 1j * np.random.randn(NX, NY))
Psi_history = []  # Pour stocker l'évolution

# Opérateur de vibration V_hat (similaire à un laplacien avec dispersion)
def V_hat(Psi):
    """Applique l'opérateur de vibration (dispersion/expansion)"""
    # Laplacien pour la dispersion
    lap = laplace(np.real(Psi)) + 1j * laplace(np.imag(Psi))
    # Terme non-linéaire (type self-interaction)
    nonlinear = -0.1 * Psi * np.abs(Psi)**2
    return lap * 0.1 + nonlinear

# Opérateur de rotation B_hat (introduit une chiralité/torsion)
def B_hat(Psi):
    """Applique l'opérateur de rotation (torsion/structure)"""
    # Création d'un champ de rotation avec une faible cohérence à grande échelle
    kx, ky = np.meshgrid(np.fft.fftfreq(NX), np.fft.fftfreq(NY), indexing='ij')
    
    # Rotation avec une direction préférentielle faible
    omega_x = 0.3 * np.exp(-(kx**2 + ky**2) / 0.1)  # Filtre passe-bas
    omega_y = 0.3 * np.exp(-(kx**2 + ky**2) / 0.1)
    
    # Ajout de bruit pour les petites échelles
    noise_scale = 0.05
    omega_x += noise_scale * np.random.randn(NX, NY)
    omega_y += noise_scale * np.random.randn(NX, NY)
    
    # Application de la rotation (dérivée directionnelle)
    dPsi_dx = np.gradient(np.real(Psi), DX, axis=0) + 1j * np.gradient(np.imag(Psi), DX, axis=0)
    dPsi_dy = np.gradient(np.real(Psi), DX, axis=1) + 1j * np.gradient(np.imag(Psi), DX, axis=1)
    
    return 1j * (omega_x * dPsi_dx + omega_y * dPsi_dy)

# Fonction pour calculer le bulk flow (flux moyen)
def compute_bulk_flow(Psi):
    """Calcule l'écoulement moyen (bulk flow) à partir du champ"""
    grad_x = np.gradient(np.real(Psi), DX, axis=0)
    grad_y = np.gradient(np.real(Psi), DX, axis=1)
    
    # Le flux est proportionnel au gradient de phase
    phase = np.angle(Psi)
    flow_x = np.gradient(phase, DX, axis=0) * np.abs(Psi)**2
    flow_y = np.gradient(phase, DX, axis=1) * np.abs(Psi)**2
    
    # Moyenne sur toute la grille
    return np.array([np.mean(flow_x), np.mean(flow_y)])

# Boucle principale de simulation
bulk_flows = []
print("Début de la simulation...")

for step in range(STEPS):
    # Sauvegarde périodique
    if step % 50 == 0:
        Psi_history.append(Psi.copy())
        print(f"Step {step}/{STEPS}")
    
    # Calcul des opérateurs
    V_Psi = V_hat(Psi)
    B_Psi = B_hat(Psi)
    
    # Calcul du commutateur [B, V]Psi (moteur de la dynamique)
    BV_Psi = B_hat(V_Psi)
    VB_Psi = V_hat(B_Psi)
    commutator = BV_Psi - VB_Psi
    
    # Évolution temporelle avec régularisation
    denominator = np.sqrt(np.abs(Psi)**2 + EPS**2)
    dPsi_dt = V_Psi + B_Psi + 0.1 * commutator / denominator
    
    # Mise à jour (Euler explicite)
    Psi = Psi + DT * dPsi_dt
    
    # Calcul et stockage du bulk flow
    bulk_flows.append(compute_bulk_flow(Psi))

print("Simulation terminée !")

# Conversion en array pour analyse
bulk_flows = np.array(bulk_flows)

# Visualisation des résultats
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# État initial et final
ax1 = axes[0, 0]
im1 = ax1.imshow(np.abs(Psi_history[0])**2, cmap='viridis')
ax1.set_title('Densité initiale')
plt.colorbar(im1, ax=ax1)

ax2 = axes[0, 1]
im2 = ax2.imshow(np.abs(Psi_history[-1])**2, cmap='viridis')
ax2.set_title('Densité finale')
plt.colorbar(im2, ax=ax2)

# Phase finale
ax3 = axes[0, 2]
im3 = ax3.imshow(np.angle(Psi_history[-1]), cmap='twilight')
ax3.set_title('Phase finale')
plt.colorbar(im3, ax=ax3)

# Évolution du bulk flow
ax4 = axes[1, 0]
ax4.plot(bulk_flows[:, 0], label='Flow X')
ax4.plot(bulk_flows[:, 1], label='Flow Y')
ax4.set_xlabel('Temps')
ax4.set_ylabel('Intensité du bulk flow')
ax4.set_title('Évolution du bulk flow')
ax4.legend()
ax4.grid(True)

# Trajectoire 2D du bulk flow
ax5 = axes[1, 1]
ax5.plot(bulk_flows[:, 0], bulk_flows[:, 1], 'b-', alpha=0.7)
ax5.plot(bulk_flows[0, 0], bulk_flows[0, 1], 'go', label='Début')
ax5.plot(bulk_flows[-1, 0], bulk_flows[-1, 1], 'ro', label='Fin')
ax5.set_xlabel('Flow X')
ax5.set_ylabel('Flow Y')
ax5.set_title('Trajectoire du bulk flow')
ax5.legend()
ax5.grid(True)
ax5.axis('equal')

# Spectre de puissance du bulk flow
ax6 = axes[1, 2]
fft_x = np.fft.fft(bulk_flows[:, 0])
freqs = np.fft.fftfreq(len(bulk_flows))
ax6.plot(freqs[1:len(freqs)//2], np.abs(fft_x[1:len(freqs)//2])**2)
ax6.set_xlabel('Fréquence')
ax6.set_ylabel('Puissance')
ax6.set_title('Spectre du bulk flow (composante X)')
ax6.grid(True)
ax6.set_xscale('log')
ax6.set_yscale('log')

plt.tight_layout()
plt.savefig('simulation_EPT_results.png', dpi=150)
plt.show()

# Analyse quantitative finale
print("\n=== ANALYSE QUANTITATIVE ===")
print(f"Bulk flow moyen final: X = {np.mean(bulk_flows[-100:, 0]):.4f}, Y = {np.mean(bulk_flows[-100:, 1]):.4f}")
print(f"Amplitude du bulk flow: {np.sqrt(np.mean(bulk_flows[-100:, 0]**2 + bulk_flows[-100:, 1]**2)):.4f}")
print(f"Fluctuations (écart-type): X = {np.std(bulk_flows[-100:, 0]):.4f}, Y = {np.std(bulk_flows[-100:, 1]):.4f}")

# Calcul de la cohérence spatiale
density_final = np.abs(Psi)**2
correlation = np.correlate(density_final.flatten(), density_final.flatten(), mode='full')
correlation = correlation / correlation.max()
print(f"Longueur de corrélation (estimée): {np.sum(correlation > 0.5):.0f} pixels")

# Vérification de la régularisation (zéro n'existe pas)
min_density = np.min(np.abs(Psi)**2)
print(f"Densité minimale (doit être > 0): {min_density:.6f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace, gaussian_filter
from scipy.fft import fft2, ifft2, fftfreq
from scipy.interpolate import interp1d
from scipy.stats import binned_statistic
import warnings
warnings.filterwarnings('ignore')

# ============================================================
# 1. PARAMÈTRES DE LA SIMULATION (HAUTE RÉSOLUTION)
# ============================================================
NX, NY = 256, 256  # Grille 256x256 (16x plus de points que la précédente)
DX = 1.0            # Pas d'espace
DT = 0.002          # Pas de temps réduit pour la stabilité
STEPS = 2000        # Plus d'itérations pour voir les échelles multiples
EPS = 0.005         # Paramètre de régularisation (principe "zéro n'existe pas")

# Paramètres physiques
G_NEWTON = 1.0      # Constante gravitationnelle (unités simu)
C_LIGHT = 10.0      # Vitesse de la lumière (réduite pour la stabilité)

# ============================================================
# 2. CHAMPS SIMULÉS
# ============================================================
print("Initialisation des champs...")

# 2.1 Champ EPT Psi (information primordiale)
np.random.seed(42)
Psi = 0.05 * (np.random.randn(NX, NY) + 1j * np.random.randn(NX, NY))

# 2.2 Champ de matière (baryonique + noire)
# Initialisé avec des fluctuations suivant un spectre de type cosmologique
k = 2 * np.pi * fftfreq(NX, DX)
kx, ky = np.meshgrid(k, k, indexing='ij')
k_mag = np.sqrt(kx**2 + ky**2)
k_mag[0, 0] = 1.0  # Éviter division par zéro

# Spectre de puissance initial de type matière noire (approximation)
Pk_initial = k_mag**(-2) * np.exp(-k_mag**2 / 0.1)

# Génération de phases aléatoires
phases = np.exp(2j * np.pi * np.random.rand(NX, NY))
matter_k = np.sqrt(Pk_initial) * phases
matter = np.real(ifft2(matter_k))

# Normalisation
matter = matter / np.std(matter) * 0.1
matter_initial = matter.copy()

# 2.3 Tenseur métrique (gravité émergente)
# Initialisé comme plat (Minkowski)
g_xx = np.ones((NX, NY))
g_yy = np.ones((NX, NY))
g_xy = np.zeros((NX, NY))

# Potentiel gravitationnel
phi_grav = np.zeros((NX, NY))

# Historiques pour l'analyse
Psi_history = []
matter_history = []
g_xx_history = []
bulk_flows = []

# ============================================================
# 3. OPÉRATEURS QFUNITY
# ============================================================

def V_hat(Psi, matter, g_xx, g_yy):
    """
    Opérateur de vibration V̂_ε.
    Inclut l'influence de la matière et de la métrique.
    """
    # Laplacien avec métrique
    lap = (g_xx * laplace(np.real(Psi), mode='constant') + 
           1j * g_xx * laplace(np.imag(Psi), mode='constant'))
    
    # Couplage à la matière (la matière courbe l'EPT)
    matter_coupling = -0.05 * matter * Psi
    
    # Terme non-linéaire EPT
    nonlinear = -0.02 * Psi * np.abs(Psi)**2
    
    # Terme de régularisation
    reg = EPS**2 * Psi / (np.abs(Psi)**2 + EPS**2)
    
    return lap * 0.1 + matter_coupling + nonlinear + reg

def B_hat(Psi, matter, g_xx, g_yy, g_xy):
    """
    Opérateur de rotation/torsion B̂_ε.
    Génère des structures cohérentes (filaments, vortex).
    """
    # Gradients avec métrique
    dPsi_dx = (g_xx * np.gradient(np.real(Psi), DX, axis=0) + 
               1j * g_xx * np.gradient(np.imag(Psi), DX, axis=0))
    dPsi_dy = (g_yy * np.gradient(np.real(Psi), DX, axis=1) + 
               1j * g_yy * np.gradient(np.imag(Psi), DX, axis=1))
    
    # Rotation avec cohérence à grande échelle
    # Construction d'un champ de rotation corrélé
    kx_2d, ky_2d = np.meshgrid(fftfreq(NX), fftfreq(NY), indexing='ij')
    k_mag_rot = np.sqrt(kx_2d**2 + ky_2d**2)
    
    # Filtre passe-bas pour la cohérence à grande échelle
    rot_amp = 0.15 * np.exp(-(k_mag_rot**2) / 0.05)
    
    # Phases aléatoires pour les petites échelles
    omega_x = np.real(ifft2(rot_amp * np.exp(2j * np.pi * np.random.rand(NX, NY))))
    omega_y = np.real(ifft2(rot_amp * np.exp(2j * np.pi * np.random.rand(NX, NY))))
    
    # Application de la rotation (terme de torsion)
    rotation = 1j * (omega_x * dPsi_dx + omega_y * dPsi_dy)
    
    # Couplage torsion-matière
    matter_torsion = 0.03 * matter * Psi * (omega_x + omega_y)
    
    return rotation + matter_torsion

def compute_commutator(Psi, matter, g_xx, g_yy, g_xy):
    """
    Calcule [B̂, V̂]Ψ, le moteur de la dynamique.
    """
    V_Psi = V_hat(Psi, matter, g_xx, g_yy)
    B_Psi = B_hat(Psi, matter, g_xx, g_yy, g_xy)
    
    BV_Psi = B_hat(V_Psi, matter, g_xx, g_yy, g_xy)
    VB_Psi = V_hat(B_Psi, matter, g_xx, g_yy)
    
    return BV_Psi - VB_Psi

def evolve_matter(matter, Psi, phi_grav, g_xx, g_yy):
    """
    Évolution du champ de matière sous l'influence de l'EPT et de la gravité.
    Équation de continuité + force de torsion.
    """
    # Gradients de matière
    grad_matter_x = np.gradient(matter, DX, axis=0)
    grad_matter_y = np.gradient(matter, DX, axis=1)
    
    # Force de torsion (couplage à l'EPT)
    torsion_force_x = 0.02 * np.real(Psi) * grad_matter_x
    torsion_force_y = 0.02 * np.real(Psi) * grad_matter_y
    
    # Force gravitationnelle
    grav_force_x = -G_NEWTON * np.gradient(phi_grav, DX, axis=0) * matter
    grav_force_y = -G_NEWTON * np.gradient(phi_grav, DX, axis=1) * matter
    
    # Courants
    J_x = matter * (torsion_force_x + grav_force_x)
    J_y = matter * (torsion_force_y + grav_force_y)
    
    # Équation de continuité: ∂ρ/∂t = -∇·J
    dmatter_dt = -(np.gradient(J_x, DX, axis=0) + np.gradient(J_y, DX, axis=1))
    
    # Diffusion minimale pour stabilité numérique
    dmatter_dt += 0.001 * laplace(matter, mode='constant')
    
    return dmatter_dt

def compute_gravity(Psi, matter, g_xx, g_yy, g_xy):
    """
    Calcule le potentiel gravitationnel à partir des sources.
    Version simplifiée de l'équation de Poisson.
    """
    # Source: matière + énergie de l'EPT + terme de torsion
    source = (matter + 
              0.1 * np.abs(Psi)**2 + 
              0.05 * np.abs(compute_commutator(Psi, matter, g_xx, g_yy, g_xy))**2)
    
    # Résolution de l'équation de Poisson dans l'espace de Fourier
    source_k = fft2(source)
    kx_2d, ky_2d = np.meshgrid(fftfreq(NX), fftfreq(NY), indexing='ij')
    k2 = (2 * np.pi * kx_2d)**2 + (2 * np.pi * ky_2d)**2
    k2[0, 0] = 1.0  # Éviter division par zéro
    
    phi_k = -4 * np.pi * G_NEWTON * source_k / k2
    phi_k[0, 0] = 0.0  # Mode zéro
    
    return np.real(ifft2(phi_k))

def compute_metric(phi_grav):
    """
    Calcule les composantes métriques à partir du potentiel gravitationnel.
    Approximation linéaire: g_ij = (1 - 2φ) δ_ij pour la partie spatiale.
    """
    g_xx_new = 1.0 - 2 * phi_grav / C_LIGHT**2
    g_yy_new = 1.0 - 2 * phi_grav / C_LIGHT**2
    g_xy_new = np.zeros_like(phi_grav)
    
    return g_xx_new, g_yy_new, g_xy_new

def compute_bulk_flow(Psi, matter):
    """
    Calcule l'écoulement moyen (bulk flow) à partir du champ Psi et de la matière.
    """
    # Phase de Psi (information directionnelle)
    phase = np.angle(Psi)
    grad_phase_x = np.gradient(phase, DX, axis=0)
    grad_phase_y = np.gradient(phase, DX, axis=1)
    
    # Flux pondéré par la matière
    flow_x = np.mean(grad_phase_x * matter)
    flow_y = np.mean(grad_phase_y * matter)
    
    return np.array([flow_x, flow_y])

def compute_power_spectrum(field, DX):
    """
    Calcule le spectre de puissance 1D d'un champ 2D.
    """
    field_k = fft2(field)
    psd2d = np.abs(field_k)**2
    
    kx = fftfreq(NX, DX)
    ky = fftfreq(NY, DX)
    kx, ky = np.meshgrid(kx, ky, indexing='ij')
    k_mag = np.sqrt(kx**2 + ky**2)
    k_mag_flat = k_mag.flatten()
    psd_flat = psd2d.flatten()
    
    # Binning en k
    k_bins = np.linspace(0, np.max(k_mag)/2, 50)
    k_centers = (k_bins[1:] + k_bins[:-1]) / 2
    psd_1d, _, _ = binned_statistic(k_mag_flat, psd_flat, statistic='mean', bins=k_bins)
    
    return k_centers[psd_1d > 0], psd_1d[psd_1d > 0]

# ============================================================
# 4. BOUCLE PRINCIPALE
# ============================================================
print("Début de la simulation haute résolution...")

for step in range(STEPS):
    # Sauvegarde périodique
    if step % 200 == 0:
        Psi_history.append(Psi.copy())
        matter_history.append(matter.copy())
        g_xx_history.append(g_xx.copy())
        print(f"Step {step}/{STEPS} - Bulk flow: ({bulk_flows[-1][0]:.4f}, {bulk_flows[-1][1]:.4f})" if bulk_flows else f"Step {step}/{STEPS}")
    
    # Calcul du commutateur (moteur de la dynamique)
    commutator = compute_commutator(Psi, matter, g_xx, g_yy, g_xy)
    
    # Évolution de Psi
    V_Psi = V_hat(Psi, matter, g_xx, g_yy)
    B_Psi = B_hat(Psi, matter, g_xx, g_yy, g_xy)
    
    denominator = np.sqrt(np.abs(Psi)**2 + EPS**2)
    dPsi_dt = V_Psi + B_Psi + 0.05 * commutator / denominator
    
    Psi = Psi + DT * dPsi_dt
    
    # Évolution de la matière
    dmatter_dt = evolve_matter(matter, Psi, phi_grav, g_xx, g_yy)
    matter = matter + DT * dmatter_dt
    matter = np.maximum(matter, 0)  # Pas de densité négative
    
    # Mise à jour de la gravité
    phi_grav = compute_gravity(Psi, matter, g_xx, g_yy, g_xy)
    g_xx, g_yy, g_xy = compute_metric(phi_grav)
    
    # Calcul du bulk flow
    bulk_flows.append(compute_bulk_flow(Psi, matter))

print("Simulation terminée !")

# Conversion en array
bulk_flows = np.array(bulk_flows)

# ============================================================
# 5. ANALYSE DES RÉSULTATS
# ============================================================
print("\n" + "="*60)
print("ANALYSE QUANTITATIVE DES RÉSULTATS")
print("="*60)

# 5.1 Statistiques de base
print(f"\n--- Statistiques générales ---")
print(f"Densité de matière finale: moyenne = {np.mean(matter):.4f}, écart-type = {np.std(matter):.4f}")
print(f"Champ Psi final: |Ψ| moyen = {np.mean(np.abs(Psi)):.4f}")
print(f"Densité minimale (principe zéro): {np.min(np.abs(Psi)**2):.6f} (doit être > 0)")
print(f"Bulk flow final moyen (derniers 200 pas): X = {np.mean(bulk_flows[-200:, 0]):.4f}, Y = {np.mean(bulk_flows[-200:, 1]):.4f}")
print(f"Amplitude du bulk flow: {np.sqrt(np.mean(bulk_flows[-200:, 0]**2 + bulk_flows[-200:, 1]**2)):.4f}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace, gaussian_filter
from scipy.fft import fft2, ifft2, fftfreq
from scipy.interpolate import interp1d
from scipy.stats import binned_statistic
from scipy.signal import welch
import warnings
warnings.filterwarnings('ignore')

# ============================================================
# 1. PARAMÈTRES DE LA SIMULATION
# ============================================================
NX, NY = 256, 256
DX = 1.0
DT = 0.002
STEPS = 2000
EPS = 0.005

# Paramètres physiques
G_NEWTON = 1.0
C_LIGHT = 10.0

# Paramètres QFunity avancés
GLOBAL_TORSION_AMP = 0.15  # Amplitude de la torsion globale (Ω_EPT)
GLOBAL_TORSION_FREQ = 0.05  # Fréquence de la respiration cosmique
PHOTON_ABSORPTION_RATE = 0.001  # Taux de perte d'énergie photonique par unité de distance

# ============================================================
# 2. INITIALISATION DES CHAMPS
# ============================================================
print("Initialisation des champs avec torsion globale...")

np.random.seed(42)
Psi = 0.05 * (np.random.randn(NX, NY) + 1j * np.random.randn(NX, NY))

# Champ de matière avec spectre cosmologique
k = 2 * np.pi * fftfreq(NX, DX)
kx, ky = np.meshgrid(k, k, indexing='ij')
k_mag = np.sqrt(kx**2 + ky**2)
k_mag[0, 0] = 1.0
Pk_initial = k_mag**(-1.5) * np.exp(-k_mag**2 / 0.1)  # Pente -1.5 typique
phases = np.exp(2j * np.pi * np.random.rand(NX, NY))
matter_k = np.sqrt(Pk_initial) * phases
matter = np.real(ifft2(matter_k))
matter = matter / np.std(matter) * 0.1
matter_initial = matter.copy()

# Champs gravitationnels
g_xx = np.ones((NX, NY))
g_yy = np.ones((NX, NY))
g_xy = np.zeros((NX, NY))
phi_grav = np.zeros((NX, NY))

# Historiques
Psi_history = []
matter_history = []
bulk_flows = []
pk_history = []  # Historique des spectres de puissance
photon_energy_history = []  # Historique de l'énergie photonique

# ============================================================
# 3. OPÉRATEURS QFUNITY AVEC TORSION GLOBALE
# ============================================================

def V_hat(Psi, matter, g_xx, g_yy):
    """Opérateur de vibration V̂_ε."""
    lap = (g_xx * laplace(np.real(Psi), mode='constant') + 
           1j * g_xx * laplace(np.imag(Psi), mode='constant'))
    matter_coupling = -0.05 * matter * Psi
    nonlinear = -0.02 * Psi * np.abs(Psi)**2
    reg = EPS**2 * Psi / (np.abs(Psi)**2 + EPS**2)
    return lap * 0.1 + matter_coupling + nonlinear + reg

def B_hat(Psi, matter, g_xx, g_yy, g_xy, t):
    """
    Opérateur de rotation/torsion B̂_ε avec injection de torsion globale.
    La torsion globale simule Ω_EPT, la respiration cosmique.
    """
    # Gradients standards
    dPsi_dx = (g_xx * np.gradient(np.real(Psi), DX, axis=0) + 
               1j * g_xx * np.gradient(np.imag(Psi), DX, axis=0))
    dPsi_dy = (g_yy * np.gradient(np.real(Psi), DX, axis=1) + 
               1j * g_yy * np.gradient(np.imag(Psi), DX, axis=1))
    
    # TORSION GLOBALE INJECTÉE (Ω_EPT)
    # Elle simule la respiration cosmique : une rotation cohérente à grande échelle
    omega_global_x = GLOBAL_TORSION_AMP * np.ones((NX, NY))
    omega_global_y = GLOBAL_TORSION_AMP * np.ones((NX, NY))
    
    # Modulation temporelle (respiration)
    respiration = np.sin(2 * np.pi * GLOBAL_TORSION_FREQ * t * DT)
    omega_global_x *= respiration
    omega_global_y *= respiration
    
    # Torsion locale (fluctuations)
    kx_2d, ky_2d = np.meshgrid(fftfreq(NX), fftfreq(NY), indexing='ij')
    k_mag_rot = np.sqrt(kx_2d**2 + ky_2d**2)
    rot_amp = 0.1 * np.exp(-(k_mag_rot**2) / 0.1)
    omega_local_x = np.real(ifft2(rot_amp * np.exp(2j * np.pi * np.random.rand(NX, NY))))
    omega_local_y = np.real(ifft2(rot_amp * np.exp(2j * np.pi * np.random.rand(NX, NY))))
    
    # Torsion totale = globale (respiration) + locale (fluctuations)
    omega_x = omega_global_x + omega_local_x
    omega_y = omega_global_y + omega_local_y
    
    # Application de la rotation
    rotation = 1j * (omega_x * dPsi_dx + omega_y * dPsi_dy)
    matter_torsion = 0.03 * matter * Psi * (omega_x + omega_y)
    
    return rotation + matter_torsion

def compute_commutator(Psi, matter, g_xx, g_yy, g_xy, t):
    """Calcule [B̂, V̂]Ψ avec dépendance temporelle."""
    V_Psi = V_hat(Psi, matter, g_xx, g_yy)
    B_Psi = B_hat(Psi, matter, g_xx, g_yy, g_xy, t)
    
    BV_Psi = B_hat(V_Psi, matter, g_xx, g_yy, g_xy, t)
    VB_Psi = V_hat(B_Psi, matter, g_xx, g_yy)
    
    return BV_Psi - VB_Psi

def evolve_matter(matter, Psi, phi_grav, g_xx, g_yy, t):
    """Évolution de la matière avec couplage à la torsion globale."""
    grad_matter_x = np.gradient(matter, DX, axis=0)
    grad_matter_y = np.gradient(matter, DX, axis=1)
    
    # Force de torsion incluant la composante globale
    torsion_force_x = 0.02 * np.real(Psi) * grad_matter_x * (1 + GLOBAL_TORSION_AMP * np.sin(2 * np.pi * GLOBAL_TORSION_FREQ * t * DT))
    torsion_force_y = 0.02 * np.real(Psi) * grad_matter_y * (1 + GLOBAL_TORSION_AMP * np.sin(2 * np.pi * GLOBAL_TORSION_FREQ * t * DT))
    
    # Force gravitationnelle
    grav_force_x = -G_NEWTON * np.gradient(phi_grav, DX, axis=0) * matter
    grav_force_y = -G_NEWTON * np.gradient(phi_grav, DX, axis=1) * matter
    
    # Courants
    J_x = matter * (torsion_force_x + grav_force_x)
    J_y = matter * (torsion_force_y + grav_force_y)
    
    # Équation de continuité
    dmatter_dt = -(np.gradient(J_x, DX, axis=0) + np.gradient(J_y, DX, axis=1))
    dmatter_dt += 0.001 * laplace(matter, mode='constant')
    
    return dmatter_dt

def compute_gravity(Psi, matter, g_xx, g_yy, g_xy, t):
    """Calcule le potentiel gravitationnel."""
    # Source avec modulation temporelle due à la respiration
    commutator = compute_commutator(Psi, matter, g_xx, g_yy, g_xy, t)
    source = (matter + 
              0.1 * np.abs(Psi)**2 + 
              0.05 * np.abs(commutator)**2)
    
    # Résolution de Poisson
    source_k = fft2(source)
    kx_2d, ky_2d = np.meshgrid(fftfreq(NX), fftfreq(NY), indexing='ij')
    k2 = (2 * np.pi * kx_2d)**2 + (2 * np.pi * ky_2d)**2
    k2[0, 0] = 1.0
    phi_k = -4 * np.pi * G_NEWTON * source_k / k2
    phi_k[0, 0] = 0.0
    return np.real(ifft2(phi_k))

def compute_metric(phi_grav):
    """Calcule les composantes métriques."""
    g_xx_new = 1.0 - 2 * phi_grav / C_LIGHT**2
    g_yy_new = 1.0 - 2 * phi_grav / C_LIGHT**2
    g_xy_new = np.zeros_like(phi_grav)
    return g_xx_new, g_yy_new, g_xy_new

def compute_bulk_flow(Psi, matter, t):
    """Calcule l'écoulement moyen avec modulation temporelle."""
    phase = np.angle(Psi)
    grad_phase_x = np.gradient(phase, DX, axis=0)
    grad_phase_y = np.gradient(phase, DX, axis=1)
    
    # Le bulk flow est modulé par la respiration cosmique
    respiration_factor = 1 + 0.5 * np.sin(2 * np.pi * GLOBAL_TORSION_FREQ * t * DT)
    
    flow_x = np.mean(grad_phase_x * matter) * respiration_factor
    flow_y = np.mean(grad_phase_y * matter) * respiration_factor
    
    return np.array([flow_x, flow_y])

def simulate_photon_propagation(Psi, matter, phi_grav, steps):
    """
    Simule la propagation d'un photon à travers l'EPT.
    Calcule la perte d'énergie cumulative due aux fluctuations O(ε).
    """
    # Création d'une ligne de visée à travers la grille
    x_path = np.linspace(0, NX-1, steps)
    y_path = NY // 2 * np.ones_like(x_path)
    
    # Échantillonnage des champs le long du chemin
    Psi_path = np.zeros(steps)
    matter_path = np.zeros(steps)
    phi_path = np.zeros(steps)
    
    for i, (x, y) in enumerate(zip(x_path, y_path)):
        xi, yi = int(x), int(y)
        Psi_path[i] = np.abs(Psi[xi, yi])
        matter_path[i] = matter[xi, yi]
        phi_path[i] = phi_grav[xi, yi]
    
    # Calcul de la perte d'énergie
    # dE/dx = -α * |[B̂,V̂]Ψ| * E  (modèle phénoménologique)
    commutator_strength = np.abs(Psi_path) * matter_path
    energy_loss_rate = PHOTON_ABSORPTION_RATE * commutator_strength
    
    # Énergie initiale (normalisée à 1)
    photon_energy = np.ones(steps)
    for i in range(1, steps):
        photon_energy[i] = photon_energy[i-1] * (1 - energy_loss_rate[i-1] * DX)
    
    return x_path, photon_energy, energy_loss_rate

def compute_power_spectrum(field, DX):
    """Calcule le spectre de puissance 1D."""
    field_k = fft2(field)
    psd2d = np.abs(field_k)**2
    
    kx = fftfreq(NX, DX)
    ky = fftfreq(NY, DX)
    kx, ky = np.meshgrid(kx, ky, indexing='ij')
    k_mag = np.sqrt(kx**2 + ky**2)
    k_mag_flat = k_mag.flatten()
    psd_flat = psd2d.flatten()
    
    k_bins = np.linspace(0, np.max(k_mag)/2, 50)
    psd_1d, _, _ = binned_statistic(k_mag_flat, psd_flat, statistic='mean', bins=k_bins)
    
    k_centers = (k_bins[1:] + k_bins[:-1]) / 2
    valid = psd_1d > 0
    return k_centers[valid], psd_1d[valid]

# ============================================================
# 4. BOUCLE PRINCIPALE
# ============================================================
print("Début de la simulation avec torsion globale...")

for step in range(STEPS):
    # Sauvegarde périodique
    if step % 200 == 0:
        Psi_history.append(Psi.copy())
        matter_history.append(matter.copy())
        # Calcul et stockage du spectre de puissance
        k_sim, pk_sim = compute_power_spectrum(matter, DX)
        pk_history.append(pk_sim)
        print(f"Step {step}/{STEPS}")
    
    # Calcul du commutateur (dépend du temps t pour la torsion globale)
    commutator = compute_commutator(Psi, matter, g_xx, g_yy, g_xy, step)
    
    # Évolution de Psi
    V_Psi = V_hat(Psi, matter, g_xx, g_yy)
    B_Psi = B_hat(Psi, matter, g_xx, g_yy, g_xy, step)
    
    denominator = np.sqrt(np.abs(Psi)**2 + EPS**2)
    dPsi_dt = V_Psi + B_Psi + 0.05 * commutator / denominator
    
    Psi = Psi + DT * dPsi_dt
    
    # Évolution de la matière
    dmatter_dt = evolve_matter(matter, Psi, phi_grav, g_xx, g_yy, step)
    matter = matter + DT * dmatter_dt
    matter = np.maximum(matter, 0)
    
    # Mise à jour de la gravité
    phi_grav = compute_gravity(Psi, matter, g_xx, g_yy, g_xy, step)
    g_xx, g_yy, g_xy = compute_metric(phi_grav)
    
    # Calcul du bulk flow
    bulk_flows.append(compute_bulk_flow(Psi, matter, step))

# Simulation de propagation photonique
print("\nSimulation de la propagation photonique...")
x_path, photon_energy, loss_rate = simulate_photon_propagation(Psi, matter, phi_grav, 1000)

print("Simulation terminée !")

# Conversion en arrays
bulk_flows = np.array(bulk_flows)
pk_history = np.array(pk_history)

# ============================================================
# 5. ANALYSE DES RÉSULTATS
# ============================================================
print("\n" + "="*60)
print("ANALYSE QUANTITATIVE DES RÉSULTATS")
print("="*60)

# 5.1 Statistiques de base
print(f"\n--- Statistiques générales ---")
print(f"Densité de matière finale: moyenne = {np.mean(matter):.4f}, écart-type = {np.std(matter):.4f}")
print(f"Champ Psi final: |Ψ| moyen = {np.mean(np.abs(Psi)):.4f}")
print(f"Densité minimale (principe zéro): {np.min(np.abs(Psi)**2):.6f}")

# 5.2 Analyse du bulk flow temporel
print(f"\n--- Analyse du bulk flow temporel ---")
print(f"Bulk flow moyen (final): X = {np.mean(bulk_flows[-500:, 0]):.4f}, Y = {np.mean(bulk_flows[-500:, 1]):.4f}")
print(f"Amplitude moyenne: {np.sqrt(np.mean(bulk_flows[-500:, 0]**2 + bulk_flows[-500:, 1]**2)):.4f}")

# Analyse spectrale du bulk flow
freqs, psd_x = welch(bulk_flows[:, 0], fs=1/DT, nperseg=256)
print(f"Fréquence dominante du bulk flow: {freqs[np.argmax(psd_x[1:])+1]:.4f} (attendue: {GLOBAL_TORSION_FREQ:.4f})")

# 5.3 Analyse de l'évolution temporelle de P(k)
print(f"\n--- Évolution temporelle du spectre de puissance ---")
pk_initial = pk_history[0]
pk_final = pk_history[-1]
pk_ratio = pk_final / pk_initial
print(f"Rapport d'amplitude P(k) final/initial: {np.mean(pk_ratio):.2f}")

# 5.4 Analyse de la perte photonique
print(f"\n--- Simulation de perte d'énergie photonique ---")
print(f"Perte d'énergie totale sur le trajet: {(1 - photon_energy[-1])*100:.2f}%")
print(f"Taux de perte moyen: {np.mean(loss_rate):.6f} par unité de distance")