"""
Spacetime Bending Simulation and Visualization Tool (Version 1.0)
This script simulates the bending of spacetime around a central mass using a simplified
Newtonian gravitational potential model (V(r) = -GM/r). It generates a 2D grid representing
spacetime curvature, visualizes it as a 3D surface, simulates particle/light trajectories,
and exports data for further analysis. Designed for educational and illustrative purposes.

Features:
- Customizable mass, grid size, and scaling via command-line or interactive input
- 3D visualization with Matplotlib, including interactive rotation
- Particle trajectory simulation using basic numerical integration
- Data export to CSV for grid and trajectory data
- Error handling and logging for robust execution

Dependencies: numpy, matplotlib
Usage: Run directly for interactive mode or use command-line arguments for batch processing.
Example: python spacetime_bend.py --mass 1.989e30 --grid-size 100 --scale 1e-10
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import argparse
import logging
import csv
import sys
from typing import Tuple, List, Optional
import warnings

# Set up logging for debugging and tracking execution
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(levelname)s - %(message)s',
    filename='spacetime_bend.log',
    filemode='w'
)
logger = logging.getLogger(__name__)

class SpacetimeSimulator:
    """
    A class to simulate and visualize spacetime bending around a central mass.
    Handles grid creation, potential calculation, trajectory simulation, and data export.
    """
    def __init__(self, mass: float = 1.989e30, grid_size: int = 100, scale: float = 1e-10):
        """
        Initialize the simulator with physical and computational parameters.
        
        Args:
            mass (float): Central mass in kg (default: Sun's mass, 1.989e30 kg)
            grid_size (int): Number of points per axis for the grid (default: 100)
            scale (float): Scaling factor to amplify curvature for visibility (default: 1e-10)
        """
        self.G = 6.67430e-11  # Gravitational constant (m^3 kg^-1 s^-2)
        self.mass = mass
        self.grid_size = grid_size
        self.scale = scale
        self.range = 50.0  # Grid range in arbitrary units (-range to +range)
        self.X = None
        self.Y = None
        self.R = None
        self.Z = None
        logger.info(f"Initialized simulator with mass={self.mass} kg, grid_size={self.grid_size}, scale={self.scale}")

    def create_grid(self) -> None:
        """Create a 2D grid representing a slice of spacetime."""
        try:
            x = np.linspace(-self.range, self.range, self.grid_size)
            y = np.linspace(-self.range, self.range, self.grid_size)
            self.X, self.Y = np.meshgrid(x, y)
            self.R = np.sqrt(self.X**2 + self.Y**2)
            # Avoid division by zero at center
            self.R[self.R == 0] = 1e-10
            logger.info("Spacetime grid created successfully.")
        except MemoryError:
            logger.error("Memory error while creating grid. Reduce grid_size.")
            raise
        except Exception as e:
            logger.error(f"Error creating grid: {str(e)}")
            raise

    def calculate_potential(self) -> None:
        """Calculate gravitational potential V(r) = -GM/(r*scale) across the grid."""
        if self.R is None:
            raise ValueError("Grid not created. Call create_grid() first.")
        try:
            self.Z = - (self.G * self.mass) / (self.R * self.scale)
            logger.info("Gravitational potential calculated across grid.")
        except Exception as e:
            logger.error(f"Error calculating potential: {str(e)}")
            raise

    def visualize_3d(self, colormap: str = 'viridis', save_plot: bool = False) -> None:
        """
        Visualize spacetime curvature as a 3D surface plot.
        
        Args:
            colormap (str): Matplotlib colormap for surface (default: 'viridis')
            save_plot (bool): Save the plot to a file if True (default: False)
        """
        if self.Z is None:
            raise ValueError("Potential not calculated. Call calculate_potential() first.")
        try:
            fig = plt.figure(figsize=(12, 10))
            ax = fig.add_subplot(111, projection='3d')
            surf = ax.plot_surface(self.X, self.Y, self.Z, cmap=colormap, edgecolor='none')
            fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5, label='Spacetime Curvature (Scaled Potential)')
            ax.set_xlabel('X (arbitrary units)')
            ax.set_ylabel('Y (arbitrary units)')
            ax.set_zlabel('Curvature Depth (Scaled)')
            ax.set_title(f'Spacetime Bending by Mass {self.mass:.2e} kg')
            if save_plot:
                plt.savefig('spacetime_curvature.png', dpi=300, bbox_inches='tight')
                logger.info("Plot saved as 'spacetime_curvature.png'")
            plt.show()
            logger.info("3D visualization displayed.")
        except Exception as e:
            logger.error(f"Error in visualization: {str(e)}")
            raise

    def simulate_trajectory(self, initial_pos: Tuple[float, float] = (40.0, 0.0),
                            initial_vel: Tuple[float, float] = (0.0, 0.5),
                            steps: int = 1000, dt: float = 0.05,
                            save_data: bool = False) -> Tuple[np.ndarray, np.ndarray]:
        """
        Simulate a test particle trajectory in the gravitational well using Euler integration.
        This is a simplified model for illustration, not fully relativistic.
        
        Args:
            initial_pos (tuple): Starting (x, y) position (default: (40.0, 0.0))
            initial_vel (tuple): Starting (vx, vy) velocity (default: (0.0, 0.5))
            steps (int): Number of time steps (default: 1000)
            dt (float): Time step size (default: 0.05)
            save_data (bool): Save trajectory to CSV if True (default: False)
        
        Returns:
            Tuple of numpy arrays (x_traj, y_traj) representing the particle path
        """
        try:
            x, y = initial_pos
            vx, vy = initial_vel
            x_traj, y_traj = [x], [y]
            scale_force = self.scale * 1e-2  # Adjust force scaling for simulation
            for _ in range(steps):
                r = np.sqrt(x**2 + y**2)
                if r < 1e-5:  # Prevent extreme forces near center
                    break
                # Gravitational acceleration (simplified: a = -GM/r^2)
                accel = - (self.G * self.mass) / (r**3 * scale_force)
                ax = accel * x
                ay = accel * y
                # Update velocity and position (Euler method)
                vx += ax * dt
                vy += ay * dt
                x += vx * dt
                y += vy * dt
                x_traj.append(x)
                y_traj.append(y)
                # Stop if particle escapes grid
                if abs(x) > self.range * 1.5 or abs(y) > self.range * 1.5:
                    break
            x_traj = np.array(x_traj)
            y_traj = np.array(y_traj)
            if save_data:
                with open('trajectory_data.csv', 'w', newline='') as f:
                    writer = csv.writer(f)
                    writer.writerow(['X', 'Y'])
                    writer.writerows(zip(x_traj, y_traj))
                logger.info("Trajectory data saved to 'trajectory_data.csv'")
            logger.info("Particle trajectory simulation completed.")
            return x_traj, y_traj
        except Exception as e:
            logger.error(f"Error in trajectory simulation: {str(e)}")
            raise

    def visualize_trajectory(self, x_traj: np.ndarray, y_traj: np.ndarray) -> None:
        """
        Plot the particle trajectory on a 2D projection of the spacetime grid.
        
        Args:
            x_traj (np.ndarray): X coordinates of trajectory
            y_traj (np.ndarray): Y coordinates of trajectory
        """
        try:
            fig, ax = plt.subplots(figsize=(8, 8))
            # Plot potential as a contour map in the background
            contour = ax.contourf(self.X, self.Y, self.Z, cmap='viridis', alpha=0.5)
            plt.colorbar(contour, ax=ax, label='Spacetime Curvature (Scaled Potential)')
            # Overlay trajectory
            ax.plot(x_traj, y_traj, 'r-', linewidth=2, label='Particle Trajectory')
            ax.set_xlabel('X (arbitrary units)')
            ax.set_ylabel('Y (arbitrary units)')
            ax.set_title('Particle Trajectory in Gravitational Well')
            ax.legend()
            ax.grid(True)
            plt.show()
            logger.info("Trajectory visualization displayed.")
        except Exception as e:
            logger.error(f"Error in trajectory visualization: {str(e)}")
            raise

    def export_grid_data(self, filename: str = 'spacetime_grid.csv') -> None:
        """
        Export the grid and potential data to a CSV file.
        
        Args:
            filename (str): Output CSV file name (default: 'spacetime_grid.csv')
        """
        if self.Z is None:
            raise ValueError("Potential not calculated. Call calculate_potential() first.")
        try:
            with open(filename, 'w', newline='') as f:
                writer = csv.writer(f)
                writer.writerow(['X', 'Y', 'Potential'])
                for i in range(self.grid_size):
                    for j in range(self.grid_size):
                        writer.writerow([self.X[i, j], self.Y[i, j], self.Z[i, j]])
            logger.info(f"Grid data exported to '{filename}'")
        except Exception as e:
            logger.error(f"Error exporting grid data: {str(e)}")
            raise

def parse_arguments() -> argparse.Namespace:
    """Parse command-line arguments for simulation parameters."""
    parser = argparse.ArgumentParser(description='Spacetime Bending Simulation Tool')
    parser.add_argument('--mass', type=float, default=1.989e30,
                        help='Central mass in kg (default: Sun\'s mass, 1.989e30)')
    parser.add_argument('--grid-size', type=int, default=100,
                        help='Grid points per axis (default: 100)')
    parser.add_argument('--scale', type=float, default=1e-10,
                        help='Scaling factor for potential visibility (default: 1e-10)')
    parser.add_argument('--no-interactive', action='store_true',
                        help='Disable interactive mode and use default or CLI parameters')
    parser.add_argument('--save-plot', action='store_true',
                        help='Save the 3D plot to a file')
    parser.add_argument('--simulate-trajectory', action='store_true',
                        help='Run particle trajectory simulation')
    parser.add_argument('--save-data', action='store_true',
                        help='Save trajectory and grid data to CSV files')
    return parser.parse_args()

def get_user_input(args: argparse.Namespace) -> Tuple[float, int, float]:
    """Prompt user for simulation parameters if interactive mode is enabled."""
    if args.no_interactive:
        return args.mass, args.grid_size, args.scale
    print("Spacetime Bending Simulation Configuration")
    try:
        mass = float(input(f"Enter central mass in kg (default: {args.mass:.2e}): ") or args.mass)
        grid_size = int(input(f"Enter grid size per axis (default: {args.grid_size}): ") or args.grid_size)
        scale = float(input(f"Enter scaling factor for curvature (default: {args.scale:.2e}): ") or args.scale)
        return mass, grid_size, scale
    except ValueError as e:
        logger.warning(f"Invalid input received: {str(e)}. Using default values.")
        print("Invalid input. Using default or command-line values.")
        return args.mass, args.grid_size, args.scale

def main():
    """Main function to orchestrate the spacetime bending simulation."""
    args = parse_arguments()
    mass, grid_size, scale = get_user_input(args)
    
    # Validate inputs
    if mass <= 0:
        logger.error("Mass must be positive.")
        raise ValueError("Mass must be positive.")
    if grid_size < 10:
        logger.error("Grid size too small.")
        raise ValueError("Grid size must be at least 10.")
    if scale <= 0:
        logger.error("Scale factor must be positive.")
        raise ValueError("Scale factor must be positive.")
    
    # Initialize simulator
    try:
        simulator = SpacetimeSimulator(mass=mass, grid_size=grid_size, scale=scale)
        simulator.create_grid()
        simulator.calculate_potential()
        simulator.visualize_3d(colormap='viridis', save_plot=args.save_plot)
        if args.save_data:
            simulator.export_grid_data()
        
        if args.simulate_trajectory:
            x_traj, y_traj = simulator.simulate_trajectory(save_data=args.save_data)
            simulator.visualize_trajectory(x_traj, y_traj)
    except Exception as e:
        logger.critical(f"Simulation failed: {str(e)}")
        print(f"An error occurred: {str(e)}")
        sys.exit(1)

if __name__ == "__main__":
    logger.info("Starting Spacetime Bending Simulation.")
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")  # Suppress minor numerical warnings
        main()
    logger.info("Simulation completed successfully.")