trinity_sim.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# --- CONFIGURATION (Searching for Stable Particles) ---
GRID_SIZE = 50          # Size of the lattice (50x50 trinities)
STEPS = 500             # How many time steps to run
COUPLING_STRENGTH = 8.0 # Reduced slightly
COLLAPSE_SENSITIVITY = 40.0 # <<-- DECREASED: Loosening confinement to allow stable non-uniformity

# --- INITIALIZATION ---
# 1. Initialize Curvature (R) for each trinity
# Random values representing different expansion rates (R)
curvature_grid = np.random.rand(GRID_SIZE, GRID_SIZE)

# 2. Initialize Selector Phase (phi)
# Random starting phases between 0 and 2*pi
phase_grid = np.random.uniform(0, 2*np.pi, (GRID_SIZE, GRID_SIZE))

def get_neighbors(i, j, size):
    """Finds top, bottom, left, right neighbors (wrapping around edges)"""
    neighbors = []
    # (row, col) indices for 4 neighbors
    candidates = [
        ((i+1)%size, j), ((i-1)%size, j),
        (i, (j+1)%size), (i, (j-1)%size)
    ]
    return candidates

def update_step(frame_num):
    """The core logic of the TGM model, running once per time step."""
    global curvature_grid, phase_grid
    
    # Create copies so we don't overwrite data while reading it
    new_curvature = curvature_grid.copy()
    new_phase = phase_grid.copy()
    
    size = GRID_SIZE
    
    for i in range(size):
        for j in range(size):
            # Get current trinity state
            R_current = curvature_grid[i, j]
            phi_current = phase_grid[i, j]
            
            neighbors = get_neighbors(i, j, size)
            
            # --- LOGIC FROM SECTION 7.2: Curvature Mismatch Collapse (Particle Formation) ---
            # Eq 48: P_collapse depends on Delta_R. High mismatch -> Low P_collapse (stabilizes difference)
            
            for ni, nj in neighbors:
                R_neighbor = curvature_grid[ni, nj]
                delta_R = abs(R_current - R_neighbor)
                
                # Calculate Probability (simplified Eq 48)
                # If mismatch is HIGH, probability is LOW. This keeps particles separate.
                prob_collapse = np.exp(-delta_R * COLLAPSE_SENSITIVITY)
                
                # Roll the dice: If P is high, they collapse and average their curvature
                if np.random.random() < prob_collapse:
                    avg_R = (R_current + R_neighbor) / 2.0
                    new_curvature[i, j] = avg_R
                    new_curvature[ni, nj] = avg_R

            # --- LOGIC FROM SECTION 7.3: Phase Synchronization (Coherence) ---
            # Eq 50: Phases couple like magnets
            
            phase_pull = 0.0
            for ni, nj in neighbors:
                phi_neighbor = phase_grid[ni, nj]
                phase_pull += np.sin(phi_neighbor - phi_current)
            
            # Update phase: intrinsic rotation (R) + neighbor pull
            intrinsic_rate = R_current * 0.1 
            new_phase[i, j] = phi_current + intrinsic_rate + (COUPLING_STRENGTH * phase_pull)

    # Update the main grids
    curvature_grid = new_curvature
    phase_grid = new_phase
    
    # Update the plot
    im.set_array(curvature_grid)
    return [im]

# --- VISUALIZATION ---
fig, ax = plt.subplots()
plt.title("TGM Simulation: Curvature Mismatch Dynamics")
# We visualize Curvature (R). Dark spots = Low Curvature, Bright spots = High Curvature
im = ax.imshow(curvature_grid, cmap='inferno', animated=True)
plt.colorbar(im, label="Curvature Intensity")

print("Starting Simulation... Close the window to stop.")
# 50 millisecond interval between frames
ani = animation.FuncAnimation(fig, update_step, frames=STEPS, interval=50, blit=True)
plt.show()