import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import minimize_scalar
import time
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# Physical constants
G = 6.67430e-11  # Gravitational constant (m³ kg⁻¹ s⁻²)
c = 299792458    # Speed of light (m/s)
M_sun = 1.989e30 # Solar mass (kg)
R_sun = 6.96e8   # Solar radius (m)

class SupernovaSimulation:
    """
    1D Spherically symmetric supernova core collapse simulation
    using the Tolman-Oppenheimer-Volkoff (TOV) equations
    """
    
    def __init__(self, M_core=1.4, R_core=1e4, grid_points=100):
        """
        Initialize the simulation parameters
        
        Parameters:
        M_core: Core mass in solar masses
        R_core: Initial core radius in km
        grid_points: Number of radial grid points
        """
        self.M_core = M_core * M_sun  # Convert to kg
        self.R_core = R_core * 1000   # Convert to meters
        self.N = grid_points
        self.r = np.linspace(1e3, self.R_core, self.N)  # Radial grid (avoid r=0)
        
        # Initialize density profile (polytropic model)
        self.rho_central = 1e15  # kg/m³ (nuclear density)
        self.gamma = 2.0  # Polytropic index
        
    def initial_density_profile(self):
        """
        Set up initial density profile using Lane-Emden equation solution
        For simplicity, we use a polytropic approximation:
        ρ(r) = ρ_c * (1 - (r/R)²)^n
        """
        n = 1.5  # Polytropic index
        density = self.rho_central * np.maximum(0, (1 - (self.r/self.R_core)**2)**n)
        return density
    
    def equation_of_state(self, rho):
        """
        Simplified equation of state: P = K * ρ^γ
        This represents degenerate electron pressure
        """
        K = 1e12  # Polytropic constant (Pa m³ᵞ kg⁻ᵞ)
        return K * rho**self.gamma
    
    def tov_equations(self, y, r):
        """
        Tolman-Oppenheimer-Volkoff equations for stellar structure
        
        dy/dr = [dP/dr, dm/dr]
        
        Where:
        dP/dr = -G*m*ρ/r² * (1 + P/(ρc²)) * (1 + 4πr³P/(mc²)) / (1 - 2Gm/(rc²))
        dm/dr = 4πr²ρ
        """
        P, m = y
        
        # Avoid division by zero
        if r < 1e-10:
            r = 1e-10
            
        # Calculate density from pressure (inverse EOS)
        if P <= 0:
            rho = 0
        else:
            rho = (P / 1e12)**(1/self.gamma)
        
        # TOV equations (simplified for non-relativistic case)
        if m > 0 and rho > 0:
            dP_dr = -G * m * rho / r**2
            dm_dr = 4 * np.pi * r**2 * rho
        else:
            dP_dr = 0
            dm_dr = 0
            
        return [dP_dr, dm_dr]
    
    def solve_structure(self):
        """
        Solve the stellar structure equations
        """
        # Initial conditions
        P_central = self.equation_of_state(self.rho_central)
        m_initial = 4/3 * np.pi * self.r[0]**3 * self.rho_central
        
        y0 = [P_central, m_initial]
        
        # Solve ODEs
        solution = odeint(self.tov_equations, y0, self.r)
        
        self.pressure = solution[:, 0]
        self.mass = solution[:, 1]
        
        # Calculate density from pressure
        self.density = np.zeros_like(self.pressure)
        for i, P in enumerate(self.pressure):
            if P > 0:
                self.density[i] = (P / 1e12)**(1/self.gamma)
        
        return self.r, self.density, self.pressure, self.mass
    
    def collapse_dynamics(self, t_span, dt=1e-6):
        """
        Simulate the collapse dynamics using free-fall approximation
        
        The collapse time scale is approximately:
        t_ff ≈ √(3π/(32Gρ))
        """
        t_ff = np.sqrt(3 * np.pi / (32 * G * self.rho_central))
        t = np.arange(0, t_span, dt)
        
        # Free-fall collapse: R(t) = R₀ * √(1 - t²/t_ff²)
        R_t = self.R_core * np.sqrt(np.maximum(0, 1 - (t/t_ff)**2))
        
        # Central density evolution: ρ(t) = ρ₀ / (1 - t²/t_ff²)^(3/2)
        rho_t = self.rho_central / np.maximum(1e-10, (1 - (t/t_ff)**2)**(3/2))
        
        return t, R_t, rho_t, t_ff

def compute_simulation_accuracy(N_grid):
    """
    Compute simulation accuracy as a function of grid resolution
    """
    sim = SupernovaSimulation(grid_points=N_grid)
    r, density, pressure, mass = sim.solve_structure()
    
    # Calculate total mass and central pressure as accuracy metrics
    total_mass = mass[-1]
    central_pressure = pressure[0]
    
    # Theoretical values for comparison (analytical estimates)
    M_theoretical = sim.M_core
    P_theoretical = sim.equation_of_state(sim.rho_central)
    
    # Relative errors
    mass_error = abs(total_mass - M_theoretical) / M_theoretical
    pressure_error = abs(central_pressure - P_theoretical) / P_theoretical
    
    # Combined accuracy metric (lower is better)
    accuracy = np.sqrt(mass_error**2 + pressure_error**2)
    
    return accuracy, total_mass, central_pressure

def compute_computational_cost(N_grid):
    """
    Estimate computational cost (execution time) for given grid resolution
    """
    start_time = time.time()
    
    sim = SupernovaSimulation(grid_points=N_grid)
    sim.solve_structure()
    
    end_time = time.time()
    
    return end_time - start_time

def optimize_grid_resolution():
    """
    Find optimal grid resolution that balances accuracy and computational cost
    
    Objective function: f(N) = α * accuracy(N) + β * cost(N)
    where α and β are weighting factors
    """
    N_values = np.arange(50, 500, 25)
    accuracies = []
    costs = []
    
    print("Optimizing grid resolution...")
    print("N_grid\tAccuracy\tCost(s)\tObjective")
    print("-" * 45)
    
    for N in N_values:
        accuracy, _, _ = compute_simulation_accuracy(N)
        cost = compute_computational_cost(N)
        
        accuracies.append(accuracy)
        costs.append(cost)
        
        # Normalize and combine metrics
        # α = 1.0 (accuracy weight), β = 0.1 (cost weight)
        objective = accuracy + 0.1 * cost
        
        print(f"{N}\t{accuracy:.6f}\t{cost:.4f}\t{objective:.6f}")
    
    accuracies = np.array(accuracies)
    costs = np.array(costs)
    
    # Find optimal N
    objectives = accuracies + 0.1 * costs
    optimal_idx = np.argmin(objectives)
    optimal_N = N_values[optimal_idx]
    
    print(f"\nOptimal grid resolution: N = {optimal_N}")
    print(f"Accuracy: {accuracies[optimal_idx]:.6f}")
    print(f"Cost: {costs[optimal_idx]:.4f} seconds")
    
    return N_values, accuracies, costs, optimal_N

# Run optimization
N_values, accuracies, costs, optimal_N = optimize_grid_resolution()

# Create detailed simulation with optimal parameters
print(f"\nRunning detailed simulation with N = {optimal_N}...")
sim_optimal = SupernovaSimulation(grid_points=optimal_N)
r, density, pressure, mass = sim_optimal.solve_structure()

# Simulate collapse dynamics
t, R_t, rho_t, t_ff = sim_optimal.collapse_dynamics(t_span=0.001, dt=1e-6)

print(f"Free-fall collapse time: {t_ff:.6f} seconds ({t_ff*1000:.3f} ms)")

# Create comprehensive visualization
fig = plt.figure(figsize=(16, 12))

# 1. Optimization Results
ax1 = plt.subplot(2, 3, 1)
plt.plot(N_values, accuracies, 'b-o', label='Accuracy', linewidth=2, markersize=6)
plt.axvline(optimal_N, color='red', linestyle='--', alpha=0.7, label=f'Optimal N={optimal_N}')
plt.xlabel('Grid Points (N)')
plt.ylabel('Accuracy Metric')
plt.title('Simulation Accuracy vs Grid Resolution')
plt.legend()
plt.grid(True, alpha=0.3)

ax2 = plt.subplot(2, 3, 2)
plt.plot(N_values, costs, 'g-s', label='Computational Cost', linewidth=2, markersize=6)
plt.axvline(optimal_N, color='red', linestyle='--', alpha=0.7, label=f'Optimal N={optimal_N}')
plt.xlabel('Grid Points (N)')
plt.ylabel('Execution Time (seconds)')
plt.title('Computational Cost vs Grid Resolution')
plt.legend()
plt.grid(True, alpha=0.3)

# 2. Stellar Structure
ax3 = plt.subplot(2, 3, 3)
plt.loglog(r/1000, density/1e15, 'b-', linewidth=3, label='Density Profile')
plt.xlabel('Radius (km)')
plt.ylabel('Density (10¹⁵ kg/m³)')
plt.title('Initial Stellar Structure')
plt.grid(True, alpha=0.3)
plt.legend()

ax4 = plt.subplot(2, 3, 4)
plt.semilogy(r/1000, pressure/1e30, 'r-', linewidth=3, label='Pressure Profile')
plt.xlabel('Radius (km)')
plt.ylabel('Pressure (10³⁰ Pa)')
plt.title('Pressure Distribution')
plt.grid(True, alpha=0.3)
plt.legend()

# 3. Collapse Dynamics
ax5 = plt.subplot(2, 3, 5)
plt.plot(t*1000, R_t/1000, 'purple', linewidth=3, label='Core Radius')
plt.xlabel('Time (ms)')
plt.ylabel('Core Radius (km)')
plt.title('Core Collapse Evolution')
plt.grid(True, alpha=0.3)
plt.legend()

ax6 = plt.subplot(2, 3, 6)
plt.semilogy(t*1000, rho_t/1e15, 'orange', linewidth=3, label='Central Density')
plt.xlabel('Time (ms)')
plt.ylabel('Central Density (10¹⁵ kg/m³)')
plt.title('Density Evolution During Collapse')
plt.grid(True, alpha=0.3)
plt.legend()

plt.tight_layout()
plt.show()

# Additional analysis: Mass-Radius relation
print(f"\nStellar Structure Analysis:")
print(f"Total mass: {mass[-1]/M_sun:.3f} M☉")
print(f"Core radius: {r[-1]/1000:.1f} km")
print(f"Central density: {density[0]/1e15:.2f} × 10¹⁵ kg/m³")
print(f"Central pressure: {pressure[0]/1e30:.2f} × 10³⁰ Pa")

# Performance comparison
print(f"\nPerformance Analysis:")
print(f"Optimal configuration saves {((costs[-1] - costs[optimal_N//25-2])/costs[-1]*100):.1f}% computation time")
print(f"while maintaining accuracy within {(accuracies[optimal_N//25-2]*100):.3f}% relative error")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import math

# Planetary and solar system constants
class CelestialBody:
    def __init__(self, name, mu, radius, orbital_radius, orbital_period):
        self.name = name
        self.mu = mu  # Gravitational parameter (km^3/s^2)
        self.radius = radius  # Planet radius (km)
        self.orbital_radius = orbital_radius  # Distance from Sun (AU)
        self.orbital_period = orbital_period  # Orbital period (years)
        # Calculate orbital velocity (avoid division by zero for Sun)
        if orbital_period > 0:
            self.orbital_velocity = 2 * np.pi * orbital_radius * 149597870.7 / (orbital_period * 365.25 * 24 * 3600)  # km/s
        else:
            self.orbital_velocity = 0  # Sun doesn't orbit anything

# Define celestial bodies
sun = CelestialBody("Sun", 1.327e11, 696340, 0, 0)
earth = CelestialBody("Earth", 3.986e5, 6371, 1.0, 1.0)
jupiter = CelestialBody("Jupiter", 1.267e8, 69911, 5.204, 11.86)
saturn = CelestialBody("Saturn", 3.794e7, 58232, 9.537, 29.46)

AU = 149597870.7  # km

def hohmann_transfer_dv(r1, r2, mu_central):
    """Calculate delta-v for Hohmann transfer between two circular orbits"""
    a_transfer = (r1 + r2) / 2
    v1 = np.sqrt(mu_central / r1)
    v2 = np.sqrt(mu_central / r2)
    v_transfer_1 = np.sqrt(mu_central * (2/r1 - 1/a_transfer))
    v_transfer_2 = np.sqrt(mu_central * (2/r2 - 1/a_transfer))
    
    dv1 = abs(v_transfer_1 - v1)
    dv2 = abs(v2 - v_transfer_2)
    transfer_time = np.pi * np.sqrt(a_transfer**3 / mu_central)
    
    return dv1 + dv2, transfer_time

def gravity_assist_deflection(v_inf, rp, mu_planet):
    """Calculate deflection angle for gravity assist"""
    # Deflection angle in radians
    if v_inf == 0:
        return 0
    
    delta = 2 * np.arcsin(1 / (1 + rp * v_inf**2 / mu_planet))
    return delta

def velocity_change_magnitude(v_inf, delta):
    """Calculate magnitude of velocity change from gravity assist"""
    return 2 * v_inf * np.sin(delta / 2)

def calculate_trajectory_with_gravity_assist(departure_dv, jupiter_rp, arrival_dv):
    """
    Calculate complete Earth-Jupiter-Saturn trajectory with gravity assist
    Returns total delta-v and transfer times
    """
    
    # Earth to Jupiter transfer
    r_earth = earth.orbital_radius * AU
    r_jupiter = jupiter.orbital_radius * AU
    
    # Calculate hyperbolic excess velocity at Jupiter
    a_ej = (r_earth + r_jupiter) / 2
    v_ej_earth = np.sqrt(sun.mu * (2/r_earth - 1/a_ej))
    v_earth_orbital = earth.orbital_velocity
    
    # Excess velocity relative to Earth
    v_inf_earth = abs(v_ej_earth - v_earth_orbital)
    
    # At Jupiter arrival
    v_ej_jupiter = np.sqrt(sun.mu * (2/r_jupiter - 1/a_ej))
    v_jupiter_orbital = jupiter.orbital_velocity
    v_inf_jupiter = abs(v_ej_jupiter - v_jupiter_orbital)
    
    # Gravity assist at Jupiter
    delta_jupiter = gravity_assist_deflection(v_inf_jupiter, jupiter_rp, jupiter.mu)
    dv_jupiter = velocity_change_magnitude(v_inf_jupiter, delta_jupiter)
    
    # Assume optimal deflection gives us the desired velocity for Jupiter-Saturn transfer
    r_saturn = saturn.orbital_radius * AU
    a_js = (r_jupiter + r_saturn) / 2
    v_js_jupiter = np.sqrt(sun.mu * (2/r_jupiter - 1/a_js))
    
    # Transfer times
    t_ej = np.pi * np.sqrt(a_ej**3 / sun.mu) / (24 * 3600)  # days
    t_js = np.pi * np.sqrt(a_js**3 / sun.mu) / (24 * 3600)  # days
    
    total_dv = departure_dv + arrival_dv  # Simplified - gravity assist provides the Jupiter departure velocity
    total_time = t_ej + t_js
    
    return total_dv, total_time, t_ej, t_js, v_inf_jupiter, delta_jupiter, dv_jupiter

def optimize_gravity_assist():
    """Optimize gravity assist parameters"""
    
    def objective(params):
        departure_dv, jupiter_rp_factor, arrival_dv = params
        jupiter_rp = jupiter.radius * jupiter_rp_factor  # Minimum safe distance factor
        
        if jupiter_rp_factor < 1.1:  # Safety constraint
            return 1e10
            
        try:
            total_dv, _, _, _, _, _, _ = calculate_trajectory_with_gravity_assist(
                departure_dv, jupiter_rp, arrival_dv)
            return total_dv
        except:
            return 1e10
    
    # Initial guess
    x0 = [3.0, 2.0, 3.0]  # departure_dv, rp_factor, arrival_dv
    bounds = [(0.5, 10), (1.1, 5), (0.5, 10)]
    
    result = minimize(objective, x0, bounds=bounds, method='L-BFGS-B')
    return result

# Calculate direct Earth-Saturn transfer (comparison)
r_earth = earth.orbital_radius * AU
r_saturn = saturn.orbital_radius * AU
direct_dv, direct_time = hohmann_transfer_dv(r_earth, r_saturn, sun.mu)
direct_time_days = direct_time / (24 * 3600)

print("=== GRAVITY ASSIST TRAJECTORY OPTIMIZATION ===")
print()
print("Direct Earth-Saturn Transfer (Hohmann):")
print(f"  Total ΔV: {direct_dv:.2f} km/s")
print(f"  Transfer time: {direct_time_days:.1f} days ({direct_time_days/365.25:.2f} years)")
print()

# Optimize gravity assist trajectory
print("Optimizing gravity assist trajectory...")
opt_result = optimize_gravity_assist()

if opt_result.success:
    opt_departure_dv, opt_rp_factor, opt_arrival_dv = opt_result.x
    opt_jupiter_rp = jupiter.radius * opt_rp_factor
    
    total_dv, total_time, t_ej, t_js, v_inf_jupiter, delta_jupiter, dv_jupiter = \
        calculate_trajectory_with_gravity_assist(opt_departure_dv, opt_jupiter_rp, opt_arrival_dv)
    
    print("Optimized Gravity Assist Trajectory:")
    print(f"  Departure ΔV (Earth): {opt_departure_dv:.2f} km/s")
    print(f"  Jupiter periapsis: {opt_jupiter_rp:.0f} km ({opt_rp_factor:.2f} × R_Jupiter)")
    print(f"  Jupiter deflection angle: {np.degrees(delta_jupiter):.1f} degrees")
    print(f"  Velocity change at Jupiter: {dv_jupiter:.2f} km/s")
    print(f"  Arrival ΔV (Saturn): {opt_arrival_dv:.2f} km/s")
    print(f"  Total ΔV: {total_dv:.2f} km/s")
    print(f"  Earth-Jupiter time: {t_ej:.1f} days ({t_ej/365.25:.2f} years)")
    print(f"  Jupiter-Saturn time: {t_js:.1f} days ({t_js/365.25:.2f} years)")
    print(f"  Total transfer time: {total_time:.1f} days ({total_time/365.25:.2f} years)")
    print()
    print(f"ΔV Savings: {direct_dv - total_dv:.2f} km/s ({100*(direct_dv - total_dv)/direct_dv:.1f}%)")
else:
    print("Optimization failed!")

# Detailed gravity assist analysis
print("\n=== GRAVITY ASSIST PHYSICS ANALYSIS ===")
print()

# Calculate different periapsis distances and their effects
rp_factors = np.linspace(1.1, 4.0, 20)
deflection_angles = []
velocity_changes = []
rp_distances = []

v_inf_analysis = 8.0  # km/s - typical hyperbolic excess velocity

for rp_factor in rp_factors:
    rp = jupiter.radius * rp_factor
    delta = gravity_assist_deflection(v_inf_analysis, rp, jupiter.mu)
    dv = velocity_change_magnitude(v_inf_analysis, delta)
    
    rp_distances.append(rp)
    deflection_angles.append(np.degrees(delta))
    velocity_changes.append(dv)

print(f"Analysis for v∞ = {v_inf_analysis} km/s at Jupiter:")
print(f"  Minimum safe periapsis: {jupiter.radius * 1.1:.0f} km")
print(f"  Maximum deflection: {max(deflection_angles):.1f} degrees")
print(f"  Maximum velocity change: {max(velocity_changes):.2f} km/s")

# Create visualization
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# Plot 1: Deflection angle vs periapsis distance
ax1.plot(np.array(rp_distances)/jupiter.radius, deflection_angles, 'b-', linewidth=2)
ax1.set_xlabel('Periapsis Distance (Jupiter Radii)')
ax1.set_ylabel('Deflection Angle (degrees)')
ax1.set_title('Gravity Assist Deflection vs Periapsis Distance')
ax1.grid(True, alpha=0.3)
ax1.axvline(x=1.1, color='r', linestyle='--', alpha=0.7, label='Minimum Safe Distance')
ax1.legend()

# Plot 2: Velocity change vs periapsis distance
ax2.plot(np.array(rp_distances)/jupiter.radius, velocity_changes, 'g-', linewidth=2)
ax2.set_xlabel('Periapsis Distance (Jupiter Radii)')
ax2.set_ylabel('Velocity Change (km/s)')
ax2.set_title('Velocity Change from Gravity Assist')
ax2.grid(True, alpha=0.3)
ax2.axvline(x=1.1, color='r', linestyle='--', alpha=0.7, label='Minimum Safe Distance')
ax2.legend()

# Plot 3: Trajectory comparison
labels = ['Direct Transfer', 'Gravity Assist']
delta_vs = [direct_dv, total_dv if opt_result.success else 0]
times = [direct_time_days/365.25, total_time/365.25 if opt_result.success else 0]

x = np.arange(len(labels))
width = 0.35

ax3.bar(x - width/2, delta_vs, width, label='Total ΔV (km/s)', color='skyblue')
ax3.bar(x + width/2, times, width, label='Transfer Time (years)', color='lightcoral')
ax3.set_xlabel('Transfer Method')
ax3.set_ylabel('Value')
ax3.set_title('Direct vs Gravity Assist Comparison')
ax3.set_xticks(x)
ax3.set_xticklabels(labels)
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Orbital diagram (simplified)
theta = np.linspace(0, 2*np.pi, 100)

# Orbital paths
earth_orbit = earth.orbital_radius
jupiter_orbit = jupiter.orbital_radius
saturn_orbit = saturn.orbital_radius

ax4.plot(earth_orbit * np.cos(theta), earth_orbit * np.sin(theta), 'b-', label='Earth Orbit', linewidth=2)
ax4.plot(jupiter_orbit * np.cos(theta), jupiter_orbit * np.sin(theta), 'orange', label='Jupiter Orbit', linewidth=2)
ax4.plot(saturn_orbit * np.cos(theta), saturn_orbit * np.sin(theta), 'brown', label='Saturn Orbit', linewidth=2)

# Planets
ax4.plot(earth_orbit, 0, 'bo', markersize=8, label='Earth')
ax4.plot(jupiter_orbit, 0, 'o', color='orange', markersize=12, label='Jupiter')
ax4.plot(-saturn_orbit, 0, 'o', color='brown', markersize=10, label='Saturn')
ax4.plot(0, 0, 'yo', markersize=15, label='Sun')

# Transfer trajectory (simplified)
transfer1_x = np.linspace(earth_orbit, jupiter_orbit, 50)
transfer1_y = 0.3 * np.sin(np.pi * (transfer1_x - earth_orbit) / (jupiter_orbit - earth_orbit))
transfer2_x = np.linspace(jupiter_orbit, -saturn_orbit, 50)
transfer2_y = -0.2 * np.sin(np.pi * (transfer2_x - jupiter_orbit) / (-saturn_orbit - jupiter_orbit))

ax4.plot(transfer1_x, transfer1_y, 'r--', linewidth=2, label='Transfer Trajectory')
ax4.plot(transfer2_x, transfer2_y, 'r--', linewidth=2)

ax4.set_xlabel('Distance (AU)')
ax4.set_ylabel('Distance (AU)')
ax4.set_title('Earth-Jupiter-Saturn Transfer Trajectory')
ax4.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
ax4.grid(True, alpha=0.3)
ax4.set_aspect('equal')
ax4.set_xlim(-11, 7)
ax4.set_ylim(-6, 6)

plt.tight_layout()
plt.show()

# Additional analysis: Effect of hyperbolic excess velocity
print("\n=== HYPERBOLIC EXCESS VELOCITY ANALYSIS ===")
v_inf_range = np.linspace(2, 15, 8)
rp_fixed = jupiter.radius * 1.5  # Fixed periapsis at 1.5 Jupiter radii

print(f"Analysis at fixed periapsis = {rp_fixed:.0f} km ({rp_fixed/jupiter.radius:.1f} R_Jupiter):")
print("v∞ (km/s) | Deflection (°) | ΔV (km/s)")
print("-" * 40)

for v_inf in v_inf_range:
    delta = gravity_assist_deflection(v_inf, rp_fixed, jupiter.mu)
    dv = velocity_change_magnitude(v_inf, delta)
    print(f"   {v_inf:5.1f}   |    {np.degrees(delta):6.1f}    |   {dv:5.2f}")

def gravity_assist_deflection(v_inf, rp, mu_planet):
    delta = 2 * np.arcsin(1 / (1 + rp * v_inf**2 / mu_planet))
    return delta

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

# Constants
MU_EARTH = 3.986004418e14  # Earth's gravitational parameter (m^3/s^2)
R_EARTH = 6.371e6  # Earth's radius (m)
LEO_ALTITUDE = 300e3  # LEO altitude (m)
GEO_ALTITUDE = 35786e3  # GEO altitude (m)

# Orbital radii
r_leo = R_EARTH + LEO_ALTITUDE
r_geo = R_EARTH + GEO_ALTITUDE

print(f"LEO radius: {r_leo/1000:.0f} km")
print(f"GEO radius: {r_geo/1000:.0f} km")

def circular_velocity(r):
    """Calculate circular orbital velocity at radius r"""
    return np.sqrt(MU_EARTH / r)

def hohmann_delta_v(r1, r2):
    """
    Calculate delta-V for Hohmann transfer between two circular orbits
    
    Parameters:
    r1: radius of initial orbit (m)
    r2: radius of target orbit (m)
    
    Returns:
    dv1: delta-V at initial orbit (m/s)
    dv2: delta-V at target orbit (m/s)
    total_dv: total delta-V (m/s)
    """
    # Velocities in circular orbits
    v1_circular = circular_velocity(r1)
    v2_circular = circular_velocity(r2)
    
    # Semi-major axis of transfer ellipse
    a_transfer = (r1 + r2) / 2
    
    # Velocities at perigee and apogee of transfer orbit
    v1_transfer = np.sqrt(MU_EARTH * (2/r1 - 1/a_transfer))
    v2_transfer = np.sqrt(MU_EARTH * (2/r2 - 1/a_transfer))
    
    # Delta-V calculations
    dv1 = v1_transfer - v1_circular  # Acceleration at perigee
    dv2 = v2_circular - v2_transfer  # Acceleration at apogee
    
    total_dv = abs(dv1) + abs(dv2)
    
    return dv1, dv2, total_dv

# Calculate Hohmann transfer for LEO to GEO
dv1_hohmann, dv2_hohmann, total_dv_hohmann = hohmann_delta_v(r_leo, r_geo)

print(f"\nHohmann Transfer LEO to GEO:")
print(f"First burn (ΔV₁): {dv1_hohmann:.2f} m/s")
print(f"Second burn (ΔV₂): {dv2_hohmann:.2f} m/s")
print(f"Total ΔV: {total_dv_hohmann:.2f} m/s")

# Now let's consider a more complex optimization problem:
# Multi-impulse transfer with intermediate orbit

def multi_impulse_transfer(params):
    """
    Calculate total delta-V for a multi-impulse transfer
    params: [r_intermediate] - radius of intermediate orbit
    """
    r_int = params[0]
    
    # Ensure intermediate orbit is between LEO and GEO
    if r_int <= r_leo or r_int >= r_geo:
        return 1e10  # Penalty for invalid orbit
    
    # First transfer: LEO to intermediate orbit
    _, _, dv_total_1 = hohmann_delta_v(r_leo, r_int)
    
    # Second transfer: intermediate orbit to GEO
    _, _, dv_total_2 = hohmann_delta_v(r_int, r_geo)
    
    return dv_total_1 + dv_total_2

# Optimize intermediate orbit radius
bounds = [(r_leo + 1000, r_geo - 1000)]  # Bounds for intermediate orbit radius
result = minimize(multi_impulse_transfer, 
                 x0=[(r_leo + r_geo) / 2], 
                 bounds=bounds, 
                 method='L-BFGS-B')

r_optimal = result.x[0]
optimal_dv = result.fun

print(f"\nOptimal Multi-Impulse Transfer:")
print(f"Optimal intermediate orbit radius: {r_optimal/1000:.0f} km")
print(f"Optimal intermediate orbit altitude: {(r_optimal - R_EARTH)/1000:.0f} km")
print(f"Total ΔV: {optimal_dv:.2f} m/s")
print(f"Fuel savings vs Hohmann: {total_dv_hohmann - optimal_dv:.2f} m/s ({((total_dv_hohmann - optimal_dv)/total_dv_hohmann)*100:.2f}%)")

# Let's analyze the transfer in more detail
def analyze_transfer_sequence(r1, r_int, r3):
    """Analyze a complete transfer sequence"""
    dv1_1, dv1_2, total_dv_1 = hohmann_delta_v(r1, r_int)
    dv2_1, dv2_2, total_dv_2 = hohmann_delta_v(r_int, r3)
    
    return {
        'first_transfer': {'dv1': dv1_1, 'dv2': dv1_2, 'total': total_dv_1},
        'second_transfer': {'dv1': dv2_1, 'dv2': dv2_2, 'total': total_dv_2},
        'total_dv': total_dv_1 + total_dv_2
    }

optimal_analysis = analyze_transfer_sequence(r_leo, r_optimal, r_geo)
hohmann_analysis = analyze_transfer_sequence(r_leo, r_geo, r_geo)  # Direct transfer

print(f"\nDetailed Analysis:")
print(f"Optimal transfer:")
print(f"  LEO to intermediate: {optimal_analysis['first_transfer']['total']:.2f} m/s")
print(f"  Intermediate to GEO: {optimal_analysis['second_transfer']['total']:.2f} m/s")
print(f"Direct Hohmann transfer: {total_dv_hohmann:.2f} m/s")

# Create visualization of different transfer options
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))

# Plot 1: Orbital geometry
theta = np.linspace(0, 2*np.pi, 100)

# Earth
earth_x = R_EARTH * np.cos(theta)
earth_y = R_EARTH * np.sin(theta)

# Orbits
leo_x = r_leo * np.cos(theta)
leo_y = r_leo * np.sin(theta)
geo_x = r_geo * np.cos(theta)
geo_y = r_geo * np.sin(theta)
int_x = r_optimal * np.cos(theta)
int_y = r_optimal * np.sin(theta)

ax1.fill(earth_x/1000, earth_y/1000, color='blue', alpha=0.7, label='Earth')
ax1.plot(leo_x/1000, leo_y/1000, 'g-', linewidth=2, label='LEO')
ax1.plot(int_x/1000, int_y/1000, 'orange', linewidth=2, label='Intermediate Orbit')
ax1.plot(geo_x/1000, geo_y/1000, 'r-', linewidth=2, label='GEO')

# Transfer ellipses
# Hohmann transfer ellipse
a_hohmann = (r_leo + r_geo) / 2
b_hohmann = np.sqrt(r_leo * r_geo)
ellipse_hohmann = patches.Ellipse((0, 0), 2*a_hohmann/1000, 2*b_hohmann/1000, 
                                 fill=False, linestyle='--', color='purple', linewidth=2)
ax1.add_patch(ellipse_hohmann)

# Optimal transfer ellipses
a1_opt = (r_leo + r_optimal) / 2
b1_opt = np.sqrt(r_leo * r_optimal)
ellipse1_opt = patches.Ellipse((0, 0), 2*a1_opt/1000, 2*b1_opt/1000, 
                              fill=False, linestyle=':', color='cyan', linewidth=2)
ax1.add_patch(ellipse1_opt)

a2_opt = (r_optimal + r_geo) / 2
b2_opt = np.sqrt(r_optimal * r_geo)
ellipse2_opt = patches.Ellipse((0, 0), 2*a2_opt/1000, 2*b2_opt/1000, 
                              fill=False, linestyle=':', color='magenta', linewidth=2)
ax1.add_patch(ellipse2_opt)

ax1.set_xlim(-50000, 50000)
ax1.set_ylim(-50000, 50000)
ax1.set_xlabel('Distance (km)')
ax1.set_ylabel('Distance (km)')
ax1.set_title('Orbital Transfer Geometry')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_aspect('equal')

# Plot 2: Delta-V comparison
transfers = ['Direct Hohmann', 'Optimal Multi-Impulse']
delta_vs = [total_dv_hohmann, optimal_dv]
colors = ['red', 'green']

bars = ax2.bar(transfers, delta_vs, color=colors, alpha=0.7)
ax2.set_ylabel('Total ΔV (m/s)')
ax2.set_title('Fuel Consumption Comparison')
ax2.grid(True, alpha=0.3)

# Add value labels on bars
for bar, dv in zip(bars, delta_vs):
    height = bar.get_height()
    ax2.text(bar.get_x() + bar.get_width()/2., height + 20,
             f'{dv:.0f} m/s', ha='center', va='bottom', fontweight='bold')

# Plot 3: Optimization landscape
r_range = np.linspace(r_leo + 1000, r_geo - 1000, 200)
dv_values = [multi_impulse_transfer([r]) for r in r_range]

ax3.plot(r_range/1000, dv_values, 'b-', linewidth=2)
ax3.axvline(r_optimal/1000, color='red', linestyle='--', linewidth=2, 
           label=f'Optimal: {r_optimal/1000:.0f} km')
ax3.axhline(total_dv_hohmann, color='orange', linestyle='--', linewidth=2, 
           label=f'Direct Hohmann: {total_dv_hohmann:.0f} m/s')
ax3.set_xlabel('Intermediate Orbit Radius (km)')
ax3.set_ylabel('Total ΔV (m/s)')
ax3.set_title('Optimization Landscape')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Transfer timeline and velocity profile
# Calculate orbital periods
def orbital_period(r):
    return 2 * np.pi * np.sqrt(r**3 / MU_EARTH)

T_leo = orbital_period(r_leo) / 3600  # Convert to hours
T_int = orbital_period(r_optimal) / 3600
T_geo = orbital_period(r_geo) / 3600
T_transfer1 = orbital_period((r_leo + r_optimal)/2) / 2 / 3600  # Half period of transfer orbit
T_transfer2 = orbital_period((r_optimal + r_geo)/2) / 2 / 3600

print(f"\nOrbital Periods:")
print(f"LEO period: {T_leo:.2f} hours")
print(f"Intermediate orbit period: {T_int:.2f} hours")
print(f"GEO period: {T_geo:.2f} hours")
print(f"First transfer time: {T_transfer1:.2f} hours")
print(f"Second transfer time: {T_transfer2:.2f} hours")

# Velocity profile
altitudes = np.array([LEO_ALTITUDE, (r_optimal - R_EARTH), GEO_ALTITUDE]) / 1000
velocities = np.array([circular_velocity(r_leo), circular_velocity(r_optimal), 
                      circular_velocity(r_geo)]) / 1000

ax4.plot(altitudes, velocities, 'bo-', linewidth=2, markersize=8, label='Circular Velocities')
ax4.set_xlabel('Altitude (km)')
ax4.set_ylabel('Orbital Velocity (km/s)')
ax4.set_title('Velocity Profile vs Altitude')
ax4.grid(True, alpha=0.3)
ax4.legend()

# Add annotations
for i, (alt, vel) in enumerate(zip(altitudes, velocities)):
    orbit_names = ['LEO', 'Intermediate', 'GEO']
    ax4.annotate(f'{orbit_names[i]}\n({alt:.0f} km, {vel:.2f} km/s)', 
                xy=(alt, vel), xytext=(10, 10), textcoords='offset points',
                bbox=dict(boxstyle='round,pad=0.3', facecolor='yellow', alpha=0.7))

plt.tight_layout()
plt.show()

# Additional analysis: Sensitivity to intermediate orbit selection
print(f"\nSensitivity Analysis:")
r_test_values = np.linspace(r_leo + 1000, r_geo - 1000, 20)
dv_sensitivity = []

for r_test in r_test_values:
    dv_test = multi_impulse_transfer([r_test])
    dv_sensitivity.append(dv_test)
    if abs(r_test - r_optimal) < 1000000:  # Within 1000 km of optimal
        print(f"  Radius: {r_test/1000:.0f} km, ΔV: {dv_test:.2f} m/s, "
              f"Penalty: {dv_test - optimal_dv:.2f} m/s")

# Create final summary plot
fig, ax = plt.subplots(figsize=(12, 8))

# Plot the sensitivity curve
ax.plot(r_test_values/1000, dv_sensitivity, 'b-', linewidth=3, label='Multi-impulse transfer')
ax.axhline(total_dv_hohmann, color='red', linestyle='--', linewidth=2, 
          label=f'Direct Hohmann ({total_dv_hohmann:.0f} m/s)')
ax.axvline(r_optimal/1000, color='green', linestyle=':', linewidth=2, 
          label=f'Optimal intermediate orbit')

# Mark the optimal point
ax.plot(r_optimal/1000, optimal_dv, 'go', markersize=10, label=f'Optimal point ({optimal_dv:.0f} m/s)')

ax.set_xlabel('Intermediate Orbit Radius (km)', fontsize=12)
ax.set_ylabel('Total ΔV (m/s)', fontsize=12)
ax.set_title('Fuel Optimization for LEO to GEO Transfer', fontsize=14, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)

# Add savings annotation
savings = total_dv_hohmann - optimal_dv
ax.annotate(f'Fuel Savings: {savings:.0f} m/s ({(savings/total_dv_hohmann)*100:.1f}%)', 
           xy=(r_optimal/1000, optimal_dv), xytext=(30000, optimal_dv + 200),
           arrowprops=dict(arrowstyle='->', color='red', lw=2),
           bbox=dict(boxstyle='round,pad=0.5', facecolor='lightgreen', alpha=0.8),
           fontsize=12, fontweight='bold')

plt.tight_layout()
plt.show()

print(f"\n" + "="*60)
print(f"OPTIMIZATION RESULTS SUMMARY")
print(f"="*60)
print(f"Problem: Minimize fuel consumption for LEO to GEO transfer")
print(f"Initial orbit (LEO): {LEO_ALTITUDE/1000:.0f} km altitude")
print(f"Target orbit (GEO): {GEO_ALTITUDE/1000:.0f} km altitude")
print(f"")
print(f"Direct Hohmann Transfer:")
print(f"  Total ΔV: {total_dv_hohmann:.2f} m/s")
print(f"")
print(f"Optimized Multi-Impulse Transfer:")
print(f"  Optimal intermediate altitude: {(r_optimal-R_EARTH)/1000:.0f} km")
print(f"  Total ΔV: {optimal_dv:.2f} m/s")
print(f"  Fuel savings: {savings:.2f} m/s ({(savings/total_dv_hohmann)*100:.2f}%)")
print(f"="*60)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import optimize
from scipy.optimize import minimize
from scipy.stats import norm, multivariate_normal
import warnings
warnings.filterwarnings('ignore')

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

class NonlinearPortfolioOptimizer:
    def __init__(self, n_assets=5, n_simulations=10000):
        self.n_assets = n_assets
        self.n_simulations = n_simulations
        self.setup_market_data()
        
    def setup_market_data(self):
        """Generate synthetic market data with nonlinear correlations"""
        # Expected returns (annualized)
        self.mu = np.array([0.08, 0.12, 0.15, 0.10, 0.06])
        
        # Asset names
        self.asset_names = ['Stock A', 'Stock B', 'Stock C', 'Bond A', 'Bond B']
        
        # Base correlation matrix
        base_corr = np.array([
            [1.0, 0.3, 0.5, -0.2, -0.1],
            [0.3, 1.0, 0.7, -0.1, 0.0],
            [0.5, 0.7, 1.0, 0.1, 0.2],
            [-0.2, -0.1, 0.1, 1.0, 0.8],
            [-0.1, 0.0, 0.2, 0.8, 1.0]
        ])
        
        # Volatilities (annualized)
        self.sigma = np.array([0.20, 0.25, 0.30, 0.08, 0.06])
        
        # Create covariance matrix
        self.cov_matrix = np.outer(self.sigma, self.sigma) * base_corr
        
        # Transaction cost coefficients
        self.transaction_costs = np.array([0.002, 0.003, 0.004, 0.001, 0.001])
        
        # Generate nonlinear correlation scenario returns
        self.generate_nonlinear_returns()
        
    def generate_nonlinear_returns(self):
        """Generate returns with nonlinear correlation structure"""
        # Generate base multivariate normal returns
        base_returns = multivariate_normal.rvs(
            mean=self.mu/252,  # Daily returns
            cov=self.cov_matrix/252,
            size=self.n_simulations
        )
        
        # Add nonlinear correlation effects
        # Regime-dependent correlations based on market stress
        market_factor = base_returns[:, :3].mean(axis=1)  # Average stock performance
        stress_indicator = market_factor < np.percentile(market_factor, 25)
        
        # Increase correlations during stress periods
        stressed_returns = base_returns.copy()
        for i in range(self.n_simulations):
            if stress_indicator[i]:
                # Increase correlations between stocks during stress
                factor = 1.5
                stressed_returns[i, :3] = (
                    base_returns[i, :3] * 0.7 + 
                    market_factor[i] * 0.3 * factor
                )
        
        self.return_scenarios = stressed_returns
        
    def nonlinear_transaction_cost(self, weights):
        """Calculate nonlinear transaction costs"""
        return np.sum(self.transaction_costs * np.power(weights, 1.5))
    
    def portfolio_var(self, weights, confidence_level=0.05):
        """Calculate portfolio VaR using Monte Carlo simulation"""
        portfolio_returns = np.dot(self.return_scenarios, weights)
        var = -np.percentile(portfolio_returns, confidence_level * 100)
        return var
    
    def portfolio_expected_return(self, weights):
        """Calculate expected portfolio return"""
        return np.dot(weights, self.mu)
    
    def objective_function(self, weights):
        """Objective function: maximize expected return minus transaction costs"""
        expected_return = self.portfolio_expected_return(weights)
        transaction_cost = self.nonlinear_transaction_cost(weights)
        return -(expected_return - transaction_cost)  # Negative for minimization
    
    def var_constraint(self, weights, max_var):
        """VaR constraint function"""
        return max_var - self.portfolio_var(weights)
    
    def optimize_portfolio(self, max_var=0.02):
        """Solve the portfolio optimization problem"""
        # Initial guess (equal weights)
        x0 = np.ones(self.n_assets) / self.n_assets
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda x: np.sum(x) - 1.0},  # Sum to 1
            {'type': 'ineq', 'fun': lambda x: self.var_constraint(x, max_var)}  # VaR constraint
        ]
        
        # Bounds (long-only)
        bounds = [(0, 1) for _ in range(self.n_assets)]
        
        # Solve optimization
        result = optimize.minimize(
            self.objective_function,
            x0,
            method='SLSQP',
            bounds=bounds,
            constraints=constraints,
            options={'disp': True, 'maxiter': 1000}
        )
        
        return result
    
    def efficient_frontier(self, var_levels):
        """Generate efficient frontier for different VaR levels"""
        results = []
        for var_level in var_levels:
            try:
                result = self.optimize_portfolio(max_var=var_level)
                if result.success:
                    weights = result.x
                    expected_return = self.portfolio_expected_return(weights)
                    actual_var = self.portfolio_var(weights)
                    transaction_cost = self.nonlinear_transaction_cost(weights)
                    
                    results.append({
                        'var_limit': var_level,
                        'expected_return': expected_return,
                        'actual_var': actual_var,
                        'transaction_cost': transaction_cost,
                        'net_return': expected_return - transaction_cost,
                        'weights': weights
                    })
            except:
                continue
        
        return pd.DataFrame(results)

# Initialize optimizer
optimizer = NonlinearPortfolioOptimizer()

# Solve for a specific VaR constraint
print("=== Portfolio Optimization with Nonlinear Correlations ===")
print(f"Assets: {optimizer.asset_names}")
print(f"Expected Returns: {optimizer.mu}")
print(f"Volatilities: {optimizer.sigma}")
print(f"Transaction Cost Coefficients: {optimizer.transaction_costs}")

# Optimize with 2% VaR limit
result = optimizer.optimize_portfolio(max_var=0.02)

if result.success:
    optimal_weights = result.x
    print(f"\n=== Optimization Results ===")
    print(f"Optimization successful: {result.success}")
    print(f"Optimal weights:")
    for i, (asset, weight) in enumerate(zip(optimizer.asset_names, optimal_weights)):
        print(f"  {asset}: {weight:.4f} ({weight*100:.2f}%)")
    
    # Calculate portfolio metrics
    expected_return = optimizer.portfolio_expected_return(optimal_weights)
    portfolio_var = optimizer.portfolio_var(optimal_weights)
    transaction_cost = optimizer.nonlinear_transaction_cost(optimal_weights)
    net_return = expected_return - transaction_cost
    
    print(f"\n=== Portfolio Metrics ===")
    print(f"Expected Return: {expected_return:.4f} ({expected_return*100:.2f}%)")
    print(f"Transaction Cost: {transaction_cost:.4f} ({transaction_cost*100:.2f}%)")
    print(f"Net Expected Return: {net_return:.4f} ({net_return*100:.2f}%)")
    print(f"Portfolio VaR (5%): {portfolio_var:.4f} ({portfolio_var*100:.2f}%)")
else:
    print("Optimization failed!")
    print(result.message)

# Generate efficient frontier
print("\n=== Generating Efficient Frontier ===")
var_levels = np.linspace(0.01, 0.05, 20)
efficient_frontier_data = optimizer.efficient_frontier(var_levels)

print(f"Successfully computed {len(efficient_frontier_data)} frontier points")

# Create comprehensive visualizations
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Advanced Portfolio Optimization Analysis', fontsize=16, fontweight='bold')

# 1. Efficient Frontier
ax1 = axes[0, 0]
if len(efficient_frontier_data) > 0:
    ax1.scatter(efficient_frontier_data['actual_var'] * 100, 
               efficient_frontier_data['net_return'] * 100,
               c='blue', s=50, alpha=0.7, label='Efficient Frontier')
    ax1.scatter(portfolio_var * 100, net_return * 100, 
               c='red', s=100, marker='*', label='Optimal Portfolio')
ax1.set_xlabel('Portfolio VaR (%)')
ax1.set_ylabel('Net Expected Return (%)')
ax1.set_title('Efficient Frontier')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Optimal Portfolio Allocation
ax2 = axes[0, 1]
if result.success:
    colors = plt.cm.Set3(np.linspace(0, 1, len(optimizer.asset_names)))
    wedges, texts, autotexts = ax2.pie(optimal_weights, labels=optimizer.asset_names, 
                                      autopct='%1.1f%%', colors=colors, startangle=90)
    ax2.set_title('Optimal Portfolio Allocation')

# 3. Asset Returns vs Transaction Costs
ax3 = axes[0, 2]
x_pos = np.arange(len(optimizer.asset_names))
width = 0.35
bars1 = ax3.bar(x_pos - width/2, optimizer.mu * 100, width, 
                label='Expected Return (%)', alpha=0.8, color='green')
bars2 = ax3.bar(x_pos + width/2, optimizer.transaction_costs * 100, width,
                label='Transaction Cost Coeff (%)', alpha=0.8, color='red')
ax3.set_xlabel('Assets')
ax3.set_ylabel('Percentage (%)')
ax3.set_title('Expected Returns vs Transaction Costs')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(optimizer.asset_names, rotation=45)
ax3.legend()
ax3.grid(True, alpha=0.3)

# 4. Correlation Heatmap
ax4 = axes[1, 0]
correlation_matrix = optimizer.cov_matrix / np.outer(optimizer.sigma, optimizer.sigma)
im = ax4.imshow(correlation_matrix, cmap='RdBu_r', aspect='auto', vmin=-1, vmax=1)
ax4.set_xticks(range(len(optimizer.asset_names)))
ax4.set_yticks(range(len(optimizer.asset_names)))
ax4.set_xticklabels(optimizer.asset_names, rotation=45)
ax4.set_yticklabels(optimizer.asset_names)
ax4.set_title('Asset Correlation Matrix')

# Add correlation values to heatmap
for i in range(len(optimizer.asset_names)):
    for j in range(len(optimizer.asset_names)):
        text = ax4.text(j, i, f'{correlation_matrix[i, j]:.2f}',
                       ha="center", va="center", color="black", fontsize=8)

plt.colorbar(im, ax=ax4, shrink=0.8)

# 5. VaR vs Expected Return Trade-off
ax5 = axes[1, 1]
if len(efficient_frontier_data) > 0:
    ax5.plot(efficient_frontier_data['actual_var'] * 100, 
             efficient_frontier_data['expected_return'] * 100,
             'b-o', label='Gross Return', markersize=4)
    ax5.plot(efficient_frontier_data['actual_var'] * 100, 
             efficient_frontier_data['net_return'] * 100,
             'r--s', label='Net Return (after costs)', markersize=4)
    ax5.scatter(portfolio_var * 100, expected_return * 100, 
               c='blue', s=100, marker='*', label='Optimal (Gross)')
    ax5.scatter(portfolio_var * 100, net_return * 100, 
               c='red', s=100, marker='*', label='Optimal (Net)')
ax5.set_xlabel('Portfolio VaR (%)')
ax5.set_ylabel('Expected Return (%)')
ax5.set_title('Return vs Risk Trade-off')
ax5.legend(fontsize=8)
ax5.grid(True, alpha=0.3)

# 6. Portfolio Return Distribution
ax6 = axes[1, 2]
if result.success:
    portfolio_returns = np.dot(optimizer.return_scenarios, optimal_weights)
    ax6.hist(portfolio_returns * 100, bins=50, density=True, alpha=0.7, 
             color='skyblue', edgecolor='black')
    ax6.axvline(-portfolio_var * 100, color='red', linestyle='--', 
                label=f'VaR (5%): {portfolio_var*100:.2f}%')
    ax6.axvline(np.mean(portfolio_returns) * 100, color='green', linestyle='-', 
                label=f'Mean: {np.mean(portfolio_returns)*100:.2f}%')
    ax6.set_xlabel('Daily Return (%)')
    ax6.set_ylabel('Density')
    ax6.set_title('Portfolio Return Distribution')
    ax6.legend()
    ax6.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Summary statistics
print("\n=== Summary Analysis ===")
if len(efficient_frontier_data) > 0:
    print(f"Efficient Frontier Statistics:")
    print(f"  Number of efficient portfolios: {len(efficient_frontier_data)}")
    print(f"  VaR range: {efficient_frontier_data['actual_var'].min()*100:.2f}% - {efficient_frontier_data['actual_var'].max()*100:.2f}%")
    print(f"  Net return range: {efficient_frontier_data['net_return'].min()*100:.2f}% - {efficient_frontier_data['net_return'].max()*100:.2f}%")
    print(f"  Average transaction cost: {efficient_frontier_data['transaction_cost'].mean()*100:.3f}%")

if result.success:
    portfolio_returns = np.dot(optimizer.return_scenarios, optimal_weights)
    print(f"\nOptimal Portfolio Risk Metrics:")
    print(f"  VaR (1%): {-np.percentile(portfolio_returns, 1)*100:.3f}%")
    print(f"  VaR (5%): {-np.percentile(portfolio_returns, 5)*100:.3f}%")
    print(f"  VaR (10%): {-np.percentile(portfolio_returns, 10)*100:.3f}%")
    print(f"  Expected Shortfall (5%): {-np.mean(portfolio_returns[portfolio_returns <= np.percentile(portfolio_returns, 5)])*100:.3f}%")
    print(f"  Portfolio Volatility: {np.std(portfolio_returns)*np.sqrt(252)*100:.2f}% (annualized)")
    print(f"  Sharpe Ratio: {(net_return*252)/(np.std(portfolio_returns)*np.sqrt(252)):.3f}")

def setup_market_data(self):

def generate_nonlinear_returns(self):

def nonlinear_transaction_cost(self, weights):
    return np.sum(self.transaction_costs * np.power(weights, 1.5))

def portfolio_var(self, weights, confidence_level=0.05):

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.integrate import solve_ivp
import warnings
warnings.filterwarnings('ignore')

class TwoLinkRobot:
    """
    2-DOF planar robotic arm for trajectory optimization
    """
    def __init__(self, L1=1.0, L2=1.0, m1=1.0, m2=1.0):
        """
        Initialize robot parameters
        L1, L2: link lengths
        m1, m2: link masses
        """
        self.L1 = L1
        self.L2 = L2
        self.m1 = m1
        self.m2 = m2
        self.g = 9.81  # gravity
        
    def forward_kinematics(self, q):
        """
        Compute end-effector position given joint angles
        q: [q1, q2] joint angles
        Returns: [x, y] end-effector position
        """
        q1, q2 = q[0], q[1]
        x = self.L1 * np.cos(q1) + self.L2 * np.cos(q1 + q2)
        y = self.L1 * np.sin(q1) + self.L2 * np.sin(q1 + q2)
        return np.array([x, y])
    
    def mass_matrix(self, q):
        """
        Compute mass matrix M(q)
        """
        q1, q2 = q[0], q[1]
        
        # Mass matrix elements
        M11 = (self.m1 + self.m2) * self.L1**2 + self.m2 * self.L2**2 + 2 * self.m2 * self.L1 * self.L2 * np.cos(q2)
        M12 = self.m2 * self.L2**2 + self.m2 * self.L1 * self.L2 * np.cos(q2)
        M21 = M12
        M22 = self.m2 * self.L2**2
        
        return np.array([[M11, M12], [M21, M22]])
    
    def coriolis_matrix(self, q, dq):
        """
        Compute Coriolis matrix C(q, dq)
        """
        q1, q2 = q[0], q[1]
        dq1, dq2 = dq[0], dq[1]
        
        # Coriolis terms
        h = -self.m2 * self.L1 * self.L2 * np.sin(q2)
        C11 = h * dq2
        C12 = h * (dq1 + dq2)
        C21 = -h * dq1
        C22 = 0
        
        return np.array([[C11, C12], [C21, C22]])
    
    def gravity_vector(self, q):
        """
        Compute gravity vector G(q)
        """
        q1, q2 = q[0], q[1]
        
        G1 = (self.m1 + self.m2) * self.g * self.L1 * np.cos(q1) + self.m2 * self.g * self.L2 * np.cos(q1 + q2)
        G2 = self.m2 * self.g * self.L2 * np.cos(q1 + q2)
        
        return np.array([G1, G2])

class TrajectoryOptimizer:
    """
    Trajectory optimization for 2-DOF robot
    """
    def __init__(self, robot, T=2.0, N=50):
        """
        Initialize optimizer
        robot: TwoLinkRobot instance
        T: time horizon
        N: number of discretization points
        """
        self.robot = robot
        self.T = T
        self.N = N
        self.dt = T / (N - 1)
        
        # Obstacle parameters (circular obstacle)
        self.obstacle_center = np.array([1.0, 0.5])
        self.obstacle_radius = 0.3
        
    def interpolate_trajectory(self, q_waypoints):
        """
        Create smooth trajectory using polynomial interpolation
        """
        t = np.linspace(0, self.T, self.N)
        
        # Use cubic polynomial for smooth trajectory
        q_traj = np.zeros((self.N, 2))
        dq_traj = np.zeros((self.N, 2))
        ddq_traj = np.zeros((self.N, 2))
        
        for i in range(2):  # for each joint
            # Cubic polynomial coefficients
            # q(t) = a0 + a1*t + a2*t^2 + a3*t^3
            # Boundary conditions: q(0) = q_start, q(T) = q_end, dq(0) = 0, dq(T) = 0
            
            q_start = q_waypoints[0, i]
            q_end = q_waypoints[-1, i]
            
            # Solve for coefficients
            A = np.array([[0, 0, 0, 1],
                         [self.T**3, self.T**2, self.T, 1],
                         [0, 0, 1, 0],
                         [3*self.T**2, 2*self.T, 1, 0]])
            b = np.array([q_start, q_end, 0, 0])
            coeffs = np.linalg.solve(A, b)
            
            # Generate trajectory
            q_traj[:, i] = coeffs[3] + coeffs[2]*t + coeffs[1]*t**2 + coeffs[0]*t**3
            dq_traj[:, i] = coeffs[2] + 2*coeffs[1]*t + 3*coeffs[0]*t**2
            ddq_traj[:, i] = 2*coeffs[1] + 6*coeffs[0]*t
            
        return q_traj, dq_traj, ddq_traj
    
    def compute_torques(self, q_traj, dq_traj, ddq_traj):
        """
        Compute required torques using inverse dynamics
        """
        torques = np.zeros((self.N, 2))
        
        for i in range(self.N):
            q = q_traj[i]
            dq = dq_traj[i]
            ddq = ddq_traj[i]
            
            M = self.robot.mass_matrix(q)
            C = self.robot.coriolis_matrix(q, dq)
            G = self.robot.gravity_vector(q)
            
            # Inverse dynamics: tau = M*ddq + C*dq + G
            torques[i] = M @ ddq + C @ dq + G
            
        return torques
    
    def objective_function(self, q_waypoints_flat):
        """
        Objective function for optimization
        """
        # Reshape waypoints
        q_waypoints = q_waypoints_flat.reshape(-1, 2)
        
        # Generate trajectory
        q_traj, dq_traj, ddq_traj = self.interpolate_trajectory(q_waypoints)
        
        # Compute torques
        torques = self.compute_torques(q_traj, dq_traj, ddq_traj)
        
        # Energy cost (integral of torque squared)
        energy_cost = np.sum(torques**2) * self.dt
        
        # Smoothness cost (integral of acceleration squared)
        smoothness_cost = np.sum(ddq_traj**2) * self.dt
        
        # Total cost
        total_cost = energy_cost + 0.1 * smoothness_cost
        
        return total_cost
    
    def obstacle_constraint(self, q_waypoints_flat):
        """
        Obstacle avoidance constraint
        """
        q_waypoints = q_waypoints_flat.reshape(-1, 2)
        q_traj, _, _ = self.interpolate_trajectory(q_waypoints)
        
        constraints = []
        for i in range(self.N):
            # End-effector position
            ee_pos = self.robot.forward_kinematics(q_traj[i])
            
            # Distance to obstacle center
            dist = np.linalg.norm(ee_pos - self.obstacle_center)
            
            # Constraint: distance should be greater than obstacle radius
            constraints.append(dist - self.obstacle_radius - 0.1)  # 0.1m safety margin
            
        return np.array(constraints)
    
    def optimize_trajectory(self, q_start, q_end, n_waypoints=5):
        """
        Optimize trajectory from start to end configuration
        """
        # Initialize waypoints with linear interpolation
        t_waypoints = np.linspace(0, 1, n_waypoints)
        q_waypoints_init = np.zeros((n_waypoints, 2))
        
        for i in range(2):
            q_waypoints_init[:, i] = q_start[i] + t_waypoints * (q_end[i] - q_start[i])
        
        # Fix start and end points
        q_waypoints_init[0] = q_start
        q_waypoints_init[-1] = q_end
        
        # Flatten for optimization
        x0 = q_waypoints_init.flatten()
        
        # Bounds (joint limits)
        bounds = []
        for i in range(n_waypoints * 2):
            if i < 2 or i >= (n_waypoints - 1) * 2:  # Fix start and end points
                bounds.append((x0[i], x0[i]))
            else:
                bounds.append((-np.pi, np.pi))
        
        # Constraints
        constraints = {'type': 'ineq', 'fun': self.obstacle_constraint}
        
        # Optimize
        result = minimize(self.objective_function, x0, method='SLSQP', 
                         bounds=bounds, constraints=constraints,
                         options={'maxiter': 1000, 'ftol': 1e-6})
        
        # Extract optimized waypoints
        q_waypoints_opt = result.x.reshape(-1, 2)
        
        return q_waypoints_opt, result

# Initialize robot and optimizer
robot = TwoLinkRobot(L1=1.0, L2=1.0, m1=1.0, m2=1.0)
optimizer = TrajectoryOptimizer(robot, T=3.0, N=100)

# Define start and end configurations
q_start = np.array([0.0, 0.0])  # Straight configuration
q_end = np.array([np.pi/2, -np.pi/4])  # Target configuration

print("Starting trajectory optimization...")
print(f"Start configuration: q1={q_start[0]:.3f}, q2={q_start[1]:.3f}")
print(f"End configuration: q1={q_end[0]:.3f}, q2={q_end[1]:.3f}")

# Optimize trajectory
q_waypoints_opt, opt_result = optimizer.optimize_trajectory(q_start, q_end, n_waypoints=8)

print(f"\nOptimization completed!")
print(f"Success: {opt_result.success}")
print(f"Final cost: {opt_result.fun:.6f}")
print(f"Number of iterations: {opt_result.nit}")

# Generate final optimized trajectory
q_traj_opt, dq_traj_opt, ddq_traj_opt = optimizer.interpolate_trajectory(q_waypoints_opt)
torques_opt = optimizer.compute_torques(q_traj_opt, dq_traj_opt, ddq_traj_opt)

# Generate comparison trajectory (straight line in joint space)
q_waypoints_straight = np.array([q_start, q_end])
q_traj_straight, dq_traj_straight, ddq_traj_straight = optimizer.interpolate_trajectory(q_waypoints_straight)
torques_straight = optimizer.compute_torques(q_traj_straight, dq_traj_straight, ddq_traj_straight)

print("\nTrajectory generation completed!")
print("Ready for visualization...")

# Visualization
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Robot Trajectory Optimization Results', fontsize=16, fontweight='bold')

# Time vector
t = np.linspace(0, optimizer.T, optimizer.N)

# Plot 1: Joint trajectories
axes[0, 0].plot(t, q_traj_straight[:, 0], 'r--', label='Straight - Joint 1', linewidth=2)
axes[0, 0].plot(t, q_traj_straight[:, 1], 'b--', label='Straight - Joint 2', linewidth=2)
axes[0, 0].plot(t, q_traj_opt[:, 0], 'r-', label='Optimized - Joint 1', linewidth=2)
axes[0, 0].plot(t, q_traj_opt[:, 1], 'b-', label='Optimized - Joint 2', linewidth=2)
axes[0, 0].set_xlabel('Time [s]')
axes[0, 0].set_ylabel('Joint Angle [rad]')
axes[0, 0].set_title('Joint Trajectories')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Plot 2: Joint velocities
axes[0, 1].plot(t, dq_traj_straight[:, 0], 'r--', label='Straight - Joint 1', linewidth=2)
axes[0, 1].plot(t, dq_traj_straight[:, 1], 'b--', label='Straight - Joint 2', linewidth=2)
axes[0, 1].plot(t, dq_traj_opt[:, 0], 'r-', label='Optimized - Joint 1', linewidth=2)
axes[0, 1].plot(t, dq_traj_opt[:, 1], 'b-', label='Optimized - Joint 2', linewidth=2)
axes[0, 1].set_xlabel('Time [s]')
axes[0, 1].set_ylabel('Joint Velocity [rad/s]')
axes[0, 1].set_title('Joint Velocities')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Plot 3: Required torques
axes[0, 2].plot(t, torques_straight[:, 0], 'r--', label='Straight - Joint 1', linewidth=2)
axes[0, 2].plot(t, torques_straight[:, 1], 'b--', label='Straight - Joint 2', linewidth=2)
axes[0, 2].plot(t, torques_opt[:, 0], 'r-', label='Optimized - Joint 1', linewidth=2)
axes[0, 2].plot(t, torques_opt[:, 1], 'b-', label='Optimized - Joint 2', linewidth=2)
axes[0, 2].set_xlabel('Time [s]')
axes[0, 2].set_ylabel('Torque [Nm]')
axes[0, 2].set_title('Required Torques')
axes[0, 2].legend()
axes[0, 2].grid(True, alpha=0.3)

# Plot 4: End-effector trajectory in Cartesian space
ee_traj_straight = np.array([robot.forward_kinematics(q) for q in q_traj_straight])
ee_traj_opt = np.array([robot.forward_kinematics(q) for q in q_traj_opt])

axes[1, 0].plot(ee_traj_straight[:, 0], ee_traj_straight[:, 1], 'r--', label='Straight', linewidth=3)
axes[1, 0].plot(ee_traj_opt[:, 0], ee_traj_opt[:, 1], 'g-', label='Optimized', linewidth=3)

# Plot obstacle
circle = plt.Circle(optimizer.obstacle_center, optimizer.obstacle_radius, 
                   color='red', alpha=0.3, label='Obstacle')
axes[1, 0].add_patch(circle)

# Mark start and end points
axes[1, 0].plot(ee_traj_straight[0, 0], ee_traj_straight[0, 1], 'ko', markersize=8, label='Start')
axes[1, 0].plot(ee_traj_straight[-1, 0], ee_traj_straight[-1, 1], 'ks', markersize=8, label='End')

axes[1, 0].set_xlabel('X [m]')
axes[1, 0].set_ylabel('Y [m]')
axes[1, 0].set_title('End-Effector Trajectory')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_aspect('equal')

# Plot 5: Energy consumption comparison
energy_straight = np.cumsum(np.sum(torques_straight**2, axis=1)) * optimizer.dt
energy_opt = np.cumsum(np.sum(torques_opt**2, axis=1)) * optimizer.dt

axes[1, 1].plot(t, energy_straight, 'r--', label='Straight', linewidth=3)
axes[1, 1].plot(t, energy_opt, 'g-', label='Optimized', linewidth=3)
axes[1, 1].set_xlabel('Time [s]')
axes[1, 1].set_ylabel('Cumulative Energy [J]')
axes[1, 1].set_title('Energy Consumption')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# Plot 6: Robot arm animation (final positions)
def plot_robot_arm(ax, q, color, label, alpha=1.0):
    # Joint positions
    x0, y0 = 0, 0  # Base
    x1 = robot.L1 * np.cos(q[0])
    y1 = robot.L1 * np.sin(q[0])
    x2 = x1 + robot.L2 * np.cos(q[0] + q[1])
    y2 = y1 + robot.L2 * np.sin(q[0] + q[1])
    
    # Plot links
    ax.plot([x0, x1], [y0, y1], color=color, linewidth=4, alpha=alpha, label=f'{label} Link 1')
    ax.plot([x1, x2], [y1, y2], color=color, linewidth=4, alpha=alpha, label=f'{label} Link 2')
    
    # Plot joints
    ax.plot([x0, x1, x2], [y0, y1, y2], 'o', color=color, markersize=8, alpha=alpha)
    
    return x2, y2

# Plot initial and final configurations
plot_robot_arm(axes[1, 2], q_start, 'blue', 'Initial', alpha=0.5)
plot_robot_arm(axes[1, 2], q_end, 'red', 'Final', alpha=0.8)

# Plot obstacle
circle2 = plt.Circle(optimizer.obstacle_center, optimizer.obstacle_radius, 
                    color='red', alpha=0.3)
axes[1, 2].add_patch(circle2)

axes[1, 2].set_xlabel('X [m]')
axes[1, 2].set_ylabel('Y [m]')
axes[1, 2].set_title('Robot Arm Configurations')
axes[1, 2].grid(True, alpha=0.3)
axes[1, 2].set_aspect('equal')
axes[1, 2].set_xlim(-0.5, 2.5)
axes[1, 2].set_ylim(-0.5, 2.5)

plt.tight_layout()
plt.show()

# Print performance comparison
total_energy_straight = energy_straight[-1]
total_energy_opt = energy_opt[-1]
energy_saving = (total_energy_straight - total_energy_opt) / total_energy_straight * 100

print(f"\n=== Performance Comparison ===")
print(f"Total Energy - Straight trajectory: {total_energy_straight:.4f} J")
print(f"Total Energy - Optimized trajectory: {total_energy_opt:.4f} J")
print(f"Energy saving: {energy_saving:.2f}%")

# Check obstacle avoidance
min_distance_straight = min([np.linalg.norm(robot.forward_kinematics(q) - optimizer.obstacle_center) 
                           for q in q_traj_straight])
min_distance_opt = min([np.linalg.norm(robot.forward_kinematics(q) - optimizer.obstacle_center) 
                       for q in q_traj_opt])

print(f"\nMinimum distance to obstacle:")
print(f"Straight trajectory: {min_distance_straight:.4f} m")
print(f"Optimized trajectory: {min_distance_opt:.4f} m")
print(f"Obstacle radius: {optimizer.obstacle_radius:.4f} m")

if min_distance_straight < optimizer.obstacle_radius:
    print("⚠️  Straight trajectory collides with obstacle!")
else:
    print("✅ Straight trajectory avoids obstacle")

if min_distance_opt < optimizer.obstacle_radius:
    print("⚠️  Optimized trajectory collides with obstacle!")
else:
    print("✅ Optimized trajectory avoids obstacle")