"""
Anti-matter Containment Simulation using Electromagnetic Field Dynamics

This module implements a deterministic simulation for charged particle confinement
in electromagnetic fields. The system uses multi-layer electromagnetic containment
with fractal-based predictive corrections to maintain stable particle trajectories.

The implementation includes:
- Lorentz force calculations for charged particle dynamics
- Boris algorithm and Runge-Kutta 4th order numerical integration
- Energy and angular momentum conservation tracking
- Multi-layer electromagnetic containment system
- Adaptive time-stepping based on physical stability criteria
- Realistic environmental perturbations

Physical Constants and Units:
- All calculations performed in SI units (meters, seconds, kilograms, etc.)
- Follows standard electromagnetic theory and relativistic corrections
- Implements conservation laws for validation

Author: [Original Author]
Date: [Date]
Version: Enhanced with scientific rigor and documentation
"""

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
import os
import datetime
import sys
import traceback

# Configure matplotlib backend for headless operation when needed
import matplotlib
if not sys.stdout.isatty():  # Detect non-interactive environment
    matplotlib.use('Agg')  # Use non-interactive backend for file output

class PhysicallyAccurateAntimatterSimulator:
    """
    Electromagnetic Containment Simulator for Charged Particles
    
    This class simulates the dynamics of charged particles (specifically antiprotons)
    confined within a multi-layer electromagnetic field system. The simulation
    incorporates realistic physical models including:
    
    - Lorentz force dynamics in combined electric and magnetic fields
    - Relativistic corrections for high-velocity particles
    - Conservation of energy and angular momentum
    - Multi-layer containment field architecture
    - Adaptive numerical integration methods
    - Environmental perturbation modeling
    
    The containment system uses a fractal-based predictive algorithm to anticipate
    and correct for instabilities before they lead to containment failure.
    
    Attributes:
        Physical constants (elementary_charge, particle masses, etc.)
        Field parameters (magnetic_field_strength, electric_field_gradient)
        Numerical integration settings (time steps, integration methods)
        Conservation tracking arrays (energy, angular momentum histories)
        Particles list containing all simulated particles
    """
    def __init__(self, output_dir=None):
        """
        Initialize the electromagnetic containment simulation.
        
        Sets up all physical constants, field parameters, numerical integration
        settings, and data tracking structures needed for the simulation.
        
        Args:
            output_dir (str, optional): Directory path for saving simulation results.
                                      If None, results are not saved to files.
        
        The initialization configures:
        - Fundamental physical constants (charge, mass, permeability, etc.)
        - Electromagnetic field parameters for containment
        - Numerical stability and integration parameters
        - Data structures for tracking particle histories and conservation laws
        """
        # Physical constants
        self.elementary_charge = 1.602e-19  # C
        self.proton_mass = 1.673e-27  # kg
        self.antiproton_mass = self.proton_mass  # kg
        self.vacuum_permeability = 4 * np.pi * 1e-7  # H/m (magnetic constant)
        self.vacuum_permittivity = 8.85e-12  # F/m (electric constant)
        self.speed_of_light = 299792458  # m/s
        
        # Confinement parameters - CORE SETTINGS
        self.magnetic_field_strength = 5.0  # Tesla
        self.electric_field_gradient = 1e4  # V/m^2 - Realistic value for antimatter containment
        self.containment_radius = 0.1  # meters
        
        # Advanced physical model parameters - NEW
        self.max_field_strength_physical = 45.0  # Tesla (world record for sustained fields)
        self.relativistic_corrections = True
        self.QED_effects_threshold = 1e18  # V/m (Schwinger limit field strength)
        self.plasma_effects_threshold = 1e15  # particles/m^3
        
        # Fractal controller parameters - CORE SETTINGS (PRESERVED)
        self.fractal_correction_strength = 0.2
        self.stability_threshold = 0.5
        self.max_correction_force = 1e-15  # Newtons
        self.fractal_memory_length = 10
        
        # Damping parameters - CORE SETTINGS (PRESERVED)
        self.damping_coefficient = 1e-12  # velocity damping
        self.adaptive_damping_threshold = 0.7
        self.adaptive_damping_factor = 3.0
        
        # Emergency parameters - CORE SETTINGS (PRESERVED)
        self.emergency_containment_threshold = 1.0
        self.pre_emergency_threshold = 0.5
        self.emergency_correction_multiplier = 15.0
        self.emergency_mode_active = False
        self.emergency_mode_duration = 100
        self.emergency_field_inversion = True
        self.field_inversion_duration = 20
        self.field_inversion_active = False
        self.field_inversion_count = 0
        self.pre_emergency_mode_active = False
        self.pre_emergency_duration = 50
        
        # Time stepping - CORE SETTINGS
        self.base_dt = 1e-10  # seconds
        self.min_dt = 1e-14  # minimum allowed time step
        self.instability_time_reduction_factor = 10.0
        self.adaptive_time_stepping = True
        self.courant_factor = 0.2  # For physical stability (Courant-Friedrichs-Lewy condition)
        
        # Enhanced settings - CAREFULLY TUNED (PRESERVED)
        self.harmonic_damping_enabled = True
        self.harmonic_frequencies = [1e3, 1e4, 1e5, 1e6, 1e7]
        self.harmonic_damping_factors = [0.3, 0.5, 0.7, 0.9, 1.1]
        self.resonance_monitoring_window = 20
        
        # Enhanced stability algorithms - NEW
        self.runge_kutta_order = 4  # Higher order integration
        self.symplectic_integration = True  # Better energy conservation
        self.boris_algorithm = True  # Specialized for charged particles in EM fields
        
        # Dynamic field adaptation - CORE SETTINGS
        self.dynamic_field_adaptation = True
        self.field_response_rate = 0.4
        self.field_memory_length = 20
        self.field_history = []
        self.field_adaptation_factor = 1.0
        self.field_adaptation_minimum = 1.3
        self.max_adaptation_factor = 2.0
        
        # Multi-layer containment - CORE SETTINGS
        self.multi_layer_containment = True
        self.containment_layers = [
            {'radius': 0.6, 'strength_factor': 2.0, 'active': True},
            {'radius': 0.8, 'strength_factor': 4.0, 'active': True},
            {'radius': 0.9, 'strength_factor': 8.0, 'active': True},
            {'radius': 1.0, 'strength_factor': 16.0, 'active': True}
        ]
        self.layer_independence = 0.7
        self.layer_cross_coupling = True
        self.layer_status_history = []
        
        # Enhanced stability guards - CORE SETTINGS (PRESERVED)
        self.numerical_stability_guard = True  # Essential for preventing numerical explosion
        self.force_ceiling = 1e9             # Maximum force magnitude allowed
        self.acceleration_ceiling = 1e11     # Maximum acceleration magnitude
        self.velocity_ceiling = 1e7          # Maximum velocity magnitude
        self.instability_ceiling = 1e10      # Maximum instability value
        
        # Cascaded control - CORE SETTINGS (PRESERVED)
        self.cascaded_control_enabled = True
        self.stability_thresholds = [0.05, 0.1, 0.3, 0.5, 0.8, 1.2, 2.0, 4.0]
        self.correction_strengths = [0.1, 0.2, 0.5, 1.0, 2.0, 4.0, 7.0, 10.0]
        self.exponential_response = True
        self.exponential_base = 1.5
        self.immediate_response = True
        
        # Instability metrics - CORE SETTINGS (PRESERVED)
        self.instability_weights = {
            'proximity': 0.15,
            'velocity': 0.4,
            'energy': 0.4,
            'fractal': 0.05,
            'acceleration': 0.2
        }
        
        # Physical metrics and conservation tracking - NEW
        self.energy_conservation = True
        self.total_energy_initial = 0.0
        self.total_energy_history = []
        self.energy_conservation_error = []
        self.angular_momentum_tracking = True
        self.total_angular_momentum_initial = np.zeros(3)
        self.angular_momentum_history = []
        self.angular_momentum_error = []
        self.kinetic_energy_history = []
        self.potential_energy_history = []
        self.field_energy_history = []
        
        # Realistic perturbation models - NEW
        self.realistic_perturbations = True
        self.earth_magnetic_field_fluctuation = 1e-9  # Tesla (natural background variations)
        self.lab_environment_fluctuation = 1e-6  # Tesla (typical lab environment)
        self.power_grid_fluctuation = 5e-5  # Tesla (from nearby power systems)
        self.thermal_noise_factor = 1.38e-23 * 300  # Boltzmann constant * room temperature
        self.cosmic_ray_flux = 1e-2  # particles/m^2/s
        
        # Simulation state
        self.particles = []
        self.state_history = []
        self.time_history = []
        self.energy_history = []
        self.stability_metric_history = []
        self.dt_history = []
        self.emergency_activations = []
        self.field_inversions = []
        self.pre_emergency_activations = []
        
        # Output directory
        self.output_dir = output_dir
        
        # Perturbation settings (PRESERVED)
        self.perturbation_mode = False
        self.perturbation_response_duration = 0
        self.perturbation_detection_sensitivity = 0.2
        self.perturbation_response_active = False
        self.perturbation_anticipation = True
        self.perturbation_history = []
        self.perturbation_response_strength = 1.2
        
    def add_antiparticle(self, position, velocity):
        """
        Add a charged particle (antiproton) to the simulation.
        
        Creates a new particle with specified initial conditions and adds it to
        the simulation. The particle is initialized with proper physical properties
        including mass, charge, and initial energy/momentum calculations.
        
        Args:
            position (array-like): 3D position vector [x, y, z] in meters
            velocity (array-like): 3D velocity vector [vx, vy, vz] in m/s
        
        The function enforces velocity limits for numerical stability while
        maintaining physical realism (limited to 10% speed of light).
        Each particle tracks its complete history for trajectory analysis.
        """
        """Add an antiproton to the simulation with initial position and velocity."""
        # Ensure position and velocity are within reasonable bounds
        position = np.array(position)
        velocity = np.array(velocity)
        
        # Apply initial velocity limit for stability - use speed of light as physical ceiling
        velocity_mag = np.linalg.norm(velocity)
        physical_limit = 0.1 * self.speed_of_light  # Cap at 10% speed of light
        artificial_limit = 1e5  # Original artificial limit for stability
        
        velocity_limit = min(physical_limit, artificial_limit)  # Use the more restrictive limit
        
        if velocity_mag > velocity_limit:
            velocity = velocity * (velocity_limit / velocity_mag)
        
        # Create particle with all required properties
        self.particles.append({
            'position': position,
            'velocity': velocity,
            'mass': self.antiproton_mass,
            'charge': -self.elementary_charge,
            'history': [],
            'instability_history': [],
            'acceleration_history': [],
            'last_acceleration': np.zeros(3),
            'perturbed': False,
            # Physical conservation properties
            'kinetic_energy': 0.5 * self.antiproton_mass * np.sum(velocity**2),
            'potential_energy': 0,
            'angular_momentum': np.cross(position, self.antiproton_mass * velocity)
        })
        
    def calculate_lorentz_force(self, particle):
        """
        Calculate the electromagnetic force on a charged particle using Lorentz force law.
        
        Implements F = q(E + v × B) with optional relativistic corrections.
        The electric field uses a quadrupole configuration for stable confinement,
        while the magnetic field provides cyclotron motion stabilization.
        
        Args:
            particle (dict): Particle data containing position, velocity, mass, charge
        
        Returns:
            numpy.ndarray: 3D force vector in Newtons
        
        Physical implementation:
        - Electric field: Quadrupole configuration E = k[x, y, -2z]
        - Magnetic field: Uniform field in z-direction
        - Multi-layer containment with field adaptation
        - Relativistic corrections for v > 0.01c
        - Force limiting for numerical stability
        """
        """Calculate the Lorentz force on a charged particle with both physical models and safety limits."""
        position = particle['position']
        velocity = particle['velocity']
        charge = particle['charge']
        
        # Base magnetic field (uniform in z direction)
        B_base = np.array([0, 0, self.magnetic_field_strength])
        
        # Base electric field (quadrupole configuration)
        # E = k[x, y, -2z] provides stable trapping in all three dimensions
        field_gradient = self.electric_field_gradient
        E_base = np.array([
            field_gradient * position[0],
            field_gradient * position[1],
            -2 * field_gradient * position[2]
        ])
        
        # Calculate forces from multi-layer containment if enabled
        if self.multi_layer_containment:
            return self.calculate_multi_layer_forces(particle, B_base, E_base)
        
        # Standard single-layer approach with field adaptation
        # Dynamic field adjustment if enabled
        if self.dynamic_field_adaptation and len(self.field_history) > 0:
            # Apply adaptation factor, ensuring it never goes below minimum
            adaptation_factor = max(self.field_adaptation_minimum, self.field_adaptation_factor)
            # Ensure adaptation factor is within reasonable bounds
            adaptation_factor = min(adaptation_factor, self.max_adaptation_factor)
            B = B_base * adaptation_factor
            E = E_base * adaptation_factor
            
            # NEW: Apply physical limit on magnetic field strength
            B_magnitude = np.linalg.norm(B)
            if B_magnitude > self.max_field_strength_physical:
                B = B * (self.max_field_strength_physical / B_magnitude)
        else:
            B = B_base
            E = E_base
        
        # Emergency field inversion if active
        if self.field_inversion_active:
            B = -B * 0.8  # Invert field but at reduced strength
            E = -E * 0.8
        
        # Apply relativistic corrections for high-velocity particles (v > 0.01c)
        if self.relativistic_corrections and np.linalg.norm(velocity) > 0.01 * self.speed_of_light:
            # Calculate relativistic parameters
            beta = np.linalg.norm(velocity) / self.speed_of_light
            # Note: gamma factor calculated but not used in current simplified model
            
            # Relativistic Lorentz force: F = q(E + v × B) with relativistic corrections
            magnetic_force = charge * np.cross(velocity, B)
            electric_force = charge * E
            
            # Additional relativistic magnetic force component (simplified model)
            v_cross_E = np.cross(velocity, E)
            relativistic_correction = charge * v_cross_E * beta**2 / self.speed_of_light
            
            total_force = magnetic_force + electric_force + relativistic_correction
        else:
            # Standard Lorentz force: F = q(E + v × B)
            magnetic_force = charge * np.cross(velocity, B)
            electric_force = charge * E
            total_force = magnetic_force + electric_force
        
        # Apply force ceiling for stability (PRESERVED)
        force_magnitude = np.linalg.norm(total_force)
        if force_magnitude > self.force_ceiling:
            total_force = total_force * (self.force_ceiling / force_magnitude)
        
        # Store potential energy for conservation tracking (NEW)
        particle['potential_energy'] = charge * np.sum(position * E)
        
        return total_force
    
    def calculate_multi_layer_forces(self, particle, B_base, E_base):
        """
        Calculate electromagnetic forces using multi-layer containment architecture.
        
        Implements a sophisticated containment system with multiple concentric
        field layers, each with different strength factors and activation zones.
        This provides graduated response to particle motion and improved stability.
        
        Args:
            particle (dict): Particle data structure
            B_base (numpy.ndarray): Base magnetic field vector
            E_base (numpy.ndarray): Base electric field vector
        
        Returns:
            numpy.ndarray: Combined force vector from all active layers
        
        Multi-layer system:
        - Inner layers: Minimal activation for stable particles
        - Outer layers: Increased activation near containment boundaries
        - Layer independence: Each layer contributes independently
        - Cross-coupling: Layers can influence each other's strength
        """
        """Calculate forces with multi-layer containment architecture and numerical stability safeguards."""
        position = particle['position']
        velocity = particle['velocity']
        charge = particle['charge']
        
        # Get distance from center
        r = np.linalg.norm(position)
        if r < 1e-15:  # Avoid division by zero
            r = 1e-15
            position = np.array([1e-15, 0, 0])  # Assign a tiny offset
        
        # Initialize combined force
        combined_magnetic_force = np.zeros(3)
        combined_electric_force = np.zeros(3)
        
        # Process each layer
        layer_statuses = []
        for i, layer in enumerate(self.containment_layers):
            # Check if this layer is active
            if not layer['active']:
                layer_statuses.append(0)  # Inactive
                continue
                
            # Normalized distance relative to this layer's containment radius
            layer_radius = layer['radius'] * self.containment_radius
            relative_distance = r / layer_radius
            
            # Layer strength factor
            strength_factor = layer['strength_factor']
            
            # Layer activation status (0=inactive, 1=minimal, 2=active, 3=full)
            if relative_distance < 0.3:
                layer_status = 1  # Minimal activation when particle is deep inside
            elif relative_distance < 0.7:
                layer_status = 2  # Normal activation
            else:
                layer_status = 3  # Full activation when near boundary
                
            layer_statuses.append(layer_status)
            
            # Calculate smoother transitions between activation states
            if layer_status == 1 and relative_distance > 0.25:
                # Blend between minimal and normal activation
                blend_factor = (relative_distance - 0.25) / 0.05
                blend_factor = min(1.0, max(0.0, blend_factor))
                effective_status = 1.0 + blend_factor
            elif layer_status == 2 and relative_distance > 0.65:
                # Blend between normal and full activation
                blend_factor = (relative_distance - 0.65) / 0.05
                blend_factor = min(1.0, max(0.0, blend_factor))
                effective_status = 2.0 + blend_factor
            else:
                effective_status = float(layer_status)
                
            # Calculate layer-specific adaptation
            if self.dynamic_field_adaptation and len(self.field_history) > 0:
                # Base adaptation with minimum enforcement
                adaptation = max(self.field_adaptation_minimum, self.field_adaptation_factor)
                # Cap adaptation factor
                adaptation = min(adaptation, self.max_adaptation_factor)
                
                # More gradual adaptation scaling based on activation level
                adaptation_scaling = 1.0 + 0.2 * (effective_status - 1.0)
                adaptation *= adaptation_scaling
            else:
                adaptation = 1.0
                
            # Apply field inversion if active
            inversion_factor = -0.5 if self.field_inversion_active else 1.0
                
            # Calculate specific layer fields with adaptation and layer strength
            B_layer = B_base * adaptation * strength_factor * inversion_factor
            E_layer = E_base * adaptation * strength_factor * inversion_factor
            
            # NEW: Apply physical maximum field strength limit
            B_norm = np.linalg.norm(B_layer)
            if B_norm > self.max_field_strength_physical:
                B_layer = B_layer * (self.max_field_strength_physical / B_norm)
            
            # Calculate electromagnetic forces for this containment layer
            if self.relativistic_corrections and np.linalg.norm(velocity) > 0.01 * self.speed_of_light:
                # Apply relativistic corrections for high-velocity particles
                beta = np.linalg.norm(velocity) / self.speed_of_light
                # Note: gamma factor not used in current simplified relativistic model
                
                # Standard Lorentz force components
                magnetic_force = charge * np.cross(velocity, B_layer)
                electric_force = charge * E_layer
                
                # Add relativistic correction term for high velocities
                velocity_cross_electric = np.cross(velocity, E_layer)
                relativistic_correction = charge * velocity_cross_electric * beta**2 / self.speed_of_light
                
                magnetic_force += relativistic_correction
            else:
                # Standard non-relativistic Lorentz force calculation
                magnetic_force = charge * np.cross(velocity, B_layer)
                electric_force = charge * E_layer
            
            # Apply numerical stability protection - VERY IMPORTANT
            if self.numerical_stability_guard:
                # Check for extremely large force components
                mag_force_mag = np.linalg.norm(magnetic_force)
                elec_force_mag = np.linalg.norm(electric_force)
                
                # Cap forces if they exceed reasonable thresholds
                if mag_force_mag > self.force_ceiling / 10:
                    magnetic_force = magnetic_force * (self.force_ceiling / 10 / mag_force_mag)
                if elec_force_mag > self.force_ceiling / 10:
                    electric_force = electric_force * (self.force_ceiling / 10 / elec_force_mag)
            
            # Add to combined forces with independence factor
            independence = self.layer_independence
            if i == 0:  # Always include the first layer fully
                combined_magnetic_force += magnetic_force
                combined_electric_force += electric_force
            else:
                # Apply cross-coupling if enabled
                if self.layer_cross_coupling and i > 0 and layer_statuses[i-1] > 0:
                    # Gentler coupling factor
                    coupling_factor = 1.0 + 0.1 * (layer_statuses[i-1] - 1)
                    magnetic_force *= coupling_factor
                    electric_force *= coupling_factor
                
                # Add this layer's forces with independence factor
                combined_magnetic_force = combined_magnetic_force * (1-independence) + magnetic_force * independence
                combined_electric_force = combined_electric_force * (1-independence) + electric_force * independence
        
        # Record layer status history
        self.layer_status_history.append(layer_statuses)
        if len(self.layer_status_history) > 50:
            self.layer_status_history.pop(0)
            
        # Calculate potential energy for conservation tracking
        particle['potential_energy'] = charge * np.sum(position * combined_electric_force)
            
        # Final force limiting for stability
        total_force = combined_magnetic_force + combined_electric_force
        force_magnitude = np.linalg.norm(total_force)
        if force_magnitude > self.force_ceiling:
            total_force = total_force * (self.force_ceiling / force_magnitude)
            
        return total_force
    
    def apply_harmonic_damping(self, particle, time_step):
        """
        Apply velocity-dependent damping to prevent resonance oscillations.
        
        Analyzes particle velocity patterns to detect and suppress harmful
        resonances that could lead to instability or containment failure.
        Uses frequency-targeted damping based on velocity magnitude and
        oscillation detection.
        
        Args:
            particle (dict): Particle data with velocity history
            time_step (float): Current simulation time step (currently unused)
        
        Returns:
            numpy.ndarray: Damping force vector to be applied
        
        Damping strategy:
        - Velocity-magnitude dependent damping coefficients
        - Oscillation detection through sign change analysis
        - Enhanced damping for detected resonance patterns
        - Force magnitude limiting for stability
        """
        """Apply enhanced frequency-targeted damping to prevent resonance."""
        if not self.harmonic_damping_enabled or len(particle['history']) < 3:
            return np.zeros(3)
            
        # Extract recent velocity history
        if len(particle['history']) >= self.resonance_monitoring_window:
            velocities = np.array([state['velocity'] for state in particle['history'][-self.resonance_monitoring_window:]])
        else:
            velocities = np.array([state['velocity'] for state in particle['history']])
        
        # If we don't have enough history, use simple damping
        if len(velocities) < 4:  # Need minimum data for analysis
            v_mag = np.linalg.norm(particle['velocity'])
            damping_factor = self.harmonic_damping_factors[0]
            return -damping_factor * particle['velocity']
        
        # Initialize damping force
        damping_force = np.zeros(3)
        
        # Apply basic velocity-dependent damping for stability
        v_mag = np.linalg.norm(particle['velocity'])
        if v_mag < 10:  # Low velocity range
            damping_factor = self.harmonic_damping_factors[0]
        elif v_mag < 100:  # Medium velocity range
            damping_factor = self.harmonic_damping_factors[1]
        else:  # High velocity range
            damping_factor = self.harmonic_damping_factors[2]
        
        damping_force = -damping_factor * particle['velocity']
            
        # Add special resonance breaking if velocity oscillation detected
        if len(particle['history']) > 5:
            # Look for alternating velocity signs (oscillation)
            recent_v = np.array([state['velocity'] for state in particle['history'][-5:]])
            for axis in range(3):
                # Skip if recent velocity is too small
                if np.max(np.abs(recent_v[:, axis])) < 1e-10:
                    continue
                
                sign_changes = np.sum(np.diff(np.signbit(recent_v[:, axis])))
                if sign_changes >= 3:  # Multiple sign changes = oscillation
                    # Apply extra damping to break resonance
                    damping_force[axis] -= 1.5 * np.sign(particle['velocity'][axis]) * abs(particle['velocity'][axis])
        
        # Ensure damping force is within reasonable bounds
        force_mag = np.linalg.norm(damping_force)
        if force_mag > self.force_ceiling / 10:
            damping_force *= (self.force_ceiling / 10) / force_mag
            
        return damping_force
    
    def update_dynamic_field(self, system_instability):
        """
        Adapt electromagnetic field strength based on system stability.
        
        Monitors system instability metrics and adjusts containment field
        strength accordingly. Higher instability triggers stronger fields
        to maintain confinement, while stable conditions allow field relaxation.
        
        Args:
            system_instability (float): Current system instability metric
        
        The adaptation algorithm:
        - Tracks instability history over multiple time steps
        - Increases field strength when instability trends upward
        - Gradually returns to baseline when system is stable
        - Enforces minimum and maximum adaptation factors
        """
        """Update field adaptation using enhanced approach."""
        if not self.dynamic_field_adaptation:
            return
        
        # Cap instability for numerical stability
        system_instability = min(system_instability, self.instability_ceiling)
            
        # Record current instability and field state
        self.field_history.append({
            'instability': system_instability,
            'adaptation_factor': self.field_adaptation_factor
        })
        
        # Keep history within memory length
        if len(self.field_history) > self.field_memory_length:
            self.field_history = self.field_history[-self.field_memory_length:]
        
        # If we have enough history, adjust field response
        if len(self.field_history) >= 3:
            # Check if instability is increasing
            recent_instabilities = [entry['instability'] for entry in self.field_history[-3:]]
            
            # Calculate instability trend
            instability_increasing = (recent_instabilities[-1] > recent_instabilities[0])
            
            # Standard field adaptation behavior for normal operation
            if instability_increasing:
                # More gradual response for stability
                response_rate = min(0.05, self.field_response_rate * 0.1)
                self.field_adaptation_factor *= (1 + response_rate)
            else:
                # Gradual return to normal
                self.field_adaptation_factor = 1.0 + (self.field_adaptation_factor - 1.0) * 0.95
            
            # Enforce minimum adaptation factor
            self.field_adaptation_factor = max(self.field_adaptation_minimum, self.field_adaptation_factor)
                
            # Enforce maximum adaptation factor
            self.field_adaptation_factor = min(self.max_adaptation_factor, self.field_adaptation_factor)
    
    def analyze_fractal_pattern(self, history):
        """
        Analyze particle trajectory for chaotic/fractal patterns.
        
        Examines particle position and velocity history to identify
        signatures of chaotic motion that could indicate impending
        containment failure. Uses scaling variance and trajectory
        curvature analysis.
        
        Args:
            history (list): Particle position and velocity history
        
        Returns:
            float: Fractal/chaos metric (higher values indicate more chaotic motion)
        
        Analysis methods:
        - Scaling behavior of inter-point distances
        - Velocity magnitude progression analysis
        - Exponential divergence detection
        - Variance-based chaos quantification
        """
        """Enhanced fractal pattern recognition from position history."""
        if len(history) < 3:
            return 0.0
        
        try:
            # Calculate trajectory curvature pattern
            positions = np.array([h['position'] for h in history])
            velocities = np.array([h['velocity'] for h in history])
            
            # Calculate distances between consecutive points
            distances = np.sqrt(np.sum(np.diff(positions, axis=0)**2, axis=1))
            if len(distances) < 2 or np.all(distances < 1e-20):
                return 0.0
                
            # Estimate scaling behavior (higher value = more chaotic)
            log_distances = np.log(distances + 1e-20)  # avoid log(0)
            scaling_variance = np.var(log_distances)
            
            # Check for exponential divergence
            velocity_magnitudes = np.sqrt(np.sum(velocities**2, axis=1))
            velocity_increase = np.diff(velocity_magnitudes)
            acceleration_factor = np.mean(velocity_increase) if len(velocity_increase) > 0 else 0
            
            # Combined chaotic pattern metric
            fractal_metric = scaling_variance + max(0, acceleration_factor)
            
            # Cap metric for stability
            return min(fractal_metric, 100.0)
            
        except Exception as e:
            print(f"Error in fractal analysis: {e}")
            return 0.0
    
    def calculate_acceleration_factor(self, particle):
        """
        Calculate instability factor from particle acceleration patterns.
        
        Analyzes recent acceleration history to identify trends that
        indicate increasing instability, such as consistent acceleration
        increases or high jerk (rate of acceleration change).
        
        Args:
            particle (dict): Particle with acceleration_history attribute
        
        Returns:
            float: Acceleration-based instability factor
        
        Analysis includes:
        - Acceleration magnitude trends
        - Jerk (third derivative of position) calculations
        - Consistent acceleration increase detection
        - Capped return values for numerical stability
        """
        """Calculate instability factor based on acceleration patterns."""
        if len(particle['acceleration_history']) < 3:
            return 0.0
            
        try:
            # Get recent accelerations
            recent_accels = np.array(particle['acceleration_history'][-5:])
            
            # Calculate acceleration magnitude increases
            accel_magnitudes = np.linalg.norm(recent_accels, axis=1)
            
            # Check for consistent acceleration increases
            accel_increases = np.diff(accel_magnitudes)
            if len(accel_increases) > 0 and np.all(accel_increases > 0):
                # Consistent acceleration - concerning
                return min(10.0, np.mean(accel_increases) * 10.0)
            
            # Calculate jerk (rate of change of acceleration)
            jerk_magnitudes = np.linalg.norm(np.diff(recent_accels, axis=0), axis=1)
            jerk_factor = min(10.0, np.mean(jerk_magnitudes) * 5.0) if len(jerk_magnitudes) > 0 else 0.0
            
            # Return capped factor for stability
            return min(10.0, max(np.mean(accel_magnitudes), jerk_factor))
            
        except Exception as e:
            print(f"Error in acceleration factor calculation: {e}")
            return 0.0
    
    def predict_instability(self, particle):
        """
        Predict future instability using multi-factor analysis.
        
        Combines multiple instability indicators to generate a comprehensive
        prediction of containment failure risk. Uses weighted combination
        of proximity, velocity, energy, fractal, and acceleration factors.
        
        Args:
            particle (dict): Complete particle data structure
        
        Returns:
            float: Combined instability prediction metric
        
        Prediction factors:
        - Proximity: Distance from containment center
        - Velocity: Radial velocity component (outward motion)
        - Energy: Kinetic energy relative to containment scales
        - Fractal: Chaotic motion patterns from trajectory analysis
        - Acceleration: Acceleration trend analysis
        """
        """Advanced fractal-based pattern recognition to predict future instabilities."""
        position = particle['position']
        velocity = particle['velocity']
        history = particle['history']
        
        try:
            # Basic distance-based assessment
            r = np.linalg.norm(position)
            proximity_factor = (r / self.containment_radius)**2
            
            # Velocity component toward/away from center
            if r > 1e-15:  # Avoid division by zero
                v_radial = np.dot(position, velocity) / r
            else:
                v_radial = 0
                
            velocity_factor = abs(v_radial) / 1e5
            
            # Energy-based assessment
            kinetic_energy = 0.5 * particle['mass'] * np.sum(velocity**2)
            energy_factor = kinetic_energy / 1e-16  # normalize to typical containment energies
            
            # Advanced fractal pattern recognition to predict future behavior
            fractal_factor = self.analyze_fractal_pattern(history) * 5.0
            
            # Acceleration-based prediction
            acceleration_factor = self.calculate_acceleration_factor(particle)
            
            # Combined instability prediction with weighted components
            weights = self.instability_weights
            instability = (
                weights.get('proximity', 0.15) * proximity_factor + 
                weights.get('velocity', 0.4) * velocity_factor + 
                weights.get('energy', 0.4) * energy_factor + 
                weights.get('fractal', 0.05) * fractal_factor +
                weights.get('acceleration', 0.2) * acceleration_factor
            )
            
            # Cap instability for numerical stability
            return min(instability, self.instability_ceiling)
            
        except Exception as e:
            print(f"Error in instability prediction: {e}")
            return 0.0
    
    def determine_correction_level(self, instability):
        """
        Determine appropriate containment field correction level.
        
        Maps instability levels to specific containment response levels
        in a cascaded control system. Higher instability triggers more
        aggressive containment measures.
        
        Args:
            instability (float): Current instability metric
        
        Returns:
            int: Correction level index (-1 for no correction, 0+ for correction levels)
        
        Correction levels correspond to different field strength multipliers
        and response strategies defined in the stability_thresholds array.
        """
        """Determine which level of containment field reinforcement is needed."""
        if not self.cascaded_control_enabled:
            # Default to single threshold behavior
            return 0 if instability > self.stability_threshold else -1
        
        try:    
            # Find the highest threshold that is exceeded to determine appropriate field response
            for i in range(len(self.stability_thresholds) - 1, -1, -1):
                if instability > self.stability_thresholds[i]:
                    return i
                    
            # Additional safety: if any instability is detected, use level 0
            if instability > 0.01 and self.immediate_response:
                return 0
                    
            return -1  # No threshold exceeded
            
        except Exception as e:
            print(f"Error in correction level determination: {e}")
            return 0  # Default to level 0 for safety
    
    def predict_containment_reinforcement(self, particle, instability):
        """
        Predict optimal containment field reinforcement strategy.
        
        Uses fractal analysis results to determine where and how strongly
        to reinforce the containment field. Returns both direction and
        intensity recommendations for field adjustment.
        
        Args:
            particle (dict): Particle requiring containment adjustment
            instability (float): Current instability level
        
        Returns:
            tuple: (direction_vector, intensity_factor)
                direction_vector: 3D unit vector for reinforcement direction
                intensity_factor: Multiplier for reinforcement strength
        
        Strategy:
        - Generally reinforces toward containment center
        - Intensity scales with instability level and distance
        - Exponential response option for severe instabilities
        - Position-dependent scaling for distance-based reinforcement
        """
        """
        Predict where containment field needs to be reinforced based on fractal analysis.
        Returns a reinforcement direction and intensity, not a force directly.
        """
        try:
            # Determine appropriate reinforcement level based on predicted instability
            correction_level = self.determine_correction_level(instability)
            
            # IMMEDIATE preemptive reinforcement even if below threshold
            if correction_level < 0 and self.immediate_response:
                # Apply minimum correction even for tiny instability
                if instability > 0.01:  # Use ultra-low threshold
                    correction_level = 0  # Use lowest correction level
            
            if correction_level < 0:
                return None, 0  # No reinforcement needed
                
            # Calculate reinforcement direction (toward center)
            position = particle['position']
            velocity = particle['velocity']
            r = np.linalg.norm(position)
            
            if r < 1e-15:  # Avoid division by zero
                # Just provide a damping recommendation
                return -velocity, 0.1  # Direction and low intensity
                
            # Direction vector where containment needs reinforcement
            direction = -position / r
            
            # Determine appropriate field strength based on correction level
            reinforcement_intensity = self.correction_strengths[correction_level]
            
            # Apply position-based scaling - stronger reinforcement when further out
            position_factor = min(5.0, max(1.0, r / (0.5 * self.containment_radius)))
            reinforcement_intensity *= position_factor
            
            # Scale reinforcement intensity based on how far above threshold
            level_threshold = self.stability_thresholds[correction_level] if correction_level > 0 else 0.01
            
            if self.exponential_response:
                # Exponential response based on how far above threshold
                excess_instability = instability / level_threshold
                excess_instability = min(10.0, excess_instability)  # Limit for stability
                
                exponential_factor = self.exponential_base ** (excess_instability - 1.0)
                exponential_factor = min(10.0, exponential_factor)  # Limit for stability
                reinforcement_intensity *= exponential_factor
            else:
                # Linear response
                reinforcement_intensity *= (instability - level_threshold)
            
            # Apply immediate intensity boost during high instability
            if instability > 1.0 and self.immediate_response:
                boost_factor = min(5.0, instability)  # Limit boost for stability
                reinforcement_intensity *= boost_factor
            
            # Limit maximum reinforcement intensity
            reinforcement_intensity = min(reinforcement_intensity, 10.0)  # Reasonable upper limit
                
            # Return the recommended reinforcement direction and intensity
            return direction, reinforcement_intensity
                
        except Exception as e:
            print(f"Error in containment reinforcement prediction: {e}")
            position = particle['position']
            r = np.linalg.norm(position)
            if r > 1e-15:
                # Safe default: suggest reinforcing toward center
                return -position / r, 1.0
            else:
                return None, 0
    
    def adjust_containment_field(self, particle, direction, intensity):
        """
        Apply calculated field adjustments as actual electromagnetic forces.
        
        Converts the abstract reinforcement recommendations from the fractal
        algorithm into concrete electromagnetic forces that can be applied
        to the particle during numerical integration.
        
        Args:
            particle (dict): Target particle for force application
            direction (numpy.ndarray): Direction vector for force application
            intensity (float): Force intensity multiplier
        
        Returns:
            numpy.ndarray: Force vector to be added to particle dynamics
        
        Implementation:
        - Converts intensity to appropriate force magnitude
        - Adds velocity damping in reinforced regions
        - Combines centering and damping forces
        - Applies force magnitude limits for stability
        """
        """
        Adjust the containment electromagnetic field based on fractal prediction.
        This function applies the actual forces based on the fractal algorithm's recommendations.
        """
        if direction is None:
            return np.zeros(3)  # No field adjustment needed
            
        try:
            position = particle['position']
            velocity = particle['velocity']
            mass = particle['mass']
            charge = particle['charge']
            
            # Convert reinforcement intensity to appropriate force magnitude
            # This is where the containment system interprets the fractal algorithm's recommendation
            base_force_magnitude = self.max_correction_force * intensity
            
            # Add velocity damping in the reinforced region
            damping_component = -self.damping_coefficient * velocity * 10.0 
            
            # Calculate the actual electromagnetic force to apply
            # Direction comes from fractal prediction, magnitude from containment system
            centering_force = direction * base_force_magnitude
            
            # Combine the forces for the final field adjustment
            total_correction = centering_force + damping_component
            
            # Apply force magnitude limits
            force_magnitude = np.linalg.norm(total_correction)
            max_force = self.max_correction_force * 100 * intensity  # Scale max force with intensity
            
            if force_magnitude > max_force:
                total_correction = total_correction * (max_force / force_magnitude)
                
            return total_correction
            
        except Exception as e:
            print(f"Error in containment field adjustment: {e}")
            return np.zeros(3)  # Return no force as a safe default
    
    def check_emergency_containment(self):
        """
        Monitor system for emergency containment activation conditions.
        
        Continuously monitors system stability metrics to determine when
        emergency containment measures should be activated. Uses multiple
        threshold levels for graduated emergency response.
        
        Returns:
            bool: True if emergency containment is currently active
        
        Emergency levels:
        - Pre-emergency: Enhanced monitoring and mild reinforcement
        - Emergency: Strong containment forces and field inversions
        - Field inversion: Temporary reversal of field polarity
        
        Activation is based on instability exceeding critical thresholds
        with immediate response to prevent containment failure.
        """
        """Ultra-sensitive emergency containment activation with IMMEDIATE response."""
        if len(self.stability_metric_history) < 2:
            return False
            
        try:
            # SIMPLIFIED: Just look at the most recent instability - respond IMMEDIATELY
            current_instability = self.stability_metric_history[-1]
            
            # Check for pre-emergency threshold
            if current_instability > self.pre_emergency_threshold and not self.emergency_mode_active and not self.pre_emergency_mode_active:
                print(f"PRE-EMERGENCY MODE ACTIVATED: Instability = {current_instability:.2f}")
                self.pre_emergency_mode_active = True
                self.pre_emergency_duration = 50
                self.pre_emergency_activations.append(len(self.time_history) - 1)
            
            # LOWERED threshold for much earlier intervention
            if current_instability > self.emergency_containment_threshold * 0.3:  # MORE SENSITIVE
                if not self.emergency_mode_active:
                    print(f"EMERGENCY CONTAINMENT ACTIVATED: Instability = {current_instability:.2f}")
                    self.emergency_activations.append(len(self.time_history) - 1)
                    
                self.emergency_mode_active = True
                self.emergency_mode_duration = 100  # Extended emergency mode duration
                self.pre_emergency_mode_active = False  # Emergency overrides pre-emergency
                
                # IMMEDIATE field inversion for ANY emergency
                if not self.field_inversion_active:
                    print(f"FIELD INVERSION ACTIVATED: Instability = {current_instability:.2f}")
                    self.field_inversions.append(len(self.time_history) - 1)
                self.field_inversion_active = True
                self.field_inversion_count = self.field_inversion_duration
                    
                return True
            
            # More aggressive emergency response for any instability above 1.0
            if current_instability > 1.0 and not self.emergency_mode_active:
                # Activate emergency mode at lower threshold
                print(f"PRE-EMPTIVE EMERGENCY CONTAINMENT: Instability = {current_instability:.2f}")
                self.emergency_activations.append(len(self.time_history) - 1)
                self.emergency_mode_active = True
                self.emergency_mode_duration = 20  # Shorter duration for pre-emptive response
                self.pre_emergency_mode_active = False  # Emergency overrides pre-emergency
            
            # Handle pre-emergency mode duration
            if self.pre_emergency_mode_active:
                self.pre_emergency_duration -= 1
                if self.pre_emergency_duration <= 0:
                    print("Pre-emergency mode deactivated")
                    self.pre_emergency_mode_active = False
            
            # Decrement emergency mode duration if active
            if self.emergency_mode_active:
                self.emergency_mode_duration -= 1
                if self.emergency_mode_duration <= 0:
                    print("Emergency containment deactivated")
                    self.emergency_mode_active = False
            
            # Handle field inversion timing separately
            if self.field_inversion_active:
                self.field_inversion_count -= 1
                if self.field_inversion_count <= 0:
                    print("Field inversion deactivated")
                    self.field_inversion_active = False
                    
            return self.emergency_mode_active
            
        except Exception as e:
            print(f"Error in emergency containment: {e}")
            return False
    
    def apply_emergency_containment(self, particle):
        """
        Apply maximum-strength emergency containment forces.
        
        Implements the strongest available containment measures when
        emergency conditions are detected. Uses aggressive braking,
        strong centering forces, and distance-dependent scaling.
        
        Args:
            particle (dict): Particle requiring emergency containment
        
        Returns:
            numpy.ndarray: Emergency containment force vector
        
        Emergency measures:
        - Ultra-strong velocity damping (immediate braking)
        - Massive centering forces toward containment center
        - Exponential force scaling with distance from center
        - Enhanced radial braking for outward-moving particles
        - Force ceiling enforcement for numerical stability
        """
        """Apply ULTRA-STRONG emergency containment forces."""
        try:
            position = particle['position']
            velocity = particle['velocity']
            
            # ULTRA-STRONG velocity damping - full stop
            braking_force = -velocity * 5.0
            
            # ULTRA-STRONG centering force
            r = np.linalg.norm(position)
            if r < 1e-15:  # Avoid division by zero
                return braking_force  # Just apply braking
                
            # Base centering force - MUCH stronger
            centering_factor = self.max_correction_force * self.emergency_correction_multiplier * 5.0
            
            # Add extreme position-dependent scaling - exponential increase with distance
            r_factor = min(10.0, (r / self.containment_radius) ** 3)
            centering_factor *= r_factor
            
            # Apply massive centering force for distant particles
            centering_force = -position / r * centering_factor
            
            # Add strong radial velocity damping - more damping for outward motion
            radial_dir = -position / r
            radial_vel = np.dot(velocity, radial_dir)
            
            if radial_vel < 0:  # Moving outward
                # Apply massive radial braking directly proportional to distance
                radial_brake = -radial_vel * radial_dir * 10.0 * r_factor
                braking_force += radial_brake
            
            # Combined emergency force with force ceiling
            total_force = braking_force + centering_force
            force_magnitude = np.linalg.norm(total_force)
            if force_magnitude > self.force_ceiling:
                total_force = total_force * (self.force_ceiling / force_magnitude)
                
            return total_force
            
        except Exception as e:
            print(f"Error in emergency containment: {e}")
            # Safe default: just apply strong braking
            return -particle['velocity'] * 10.0
    
    def apply_pre_emergency_containment(self, particle):
        """
        Apply moderate enhancement to containment forces during pre-emergency.
        
        Provides intermediate containment response between normal operation
        and full emergency mode. Helps prevent progression to emergency
        conditions through early intervention.
        
        Args:
            particle (dict): Particle requiring enhanced containment
        
        Returns:
            numpy.ndarray: Enhanced containment force vector
        
        Pre-emergency measures:
        - Moderately enhanced velocity damping
        - Increased centering forces with position scaling
        - Radial velocity braking for outward motion
        - Force limiting for numerical stability
        """
        """Apply moderately enhanced containment forces during pre-emergency mode."""
        try:
            position = particle['position']
            velocity = particle['velocity']
            
            # Enhanced velocity damping - not as strong as full emergency
            braking_force = -velocity * 2.0
            
            # Enhanced centering force
            r = np.linalg.norm(position)
            if r < 1e-15:  # Avoid division by zero
                return braking_force  # Just apply braking
                
            # Base centering force - moderately stronger
            centering_factor = self.max_correction_force * self.emergency_correction_multiplier * 2.0
            
            # Add position-dependent scaling
            r_factor = min(5.0, (r / self.containment_radius) ** 2)
            centering_factor *= r_factor
            
            # Apply centering force
            centering_force = -position / r * centering_factor
            
            # Add moderate radial velocity damping
            radial_dir = -position / r
            radial_vel = np.dot(velocity, radial_dir)
            
            if radial_vel < 0:  # Moving outward
                # Apply radial braking
                radial_brake = -radial_vel * radial_dir * 5.0 * r_factor
                braking_force += radial_brake
            
            # Combined pre-emergency force with force ceiling
            total_force = braking_force + centering_force
            force_magnitude = np.linalg.norm(total_force)
            if force_magnitude > self.force_ceiling:
                total_force = total_force * (self.force_ceiling / force_magnitude)
                
            return total_force
            
        except Exception as e:
            print(f"Error in pre-emergency containment: {e}")
            # Safe default: just apply moderate braking
            return -particle['velocity'] * 2.0
    
    def calculate_system_energy(self):
        """
        Calculate total system energy for conservation tracking.
        
        Computes the total energy of the system including kinetic energy
        of all particles, potential energy in the electromagnetic fields,
        and the field energy itself. Used for conservation law validation.
        
        Returns:
            float: Total system energy in Joules
        
        Energy components:
        - Kinetic energy: Sum of (1/2)mv² for all particles
        - Potential energy: Particle potential energy in fields
        - Magnetic field energy: B²/(2μ₀) integrated over volume
        - Electric field energy: ε₀E²/2 integrated over volume
        
        All energy components are tracked separately for detailed analysis.
        """
        """Calculate total system energy (kinetic + potential + field)."""
        # Kinetic energy
        kinetic_energy = sum(particle['kinetic_energy'] for particle in self.particles)
        
        # Potential energy
        potential_energy = sum(particle['potential_energy'] for particle in self.particles)
        
        # Field energy (simplified model)
        # B-field energy density: u_B = B^2/(2*mu_0)
        # E-field energy density: u_E = epsilon_0*E^2/2
        # Total field energy: integrate u_B + u_E over volume
        # This is an approximation using average field values
        volume = (4/3) * np.pi * self.containment_radius**3
        b_squared = self.magnetic_field_strength**2
        e_squared = (self.electric_field_gradient * self.containment_radius)**2
        
        magnetic_field_energy = b_squared / (2 * self.vacuum_permeability) * volume
        electric_field_energy = self.vacuum_permittivity * e_squared / 2 * volume
        
        field_energy = magnetic_field_energy + electric_field_energy
        
        # Store separate components
        self.kinetic_energy_history.append(kinetic_energy)
        self.potential_energy_history.append(potential_energy)
        self.field_energy_history.append(field_energy)
        
        return kinetic_energy + potential_energy + field_energy
    
    def calculate_total_angular_momentum(self):
        """
        Calculate total angular momentum of all particles.
        
        Computes the vector sum of angular momentum (L = r × p) for all
        particles in the system. Used for angular momentum conservation
        validation in the electromagnetic confinement.
        
        Returns:
            numpy.ndarray: Total angular momentum vector [Lx, Ly, Lz] in kg⋅m²/s
        
        Angular momentum L = r × p where:
        - r is position vector from origin
        - p is momentum vector (mv)
        - × denotes vector cross product
        """
        """Calculate total angular momentum of the system."""
        total_angular_momentum = np.zeros(3)
        
        for particle in self.particles:
            total_angular_momentum += particle['angular_momentum']
            
        return total_angular_momentum
        
    def determine_time_step(self):
        """
        Calculate optimal time step based on physical stability criteria.
        
        Uses multiple physical timescales to determine a numerically stable
        time step that maintains accuracy while preventing numerical instabilities.
        
        Returns:
            float: Optimal time step in seconds
        
        Time step criteria:
        - Courant-Friedrichs-Lewy condition: dt ≤ C⋅dx/v_max
        - Cyclotron period resolution: dt ≤ T_cyclotron/20
        - Maximum velocity constraint
        - Minimum time step enforcement
        
        The algorithm selects the most restrictive criterion to ensure stability.
        """
        """Determine time step based on physical considerations including Courant condition."""
        if not self.adaptive_time_stepping:
            return self.base_dt
            
        # Find maximum particle velocity
        max_velocity = 0
        for particle in self.particles:
            max_velocity = max(max_velocity, np.linalg.norm(particle['velocity']))
            
        if max_velocity < 1e-10:  # Avoid division by zero
            max_velocity = 1e-10
            
        # Courant-Friedrichs-Lewy condition for stability:
        # dt <= courant_factor * dx / max_velocity
        # where dx is the characteristic spatial resolution
        # For our purposes, we'll use containment_radius/100 as a rough estimate
        characteristic_length = self.containment_radius / 100
        courant_dt = self.courant_factor * characteristic_length / max_velocity
        
        # Additional physical timescale: cyclotron period
        # T_cyclotron = 2π*m/(q*B)
        if self.magnetic_field_strength > 1e-10:  # Avoid division by zero
            cyclotron_period = 2 * np.pi * self.antiproton_mass / (self.elementary_charge * self.magnetic_field_strength)
            # Use a fraction of cyclotron period for accurate resolution
            cyclotron_dt = cyclotron_period / 20
        else:
            cyclotron_dt = self.base_dt
            
        # Take minimum of various physical timescales for stability
        dt = min(self.base_dt, courant_dt, cyclotron_dt)
        
        # Enforce minimum time step
        return max(dt, self.min_dt)
        
    def update_particle_boris(self, particle, dt):
        """
        Update particle dynamics using the Boris algorithm.
        
        Implements the Boris algorithm, a specialized numerical method for
        integrating charged particle motion in electromagnetic fields.
        This method is particularly effective for maintaining energy conservation
        and handling the magnetic field rotation correctly.
        
        Args:
            particle (dict): Particle data structure to update
            dt (float): Integration time step
        
        Boris algorithm steps:
        1. Half electric field acceleration
        2. Magnetic field rotation using Boris rotation matrices
        3. Second half electric field acceleration
        4. Position update using final velocity
        
        The algorithm separates electric and magnetic field effects,
        allowing accurate treatment of the Lorentz force components.
        """
        """
        Update particle position and velocity using the Boris algorithm.
        This is a specialized method for integrating charged particles in electromagnetic fields.
        """
        try:
            # Get current position, velocity, and physical properties
            position = particle['position'].copy()
            velocity = particle['velocity'].copy()
            mass = particle['mass']
            charge = particle['charge']
            
            # Calculate the electromagnetic force at the current position
            force = self.calculate_lorentz_force(particle)
            
            # Extract E and B fields from the force
            # For a charged particle: F = q(E + v × B)
            # We need to extract E and B from the total force
            
            # Extract E field (approximation)
            # First create zero velocity particle to calculate pure E field force
            zero_vel_particle = {
                'position': position.copy(),
                'velocity': np.zeros(3),
                'mass': mass,
                'charge': charge,
                'potential_energy': 0
            }
            
            # Get force with zero velocity (only E field component)
            E_force = self.calculate_lorentz_force(zero_vel_particle)
            
            # E = F_E / q
            E_field = E_force / charge
            
            # Calculate B field from the difference
            # F_B = q(v × B)
            # F = F_E + F_B, so F_B = F - F_E
            B_force = force - E_force
            
            # For B field, use approximation: B = (F_B × v) / (q|v|²)
            v_squared = np.sum(velocity**2)
            if v_squared > 1e-20:  # Avoid division by zero
                v_cross_F = np.cross(velocity, B_force)
                B_field = v_cross_F / (charge * v_squared)
            else:
                # If velocity is near zero, use default B field orientation
                B_field = np.array([0, 0, self.magnetic_field_strength])
            
            # Now implement the Boris algorithm
            
            # Half acceleration from electric field
            v_minus = velocity + (charge * E_field / mass) * (dt/2)
            
            # Calculate rotation factor for magnetic field
            q_dt_over_2m = (charge * dt) / (2 * mass)
            B_factor = B_field * q_dt_over_2m
            
            # Calculate auxiliary vector
            t = B_factor
            t_squared = np.sum(t**2)
            
            if t_squared > 1e-20:  # Check if magnetic field is significant
                s = 2*t / (1 + t_squared)
                
                # Boris rotation algorithm
                v_prime = v_minus + np.cross(v_minus, t)
                v_plus = v_minus + np.cross(v_prime, s)
            else:
                # No magnetic field, just pass through
                v_plus = v_minus
            
            # Half acceleration from electric field again
            new_velocity = v_plus + (charge * E_field / mass) * (dt/2)
            
            # Apply velocity ceiling for stability
            velocity_magnitude = np.linalg.norm(new_velocity)
            if velocity_magnitude > self.velocity_ceiling:
                new_velocity = new_velocity * (self.velocity_ceiling / velocity_magnitude)
            
            # Update position using new velocity
            new_position = position + new_velocity * dt
            
            # Record acceleration for analysis
            acceleration = (new_velocity - velocity) / dt
            particle['acceleration_history'].append(acceleration.copy())
            particle['last_acceleration'] = acceleration.copy()
            
            # Apply velocity and position updates
            particle['velocity'] = new_velocity
            particle['position'] = new_position
            
            # Update energy and momentum (for conservation tracking)
            particle['kinetic_energy'] = 0.5 * mass * np.sum(new_velocity**2)
            particle['angular_momentum'] = np.cross(new_position, mass * new_velocity)
            
            # Calculate potential energy (for conservation tracking)
            # This is handled in calculate_lorentz_force but we need to update with new position
            temp_particle = {
                'position': new_position,
                'velocity': new_velocity,
                'mass': mass,
                'charge': charge,
                'potential_energy': 0
            }
            _ = self.calculate_lorentz_force(temp_particle)
            particle['potential_energy'] = temp_particle['potential_energy']
            
        except Exception as e:
            print(f"Error in Boris algorithm: {e}")
            # Fall back to simple Euler method if Boris fails
            force = self.calculate_lorentz_force(particle)
            acceleration = force / particle['mass']
            
            # Store acceleration
            particle['acceleration_history'].append(acceleration.copy())
            particle['last_acceleration'] = acceleration.copy()
            
            # Update velocity and position (simple Euler)
            particle['velocity'] += acceleration * dt
            particle['position'] += particle['velocity'] * dt
            
            # Update energy and momentum (for conservation tracking)
            particle['kinetic_energy'] = 0.5 * particle['mass'] * np.sum(particle['velocity']**2)
            particle['angular_momentum'] = np.cross(particle['position'], particle['mass'] * particle['velocity'])
    
    def update_particle_rk4(self, particle, dt):
        """
        Update particle dynamics using 4th-order Runge-Kutta integration.
        
        Implements the classical RK4 method for solving the differential
        equation system governing charged particle motion. Provides higher
        accuracy than simple Euler integration at the cost of more force
        evaluations per time step.
        
        Args:
            particle (dict): Particle data structure to update
            dt (float): Integration time step
        
        RK4 method:
        - Evaluates derivatives at 4 points within the time step
        - Combines evaluations with weighted average: (k1 + 2k2 + 2k3 + k4)/6
        - Provides 4th-order accuracy in time step size
        - Better error control than lower-order methods
        """
        """
        Update particle position and velocity using 4th-order Runge-Kutta integration.
        This provides better accuracy than simple Euler integration.
        """
        try:
            # Get current state
            x0 = particle['position'].copy()
            v0 = particle['velocity'].copy()
            m = particle['mass']
            
            # Define function for position derivative
            def x_derivative(x, v, t):
                return v
            
            # Define function for velocity derivative
            def v_derivative(x, v, t):
                # Create temporary particle at given position/velocity
                temp_particle = {
                    'position': x.copy(),
                    'velocity': v.copy(),
                    'mass': particle['mass'],
                    'charge': particle['charge'],
                    'potential_energy': 0
                }
                # Calculate force at this state
                force = self.calculate_lorentz_force(temp_particle)
                # F = ma, so a = F/m
                return force / m
            
            # RK4 integration for position
            k1x = x_derivative(x0, v0, 0) * dt
            k1v = v_derivative(x0, v0, 0) * dt
            
            k2x = x_derivative(x0 + k1x/2, v0 + k1v/2, dt/2) * dt
            k2v = v_derivative(x0 + k1x/2, v0 + k1v/2, dt/2) * dt
            
            k3x = x_derivative(x0 + k2x/2, v0 + k2v/2, dt/2) * dt
            k3v = v_derivative(x0 + k2x/2, v0 + k2v/2, dt/2) * dt
            
            k4x = x_derivative(x0 + k3x, v0 + k3v, dt) * dt
            k4v = v_derivative(x0 + k3x, v0 + k3v, dt) * dt
            
            # Calculate new position and velocity using weighted average of intermediate steps
            new_position = x0 + (k1x + 2*k2x + 2*k3x + k4x) / 6
            new_velocity = v0 + (k1v + 2*k2v + 2*k3v + k4v) / 6
            
            # Apply velocity ceiling for stability
            velocity_magnitude = np.linalg.norm(new_velocity)
            if velocity_magnitude > self.velocity_ceiling:
                new_velocity = new_velocity * (self.velocity_ceiling / velocity_magnitude)
            
            # Calculate acceleration for analysis
            acceleration = (new_velocity - v0) / dt
            particle['acceleration_history'].append(acceleration.copy())
            particle['last_acceleration'] = acceleration.copy()
            
            # Update particle state
            particle['position'] = new_position
            particle['velocity'] = new_velocity
            
            # Update energy and momentum (for conservation tracking)
            particle['kinetic_energy'] = 0.5 * m * np.sum(new_velocity**2)
            particle['angular_momentum'] = np.cross(new_position, m * new_velocity)
            
            # Calculate potential energy (for conservation tracking)
            temp_particle = {
                'position': new_position,
                'velocity': new_velocity,
                'mass': m,
                'charge': particle['charge'],
                'potential_energy': 0
            }
            _ = self.calculate_lorentz_force(temp_particle)
            particle['potential_energy'] = temp_particle['potential_energy']
            
        except Exception as e:
            print(f"Error in RK4 integration: {e}")
            # Fall back to simple Euler method if RK4 fails
            force = self.calculate_lorentz_force(particle)
            acceleration = force / particle['mass']
            
            # Store acceleration
            particle['acceleration_history'].append(acceleration.copy())
            particle['last_acceleration'] = acceleration.copy()
            
            # Update velocity and position
            particle['velocity'] += acceleration * dt
            particle['position'] += particle['velocity'] * dt
            
            # Update energy and momentum
            particle['kinetic_energy'] = 0.5 * particle['mass'] * np.sum(particle['velocity']**2)
            particle['angular_momentum'] = np.cross(particle['position'], particle['mass'] * particle['velocity'])
    
    def apply_physical_perturbation(self, perturbation_type):
        """
        Apply realistic environmental perturbations to the system.
        
        Simulates real-world disturbances that would affect a physical
        containment system, including magnetic field fluctuations,
        thermal noise, and cosmic ray interactions.
        
        Args:
            perturbation_type (str): Type of perturbation to apply
                - "magnetic_fluctuation": Environmental magnetic field variations
                - "thermal_noise": Brownian motion from thermal energy
                - "cosmic_ray": High-energy particle interactions
        
        Perturbation modeling:
        - Magnetic fluctuations: Random field variations from power systems
        - Thermal noise: Maxwell-Boltzmann distributed velocity kicks
        - Cosmic rays: Exponentially distributed energy transfers
        
        All perturbations are scaled to realistic magnitudes based on
        laboratory and environmental conditions.
        """
        """Apply realistic physical perturbations to the containment system."""
        try:
            if perturbation_type == "magnetic_fluctuation":
                # Simulate environmental magnetic field fluctuations
                fluctuation_strength = np.random.uniform(
                    self.earth_magnetic_field_fluctuation, 
                    self.power_grid_fluctuation
                )
                
                # Random direction for perturbation
                perturbation_dir = np.random.rand(3) - 0.5
                perturbation_dir = perturbation_dir / np.linalg.norm(perturbation_dir)
                
                # Apply small random fluctuation to each particle
                for particle in self.particles:
                    # Magnetic perturbation force: F = q(v × B')
                    velocity = particle['velocity']
                    charge = particle['charge']
                    
                    perturbation_field = perturbation_dir * fluctuation_strength
                    perturbation_force = charge * np.cross(velocity, perturbation_field)
                    
                    # Apply small impulse
                    acceleration = perturbation_force / particle['mass']
                    particle['velocity'] += acceleration * self.base_dt
                    particle['perturbed'] = True
                    
            elif perturbation_type == "thermal_noise":
                # Simulate thermal noise (Brownian motion)
                # Maxwell-Boltzmann distribution for temperature
                for particle in self.particles:
                    # Calculate thermal velocity scale
                    thermal_velocity = np.sqrt(2 * self.thermal_noise_factor / particle['mass'])
                    
                    # Generate random thermal kick
                    thermal_kick = np.random.randn(3) * thermal_velocity * 0.01
                    
                    # Apply small thermal impulse
                    particle['velocity'] += thermal_kick
                    particle['perturbed'] = True
                    
            elif perturbation_type == "cosmic_ray":
                # Simulate cosmic ray interaction with random particle
                if self.particles:
                    # Randomly select one particle to be affected
                    particle_idx = np.random.randint(0, len(self.particles))
                    particle = self.particles[particle_idx]
                    
                    # Random energy transfer from cosmic ray
                    energy_transfer = np.random.exponential(1e-18)  # Joules
                    
                    # Convert to velocity impulse: K = 1/2 mv²
                    velocity_impulse = np.sqrt(2 * energy_transfer / particle['mass'])
                    
                    # Random direction
                    impulse_dir = np.random.rand(3) - 0.5
                    impulse_dir = impulse_dir / np.linalg.norm(impulse_dir)
                    
                    # Apply impulse
                    particle['velocity'] += impulse_dir * velocity_impulse
                    particle['perturbed'] = True
                    
            # Record perturbation in history
            self.perturbation_history.append({
                'time': len(self.time_history) - 1 if self.time_history else 0,
                'type': perturbation_type
            })
            
        except Exception as e:
            print(f"Error applying physical perturbation: {e}")
    
    def run_simulation(self, duration):
        """
        Execute the main simulation loop.
        
        Runs the electromagnetic containment simulation for the specified
        duration, updating all particles at each time step and tracking
        conservation laws and system stability.
        
        Args:
            duration (float): Simulation duration in seconds
        
        Simulation loop:
        1. Determine optimal time step based on physical criteria
        2. Store current particle states in history buffers
        3. Apply random environmental perturbations (low probability)
        4. Update particle positions/velocities using chosen integration method
        5. Calculate and track energy and angular momentum conservation
        6. Check for containment breach conditions
        7. Repeat until duration is reached or containment fails
        
        The simulation maintains detailed histories for post-analysis
        and conservation law validation.
        """
        """Run the simulation with physical conservation tracking."""
        time = 0 if not self.time_history else self.time_history[-1]
        end_time = time + duration
        
        # If this is the first run, add initial time point
        if not self.time_history:
            self.time_history.append(time)
            
            # Calculate initial energy
            self.total_energy_initial = self.calculate_system_energy()
            self.total_energy_history.append(self.total_energy_initial)
            
            # Calculate initial angular momentum
            self.total_angular_momentum_initial = self.calculate_total_angular_momentum()
            self.angular_momentum_history.append(self.total_angular_momentum_initial)
            
            # Initialize error tracking
            self.energy_conservation_error.append(0.0)
            self.angular_momentum_error.append(0.0)
        
        # Run simulation until duration is reached
        while time < end_time:
            # Determine appropriate time step based on physical considerations
            dt = self.determine_time_step()
            
            # Don't exceed remaining duration
            if time + dt > end_time:
                dt = end_time - time
            
            # Store current state in history buffer for each particle
            for particle in self.particles:
                particle['history'].append({
                    'position': particle['position'].copy(),
                    'velocity': particle['velocity'].copy(),
                    'time': time
                })
            
            # Apply any physical perturbations (random occurrences)
            if self.realistic_perturbations and np.random.random() < 0.001:  # Low probability of perturbation
                perturbation_types = ["magnetic_fluctuation", "thermal_noise", "cosmic_ray"]
                perturbation_type = np.random.choice(perturbation_types)
                self.apply_physical_perturbation(perturbation_type)
            
            # Update each particle with appropriate integration method
            for particle in self.particles:
                if self.boris_algorithm:
                    self.update_particle_boris(particle, dt)
                elif self.runge_kutta_order == 4:
                    self.update_particle_rk4(particle, dt)
                else:
                    # Simple Euler method as fallback
                    force = self.calculate_lorentz_force(particle)
                    acceleration = force / particle['mass']
                    
                    # Store acceleration
                    particle['acceleration_history'].append(acceleration.copy())
                    particle['last_acceleration'] = acceleration.copy()
                    
                    # Update velocity and position
                    particle['velocity'] += acceleration * dt
                    particle['position'] += particle['velocity'] * dt
                    
                    # Update energy and momentum
                    particle['kinetic_energy'] = 0.5 * particle['mass'] * np.sum(particle['velocity']**2)
                    particle['angular_momentum'] = np.cross(particle['position'], particle['mass'] * particle['velocity'])
            
            # Calculate system energy
            total_energy = self.calculate_system_energy()
            self.total_energy_history.append(total_energy)
            
            # Calculate energy conservation error
            if self.total_energy_initial != 0:
                energy_error = abs((total_energy - self.total_energy_initial) / self.total_energy_initial)
            else:
                energy_error = 0
            self.energy_conservation_error.append(energy_error)
            
            # Calculate angular momentum
            if self.angular_momentum_tracking:
                total_angular_momentum = self.calculate_total_angular_momentum()
                self.angular_momentum_history.append(total_angular_momentum)
                
                # Calculate angular momentum conservation error
                if np.linalg.norm(self.total_angular_momentum_initial) > 1e-20:
                    am_error = np.linalg.norm(total_angular_momentum - self.total_angular_momentum_initial) / np.linalg.norm(self.total_angular_momentum_initial)
                else:
                    am_error = np.linalg.norm(total_angular_momentum)
                self.angular_momentum_error.append(am_error)
            
            # Record history
            time += dt
            self.time_history.append(time)
            
            # Check for containment breach
            max_distance = 0
            for particle in self.particles:
                distance = np.linalg.norm(particle['position'])
                max_distance = max(max_distance, distance)
                
            if max_distance > self.containment_radius * 10:
                print(f"Containment breach at t={time:.4e}s")
                break
    
    def visualize_results(self, save_only=False):
        """
        Generate comprehensive visualization of simulation results.
        
        Creates multi-panel plots showing particle trajectories,
        conservation law tracking, energy distribution, and containment
        effectiveness. Supports both interactive display and file output.
        
        Args:
            save_only (bool): If True, save to file without displaying
        
        Visualization panels:
        1. 3D particle trajectories and final positions
        2. Energy conservation error over time
        3. Angular momentum conservation error
        4. Energy distribution (kinetic, potential, field)
        5. Containment effectiveness (max particle distance)
        6. Particle velocity evolution
        
        All plots include appropriate axis labels, legends, and grid lines
        for scientific presentation quality.
        """
        """Visualize the simulation results with physical validation metrics."""
        try:
            # Check if running in non-interactive environment
            if save_only or not plt.isinteractive():
                # Use non-interactive backend if needed
                plt.switch_backend('Agg')

            # Create enhanced figure with FCE diagnostics (8-panel layout)
            fig = plt.figure(figsize=(20, 24))
            
            # 3D particle trajectories
            ax1 = fig.add_subplot(421, projection='3d')
            for i, particle in enumerate(self.particles):
                x, y, z = particle['position']
                ax1.scatter(x, y, z, c='red', marker='o')
                
                # Plot trajectory if available
                if particle['history']:
                    trajectory = np.array([state['position'] for state in particle['history']])
                    ax1.plot(trajectory[:,0], trajectory[:,1], trajectory[:,2], 'r-', alpha=0.3)
                    
            ax1.set_xlabel('X')
            ax1.set_ylabel('Y')
            ax1.set_zlabel('Z')
            ax1.set_title('Antiparticle Positions and Trajectories')
            
            # Energy conservation plot
            ax2 = fig.add_subplot(422)
            if self.energy_conservation and len(self.energy_conservation_error) > 1:
                ax2.plot(self.time_history[:len(self.energy_conservation_error)], 
                         self.energy_conservation_error, label='Energy Error')
                ax2.set_xlabel('Time (s)')
                ax2.set_ylabel('Relative Error')
                ax2.set_title('Energy Conservation Error')
                ax2.set_yscale('log')
                ax2.grid(True)
            else:
                ax2.text(0.5, 0.5, "Energy conservation tracking not enabled",
                        horizontalalignment='center', verticalalignment='center')
            
            # Angular momentum conservation
            ax3 = fig.add_subplot(423)
            if self.angular_momentum_tracking and len(self.angular_momentum_error) > 1:
                ax3.plot(self.time_history[:len(self.angular_momentum_error)], 
                         self.angular_momentum_error, label='Angular Momentum Error')
                ax3.set_xlabel('Time (s)')
                ax3.set_ylabel('Relative Error')
                ax3.set_title('Angular Momentum Conservation Error')
                ax3.set_yscale('log')
                ax3.grid(True)
            else:
                ax3.text(0.5, 0.5, "Angular momentum tracking not enabled",
                        horizontalalignment='center', verticalalignment='center')
            
            # Energy distribution
            ax4 = fig.add_subplot(424)
            if len(self.kinetic_energy_history) > 0 and len(self.potential_energy_history) > 0:
                ax4.plot(self.time_history[:len(self.kinetic_energy_history)], 
                         self.kinetic_energy_history, label='Kinetic')
                ax4.plot(self.time_history[:len(self.potential_energy_history)], 
                         self.potential_energy_history, label='Potential')
                ax4.plot(self.time_history[:len(self.field_energy_history)], 
                         self.field_energy_history, label='Field')
                ax4.set_xlabel('Time (s)')
                ax4.set_ylabel('Energy (J)')
                ax4.set_title('Energy Distribution')
                ax4.legend()
                ax4.grid(True)
            else:
                ax4.text(0.5, 0.5, "Energy distribution data not available",
                        horizontalalignment='center', verticalalignment='center')
            
            # Containment effectiveness visualization
            ax5 = fig.add_subplot(425)
            max_distances = []
            for t_idx, t in enumerate(self.time_history):
                max_dist = 0
                for particle in self.particles:
                    if t_idx < len(particle['history']):
                        pos = particle['history'][t_idx]['position']
                        dist = np.linalg.norm(pos)
                        max_dist = max(max_dist, dist)
                max_distances.append(max_dist)
                
            if max_distances:
                ax5.plot(self.time_history[:len(max_distances)], max_distances)
                ax5.axhline(y=self.containment_radius, color='r', linestyle='--', label='Containment Radius')
                ax5.set_xlabel('Time (s)')
                ax5.set_ylabel('Maximum Particle Distance (m)')
                ax5.set_title('Containment Effectiveness')
                ax5.grid(True)
            else:
                ax5.text(0.5, 0.5, "Particle distance data not available",
                        horizontalalignment='center', verticalalignment='center')
            
            # Velocity distribution
            ax6 = fig.add_subplot(426)
            if self.particles and self.particles[0]['history']:
                times = np.array(self.time_history[:len(self.particles[0]['history'])])
                vels = []
                for particle in self.particles:
                    if particle['history']:
                        v = np.array([np.linalg.norm(state['velocity']) for state in particle['history']])
                        if len(v) == len(times):
                            vels.append(v)
                
                if vels:
                    mean_vel = np.mean(vels, axis=0)
                    ax6.plot(times, mean_vel)
                    ax6.set_xlabel('Time (s)')
                    ax6.set_ylabel('Average Velocity (m/s)')
                    ax6.set_title('Particle Velocity Evolution')
                    ax6.grid(True)
                else:
                    ax6.text(0.5, 0.5, "Velocity data not available",
                            horizontalalignment='center', verticalalignment='center')
            else:
                ax6.text(0.5, 0.5, "Velocity data not available",
                        horizontalalignment='center', verticalalignment='center')

            # FCE Effectiveness Analysis (New Panel 7)
            ax7 = fig.add_subplot(427)
            if self.particles and self.particles[0]['instability_history']:
                # Calculate FCE correction strength over time
                fce_corrections = []
                for t_idx, t in enumerate(self.time_history):
                    total_correction = 0
                    particle_count = 0
                    for particle in self.particles:
                        if t_idx < len(particle['instability_history']):
                            instability = particle['instability_history'][t_idx]
                            if instability > self.stability_threshold:
                                correction_force = min(self.fractal_correction_strength * instability,
                                                     self.max_correction_force)
                                total_correction += correction_force
                                particle_count += 1
                    avg_correction = total_correction / max(1, particle_count)
                    fce_corrections.append(avg_correction)

                if fce_corrections:
                    ax7.plot(self.time_history[:len(fce_corrections)],
                            np.array(fce_corrections) * 1e15, 'purple', linewidth=2)
                    ax7.set_xlabel('Time (s)')
                    ax7.set_ylabel('FCE Correction Force (fN)')
                    ax7.set_title('Fractal Correction Engine Effectiveness')
                    ax7.grid(True, alpha=0.3)
                    ax7.set_yscale('log')
                else:
                    ax7.text(0.5, 0.5, "FCE operating in stable regime\n(No corrections needed)",
                            horizontalalignment='center', verticalalignment='center',
                            bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
            else:
                ax7.text(0.5, 0.5, "FCE data not available",
                        horizontalalignment='center', verticalalignment='center')

            # Fractal Dimension & Phase Space Analysis (New Panel 8)
            ax8 = fig.add_subplot(428)
            if self.particles and self.particles[0]['history']:
                # Calculate fractal dimension using box-counting method
                try:
                    # Collect all trajectory points
                    all_points = []
                    for particle in self.particles:
                        if particle['history']:
                            trajectory = np.array([state['position'] for state in particle['history']])
                            all_points.extend(trajectory)

                    if len(all_points) > 10:
                        all_points = np.array(all_points)

                        # Simple fractal dimension estimate using box-counting
                        # This is a simplified version for demonstration
                        scales = np.logspace(-3, -1, 10)  # Box sizes from 1mm to 10cm
                        counts = []

                        for scale in scales:
                            # Count occupied boxes at this scale
                            boxes = set()
                            for point in all_points:
                                box_coords = tuple(np.floor(point / scale).astype(int))
                                boxes.add(box_coords)
                            counts.append(len(boxes))

                        # Fit log-log relationship
                        log_scales = np.log(scales)
                        log_counts = np.log(counts)

                        # Remove any invalid values
                        valid_idx = np.isfinite(log_scales) & np.isfinite(log_counts)
                        if np.sum(valid_idx) > 2:
                            log_scales = log_scales[valid_idx]
                            log_counts = log_counts[valid_idx]

                            # Linear fit to get fractal dimension
                            coeffs = np.polyfit(log_scales, log_counts, 1)
                            fractal_dim = -coeffs[0]  # Negative slope gives dimension

                            ax8.loglog(scales[valid_idx], np.array(counts)[valid_idx], 'bo-',
                                     label=f'Box Count Data')
                            ax8.loglog(scales[valid_idx], np.exp(np.polyval(coeffs, log_scales)),
                                     'r--', label=f'Fit: D = {fractal_dim:.2f}')
                            ax8.set_xlabel('Box Size (m)')
                            ax8.set_ylabel('Number of Occupied Boxes')
                            ax8.set_title(f'Phase Space Fractal Analysis\nDimension = {fractal_dim:.2f}')
                            ax8.legend()
                            ax8.grid(True, alpha=0.3)
                        else:
                            ax8.text(0.5, 0.5, "Insufficient data for fractal analysis",
                                    horizontalalignment='center', verticalalignment='center')
                    else:
                        ax8.text(0.5, 0.5, "Limited trajectory data\nFractal analysis not reliable",
                                horizontalalignment='center', verticalalignment='center')

                except Exception as e:
                    ax8.text(0.5, 0.5, f"Fractal analysis error:\n{str(e)[:50]}...",
                            horizontalalignment='center', verticalalignment='center',
                            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
            else:
                ax8.text(0.5, 0.5, "No trajectory data for fractal analysis",
                        horizontalalignment='center', verticalalignment='center')

            plt.tight_layout()
            
            # Save if output directory is specified
            if self.output_dir:
                os.makedirs(self.output_dir, exist_ok=True)  # Ensure directory exists
                plt.savefig(os.path.join(self.output_dir, "enhanced_fce_validation.png"), dpi=300)
            
            # Show unless save_only is specified
            if not save_only:
                try:
                    plt.show()
                except Exception as e:
                    print(f"Unable to display plot interactively: {e}")
            
            plt.close(fig)
                
        except Exception as e:
            print(f"Error in visualization: {e}")
            traceback.print_exc()
    
    def save_data(self):
        """
        Save simulation data to CSV files for analysis.
        
        Exports all tracking data to structured CSV files suitable for
        post-processing analysis, statistical evaluation, and external
        validation of conservation laws.
        
        Saved data files:
        - energy_conservation.csv: Time, total energy, conservation errors
        - angular_momentum.csv: Time, L components, conservation errors
        - final_positions.csv: Final particle states and properties
        - physical_config.txt: Simulation configuration parameters
        
        All CSV files include headers and use standard comma separation
        for compatibility with analysis software.
        """
        """Save simulation data to output directory with physical metrics"""
        if not self.output_dir:
            print("No output directory specified. Cannot save data.")
            return
            
        try:
            # Ensure output directory exists
            os.makedirs(self.output_dir, exist_ok=True)
            
            # Save time history and energy conservation data
            if self.energy_conservation and len(self.energy_conservation_error) > 0:
                energy_data = np.column_stack((
                    self.time_history[:len(self.energy_conservation_error)],
                    self.total_energy_history,
                    self.energy_conservation_error,
                    self.kinetic_energy_history,
                    self.potential_energy_history,
                    self.field_energy_history
                ))
                np.savetxt(os.path.join(self.output_dir, "energy_conservation.csv"), energy_data, 
                          delimiter=',', header="time,total_energy,energy_error,kinetic_energy,potential_energy,field_energy", comments='')
            
            # Save angular momentum data
            if self.angular_momentum_tracking and len(self.angular_momentum_error) > 0:
                # Convert 3D vectors to columns
                am_history = np.array(self.angular_momentum_history)
                am_data = np.column_stack((
                    self.time_history[:len(self.angular_momentum_error)],
                    am_history[:, 0],
                    am_history[:, 1],
                    am_history[:, 2],
                    self.angular_momentum_error
                ))
                np.savetxt(os.path.join(self.output_dir, "angular_momentum.csv"), am_data, 
                          delimiter=',', header="time,L_x,L_y,L_z,am_error", comments='')
            
            # Save final particle positions
            with open(os.path.join(self.output_dir, "final_positions.csv"), 'w') as f:
                f.write("particle_id,x,y,z,vx,vy,vz,kinetic_energy,lx,ly,lz\n")
                for i, particle in enumerate(self.particles):
                    pos = particle['position']
                    vel = particle['velocity']
                    ke = particle['kinetic_energy']
                    am = particle['angular_momentum']
                    f.write(f"{i},{pos[0]},{pos[1]},{pos[2]},{vel[0]},{vel[1]},{vel[2]},{ke},{am[0]},{am[1]},{am[2]}\n")
            
            # Save configuration
            with open(os.path.join(self.output_dir, "physical_config.txt"), 'w') as f:
                f.write("Physical Simulation Configuration:\n")
                f.write(f"Relativistic Corrections: {self.relativistic_corrections}\n")
                f.write(f"Symplectic Integration: {self.symplectic_integration}\n")
                f.write(f"Boris Algorithm: {self.boris_algorithm}\n")
                f.write(f"Energy Conservation Tracking: {self.energy_conservation}\n")
                f.write(f"Angular Momentum Tracking: {self.angular_momentum_tracking}\n")
                f.write(f"Adaptive Time Stepping: {self.adaptive_time_stepping}\n")
                f.write(f"Courant Factor: {self.courant_factor}\n")
                f.write(f"Multi-Layer Containment: {self.multi_layer_containment}\n")
                f.write(f"Realistic Perturbations: {self.realistic_perturbations}\n")
                
            print(f"Data saved to {self.output_dir}")
            
        except Exception as e:
            print(f"Error saving data: {e}")
            traceback.print_exc()


def create_output_directory(test_name):
    """
    Create timestamped output directory for simulation results.
    
    Generates a unique directory name using timestamp and test name
    to organize simulation outputs and prevent overwrites of previous
    results.
    
    Args:
        test_name (str): Descriptive name for the test being run
    
    Returns:
        str: Path to created output directory
    
    Directory structure:
    - Base directory: "antimatter_physical_results"
    - Subdirectories: "{timestamp}_{test_name}"
    - Fallback: Simple timestamp-based name if creation fails
    """
    """Create a unique output directory for the test results"""
    try:
        # Create base directory for all test results if it doesn't exist
        base_dir = "antimatter_physical_results"
        if not os.path.exists(base_dir):
            os.makedirs(base_dir)
        
        # Create a unique directory name with timestamp and test name
        timestamp = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")
        test_dir_name = f"{timestamp}_{test_name.replace(' ', '_')}"
        output_dir = os.path.join(base_dir, test_dir_name)
        
        # Create the directory
        os.makedirs(output_dir, exist_ok=True)
        
        return output_dir
    except Exception as e:
        print(f"Error creating output directory: {e}")
        # Fallback to a simple directory if there's an error
        fallback_dir = f"antimatter_results_{int(time.time())}"
        os.makedirs(fallback_dir, exist_ok=True)
        return fallback_dir


def run_physical_simulation(test_name, particle_count=5, duration=1e-8, 
                          containment_radius=0.1, enable_relativistic=True,
                          integration_method="boris"):
    """
    Execute a complete electromagnetic containment simulation test.
    
    Runs a full simulation with specified parameters, including particle
    initialization, physics integration, conservation tracking, and results
    analysis. Provides comprehensive output suitable for scientific validation.
    
    Args:
        test_name (str): Descriptive name for the simulation test
        particle_count (int): Number of particles to simulate (default: 5)
        duration (float): Simulation time in seconds (default: 1e-8)
        containment_radius (float): Containment boundary in meters (default: 0.1)
        enable_relativistic (bool): Enable relativistic corrections (default: True)
        integration_method (str): "boris" or "rk4" integration (default: "boris")
    
    Returns:
        tuple: (simulator_object, result_string, output_directory_path)
    
    Test procedure:
    1. Create timestamped output directory
    2. Initialize simulator with specified physics models
    3. Add particles with realistic initial conditions
    4. Execute simulation with conservation tracking
    5. Analyze results for energy/momentum conservation
    6. Generate comprehensive visualizations
    7. Save all data and configuration for reproducibility
    
    Results are validated against conservation laws with quantitative
    error analysis and containment effectiveness metrics.
    """
    """Run antimatter simulation with physical models and conservation tracking"""
    try:
        # Create output directory
        output_dir = create_output_directory(test_name)
        
        print(f"\n{'='*50}")
        print(f"STARTING PHYSICAL SIMULATION: {test_name}")
        print(f"{'='*50}")
        print(f"Parameters:")
        print(f"  - Particle Count: {particle_count}")
        print(f"  - Simulation Duration: {duration:.1e} seconds")
        print(f"  - Containment Radius: {containment_radius:.3f} meters")
        print(f"  - Relativistic Effects: {enable_relativistic}")
        print(f"  - Integration Method: {integration_method}")
        print(f"  - Output Directory: {output_dir}")
        
        # Save test configuration
        with open(os.path.join(output_dir, "test_config.txt"), 'w') as f:
            f.write(f"Test Name: {test_name}\n")
            f.write(f"Date: {datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n")
            f.write(f"Particle Count: {particle_count}\n")
            f.write(f"Simulation Duration: {duration:.1e} seconds\n")
            f.write(f"Containment Radius: {containment_radius:.3f} meters\n")
            f.write(f"Relativistic Effects: {enable_relativistic}\n")
            f.write(f"Integration Method: {integration_method}\n")
        
        # Create simulation with physical features
        sim = PhysicallyAccurateAntimatterSimulator(output_dir=output_dir)
        sim.containment_radius = containment_radius
        
        # Configure physics models
        sim.relativistic_corrections = enable_relativistic
        sim.energy_conservation = True
        sim.angular_momentum_tracking = True
        sim.realistic_perturbations = True
        
        # Set integration method
        if integration_method.lower() == "boris":
            sim.boris_algorithm = True
            sim.runge_kutta_order = 0
            sim.symplectic_integration = True
        elif integration_method.lower() == "rk4":
            sim.boris_algorithm = False
            sim.runge_kutta_order = 4
            sim.symplectic_integration = False
        
        # Add antiparticles with physically realistic initial conditions
        np.random.seed(42)  # For reproducibility
        for i in range(particle_count):
            # Random position within containment (starting closer to center)
            position = (np.random.rand(3) - 0.5) * containment_radius * 0.2
            
            # Random initial velocity - physically reasonable for contained particles
            # For antiprotons at room temperature, thermal velocity ~2000 m/s
            thermal_velocity = 2000  # m/s
            velocity = (np.random.rand(3) - 0.5) * thermal_velocity
            
            sim.add_antiparticle(position, velocity)
        
        print(f"Added {particle_count} antiparticles to the simulation")
        
        # Run simulation
        start_time = time.time()
        sim.run_simulation(duration=duration)
        sim_time = time.time() - start_time
        
        # Calculate physical metrics
        if sim.energy_conservation_error:
            avg_energy_error = np.mean(sim.energy_conservation_error)
            max_energy_error = np.max(sim.energy_conservation_error)
        else:
            avg_energy_error = 0
            max_energy_error = 0
            
        if sim.angular_momentum_error:
            avg_am_error = np.mean(sim.angular_momentum_error)
            max_am_error = np.max(sim.angular_momentum_error)
        else:
            avg_am_error = 0
            max_am_error = 0
        
        # Calculate confinement efficiency
        max_distance = 0
        for particle in sim.particles:
            distance = np.linalg.norm(particle['position'])
            max_distance = max(max_distance, distance)
        
        confinement_ratio = max_distance / containment_radius
        
        # Save results
        with open(os.path.join(output_dir, "physical_results.txt"), 'w') as f:
            f.write("Physical Simulation Results:\n")
            f.write(f"  - Runtime: {sim_time:.2f} seconds\n")
            f.write(f"  - Time Steps: {len(sim.time_history)}\n")
            f.write(f"  - Avg Energy Conservation Error: {avg_energy_error:.4e}\n")
            f.write(f"  - Max Energy Conservation Error: {max_energy_error:.4e}\n")
            f.write(f"  - Avg Angular Momentum Error: {avg_am_error:.4e}\n")
            f.write(f"  - Max Angular Momentum Error: {max_am_error:.4e}\n")
            f.write(f"  - Max Particle Distance: {max_distance:.6f} m\n")
            f.write(f"  - Confinement Ratio: {confinement_ratio:.2%} of radius\n")
        
        # Determine if test was successful based on physical metrics
        if max_energy_error < 0.01 and max_am_error < 0.01 and confinement_ratio < 0.9:
            result = "SUCCESS - Physical conservation laws maintained"
        elif max_energy_error >= 0.01:
            result = "PARTIAL SUCCESS - Energy conservation errors detected"
        elif max_am_error >= 0.01:
            result = "PARTIAL SUCCESS - Angular momentum conservation errors detected"
        elif confinement_ratio >= 0.9:
            result = "PARTIAL SUCCESS - Particles approaching containment boundary"
        
        # Update results file with test outcome
        with open(os.path.join(output_dir, "physical_results.txt"), 'a') as f:
            f.write(f"\nTEST RESULT: {result}\n")
        
        print("\nPhysical Simulation Results:")
        print(f"  - Runtime: {sim_time:.2f} seconds")
        print(f"  - Time Steps: {len(sim.time_history)}")
        print(f"  - Avg Energy Conservation Error: {avg_energy_error:.4e}")
        print(f"  - Max Energy Conservation Error: {max_energy_error:.4e}")
        print(f"  - Avg Angular Momentum Error: {avg_am_error:.4e}")
        print(f"  - Max Angular Momentum Error: {max_am_error:.4e}")
        print(f"  - Max Particle Distance: {max_distance:.6f} m")
        print(f"  - Confinement Ratio: {confinement_ratio:.2%} of radius")
        print(f"\nTEST RESULT: {result}")
        
        # Save data
        print("\nSaving simulation data...")
        sim.save_data()
        
        # Visualize results
        print("\nGenerating physical validation visualizations...")
        sim.visualize_results()
        sim.visualize_results(save_only=True)
        
        return sim, result, output_dir
        
    except Exception as e:
        print(f"Error in physical simulation: {e}")
        traceback.print_exc()
        return None, "ERROR", None


if __name__ == "__main__":
    try:
        print("Running ENHANCED PUBLICATION-READY antimatter containment simulations with FCE diagnostics...")
        
        # Test 1: Enhanced FCE Model with Boris Integration (Publication Ready)
        print("\n\n" + "="*80)
        print("TEST 1: Enhanced FCE Model with Boris Integration (Publication Ready)")
        print("="*80)
        test1_params = {
            "test_name": "Publication_FCE_Boris_Integration",
            "particle_count": 20,
            "duration": 1e-6,  # 1 microsecond for more data
            "containment_radius": 0.1,
            "enable_relativistic": True,
            "integration_method": "boris"
        }
        sim1, result1, output_dir1 = run_physical_simulation(**test1_params)
        print(f"\nTest 1 completed with result: {result1}")
        if output_dir1:
            print(f"Results saved to: {output_dir1}")
        
        # Test 2: Enhanced FCE Model with RK4 Integration (Comparison Study)
        print("\n\n" + "="*80)
        print("TEST 2: Enhanced FCE Model with RK4 Integration (Comparison Study)")
        print("="*80)
        test2_params = {
            "test_name": "Publication_FCE_RK4_Integration",
            "particle_count": 20,
            "duration": 1e-6,  # 1 microsecond for more data
            "containment_radius": 0.1,
            "enable_relativistic": True,
            "integration_method": "rk4"
        }
        sim2, result2, output_dir2 = run_physical_simulation(**test2_params)
        print(f"\nTest 2 completed with result: {result2}")
        if output_dir2:
            print(f"Results saved to: {output_dir2}")
        
        # Test 3: Enhanced FCE Model - Extended Duration Study (Research Grade)
        print("\n\n" + "="*80)
        print("TEST 3: Enhanced FCE Model - Extended Duration Study (Research Grade)")
        print("="*80)
        test3_params = {
            "test_name": "Publication_FCE_Extended_Duration",
            "particle_count": 20,
            "duration": 5e-6,  # 5 microseconds for comprehensive analysis
            "containment_radius": 0.1,
            "enable_relativistic": True,
            "integration_method": "boris"
        }
        sim3, result3, output_dir3 = run_physical_simulation(**test3_params)
        print(f"\nTest 3 completed with result: {result3}")
        if output_dir3:
            print(f"Results saved to: {output_dir3}")
        
        # Summary of all tests
        print("\n\n" + "="*80)
        print("SUMMARY OF PHYSICAL SIMULATION TESTS")
        print("="*80)
        print(f"Test 1 (Boris Integration): {result1}")
        print(f"Test 2 (RK4 Integration): {result2}")
        print(f"Test 3 (Non-Relativistic): {result3}")
            
    except Exception as e:
        print(f"Error in main execution: {e}")
        traceback.print_exc()