import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, differential_evolution
from mpl_toolkits.mplot3d import Axes3D

# Set random seed for reproducibility
np.random.seed(42)

class TelescopeOpticalSystem:
    """
    A class to model and optimize a two-mirror Cassegrain telescope system.
    
    The system aims to minimize optical aberrations (spherical, coma, astigmatism)
    by optimizing mirror shapes (conic constants) and positions.
    """
    
    def __init__(self, focal_length=10.0, diameter=2.0, field_of_view=0.5):
        """
        Initialize telescope parameters.
        
        Parameters:
        -----------
        focal_length : float
            System focal length in meters
        diameter : float
            Primary mirror diameter in meters
        field_of_view : float
            Field of view in degrees
        """
        self.F = focal_length  # System focal length (m)
        self.D = diameter      # Primary mirror diameter (m)
        self.fov = field_of_view * np.pi / 180  # Field of view (radians)
        
    def cassegrain_geometry(self, params):
        """
        Calculate Cassegrain geometry from optimization parameters.
        
        Parameters:
        -----------
        params : array-like
            [R1, R2, d, kappa1, kappa2]
            R1, R2: radii of curvature (m)
            d: mirror separation (m)
            kappa1, kappa2: conic constants
            
        Returns:
        --------
        dict : Geometric parameters of the system
        """
        R1, R2, d, kappa1, kappa2 = params
        
        # Calculate focal lengths from radii (paraxial approximation)
        f1 = R1 / 2  # Primary focal length
        f2 = R2 / 2  # Secondary focal length
        
        # Calculate magnification
        m = (d - f1) / (d + f2)
        
        # System focal length
        F_system = f1 * abs(m)
        
        return {
            'f1': f1, 'f2': f2, 'd': d, 
            'kappa1': kappa1, 'kappa2': kappa2,
            'm': m, 'F_system': F_system
        }
    
    def seidel_aberrations(self, params):
        """
        Calculate Seidel aberration coefficients.
        
        These coefficients quantify the optical aberrations:
        W040: Spherical aberration
        W131: Coma
        W222: Astigmatism
        """
        geom = self.cassegrain_geometry(params)
        f1, f2 = geom['f1'], geom['f2']
        kappa1, kappa2 = geom['kappa1'], geom['kappa2']
        
        # Ray height at primary mirror (edge of aperture)
        h1 = self.D / 2
        
        # Ray height at secondary (scaled by magnification)
        h2 = h1 * abs(geom['m'])
        
        # Field angle (paraxial)
        u_bar = self.fov
        
        # Seidel coefficients using classical formulas
        # Spherical aberration (W040)
        W040_1 = -(h1**4) / (128 * f1**3) * (1 + kappa1)**2
        W040_2 = -(h2**4) / (128 * f2**3) * (1 + kappa2)**2
        W040 = W040_1 + W040_2
        
        # Coma (W131)
        W131_1 = -(h1**3 * u_bar) / (16 * f1**2) * (1 + kappa1)
        W131_2 = -(h2**3 * u_bar) / (16 * f2**2) * (1 + kappa2)
        W131 = W131_1 + W131_2
        
        # Astigmatism (W222)
        W222_1 = -(h1**2 * u_bar**2) / (2 * f1)
        W222_2 = -(h2**2 * u_bar**2) / (2 * f2)
        W222 = W222_1 + W222_2
        
        return {
            'W040': W040,
            'W131': W131,
            'W222': W222,
            'components': {
                'W040_1': W040_1, 'W040_2': W040_2,
                'W131_1': W131_1, 'W131_2': W131_2,
                'W222_1': W222_1, 'W222_2': W222_2
            }
        }
    
    def objective_function(self, params):
        """
        Objective function to minimize: weighted sum of aberrations.
        
        We use RMS (root mean square) of aberration coefficients as our metric.
        """
        try:
            aberrations = self.seidel_aberrations(params)
            
            # Weight different aberrations
            w_spherical = 1.0
            w_coma = 1.5  # Coma is often more visually disturbing
            w_astig = 1.0
            
            # Calculate weighted RMS aberration
            rms = np.sqrt(
                w_spherical * aberrations['W040']**2 +
                w_coma * aberrations['W131']**2 +
                w_astig * aberrations['W222']**2
            )
            
            # Add penalty for unrealistic focal length deviation
            geom = self.cassegrain_geometry(params)
            focal_penalty = 100 * (geom['F_system'] - self.F)**2
            
            return rms + focal_penalty
            
        except:
            return 1e10  # Return large value if calculation fails
    
    def optimize(self):
        """
        Optimize telescope parameters using differential evolution.
        
        This global optimization method is robust for non-convex problems.
        """
        # Parameter bounds: [R1, R2, d, kappa1, kappa2]
        # R1: Primary radius of curvature (4-8 m)
        # R2: Secondary radius of curvature (1-3 m)
        # d: Mirror separation (2-5 m)
        # kappa1, kappa2: Conic constants (-2 to 0 for typical mirrors)
        bounds = [
            (4.0, 8.0),    # R1
            (1.0, 3.0),    # R2
            (2.0, 5.0),    # d
            (-2.0, 0.0),   # kappa1
            (-2.0, 0.0)    # kappa2
        ]
        
        print("Starting optimization...")
        print("=" * 60)
        
        # Use differential evolution for global optimization
        result = differential_evolution(
            self.objective_function,
            bounds,
            maxiter=300,
            popsize=15,
            tol=1e-10,
            seed=42,
            disp=True
        )
        
        return result
    
    def analyze_results(self, initial_params, optimized_params):
        """
        Comprehensive analysis comparing initial and optimized designs.
        """
        print("\n" + "=" * 60)
        print("OPTIMIZATION RESULTS")
        print("=" * 60)
        
        # Initial system
        print("\nINITIAL DESIGN:")
        print("-" * 60)
        geom_init = self.cassegrain_geometry(initial_params)
        aber_init = self.seidel_aberrations(initial_params)
        
        print(f"Primary radius of curvature (R1):  {initial_params[0]:.3f} m")
        print(f"Secondary radius of curvature (R2): {initial_params[1]:.3f} m")
        print(f"Mirror separation (d):              {initial_params[2]:.3f} m")
        print(f"Primary conic constant (κ1):        {initial_params[3]:.3f}")
        print(f"Secondary conic constant (κ2):      {initial_params[4]:.3f}")
        print(f"System focal length:                 {geom_init['F_system']:.3f} m")
        print(f"\nAberrations:")
        print(f"  Spherical (W040): {aber_init['W040']:.6e} λ")
        print(f"  Coma (W131):      {aber_init['W131']:.6e} λ")
        print(f"  Astigmatism (W222): {aber_init['W222']:.6e} λ")
        print(f"  RMS aberration:   {np.sqrt(aber_init['W040']**2 + aber_init['W131']**2 + aber_init['W222']**2):.6e} λ")
        
        # Optimized system
        print("\n" + "=" * 60)
        print("OPTIMIZED DESIGN:")
        print("-" * 60)
        geom_opt = self.cassegrain_geometry(optimized_params)
        aber_opt = self.seidel_aberrations(optimized_params)
        
        print(f"Primary radius of curvature (R1):  {optimized_params[0]:.3f} m")
        print(f"Secondary radius of curvature (R2): {optimized_params[1]:.3f} m")
        print(f"Mirror separation (d):              {optimized_params[2]:.3f} m")
        print(f"Primary conic constant (κ1):        {optimized_params[3]:.3f}")
        print(f"Secondary conic constant (κ2):      {optimized_params[4]:.3f}")
        print(f"System focal length:                 {geom_opt['F_system']:.3f} m")
        print(f"\nAberrations:")
        print(f"  Spherical (W040): {aber_opt['W040']:.6e} λ")
        print(f"  Coma (W131):      {aber_opt['W131']:.6e} λ")
        print(f"  Astigmatism (W222): {aber_opt['W222']:.6e} λ")
        print(f"  RMS aberration:   {np.sqrt(aber_opt['W040']**2 + aber_opt['W131']**2 + aber_opt['W222']**2):.6e} λ")
        
        # Improvement metrics
        print("\n" + "=" * 60)
        print("IMPROVEMENT METRICS:")
        print("-" * 60)
        rms_init = np.sqrt(aber_init['W040']**2 + aber_init['W131']**2 + aber_init['W222']**2)
        rms_opt = np.sqrt(aber_opt['W040']**2 + aber_opt['W131']**2 + aber_opt['W222']**2)
        improvement = (rms_init - rms_opt) / rms_init * 100
        
        print(f"RMS aberration improvement: {improvement:.2f}%")
        print(f"Spherical aberration reduction: {(abs(aber_init['W040']) - abs(aber_opt['W040'])) / abs(aber_init['W040']) * 100:.2f}%")
        print(f"Coma reduction: {(abs(aber_init['W131']) - abs(aber_opt['W131'])) / abs(aber_init['W131']) * 100:.2f}%")
        print(f"Astigmatism reduction: {(abs(aber_init['W222']) - abs(aber_opt['W222'])) / abs(aber_init['W222']) * 100:.2f}%")
        print("=" * 60)
        
        return geom_init, geom_opt, aber_init, aber_opt

def visualize_optimization_results(telescope, initial_params, optimized_params):
    """
    Create comprehensive visualization of optimization results.
    """
    # Calculate data for both systems
    aber_init = telescope.seidel_aberrations(initial_params)
    aber_opt = telescope.seidel_aberrations(optimized_params)
    geom_init = telescope.cassegrain_geometry(initial_params)
    geom_opt = telescope.cassegrain_geometry(optimized_params)
    
    # Create figure with subplots
    fig = plt.figure(figsize=(16, 12))
    
    # 1. Aberration Comparison (Bar Chart)
    ax1 = plt.subplot(2, 3, 1)
    aberration_types = ['Spherical\n(W040)', 'Coma\n(W131)', 'Astigmatism\n(W222)']
    initial_values = [abs(aber_init['W040']), abs(aber_init['W131']), abs(aber_init['W222'])]
    optimized_values = [abs(aber_opt['W040']), abs(aber_opt['W131']), abs(aber_opt['W222'])]
    
    x = np.arange(len(aberration_types))
    width = 0.35
    
    bars1 = ax1.bar(x - width/2, initial_values, width, label='Initial', alpha=0.8, color='#FF6B6B')
    bars2 = ax1.bar(x + width/2, optimized_values, width, label='Optimized', alpha=0.8, color='#4ECDC4')
    
    ax1.set_ylabel('Aberration Magnitude (λ)', fontsize=11, fontweight='bold')
    ax1.set_title('Aberration Comparison', fontsize=12, fontweight='bold', pad=15)
    ax1.set_xticks(x)
    ax1.set_xticklabels(aberration_types, fontsize=9)
    ax1.legend(fontsize=10)
    ax1.grid(axis='y', alpha=0.3, linestyle='--')
    ax1.set_yscale('log')
    
    # 2. Parameter Changes (Radar Chart)
    ax2 = plt.subplot(2, 3, 2, projection='polar')
    
    # Normalize parameters to 0-1 scale for radar chart
    params_labels = ['R1\n(Primary)', 'R2\n(Secondary)', 'd\n(Separation)', 'κ1', 'κ2']
    
    # Normalize to bounds used in optimization
    bounds = [(4.0, 8.0), (1.0, 3.0), (2.0, 5.0), (-2.0, 0.0), (-2.0, 0.0)]
    initial_norm = [(initial_params[i] - bounds[i][0]) / (bounds[i][1] - bounds[i][0]) for i in range(5)]
    optimized_norm = [(optimized_params[i] - bounds[i][0]) / (bounds[i][1] - bounds[i][0]) for i in range(5)]
    
    angles = np.linspace(0, 2 * np.pi, len(params_labels), endpoint=False).tolist()
    initial_norm += initial_norm[:1]
    optimized_norm += optimized_norm[:1]
    angles += angles[:1]
    
    ax2.plot(angles, initial_norm, 'o-', linewidth=2, label='Initial', color='#FF6B6B', markersize=8)
    ax2.fill(angles, initial_norm, alpha=0.15, color='#FF6B6B')
    ax2.plot(angles, optimized_norm, 'o-', linewidth=2, label='Optimized', color='#4ECDC4', markersize=8)
    ax2.fill(angles, optimized_norm, alpha=0.15, color='#4ECDC4')
    
    ax2.set_xticks(angles[:-1])
    ax2.set_xticklabels(params_labels, fontsize=9)
    ax2.set_ylim(0, 1)
    ax2.set_title('Parameter Space\n(Normalized)', fontsize=12, fontweight='bold', pad=20)
    ax2.legend(loc='upper right', bbox_to_anchor=(1.3, 1.1), fontsize=10)
    ax2.grid(True, alpha=0.3)
    
    # 3. RMS Aberration
    ax3 = plt.subplot(2, 3, 3)
    rms_init = np.sqrt(aber_init['W040']**2 + aber_init['W131']**2 + aber_init['W222']**2)
    rms_opt = np.sqrt(aber_opt['W040']**2 + aber_opt['W131']**2 + aber_opt['W222']**2)
    
    bars = ax3.bar(['Initial', 'Optimized'], [rms_init, rms_opt], 
                   color=['#FF6B6B', '#4ECDC4'], alpha=0.8, width=0.5)
    ax3.set_ylabel('RMS Aberration (λ)', fontsize=11, fontweight='bold')
    ax3.set_title('Total RMS Aberration', fontsize=12, fontweight='bold', pad=15)
    ax3.grid(axis='y', alpha=0.3, linestyle='--')
    ax3.set_yscale('log')
    
    # Add percentage improvement text
    improvement = (rms_init - rms_opt) / rms_init * 100
    ax3.text(0.5, (rms_init + rms_opt) / 2, f'{improvement:.1f}%\nimprovement', 
             ha='center', va='center', fontsize=11, fontweight='bold',
             bbox=dict(boxstyle='round', facecolor='white', alpha=0.8, edgecolor='green', linewidth=2))
    
    # 4. Optical System Schematic (Side View)
    ax4 = plt.subplot(2, 3, 4)
    
    # Draw initial system
    d_init = initial_params[2]
    R1_init = initial_params[0]
    R2_init = initial_params[1]
    
    # Primary mirror (parabolic approximation)
    y_primary = np.linspace(-telescope.D/2, telescope.D/2, 100)
    z_primary_init = -y_primary**2 / (2 * R1_init)
    
    ax4.plot(z_primary_init, y_primary, 'r-', linewidth=2.5, label='Initial Primary', alpha=0.6)
    
    # Secondary mirror position
    z_secondary_init = d_init
    y_secondary = np.linspace(-telescope.D/4, telescope.D/4, 50)
    z_secondary_curve_init = z_secondary_init + y_secondary**2 / (2 * R2_init)
    
    ax4.plot(z_secondary_curve_init, y_secondary, 'r--', linewidth=2.5, label='Initial Secondary', alpha=0.6)
    
    ax4.set_xlabel('Optical Axis (m)', fontsize=11, fontweight='bold')
    ax4.set_ylabel('Height (m)', fontsize=11, fontweight='bold')
    ax4.set_title('Optical System Layout - Initial', fontsize=12, fontweight='bold', pad=15)
    ax4.grid(True, alpha=0.3, linestyle='--')
    ax4.legend(fontsize=9)
    ax4.set_aspect('equal', adjustable='box')
    
    # 5. Optimized System Schematic
    ax5 = plt.subplot(2, 3, 5)
    
    d_opt = optimized_params[2]
    R1_opt = optimized_params[0]
    R2_opt = optimized_params[1]
    
    # Primary mirror
    z_primary_opt = -y_primary**2 / (2 * R1_opt)
    ax5.plot(z_primary_opt, y_primary, 'b-', linewidth=2.5, label='Optimized Primary', alpha=0.8)
    
    # Secondary mirror
    z_secondary_opt = d_opt
    z_secondary_curve_opt = z_secondary_opt + y_secondary**2 / (2 * R2_opt)
    ax5.plot(z_secondary_curve_opt, y_secondary, 'b--', linewidth=2.5, label='Optimized Secondary', alpha=0.8)
    
    ax5.set_xlabel('Optical Axis (m)', fontsize=11, fontweight='bold')
    ax5.set_ylabel('Height (m)', fontsize=11, fontweight='bold')
    ax5.set_title('Optical System Layout - Optimized', fontsize=12, fontweight='bold', pad=15)
    ax5.grid(True, alpha=0.3, linestyle='--')
    ax5.legend(fontsize=9)
    ax5.set_aspect('equal', adjustable='box')
    
    # 6. Aberration Components Breakdown
    ax6 = plt.subplot(2, 3, 6)
    
    components_init = [
        abs(aber_init['components']['W040_1']), 
        abs(aber_init['components']['W040_2']),
        abs(aber_init['components']['W131_1']), 
        abs(aber_init['components']['W131_2']),
        abs(aber_init['components']['W222_1']), 
        abs(aber_init['components']['W222_2'])
    ]
    
    components_opt = [
        abs(aber_opt['components']['W040_1']), 
        abs(aber_opt['components']['W040_2']),
        abs(aber_opt['components']['W131_1']), 
        abs(aber_opt['components']['W131_2']),
        abs(aber_opt['components']['W222_1']), 
        abs(aber_opt['components']['W222_2'])
    ]
    
    component_labels = ['W040\nPrimary', 'W040\nSecondary', 'W131\nPrimary', 
                       'W131\nSecondary', 'W222\nPrimary', 'W222\nSecondary']
    
    x = np.arange(len(component_labels))
    width = 0.35
    
    ax6.bar(x - width/2, components_init, width, label='Initial', alpha=0.8, color='#FF6B6B')
    ax6.bar(x + width/2, components_opt, width, label='Optimized', alpha=0.8, color='#4ECDC4')
    
    ax6.set_ylabel('Aberration Component (λ)', fontsize=11, fontweight='bold')
    ax6.set_title('Individual Aberration Components', fontsize=12, fontweight='bold', pad=15)
    ax6.set_xticks(x)
    ax6.set_xticklabels(component_labels, fontsize=8)
    ax6.legend(fontsize=10)
    ax6.grid(axis='y', alpha=0.3, linestyle='--')
    ax6.set_yscale('log')
    
    plt.tight_layout()
    plt.savefig('telescope_optimization_results.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    # Create additional 3D visualization
    fig2 = plt.figure(figsize=(16, 6))
    
    # 3D view of initial system
    ax7 = fig2.add_subplot(121, projection='3d')
    
    theta = np.linspace(0, 2*np.pi, 50)
    y_circle = np.outer(y_primary, np.cos(theta))
    x_circle = np.outer(y_primary, np.sin(theta))
    z_circle_primary_init = np.outer(z_primary_init, np.ones(50))
    
    ax7.plot_surface(x_circle, y_circle, z_circle_primary_init, alpha=0.6, cmap='Reds', edgecolor='none')
    
    y_circle_sec = np.outer(y_secondary, np.cos(theta))
    x_circle_sec = np.outer(y_secondary, np.sin(theta))
    z_circle_secondary_init = np.outer(z_secondary_curve_init, np.ones(50))
    
    ax7.plot_surface(x_circle_sec, y_circle_sec, z_circle_secondary_init, alpha=0.6, cmap='Oranges', edgecolor='none')
    
    ax7.set_xlabel('X (m)', fontsize=10, fontweight='bold')
    ax7.set_ylabel('Y (m)', fontsize=10, fontweight='bold')
    ax7.set_zlabel('Z (m)', fontsize=10, fontweight='bold')
    ax7.set_title('Initial System - 3D View', fontsize=12, fontweight='bold', pad=15)
    
    # 3D view of optimized system
    ax8 = fig2.add_subplot(122, projection='3d')
    
    z_circle_primary_opt = np.outer(z_primary_opt, np.ones(50))
    ax8.plot_surface(x_circle, y_circle, z_circle_primary_opt, alpha=0.7, cmap='Blues', edgecolor='none')
    
    z_circle_secondary_opt = np.outer(z_secondary_curve_opt, np.ones(50))
    ax8.plot_surface(x_circle_sec, y_circle_sec, z_circle_secondary_opt, alpha=0.7, cmap='Greens', edgecolor='none')
    
    ax8.set_xlabel('X (m)', fontsize=10, fontweight='bold')
    ax8.set_ylabel('Y (m)', fontsize=10, fontweight='bold')
    ax8.set_zlabel('Z (m)', fontsize=10, fontweight='bold')
    ax8.set_title('Optimized System - 3D View', fontsize=12, fontweight='bold', pad=15)
    
    plt.tight_layout()
    plt.savefig('telescope_3d_comparison.png', dpi=300, bbox_inches='tight')
    plt.show()

# Main execution
if __name__ == "__main__":
    print("=" * 60)
    print("SPACE TELESCOPE OPTICAL SYSTEM OPTIMIZATION")
    print("=" * 60)
    print("\nObjective: Minimize optical aberrations through optimization")
    print("of mirror shapes and positions in a Cassegrain telescope.\n")
    
    # Initialize telescope system
    # Specifications: 10m focal length, 2m diameter, 0.5° field of view
    telescope = TelescopeOpticalSystem(
        focal_length=10.0,  # meters
        diameter=2.0,       # meters
        field_of_view=0.5   # degrees
    )
    
    # Initial design parameters (suboptimal starting point)
    # [R1, R2, d, kappa1, kappa2]
    initial_params = np.array([6.0, 2.0, 3.5, -1.0, -1.0])
    
    print("Initial Parameters:")
    print(f"  Primary radius (R1): {initial_params[0]} m")
    print(f"  Secondary radius (R2): {initial_params[1]} m")
    print(f"  Mirror separation (d): {initial_params[2]} m")
    print(f"  Primary conic (κ1): {initial_params[3]}")
    print(f"  Secondary conic (κ2): {initial_params[4]}")
    
    # Perform optimization
    result = telescope.optimize()
    
    # Extract optimized parameters
    optimized_params = result.x
    
    # Detailed analysis
    geom_init, geom_opt, aber_init, aber_opt = telescope.analyze_results(
        initial_params, optimized_params
    )
    
    # Visualize results
    print("\nGenerating visualizations...")
    visualize_optimization_results(telescope, initial_params, optimized_params)
    
    print("\n" + "=" * 60)
    print("Optimization complete! Graphs saved.")
    print("=" * 60)

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')

# Physical constants
AU = 1.496e11  # Astronomical Unit in meters
G = 6.674e-11  # Gravitational constant
M_sun = 1.989e30  # Solar mass in kg

# Planetary data: [semi-major axis (AU), orbital period (days), mass (kg)]
planets = {
    'Earth': {'a': 1.0, 'T': 365.25, 'mass': 5.972e24, 'radius': 6.371e6},
    'Mars': {'a': 1.524, 'T': 687.0, 'mass': 6.417e23, 'radius': 3.390e6},
    'Jupiter': {'a': 5.203, 'T': 4332.6, 'mass': 1.898e27, 'radius': 6.9911e7}
}

def planetary_position(planet_name, time_days):
    """
    Calculate planetary position in heliocentric coordinates.
    Assumes circular orbits for simplicity.
    
    Parameters:
    - planet_name: Name of the planet
    - time_days: Time in days from reference epoch
    
    Returns:
    - x, y, z coordinates in meters
    """
    planet = planets[planet_name]
    a = planet['a'] * AU  # Semi-major axis
    T = planet['T']  # Orbital period
    
    # Mean anomaly
    M = 2 * np.pi * time_days / T
    
    # Position (circular orbit approximation)
    x = a * np.cos(M)
    y = a * np.sin(M)
    z = 0  # Assume coplanar orbits
    
    return np.array([x, y, z])

def velocity_at_position(planet_name, time_days):
    """
    Calculate orbital velocity of a planet.
    
    Returns:
    - vx, vy, vz velocity components in m/s
    """
    planet = planets[planet_name]
    a = planet['a'] * AU
    T = planet['T']
    
    # Angular velocity
    omega = 2 * np.pi / (T * 86400)  # rad/s
    
    # Mean anomaly
    M = 2 * np.pi * time_days / T
    
    # Velocity (circular orbit)
    v_mag = omega * a
    vx = -v_mag * np.sin(M)
    vy = v_mag * np.cos(M)
    vz = 0
    
    return np.array([vx, vy, vz])

def hohmann_delta_v(r1, r2):
    """
    Calculate delta-v for Hohmann transfer between two circular orbits.
    
    Parameters:
    - r1: Initial orbital radius (m)
    - r2: Final orbital radius (m)
    
    Returns:
    - delta_v1: Delta-v at departure (m/s)
    - delta_v2: Delta-v at arrival (m/s)
    """
    mu = G * M_sun
    
    # Velocities in circular orbits
    v1 = np.sqrt(mu / r1)
    v2 = np.sqrt(mu / r2)
    
    # Transfer orbit parameters
    a_transfer = (r1 + r2) / 2
    
    # Velocities at perihelion and aphelion of transfer orbit
    v_transfer_peri = np.sqrt(mu * (2/r1 - 1/a_transfer))
    v_transfer_apo = np.sqrt(mu * (2/r2 - 1/a_transfer))
    
    # Delta-v requirements
    delta_v1 = abs(v_transfer_peri - v1)
    delta_v2 = abs(v2 - v_transfer_apo)
    
    return delta_v1, delta_v2

def gravity_assist_delta_v(planet_name, v_inf_in, v_inf_out):
    """
    Calculate the effect of gravity assist.
    
    Parameters:
    - planet_name: Name of the planet for gravity assist
    - v_inf_in: Incoming hyperbolic excess velocity magnitude (m/s)
    - v_inf_out: Outgoing hyperbolic excess velocity magnitude (m/s)
    
    Returns:
    - delta_v_saved: Effective delta-v saved by gravity assist (m/s)
    """
    planet = planets[planet_name]
    mu_planet = G * planet['mass']
    r_p = planet['radius'] * 2  # Periapsis radius (safety factor)
    
    # Turn angle (simplified)
    delta = 2 * np.arcsin(1 / (1 + r_p * v_inf_in**2 / mu_planet))
    
    # Gravity assist benefit
    delta_v_saved = abs(v_inf_out - v_inf_in) * np.sin(delta/2)
    
    return delta_v_saved

def mission_cost(params):
    """
    Objective function for trajectory optimization.
    
    Parameters:
    - params: [t_departure, t_mars_arrival, t_jupiter_arrival] in days
    
    Returns:
    - cost: Weighted sum of fuel consumption and time
    """
    t_dep, t_mars, t_jup = params
    
    # Ensure proper time ordering
    if t_mars <= t_dep or t_jup <= t_mars:
        return 1e10  # Penalty for invalid trajectory
    
    try:
        # Earth departure
        pos_earth_dep = planetary_position('Earth', t_dep)
        vel_earth_dep = velocity_at_position('Earth', t_dep)
        r_earth = np.linalg.norm(pos_earth_dep)
        
        # Mars arrival
        pos_mars_arr = planetary_position('Mars', t_mars)
        vel_mars_arr = velocity_at_position('Mars', t_mars)
        r_mars = np.linalg.norm(pos_mars_arr)
        
        # Jupiter arrival
        pos_jup_arr = planetary_position('Jupiter', t_jup)
        vel_jup_arr = velocity_at_position('Jupiter', t_jup)
        r_jup = np.linalg.norm(pos_jup_arr)
        
        # Delta-v calculation
        # Earth to Mars leg
        dv1_em, dv2_em = hohmann_delta_v(r_earth, r_mars)
        
        # Mars to Jupiter leg with gravity assist
        dv1_mj, dv2_mj = hohmann_delta_v(r_mars, r_jup)
        
        # Simplified gravity assist benefit at Mars
        v_inf_mars = 5000  # Simplified hyperbolic excess velocity (m/s)
        ga_benefit = gravity_assist_delta_v('Mars', v_inf_mars, v_inf_mars * 1.2)
        
        # Total delta-v
        delta_v_total = dv1_em + dv2_em + dv1_mj + dv2_mj - ga_benefit
        
        # Mission duration
        duration = t_jup - t_dep
        
        # Penalties for unrealistic trajectories
        leg1_time = t_mars - t_dep
        leg2_time = t_jup - t_mars
        
        penalty = 0
        if leg1_time < 150 or leg1_time > 400:  # Earth-Mars should be 150-400 days
            penalty += 1e6
        if leg2_time < 400 or leg2_time > 800:  # Mars-Jupiter should be 400-800 days
            penalty += 1e6
        if delta_v_total < 0 or delta_v_total > 50000:  # Unrealistic delta-v
            penalty += 1e6
        
        # Weights for multi-objective optimization
        w_fuel = 1.0      # Weight for fuel (delta-v)
        w_time = 0.01     # Weight for time
        w_penalty = 1.0   # Weight for constraints
        
        cost = w_fuel * delta_v_total + w_time * duration + w_penalty * penalty
        
        return cost
        
    except:
        return 1e10

def optimize_trajectory():
    """
    Optimize the mission trajectory using differential evolution.
    
    Returns:
    - result: Optimization result
    """
    # Bounds for optimization variables (in days)
    bounds = [
        (0, 365),      # t_departure: within first year
        (200, 600),    # t_mars_arrival: 200-600 days after epoch
        (700, 1500)    # t_jupiter_arrival: 700-1500 days after epoch
    ]
    
    print("Starting trajectory optimization...")
    print("This may take a minute...\n")
    
    result = differential_evolution(
        mission_cost,
        bounds,
        strategy='best1bin',
        maxiter=300,
        popsize=20,
        tol=0.01,
        mutation=(0.5, 1),
        recombination=0.7,
        seed=42,
        disp=True
    )
    
    return result

def calculate_trajectory_details(params):
    """
    Calculate detailed trajectory information.
    """
    t_dep, t_mars, t_jup = params
    
    # Positions
    pos_earth_dep = planetary_position('Earth', t_dep)
    pos_mars_arr = planetary_position('Mars', t_mars)
    pos_jup_arr = planetary_position('Jupiter', t_jup)
    
    # Distances
    r_earth = np.linalg.norm(pos_earth_dep)
    r_mars = np.linalg.norm(pos_mars_arr)
    r_jup = np.linalg.norm(pos_jup_arr)
    
    # Delta-v calculations
    dv1_em, dv2_em = hohmann_delta_v(r_earth, r_mars)
    dv1_mj, dv2_mj = hohmann_delta_v(r_mars, r_jup)
    
    # Gravity assist benefit
    v_inf_mars = 5000
    ga_benefit = gravity_assist_delta_v('Mars', v_inf_mars, v_inf_mars * 1.2)
    
    delta_v_total = dv1_em + dv2_em + dv1_mj + dv2_mj - ga_benefit
    duration = t_jup - t_dep
    
    return {
        'positions': {
            'earth_dep': pos_earth_dep,
            'mars_arr': pos_mars_arr,
            'jupiter_arr': pos_jup_arr
        },
        'times': {
            't_dep': t_dep,
            't_mars': t_mars,
            't_jup': t_jup,
            'leg1': t_mars - t_dep,
            'leg2': t_jup - t_mars,
            'total': duration
        },
        'delta_v': {
            'earth_departure': dv1_em,
            'mars_arrival': dv2_em,
            'mars_departure': dv1_mj,
            'jupiter_arrival': dv2_mj,
            'gravity_assist_benefit': ga_benefit,
            'total': delta_v_total
        }
    }

def plot_trajectory_3d(params):
    """
    Plot 3D trajectory visualization.
    """
    t_dep, t_mars, t_jup = params
    
    fig = plt.figure(figsize=(14, 10))
    
    # 3D trajectory plot
    ax1 = fig.add_subplot(221, projection='3d')
    
    # Generate planetary orbits
    time_orbit = np.linspace(0, 4332, 500)
    
    # Earth orbit
    earth_orbit = np.array([planetary_position('Earth', t) for t in time_orbit])
    ax1.plot(earth_orbit[:, 0]/AU, earth_orbit[:, 1]/AU, earth_orbit[:, 2]/AU, 
             'b-', alpha=0.3, label='Earth Orbit')
    
    # Mars orbit
    mars_orbit = np.array([planetary_position('Mars', t) for t in time_orbit])
    ax1.plot(mars_orbit[:, 0]/AU, mars_orbit[:, 1]/AU, mars_orbit[:, 2]/AU, 
             'r-', alpha=0.3, label='Mars Orbit')
    
    # Jupiter orbit
    jupiter_orbit = np.array([planetary_position('Jupiter', t) for t in time_orbit])
    ax1.plot(jupiter_orbit[:, 0]/AU, jupiter_orbit[:, 1]/AU, jupiter_orbit[:, 2]/AU, 
             'orange', alpha=0.3, label='Jupiter Orbit')
    
    # Spacecraft trajectory (simplified)
    trajectory_times = np.linspace(t_dep, t_jup, 200)
    spacecraft_pos = []
    
    for t in trajectory_times:
        if t <= t_mars:
            # Earth to Mars leg (linear interpolation for visualization)
            alpha = (t - t_dep) / (t_mars - t_dep)
            pos_earth = planetary_position('Earth', t_dep)
            pos_mars = planetary_position('Mars', t_mars)
            pos = pos_earth + alpha * (pos_mars - pos_earth)
        else:
            # Mars to Jupiter leg
            alpha = (t - t_mars) / (t_jup - t_mars)
            pos_mars = planetary_position('Mars', t_mars)
            pos_jup = planetary_position('Jupiter', t_jup)
            pos = pos_mars + alpha * (pos_jup - pos_mars)
        spacecraft_pos.append(pos)
    
    spacecraft_pos = np.array(spacecraft_pos)
    ax1.plot(spacecraft_pos[:, 0]/AU, spacecraft_pos[:, 1]/AU, spacecraft_pos[:, 2]/AU, 
             'g-', linewidth=2, label='Spacecraft Trajectory')
    
    # Mark key positions
    pos_earth_dep = planetary_position('Earth', t_dep)
    pos_mars_arr = planetary_position('Mars', t_mars)
    pos_jup_arr = planetary_position('Jupiter', t_jup)
    
    ax1.scatter([0], [0], [0], c='yellow', s=200, marker='o', label='Sun')
    ax1.scatter([pos_earth_dep[0]/AU], [pos_earth_dep[1]/AU], [pos_earth_dep[2]/AU], 
                c='blue', s=100, marker='o', label='Earth Departure')
    ax1.scatter([pos_mars_arr[0]/AU], [pos_mars_arr[1]/AU], [pos_mars_arr[2]/AU], 
                c='red', s=100, marker='o', label='Mars Flyby')
    ax1.scatter([pos_jup_arr[0]/AU], [pos_jup_arr[1]/AU], [pos_jup_arr[2]/AU], 
                c='orange', s=100, marker='o', label='Jupiter Arrival')
    
    ax1.set_xlabel('X (AU)')
    ax1.set_ylabel('Y (AU)')
    ax1.set_zlabel('Z (AU)')
    ax1.set_title('3D Trajectory Visualization', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=8)
    ax1.grid(True, alpha=0.3)
    
    # Top view (XY plane)
    ax2 = fig.add_subplot(222)
    ax2.plot(earth_orbit[:, 0]/AU, earth_orbit[:, 1]/AU, 'b-', alpha=0.3, label='Earth')
    ax2.plot(mars_orbit[:, 0]/AU, mars_orbit[:, 1]/AU, 'r-', alpha=0.3, label='Mars')
    ax2.plot(jupiter_orbit[:, 0]/AU, jupiter_orbit[:, 1]/AU, 'orange', alpha=0.3, label='Jupiter')
    ax2.plot(spacecraft_pos[:, 0]/AU, spacecraft_pos[:, 1]/AU, 'g-', linewidth=2, label='Spacecraft')
    
    ax2.scatter([0], [0], c='yellow', s=200, marker='o')
    ax2.scatter([pos_earth_dep[0]/AU], [pos_earth_dep[1]/AU], c='blue', s=100, marker='o')
    ax2.scatter([pos_mars_arr[0]/AU], [pos_mars_arr[1]/AU], c='red', s=100, marker='o')
    ax2.scatter([pos_jup_arr[0]/AU], [pos_jup_arr[1]/AU], c='orange', s=100, marker='o')
    
    ax2.set_xlabel('X (AU)')
    ax2.set_ylabel('Y (AU)')
    ax2.set_title('Top View (XY Plane)', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=8)
    ax2.grid(True, alpha=0.3)
    ax2.axis('equal')
    
    # Delta-v breakdown
    details = calculate_trajectory_details(params)
    dv_data = details['delta_v']
    
    ax3 = fig.add_subplot(223)
    dv_components = [
        dv_data['earth_departure'],
        dv_data['mars_arrival'],
        dv_data['mars_departure'],
        dv_data['jupiter_arrival']
    ]
    dv_labels = ['Earth\nDeparture', 'Mars\nArrival', 'Mars\nDeparture', 'Jupiter\nArrival']
    colors = ['#3498db', '#e74c3c', '#2ecc71', '#f39c12']
    
    bars = ax3.bar(dv_labels, dv_components, color=colors, alpha=0.7, edgecolor='black')
    ax3.axhline(y=dv_data['gravity_assist_benefit'], color='purple', 
                linestyle='--', linewidth=2, label=f'GA Benefit: {dv_data["gravity_assist_benefit"]:.0f} m/s')
    
    ax3.set_ylabel('Delta-v (m/s)')
    ax3.set_title('Delta-v Budget Breakdown', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    
    # Add values on bars
    for bar, val in zip(bars, dv_components):
        height = bar.get_height()
        ax3.text(bar.get_x() + bar.get_width()/2., height,
                f'{val:.0f}', ha='center', va='bottom', fontweight='bold')
    
    # Mission timeline
    ax4 = fig.add_subplot(224)
    time_data = details['times']
    
    timeline_events = [0, time_data['leg1'], time_data['total']]
    timeline_labels = ['Earth\nDeparture', 'Mars\nFlyby', 'Jupiter\nArrival']
    
    ax4.plot(timeline_events, [0]*3, 'o-', markersize=15, linewidth=3, color='#2c3e50')
    
    for i, (x, label) in enumerate(zip(timeline_events, timeline_labels)):
        ax4.text(x, -0.15, label, ha='center', va='top', fontsize=10, fontweight='bold')
        ax4.text(x, 0.15, f'Day {x:.0f}', ha='center', va='bottom', fontsize=9)
    
    # Add leg durations
    ax4.text(time_data['leg1']/2, 0.25, f"{time_data['leg1']:.0f} days", 
             ha='center', fontsize=9, style='italic', color='#e74c3c')
    ax4.text((time_data['leg1'] + time_data['total'])/2, 0.25, f"{time_data['leg2']:.0f} days", 
             ha='center', fontsize=9, style='italic', color='#e74c3c')
    
    ax4.set_xlim(-50, time_data['total'] + 50)
    ax4.set_ylim(-0.3, 0.4)
    ax4.set_xlabel('Mission Time (days)')
    ax4.set_title('Mission Timeline', fontsize=12, fontweight='bold')
    ax4.set_yticks([])
    ax4.grid(True, alpha=0.3, axis='x')
    
    plt.tight_layout()
    plt.savefig('trajectory_optimization.png', dpi=300, bbox_inches='tight')
    plt.show()

def print_mission_summary(params):
    """
    Print detailed mission summary.
    """
    details = calculate_trajectory_details(params)
    
    print("\n" + "="*70)
    print(" "*20 + "MISSION SUMMARY")
    print("="*70)
    
    print("\n📅 TIMELINE:")
    print(f"  Launch Date:              Day {details['times']['t_dep']:.1f}")
    print(f"  Mars Flyby:               Day {details['times']['t_mars']:.1f}")
    print(f"  Jupiter Arrival:          Day {details['times']['t_jup']:.1f}")
    print(f"  Earth-Mars Transit:       {details['times']['leg1']:.1f} days ({details['times']['leg1']/30:.1f} months)")
    print(f"  Mars-Jupiter Transit:     {details['times']['leg2']:.1f} days ({details['times']['leg2']/30:.1f} months)")
    print(f"  Total Mission Duration:   {details['times']['total']:.1f} days ({details['times']['total']/365:.2f} years)")
    
    print("\n🚀 DELTA-V BUDGET:")
    print(f"  Earth Departure:          {details['delta_v']['earth_departure']:,.0f} m/s")
    print(f"  Mars Arrival:             {details['delta_v']['mars_arrival']:,.0f} m/s")
    print(f"  Mars Departure:           {details['delta_v']['mars_departure']:,.0f} m/s")
    print(f"  Jupiter Arrival:          {details['delta_v']['jupiter_arrival']:,.0f} m/s")
    print(f"  Gravity Assist Benefit:   {details['delta_v']['gravity_assist_benefit']:,.0f} m/s (saved)")
    print(f"  ----------------------------------------")
    print(f"  TOTAL DELTA-V:            {details['delta_v']['total']:,.0f} m/s")
    
    print("\n📍 KEY POSITIONS (in AU):")
    for key, pos in details['positions'].items():
        print(f"  {key:20s}: ({pos[0]/AU:6.3f}, {pos[1]/AU:6.3f}, {pos[2]/AU:6.3f})")
    
    print("\n💡 MISSION HIGHLIGHTS:")
    print(f"  • Utilizes Mars gravity assist to save {details['delta_v']['gravity_assist_benefit']:,.0f} m/s")
    print(f"  • Total mission time: {details['times']['total']/365:.2f} years")
    print(f"  • Fuel efficiency optimized through multi-objective optimization")
    print("="*70 + "\n")

# Main execution
if __name__ == "__main__":
    print("╔══════════════════════════════════════════════════════════════════╗")
    print("║     SPACE MISSION TRAJECTORY OPTIMIZATION                        ║")
    print("║     Earth → Mars (Gravity Assist) → Jupiter                     ║")
    print("╚══════════════════════════════════════════════════════════════════╝\n")
    
    # Run optimization
    result = optimize_trajectory()
    
    # Extract optimal parameters
    optimal_params = result.x
    
    print("\n✅ Optimization completed successfully!\n")
    
    # Print detailed mission summary
    print_mission_summary(optimal_params)
    
    # Generate visualizations
    print("📊 Generating trajectory visualizations...\n")
    plot_trajectory_3d(optimal_params)
    
    print("✨ Analysis complete! Check the generated plots above.")
    print("\nThe optimization balanced fuel consumption (delta-v) and mission")
    print("duration to find an efficient trajectory with Mars gravity assist.")

AU = 1.496e11  # Astronomical Unit in meters
G = 6.674e-11  # Gravitational constant
M_sun = 1.989e30  # Solar mass