python3 -m venv ca-env
source ca-env/bin/activate # On Windows: .\ca-env\Scripts\activate
pip install numpy matplotlib

# ca_thermo.py
import numpy as np
import matplotlib.pyplot as plt
import random
import time

# --- Configuration ---
GRID_SIZE = 60 # Grid will be GRID_SIZE x GRID_SIZE
NUM_PARTICLES = 1200 # Number of particles to simulate
ITERATIONS = 500 # Number of simulation steps
REFRESH_RATE = 0.01 # Time delay between frames for animation

# --- Initialization ---
def initialize_grid(size, num_particles):
    """
    Initializes a square grid with particles clustered in the center.
    Returns a NumPy array where 1 indicates a particle, 0 indicates empty.
    """
    grid = np.zeros((size, size), dtype=int)

    # Calculate a small square area in the center to place particles
    center_x, center_y = size // 2, size // 2
    cluster_radius = int(np.sqrt(num_particles) / 2) + 2 # Adjust radius based on particle count

    placed_particles = 0
    while placed_particles < num_particles:
        # Randomly pick a cell within the cluster radius
        x = random.randint(max(0, center_x - cluster_radius), min(size - 1, center_x + cluster_radius))
        y = random.randint(max(0, center_y - cluster_radius), min(size - 1, center_y + cluster_radius))

        if grid[x, y] == 0:
            grid[x, y] = 1
            placed_particles += 1
    
    return grid

# Example usage (not part of main loop yet)
if __name__ == '__main__':
    initial_grid = initialize_grid(GRID_SIZE, NUM_PARTICLES)
    print("Initial Grid (top-left corner):")
    print(initial_grid[:10, :10])
    
    # You can visualize this initial grid quickly:
    # plt.imshow(initial_grid, cmap='binary')
    # plt.title("Initial Particle Distribution")
    # plt.show()

python ca_thermo.py

# Add this function to ca_thermo.py
def get_empty_neighbors(grid, r, c):
    """
    Finds and returns a list of coordinates for empty 8-neighbors of (r, c).
    Applies reflective boundary conditions.
    """
    size = grid.shape[0]
    empty_neighbors = []
    
    # Define 8-directional offsets (Moore neighborhood)
    # Includes diagonal and cardinal directions
    offsets = [(-1, -1), (-1, 0), (-1, 1),
               (0, -1),           (0, 1),
               (1, -1),  (1, 0),  (1, 1)]

    for dr, dc in offsets:
        nr, nc = r + dr, c + dc

        # Apply reflective boundary conditions
        if nr < 0: nr = -nr - 1 # Reflect from top edge
        if nr >= size: nr = size - (nr - size) -1 # Reflect from bottom edge
        if nc < 0: nc = -nc - 1 # Reflect from left edge
        if nc >= size: nc = size - (nc - size) - 1 # Reflect from right edge
        
        # Ensure coordinates are within bounds after reflection (should be by design)
        nr = max(0, min(size - 1, nr))
        nc = max(0, min(size - 1, nc))
        
        if grid[nr, nc] == 0:
            empty_neighbors.append((nr, nc))
            
    return empty_neighbors

def update_grid(current_grid):
    """
    Applies the CA rules for one time step.
    Particles move to a random empty adjacent cell.
    """
    size = current_grid.shape[0]
    next_grid = np.copy(current_grid) # Start with a copy of the current state

    # Get coordinates of all particles in the current grid
    particle_coords = list(zip(*np.where(current_grid == 1)))
    random.shuffle(particle_coords) # Process particles in a random order

    # To handle multiple particles attempting to move to the same empty cell
    # and ensuring no two particles occupy the same cell in next_grid:
    # We'll keep track of proposed moves.
    proposed_moves = {} # (new_r, new_c) -> [(old_r, old_c), ...]
    
    # First pass: Collect all proposed moves
    for r, c in particle_coords:
        empty_neighbors = get_empty_neighbors(current_grid, r, c)
        
        if empty_neighbors:
            # Particle wants to move
            new_r, new_c = random.choice(empty_neighbors)
            if (new_r, new_c) not in proposed_moves:
                proposed_moves[(new_r, new_c)] = []
            proposed_moves[(new_r, new_c)].append((r, c))
        else:
            # Particle cannot move, stays in current position
            if (r, c) not in proposed_moves:
                proposed_moves[(r, c)] = []
            proposed_moves[(r, c)].append((r, c)) # Particle stays put

    # Second pass: Apply moves to next_grid
    # Reset next_grid to all zeros, then fill based on accepted moves
    next_grid.fill(0) 
    
    for (target_r, target_c), original_positions in proposed_moves.items():
        if len(original_positions) == 1:
            # Only one particle proposed this cell, it moves there
            next_grid[target_r, target_c] = 1
        else:
            # Multiple particles proposed this cell (collision).
            # Randomly select one particle to move there, others stay in their original spot.
            # OR, more simply for this model: none move, they bounce.
            # For simplicity of diffusion: if a collision occurs on the *target* cell,
            # we can decide that only one particle successfully moves there,
            # and the others remain in their original positions.
            # Let's say one randomly selected particle from the collision group moves.
            # The others effectively "bounce" and stay in their current_grid spot.
            
            # Note: This collision resolution is a simplification. 
            # A more detailed model might simulate elastic collisions.
            # Here, we prioritize one particle moving, reflecting the "diffusion" aspect.
            # We will pick one random particle to move to the target cell,
            # and all others in the 'original_positions' list for that 'target_cell'
            # will stay in their current position.
            
            # Pick one particle to successfully move
            winner_original_r, winner_original_c = random.choice(original_positions)
            next_grid[target_r, target_c] = 1 # The winner moves
            
            # The losers stay in their original cells (if those cells are still empty in next_grid)
            # This logic needs care, as original cells might already be vacated by other particles
            # or filled by a "winner".
            
            # A simpler way to handle collisions: If *any* collision occurs,
            # the colliding particles effectively don't move for this step.
            # This is simpler and still leads to equilibrium.
            # Let's adjust the logic slightly.

    # Revised collision handling for simplicity:
    # We iterate all particles, determine their intended move.
    # We then create a NEW grid. If a target cell is already taken in the NEW grid,
    # the current particle stays put. This is a common way to handle "simultaneous" updates.
    
    # Re-writing update_grid function for clarity and correctness:
    # A particle based approach often simplifies this.
    # Let's use a particle list and update the grid from it.
    
    # Get initial particle positions
    current_particle_positions = list(zip(*np.where(current_grid == 1)))
    random.shuffle(current_particle_positions) # Process in random order
    
    next_particle_positions = set() # Use a set to automatically handle unique positions
    
    # To manage conflicts: track which cells are *already claimed* for the next step.
    # This acts as a 'reservation' system for the next_grid.
    claimed_next_cells = set()
    
    for r, c in current_particle_positions:
        empty_neighbors = get_empty_neighbors(current_grid, r, c)
        
        moved = False
        if empty_neighbors:
            # Try to move to a random empty neighbor
            target_r, target_c = random.choice(empty_neighbors)
            
            # If this target cell hasn't been claimed by another particle for the next step
            if (target_r, target_c) not in claimed_next_cells:
                next_particle_positions.add((target_r, target_c))
                claimed_next_cells.add((target_r, target_c)) # Claim it!
                moved = True
        
        if not moved:
            # If no empty neighbors, or target cell was already claimed by another particle,
            # the particle stays in its current cell.
            # We must check if its current cell has NOT been vacated by another particle moving OUT
            # AND NOT claimed by another particle moving IN.
            # This specific collision model (try to move, if target taken, stay) is fine.
            if (r, c) not in claimed_next_cells:
                next_particle_positions.add((r, c))
                claimed_next_cells.add((r, c)) # Claim its current spot
            else:
                # This case means *another* particle claimed this spot.
                # This shouldn't happen with the 'claimed_next_cells' logic for particles.
                # It means a particle that wanted to move, couldn't, AND its original cell
                # was taken by another particle that moved.
                # For simplicity, if this happens, the particle is "lost" or stays put regardless.
                # But our current_grid logic ensures particles are distinct.
                # So the `else` here is mostly for robustness.
                # For a true "particle" model, just add (r,c) to next_particle_positions
                # if it couldn't move.
                next_particle_positions.add((r,c)) # Particle implicitly stays if it failed to move.

    # Convert the set of next positions back to a grid
    next_grid = np.zeros((size, size), dtype=int)
    for r, c in next_particle_positions:
        next_grid[r, c] = 1

    return next_grid

# Add this test inside your if __name__ == '__main__': block
if __name__ == '__main__':
    # ... (existing code) ...

    print("\nTesting get_empty_neighbors:")
    test_grid = np.array([
        [0, 1, 0],
        [1, 0, 1],
        [0, 1, 0]
    ])
    
    # Test a particle at (1,1) (middle)
    empty_n = get_empty_neighbors(test_grid, 1, 1)
    print(f"Empty neighbors for (1,1): {empty_n}") # Expected: any cell with 0 around (1,1)
    
    # Test a particle at (0,0) (corner with reflective boundaries)
    # Reflective boundaries: e.g., (-1,-1) reflects to (0,0), (-1,0) reflects to (0,0)
    # The `get_empty_neighbors` logic applies reflective transformation on coordinates,
    # then checks if that *transformed* coordinate on the `current_grid` is empty.
    # This means a particle can theoretically move to its own reflected image.
    # A more robust reflective boundary condition for particle movement is
    # to simply ensure `nr` and `nc` stay within `0` and `size-1`.
    # Let's simplify the reflective boundary in `get_empty_neighbors` for clarity.
    
    # Simplified reflective boundary logic for get_empty_neighbors:
    # Instead of reflecting coordinates for neighbor *checking*, 
    # we just clamp the coordinates, and a particle won't move onto itself.
    # A particle at (0,0) can't move to (-1,-1), but it can try to move to (0,0) if that was empty.
    # This is a bit of a tricky point in CAs and reflective boundaries.
    # For diffusion, just making sure moves are within bounds and into empty cells is key.
    # The current get_empty_neighbors is effectively just clamping `nr` and `nc` 
    # for checking emptiness. For actual particle movement, 
    # the reflected coordinate is where the particle *lands*.

    # Let's adjust `get_empty_neighbors` to be simpler: it just returns valid
    # empty neighbors *within* the grid bounds. The reflective bounce happens
    # when `update_grid` determines the *actual* new position.
    
    # Note: I'll stick to the provided `get_empty_neighbors` for now as it handles reflections for the *neighbors*.
    # A simpler reflection in `update_grid` might be more intuitive:
    # Clamp the next position to be within grid bounds.
    
    # Let's re-evaluate the boundary condition strategy in `update_grid`.
    # Current `get_empty_neighbors` will actually make a particle 'see' its own cell if it reflects there.
    # A cleaner approach for reflective boundaries is to apply reflection *after* calculating `nx, ny`.
    # So `get_empty_neighbors` just checks true adjacent cells, and `update_grid` handles the boundary reflection.

    # Revised `get_empty_neighbors` ( simpler, no reflection here)
    def get_empty_neighbors(grid, r, c):
        """
        Finds and returns a list of coordinates for empty 8-neighbors of (r, c) *within grid bounds*.
        """
        size = grid.shape[0]
        empty_neighbors = []
        offsets = [(-1, -1), (-1, 0), (-1, 1),
                   (0, -1),           (0, 1),
                   (1, -1),  (1, 0),  (1, 1)]

        for dr, dc in offsets:
            nr, nc = r + dr, c + dc
            if 0 <= nr < size and 0 <= nc < size and grid[nr, nc] == 0:
                empty_neighbors.append((nr, nc))
                
        return empty_neighbors

    # Now, `update_grid` needs to explicitly handle reflective boundaries when picking `target_r, target_c`.
    # No, wait, `get_empty_neighbors` already returns `(nr, nc)` that are *valid* cells in the grid.
    # If the particle tries to move *out of bounds*, that's where reflection applies.
    # My initial `get_empty_neighbors` *was* trying to incorporate reflection into what it considers a "neighbor."
    # Let's simplify. A particle just *chooses a direction*. Then we check if that direction lands it out of bounds.
    # If it lands out of bounds, we reflect it. Then, if the *final, reflected* cell is empty, it moves.

    # Let's scrap the previous `update_grid` and `get_empty_neighbors` and implement a more direct particle update.

# --- Revised Particle Update Logic ---
def update_grid_v2(current_grid):
    """
    Applies the CA rules for one time step using a particle-based update.
    Each particle attempts to move to a random adjacent cell.
    If the target cell is occupied or already claimed for the next step, it stays put.
    Reflective boundary conditions are applied to target coordinates.
    """
    size = current_grid.shape[0]
    
    # Step 1: Collect current particle positions
    current_particle_coords = list(zip(*np.where(current_grid == 1)))
    random.shuffle(current_particle_coords) # Process in a random order
    
    next_particle_positions = set() # Use a set to ensure uniqueness for the next state
    
    # Define 8-directional offsets (Moore neighborhood) including 'stay put' (0,0)
    # The 'stay put' option is implicitly handled if no move is successful.
    offsets = [(-1, -1), (-1, 0), (-1, 1),
               (0, -1),           (0, 1),
               (1, -1),  (1, 0),  (1, 1)]

    # Temporary grid to track which cells are claimed for the *next* state during this iteration
    # This prevents multiple particles from trying to move into the same cell simultaneously.
    # Initialize it to reflect current particles that *will* stay put.
    claimed_next_cells = set(current_particle_coords) # Initially, assume all particles stay where they are

    for r, c in current_particle_coords:
        possible_moves = []
        for dr, dc in offsets:
            nr, nc = r + dr, c + dc
            
            # Apply reflective boundary conditions to the new position
            if nr < 0: nr = -nr - 1
            if nr >= size: nr = size - 1 - (nr - size)
            if nc < 0: nc = -nc - 1
            if nc >= size: nc = size - 1 - (nc - size)
            
            # Ensure coordinates are still within bounds after reflection (should be by design now)
            nr = max(0, min(size - 1, nr))
            nc = max(0, min(size - 1, nc))
            
            # A particle cannot move to its own current cell, unless it's the only option.
            # But the rule is "move to empty adjacent cell".
            # For simplicity: just consider valid positions.
            possible_moves.append((nr, nc))
            
        # Try to move to a random empty-looking cell
        random.shuffle(possible_moves) # Shuffle options for this particle
        
        moved = False
        for target_r, target_c in possible_moves:
            # Check if target is empty in *current* grid and not yet claimed in *next* step.
            # Note: This is crucial. We check current_grid[target_r, target_c] == 0 for emptiness
            # and (target_r, target_c) not in next_particle_positions for a collision avoidance.
            
            # This logic is tricky for "simultaneous" update.
            # The most straightforward way to avoid simultaneous conflicts is:
            # 1. Particles propose moves.
            # 2. Conflicts are resolved (e.g., only one particle wins a cell, others stay put).
            # 3. Update the grid based on resolved moves.

            # Let's simplify the rule:
            # A particle (r,c) looks for an empty neighbor (nr, nc) in the *current_grid*.
            # If it finds one, it attempts to move there.
            # If *multiple* particles try to move to the same (nr, nc) in the *next_grid*,
            # then a conflict resolution is needed.

            # Easiest way for diffusion:
            # Create a `new_grid`. Iterate through *all cells* of `current_grid` in random order.
            # If `current_grid[r,c]` has a particle:
            #   Find empty neighbors in `current_grid`.
            #   If any, pick one `(nr, nc)`. Set `new_grid[r,c] = 0` and `new_grid[nr,nc] = 1`.
            #   This must be done carefully to prevent double-counting or overwriting.

            # A classic CA approach uses two grids: current and next.
            # Let's refine the `update_grid` based on this robust pattern.

    # --- Re-Revised update_grid for clarity and correctness ---
    # This implements the common "particle attempts to move to random empty neighbor" model.
    # It ensures no two particles occupy the same spot in the *next* grid.
def update_grid(current_grid):
    """
    Applies the CA rules for one time step.
    Each particle in current_grid attempts to move to a random empty adjacent cell.
    If multiple particles target the same cell, only one succeeds.
    """
    size = current_grid.shape[0]
    next_grid = np.zeros((size, size), dtype=int) # The new state of the grid

    # Get a list of all particle coordinates and shuffle them to randomize update order
    particle_locations = list(zip(*np.where(current_grid == 1)))
    random.shuffle(particle_locations)

    # To track which target cells have already been claimed by a successful move in this step
    claimed_target_cells = set()

    # To track particles that successfully moved or stayed put in a non-conflicting way
    successful_moves = {} # (new_r, new_c) -> (old_r, old_c)

    # First pass: Propose moves and resolve direct conflicts (two particles targeting the same *empty* cell)
    for r, c in particle_locations:
        # Get potential empty neighbors from the *current* grid state
        potential_moves = []
        for dr, dc in [(-1, -1), (-1, 0), (-1, 1),
                       (0, -1),           (0, 1),
                       (1, -1),  (1, 0),  (1, 1)]:
            
            nr, nc = r + dr, c + dc
            
            # Apply reflective boundary conditions immediately for the target cell
            if nr < 0: nr = -nr - 1
            if nr >= size: nr = size - 1 - (nr - size - 1)
            if nc < 0: nc = -nc - 1
            if nc >= size: nc = size - 1 - (nc - size - 1)
            
            # Ensure coordinates are within bounds
            nr = max(0, min(size - 1, nr))
            nc = max(0, min(size - 1, nc))
            
            # Only consider moving to an empty cell (in the current grid)
            if current_grid[nr, nc] == 0:
                potential_moves.append((nr, nc))
        
        # If there are empty neighbors, try to move
        if potential_moves:
            random.shuffle(potential_moves) # Randomize choice among empty neighbors
            for target_r, target_c in potential_moves:
                if (target_r, target_c) not in claimed_target_cells:
                    # This particle successfully claimed the target cell
                    successful_moves[(target_r, target_c)] = (r, c)
                    claimed_target_cells.add((target_r, target_c))
                    break # Particle has moved, exit this inner loop
            else:
                # No empty neighbor was available or all were claimed by other particles this step.
                # Particle stays put (at its original location (r,c)), if that spot isn't claimed.
                if (r, c) not in claimed_target_cells:
                    successful_moves[(r, c)] = (r, c)
                    claimed_target_cells.add((r, c))
                # Else: (r,c) was claimed by another particle moving into it, this particle is 'lost'
                # (This shouldn't happen with our logic: particles only move *out* of occupied cells)
        else:
            # No empty neighbors at all, particle stays put
            if (r, c) not in claimed_target_cells:
                successful_moves[(r, c)] = (r, c)
                claimed_target_cells.add((r, c))
            # Else: (r,c) was claimed, similar 'lost' scenario.

    # Second pass: Construct the next_grid based on resolved moves
    for (target_r, target_c) in successful_moves:
        next_grid[target_r, target_c] = 1

    return next_grid

# Add this visualization code to ca_thermo.py
def visualize_grid(grid, iteration, fig, ax):
    """
    Updates the matplotlib plot with the current grid state.
    """
    ax.clear() # Clear previous frame
    ax.imshow(grid, cmap='binary', origin='upper') # 'binary' for black/white, 'upper' for 0,0 at top-left
    ax.set_title(f"Thermodynamic Equilibrium Simulation (Iteration: {iteration})")
    ax.set_xticks([]) # Remove ticks
    ax.set_yticks([])
    plt.draw()
    plt.pause(REFRESH_RATE) # Pause for animation effect

# Add this to ca_thermo.py, replacing the if __name__ == '__main__': block
if __name__ == '__main__':
    current_grid = initialize_grid(GRID_SIZE, NUM_PARTICLES)

    # Set up the plot for animation
    fig, ax = plt.subplots(figsize=(6, 6))
    plt.ion() # Turn on interactive mode for live plotting

    print(f"Starting simulation for {ITERATIONS} iterations on a {GRID_SIZE}x{GRID_SIZE} grid with {NUM_PARTICLES} particles.")
    print("Initial state: Particles clustered in the center.")
    print("Goal: Observe particles diffuse and reach thermodynamic equilibrium.")
    
    for i in range(ITERATIONS):
        visualize_grid(current_grid, i + 1, fig, ax)
        current_grid = update_grid(current_grid)
        
    print("\nSimulation complete. Press Enter to close the plot.")
    plt.ioff() # Turn off interactive mode
    plt.show() # Keep the final plot open until closed manually

python ca_thermo.py