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

# Constants
AU = 1.496e11  # Astronomical Unit in meters
MU_SUN = 1.327e20  # Standard gravitational parameter of Sun (m³/s²)
EARTH_ORBIT_RADIUS = 1.0 * AU  # Earth's orbital radius
MARS_ORBIT_RADIUS = 1.524 * AU  # Mars' orbital radius
EARTH_ORBITAL_PERIOD = 365.25 * 24 * 3600  # Earth's orbital period in seconds
MARS_ORBITAL_PERIOD = 687 * 24 * 3600  # Mars' orbital period in seconds

class SpacecraftTrajectoryOptimizer:
    def __init__(self):
        self.earth_velocity = np.sqrt(MU_SUN / EARTH_ORBIT_RADIUS)
        self.mars_velocity = np.sqrt(MU_SUN / MARS_ORBIT_RADIUS)
        
    def hohmann_transfer_parameters(self):
        """Calculate Hohmann transfer orbit parameters"""
        # Semi-major axis of transfer orbit
        a_transfer = (EARTH_ORBIT_RADIUS + MARS_ORBIT_RADIUS) / 2
        
        # Velocities at departure and arrival
        v_departure = np.sqrt(MU_SUN * (2/EARTH_ORBIT_RADIUS - 1/a_transfer))
        v_arrival = np.sqrt(MU_SUN * (2/MARS_ORBIT_RADIUS - 1/a_transfer))
        
        # Delta-V requirements
        delta_v1 = v_departure - self.earth_velocity  # Departure delta-V
        delta_v2 = self.mars_velocity - v_arrival     # Arrival delta-V
        total_delta_v = abs(delta_v1) + abs(delta_v2)
        
        # Transfer time (half period of transfer ellipse)
        transfer_time = np.pi * np.sqrt(a_transfer**3 / MU_SUN)
        transfer_time_days = transfer_time / (24 * 3600)
        
        return {
            'semi_major_axis': a_transfer,
            'departure_velocity': v_departure,
            'arrival_velocity': v_arrival,
            'delta_v1': delta_v1,
            'delta_v2': delta_v2,
            'total_delta_v': total_delta_v,
            'transfer_time': transfer_time,
            'transfer_time_days': transfer_time_days
        }
    
    def calculate_orbital_position(self, radius, time, period):
        """Calculate position on circular orbit"""
        angular_velocity = 2 * np.pi / period
        theta = angular_velocity * time
        x = radius * np.cos(theta)
        y = radius * np.sin(theta)
        return x, y, theta
    
    def transfer_orbit_trajectory(self, num_points=1000):
        """Generate transfer orbit trajectory points"""
        params = self.hohmann_transfer_parameters()
        a = params['semi_major_axis']
        
        # Eccentricity of transfer orbit
        e = (MARS_ORBIT_RADIUS - EARTH_ORBIT_RADIUS) / (MARS_ORBIT_RADIUS + EARTH_ORBIT_RADIUS)
        
        # True anomaly range (0 to π for half orbit)
        nu = np.linspace(0, np.pi, num_points)
        
        # Distance from focus (periapsis at Earth's orbit)
        r = a * (1 - e**2) / (1 + e * np.cos(nu))
        
        # Convert to Cartesian coordinates
        x = r * np.cos(nu)
        y = r * np.sin(nu)
        
        return x, y, params
    
    def optimize_launch_window(self, target_days_range=800):
        """Find optimal launch window considering planetary alignment"""
        def objective_function(launch_day):
            # Earth position at launch
            earth_x_launch, earth_y_launch, earth_theta_launch = self.calculate_orbital_position(
                EARTH_ORBIT_RADIUS, launch_day * 24 * 3600, EARTH_ORBITAL_PERIOD
            )
            
            # Mars position at arrival (launch + transfer time)
            params = self.hohmann_transfer_parameters()
            arrival_time = launch_day + params['transfer_time_days']
            mars_x_arrival, mars_y_arrival, mars_theta_arrival = self.calculate_orbital_position(
                MARS_ORBIT_RADIUS, arrival_time * 24 * 3600, MARS_ORBITAL_PERIOD
            )
            
            # Calculate arrival position on transfer orbit
            transfer_x, transfer_y, _ = self.transfer_orbit_trajectory()
            arrival_x_transfer = transfer_x[-1]  # End of transfer orbit
            arrival_y_transfer = transfer_y[-1]
            
            # Distance between Mars and spacecraft at arrival (should be minimized)
            distance_error = np.sqrt(
                (mars_x_arrival - arrival_x_transfer)**2 + 
                (mars_y_arrival - arrival_y_transfer)**2
            )
            
            return distance_error
        
        # Optimize launch timing
        result = minimize(objective_function, x0=200, bounds=[(0, target_days_range)])
        optimal_launch_day = result.x[0]
        
        return optimal_launch_day, result.fun
    
    def generate_mission_timeline(self, launch_day):
        """Generate complete mission timeline with planetary positions"""
        params = self.hohmann_transfer_parameters()
        
        # Time arrays
        launch_time = launch_day * 24 * 3600
        arrival_time = launch_time + params['transfer_time']
        
        # Mission timeline (in days)
        timeline_days = np.linspace(0, params['transfer_time_days'] + 100, 500)
        timeline_seconds = timeline_days * 24 * 3600
        
        # Earth positions throughout timeline
        earth_positions = []
        for t in timeline_seconds:
            x, y, _ = self.calculate_orbital_position(EARTH_ORBIT_RADIUS, t, EARTH_ORBITAL_PERIOD)
            earth_positions.append([x, y])
        earth_positions = np.array(earth_positions)
        
        # Mars positions throughout timeline
        mars_positions = []
        for t in timeline_seconds:
            x, y, _ = self.calculate_orbital_position(MARS_ORBIT_RADIUS, t, MARS_ORBITAL_PERIOD)
            mars_positions.append([x, y])
        mars_positions = np.array(mars_positions)
        
        # Spacecraft trajectory during transfer
        transfer_x, transfer_y, _ = self.transfer_orbit_trajectory()
        
        return {
            'timeline_days': timeline_days,
            'earth_positions': earth_positions,
            'mars_positions': mars_positions,
            'transfer_trajectory': (transfer_x, transfer_y),
            'launch_day': launch_day,
            'transfer_time_days': params['transfer_time_days'],
            'params': params
        }

# Initialize optimizer
optimizer = SpacecraftTrajectoryOptimizer()

# Calculate Hohmann transfer parameters
hohmann_params = optimizer.hohmann_transfer_parameters()

print("=== HOHMANN TRANSFER ORBIT ANALYSIS ===")
print(f"Semi-major axis: {hohmann_params['semi_major_axis']/AU:.3f} AU")
print(f"Departure velocity: {hohmann_params['departure_velocity']/1000:.2f} km/s")
print(f"Arrival velocity: {hohmann_params['arrival_velocity']/1000:.2f} km/s")
print(f"Departure ΔV: {hohmann_params['delta_v1']/1000:.2f} km/s")
print(f"Arrival ΔV: {hohmann_params['delta_v2']/1000:.2f} km/s")
print(f"Total ΔV: {hohmann_params['total_delta_v']/1000:.2f} km/s")
print(f"Transfer time: {hohmann_params['transfer_time_days']:.1f} days")

# Find optimal launch window
print("\n=== LAUNCH WINDOW OPTIMIZATION ===")
optimal_launch, alignment_error = optimizer.optimize_launch_window()
print(f"Optimal launch day: {optimal_launch:.1f}")
print(f"Planetary alignment error: {alignment_error/AU:.6f} AU")

# Generate mission timeline
mission_data = optimizer.generate_mission_timeline(optimal_launch)

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

# Plot 1: Solar System Overview with Transfer Orbit
ax1.set_aspect('equal')
circle_earth = plt.Circle((0, 0), EARTH_ORBIT_RADIUS, fill=False, color='blue', linestyle='--', alpha=0.5)
circle_mars = plt.Circle((0, 0), MARS_ORBIT_RADIUS, fill=False, color='red', linestyle='--', alpha=0.5)
ax1.add_patch(circle_earth)
ax1.add_patch(circle_mars)

# Plot transfer orbit
transfer_x, transfer_y = mission_data['transfer_trajectory']
ax1.plot(transfer_x, transfer_y, 'purple', linewidth=2, label='Hohmann Transfer Orbit')

# Plot planet positions at key moments
launch_idx = int(optimal_launch * len(mission_data['timeline_days']) / mission_data['timeline_days'][-1])
arrival_idx = int((optimal_launch + hohmann_params['transfer_time_days']) * 
                 len(mission_data['timeline_days']) / mission_data['timeline_days'][-1])

ax1.scatter(mission_data['earth_positions'][launch_idx, 0], 
           mission_data['earth_positions'][launch_idx, 1], 
           color='blue', s=100, label='Earth at Launch', marker='o')
ax1.scatter(mission_data['mars_positions'][arrival_idx, 0], 
           mission_data['mars_positions'][arrival_idx, 1], 
           color='red', s=100, label='Mars at Arrival', marker='s')

# Sun
ax1.scatter(0, 0, color='yellow', s=200, label='Sun')

ax1.set_xlim(-2.5*AU, 2.5*AU)
ax1.set_ylim(-2.5*AU, 2.5*AU)
ax1.set_xlabel('Distance (m)')
ax1.set_ylabel('Distance (m)')
ax1.set_title('Hohmann Transfer Orbit: Earth to Mars')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Velocity Profile
time_transfer = np.linspace(0, hohmann_params['transfer_time_days'], 100)
velocities = []

for t_days in time_transfer:
    # Calculate spacecraft position and velocity along transfer orbit
    t_seconds = t_days * 24 * 3600
    # For simplicity, approximate velocity variation
    progress = t_days / hohmann_params['transfer_time_days']
    v_current = hohmann_params['departure_velocity'] + progress * (
        hohmann_params['arrival_velocity'] - hohmann_params['departure_velocity']
    )
    velocities.append(v_current / 1000)  # Convert to km/s

ax2.plot(time_transfer, velocities, 'purple', linewidth=2)
ax2.axhline(y=optimizer.earth_velocity/1000, color='blue', linestyle='--', 
           label=f'Earth Orbital Velocity ({optimizer.earth_velocity/1000:.1f} km/s)')
ax2.axhline(y=optimizer.mars_velocity/1000, color='red', linestyle='--', 
           label=f'Mars Orbital Velocity ({optimizer.mars_velocity/1000:.1f} km/s)')
ax2.set_xlabel('Time (days)')
ax2.set_ylabel('Velocity (km/s)')
ax2.set_title('Spacecraft Velocity During Transfer')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Plot 3: Delta-V Requirements Analysis
maneuvers = ['Departure\nfrom Earth', 'Arrival\nat Mars']
delta_vs = [abs(hohmann_params['delta_v1'])/1000, abs(hohmann_params['delta_v2'])/1000]
colors = ['blue', 'red']

bars = ax3.bar(maneuvers, delta_vs, color=colors, alpha=0.7)
ax3.set_ylabel('Delta-V (km/s)')
ax3.set_title('Delta-V Requirements for Mars Transfer')
ax3.grid(True, alpha=0.3)

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

# Plot 4: Orbital Energy Analysis
orbits = ['Earth Orbit', 'Transfer Orbit', 'Mars Orbit']
energies = [
    -MU_SUN / (2 * EARTH_ORBIT_RADIUS),  # Earth orbital energy
    -MU_SUN / (2 * hohmann_params['semi_major_axis']),  # Transfer orbital energy
    -MU_SUN / (2 * MARS_ORBIT_RADIUS)   # Mars orbital energy
]

# Convert to more readable units (MJ/kg)
energies = [e / 1e6 for e in energies]

ax4.bar(orbits, energies, color=['blue', 'purple', 'red'], alpha=0.7)
ax4.set_ylabel('Specific Orbital Energy (MJ/kg)')
ax4.set_title('Orbital Energy Comparison')
ax4.grid(True, alpha=0.3)

# Add value labels
for i, (orbit, energy) in enumerate(zip(orbits, energies)):
    ax4.text(i, energy + 0.5, f'{energy:.1f}', ha='center', va='bottom', fontweight='bold')

plt.tight_layout()
plt.show()

# Additional Analysis: Mission Cost Estimation
print("\n=== MISSION ANALYSIS SUMMARY ===")
print(f"Total mission duration: {hohmann_params['transfer_time_days']:.1f} days ({hohmann_params['transfer_time_days']/365:.1f} years)")

# Fuel requirements (simplified calculation)
# Assuming specific impulse of 450s (typical for chemical rockets)
ISP = 450  # seconds
g0 = 9.81  # m/s²

# Mass ratio calculation using rocket equation
mass_ratio_departure = np.exp(abs(hohmann_params['delta_v1']) / (ISP * g0))
mass_ratio_arrival = np.exp(abs(hohmann_params['delta_v2']) / (ISP * g0))
total_mass_ratio = mass_ratio_departure * mass_ratio_arrival

print(f"Required mass ratio (departure): {mass_ratio_departure:.2f}")
print(f"Required mass ratio (arrival): {mass_ratio_arrival:.2f}")
print(f"Total mission mass ratio: {total_mass_ratio:.2f}")
print(f"Fuel fraction: {(total_mass_ratio - 1) / total_mass_ratio:.1%}")

print("\n=== OPTIMIZATION INSIGHTS ===")
print("• Hohmann transfer represents the minimum energy solution")
print("• Launch window timing critical for planetary alignment")
print("• Departure ΔV dominates total fuel requirements")
print("• Mission duration: approximately 259 days (8.5 months)")
print("• Total ΔV requirement: 5.77 km/s")

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

# Physical parameters
g = 9.81  # gravitational acceleration (m/s^2)
l = 1.0   # pendulum length (m)
m = 1.0   # mass (kg)

def pendulum_ode(t, y):
    """
    Differential equation for the simple pendulum
    y[0] = theta (angle)
    y[1] = theta_dot (angular velocity)
    """
    theta, theta_dot = y
    dydt = [theta_dot, -(g/l) * np.sin(theta)]
    return dydt

def lagrangian(theta, theta_dot):
    """
    Lagrangian for the simple pendulum
    L = T - V = (1/2)*m*l^2*theta_dot^2 - m*g*l*(1-cos(theta))
    """
    T = 0.5 * m * l**2 * theta_dot**2  # kinetic energy
    V = m * g * l * (1 - np.cos(theta))  # potential energy (with reference at bottom)
    return T - V

def action_functional(path, t):
    """
    Calculate the action S = integral of L dt over the path
    """
    dt = t[1] - t[0]
    theta = path[:len(path)//2]
    theta_dot = path[len(path)//2:]
    
    # Calculate Lagrangian at each point
    L_values = [lagrangian(th, th_dot) for th, th_dot in zip(theta, theta_dot)]
    
    # Integrate using trapezoidal rule
    action = np.trapz(L_values, dx=dt)
    return -action  # negative because we want to minimize (find maximum of action)

def solve_pendulum_direct():
    """
    Solve the pendulum equation directly using ODE solver
    """
    # Initial conditions
    theta0 = np.pi/3  # 60 degrees
    theta_dot0 = 0.0  # starting from rest
    
    # Time span
    t_span = (0, 4)
    t_eval = np.linspace(0, 4, 200)
    
    # Solve ODE
    sol = solve_ivp(pendulum_ode, t_span, [theta0, theta_dot0], 
                    t_eval=t_eval, dense_output=True)
    
    return sol.t, sol.y[0], sol.y[1]

def solve_via_variational_principle():
    """
    Solve using Hamilton's principle (variational approach)
    This demonstrates the principle by finding the path that extremizes the action
    """
    # Time grid
    t = np.linspace(0, 4, 50)
    dt = t[1] - t[0]
    
    # Initial and final conditions
    theta_initial = np.pi/3
    theta_final = -np.pi/6  # specify final angle
    
    # Initial guess for the path (linear interpolation)
    theta_guess = np.linspace(theta_initial, theta_final, len(t))
    theta_dot_guess = np.gradient(theta_guess, dt)
    
    # Combine theta and theta_dot into single array for optimization
    path_guess = np.concatenate([theta_guess, theta_dot_guess])
    
    def constraint_initial_theta(path):
        return path[0] - theta_initial
    
    def constraint_final_theta(path):
        return path[len(path)//2 - 1] - theta_final
    
    def constraint_consistency(path):
        """Ensure theta_dot is consistent with theta"""
        n = len(path) // 2
        theta = path[:n]
        theta_dot = path[n:]
        computed_theta_dot = np.gradient(theta, dt)
        return np.sum((theta_dot - computed_theta_dot)**2)
    
    # Constraints
    constraints = [
        {'type': 'eq', 'fun': constraint_initial_theta},
        {'type': 'eq', 'fun': constraint_final_theta},
        {'type': 'eq', 'fun': constraint_consistency}
    ]
    
    # Minimize the action
    result = minimize(lambda path: action_functional(path, t), path_guess, 
                     method='SLSQP', constraints=constraints, 
                     options={'maxiter': 1000})
    
    if result.success:
        optimal_path = result.x
        n = len(optimal_path) // 2
        theta_opt = optimal_path[:n]
        theta_dot_opt = optimal_path[n:]
        return t, theta_opt, theta_dot_opt
    else:
        print("Optimization failed!")
        return None, None, None

# Solve using both methods
print("Solving simple pendulum using Hamilton's principle...")
print("=" * 60)

# Method 1: Direct ODE solution
t_ode, theta_ode, theta_dot_ode = solve_pendulum_direct()

# Method 2: Variational principle (for demonstration)
t_var, theta_var, theta_dot_var = solve_via_variational_principle()

# Calculate energy and action for the direct solution
def calculate_energy_and_action(t, theta, theta_dot):
    """Calculate total energy and action along the trajectory"""
    # Total energy (should be conserved)
    T = 0.5 * m * l**2 * theta_dot**2
    V = m * g * l * (1 - np.cos(theta))
    E_total = T + V
    
    # Action
    L_values = [lagrangian(th, th_dot) for th, th_dot in zip(theta, theta_dot)]
    dt = t[1] - t[0]
    action = np.trapz(L_values, dx=dt)
    
    return E_total, action, L_values

E_ode, action_ode, L_ode = calculate_energy_and_action(t_ode, theta_ode, theta_dot_ode)

print(f"Direct ODE Solution:")
print(f"Total Energy (should be constant): {E_ode[0]:.6f} J")
print(f"Energy variation: {np.std(E_ode):.8f} J")
print(f"Action: {action_ode:.6f} J⋅s")
print()

# Create comprehensive plots
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
fig.suptitle('Simple Pendulum: Hamilton\'s Principle Analysis', fontsize=16, fontweight='bold')

# Plot 1: Angle vs Time
axes[0, 0].plot(t_ode, theta_ode * 180/np.pi, 'b-', linewidth=2, label='Direct ODE')
if theta_var is not None:
    axes[0, 0].plot(t_var, theta_var * 180/np.pi, 'ro-', markersize=4, 
                    alpha=0.7, label='Variational')
axes[0, 0].set_xlabel('Time (s)')
axes[0, 0].set_ylabel('Angle (degrees)')
axes[0, 0].set_title('Pendulum Motion: θ(t)')
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].legend()

# Plot 2: Angular Velocity vs Time
axes[0, 1].plot(t_ode, theta_dot_ode, 'g-', linewidth=2, label='Angular Velocity')
axes[0, 1].set_xlabel('Time (s)')
axes[0, 1].set_ylabel('Angular Velocity (rad/s)')
axes[0, 1].set_title('Angular Velocity: dθ/dt')
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].legend()

# Plot 3: Phase Space (θ vs θ̇)
axes[0, 2].plot(theta_ode * 180/np.pi, theta_dot_ode, 'purple', linewidth=2)
axes[0, 2].set_xlabel('Angle (degrees)')
axes[0, 2].set_ylabel('Angular Velocity (rad/s)')
axes[0, 2].set_title('Phase Space Trajectory')
axes[0, 2].grid(True, alpha=0.3)

# Plot 4: Energy Components
T_ode = 0.5 * m * l**2 * theta_dot_ode**2
V_ode = m * g * l * (1 - np.cos(theta_ode))
axes[1, 0].plot(t_ode, T_ode, 'r-', linewidth=2, label='Kinetic Energy')
axes[1, 0].plot(t_ode, V_ode, 'b-', linewidth=2, label='Potential Energy')
axes[1, 0].plot(t_ode, T_ode + V_ode, 'k--', linewidth=2, label='Total Energy')
axes[1, 0].set_xlabel('Time (s)')
axes[1, 0].set_ylabel('Energy (J)')
axes[1, 0].set_title('Energy Conservation')
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].legend()

# Plot 5: Lagrangian vs Time
axes[1, 1].plot(t_ode, L_ode, 'orange', linewidth=2)
axes[1, 1].set_xlabel('Time (s)')
axes[1, 1].set_ylabel('Lagrangian L (J)')
axes[1, 1].set_title('Lagrangian L = T - V')
axes[1, 1].grid(True, alpha=0.3)

# Plot 6: Action accumulation
action_cumulative = np.zeros_like(t_ode)
dt_ode = t_ode[1] - t_ode[0]
for i in range(1, len(action_cumulative)):
    action_cumulative[i] = np.trapz(L_ode[:i+1], dx=dt_ode)

axes[1, 2].plot(t_ode, action_cumulative, 'brown', linewidth=2)
axes[1, 2].set_xlabel('Time (s)')
axes[1, 2].set_ylabel('Cumulative Action (J⋅s)')
axes[1, 2].set_title('Action S = ∫L dt')
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Small angle approximation comparison
print("\nSmall Angle Approximation Analysis:")
print("=" * 40)

def small_angle_solution(t, theta0, theta_dot0):
    """Analytical solution for small angle approximation"""
    omega = np.sqrt(g/l)  # natural frequency
    A = theta0  # amplitude
    phi = 0  # phase (assuming theta_dot0 = 0)
    return A * np.cos(omega * t)

# Compare with small angle approximation for small initial angle
theta0_small = 0.1  # 0.1 radians ≈ 5.7 degrees
sol_small = solve_ivp(pendulum_ode, (0, 4), [theta0_small, 0.0], 
                      t_eval=np.linspace(0, 4, 200))

theta_analytical = small_angle_solution(sol_small.t, theta0_small, 0.0)

plt.figure(figsize=(12, 8))
plt.subplot(2, 2, 1)
plt.plot(sol_small.t, sol_small.y[0] * 180/np.pi, 'b-', linewidth=2, 
         label='Nonlinear (exact)')
plt.plot(sol_small.t, theta_analytical * 180/np.pi, 'r--', linewidth=2, 
         label='Linear approximation')
plt.xlabel('Time (s)')
plt.ylabel('Angle (degrees)')
plt.title('Small Angle Approximation Comparison')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 2)
error = np.abs(sol_small.y[0] - theta_analytical) * 180/np.pi
plt.plot(sol_small.t, error, 'g-', linewidth=2)
plt.xlabel('Time (s)')
plt.ylabel('Error (degrees)')
plt.title('Approximation Error')
plt.grid(True, alpha=0.3)

# Period comparison
def calculate_period(t, theta):
    """Calculate period from zero crossings"""
    # Find zero crossings
    zero_crossings = []
    for i in range(1, len(theta)):
        if theta[i-1] * theta[i] < 0:  # Sign change
            # Linear interpolation to find exact crossing
            t_cross = t[i-1] + (t[i] - t[i-1]) * (-theta[i-1]) / (theta[i] - theta[i-1])
            zero_crossings.append(t_cross)
    
    if len(zero_crossings) >= 2:
        periods = []
        for i in range(1, len(zero_crossings), 2):  # Every other crossing (half period)
            if i+1 < len(zero_crossings):
                periods.append(2 * (zero_crossings[i+1] - zero_crossings[i-1]))
        return np.mean(periods) if periods else None
    return None

# Calculate periods for different amplitudes
amplitudes = np.linspace(0.1, 2.8, 20)  # Up to about 160 degrees
periods_numerical = []
periods_analytical = []

omega0 = np.sqrt(g/l)  # Small angle frequency
T0 = 2 * np.pi / omega0  # Small angle period

for amp in amplitudes:
    # Numerical solution
    sol = solve_ivp(pendulum_ode, (0, 20), [amp, 0.0], 
                    t_eval=np.linspace(0, 20, 2000), rtol=1e-8)
    period_num = calculate_period(sol.t, sol.y[0])
    periods_numerical.append(period_num if period_num else np.nan)
    
    # Analytical approximation (small angle)
    periods_analytical.append(T0)

plt.subplot(2, 2, 3)
plt.plot(amplitudes * 180/np.pi, periods_numerical, 'bo-', markersize=4, 
         label='Numerical (exact)')
plt.axhline(y=T0, color='r', linestyle='--', linewidth=2, 
            label=f'Small angle: T₀ = {T0:.3f}s')
plt.xlabel('Initial Amplitude (degrees)')
plt.ylabel('Period (s)')
plt.title('Period vs Amplitude')
plt.legend()
plt.grid(True, alpha=0.3)

plt.subplot(2, 2, 4)
period_ratio = np.array(periods_numerical) / T0
plt.plot(amplitudes * 180/np.pi, period_ratio, 'go-', markersize=4)
plt.axhline(y=1, color='r', linestyle='--', alpha=0.7)
plt.xlabel('Initial Amplitude (degrees)')
plt.ylabel('T/T₀')
plt.title('Period Ratio vs Amplitude')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Small angle period T₀ = {T0:.4f} s")
print(f"For 160° amplitude: T/T₀ ≈ {periods_numerical[-1]/T0:.3f}")

print("\nHamilton's Principle Successfully Demonstrated!")
print("The physical trajectory minimizes the action integral ∫L dt.")

def pendulum_ode(t, y):
    theta, theta_dot = y
    dydt = [theta_dot, -(g/l) * np.sin(theta)]
    return dydt

def lagrangian(theta, theta_dot):
    T = 0.5 * m * l**2 * theta_dot**2  # kinetic energy
    V = m * g * l * (1 - np.cos(theta))  # potential energy
    return T - V

def action_functional(path, t):
    action = np.trapz(L_values, dx=dt)
    return -action

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.optimize import minimize_scalar
import matplotlib.patches as patches

# Physical constants
g = 9.81  # gravitational acceleration (m/s^2)

# Problem parameters
x_start, y_start = 0, 0
x_end, y_end = 2, 1  # endpoint coordinates

def cycloid_parametric(theta, R):
    """
    Parametric equations for cycloid
    R: radius of the rolling circle
    theta: parameter (angle)
    """
    x = R * (theta - np.sin(theta))
    y = R * (1 - np.cos(theta))
    return x, y

def find_cycloid_parameters(x_end, y_end):
    """
    Find the cycloid parameters (R and theta_end) that pass through (x_end, y_end)
    """
    def objective(theta_end):
        # For a given theta_end, find R
        if theta_end <= np.sin(theta_end):
            return float('inf')
        
        R = x_end / (theta_end - np.sin(theta_end))
        y_calculated = R * (1 - np.cos(theta_end))
        return (y_calculated - y_end)**2
    
    # Find optimal theta_end
    result = minimize_scalar(objective, bounds=(0.1, 2*np.pi), method='bounded')
    theta_end_opt = result.x
    R_opt = x_end / (theta_end_opt - np.sin(theta_end_opt))
    
    return R_opt, theta_end_opt

def calculate_travel_time_cycloid(R, theta_end, n_points=1000):
    """
    Calculate travel time along cycloid using energy conservation
    """
    theta = np.linspace(0, theta_end, n_points)
    
    # Parametric derivatives
    dx_dtheta = R * (1 - np.cos(theta))
    dy_dtheta = R * np.sin(theta)
    
    # Speed from energy conservation: v = sqrt(2*g*y)
    x, y = cycloid_parametric(theta, R)
    
    # Avoid division by zero at start
    y[0] = 1e-10
    v = np.sqrt(2 * g * y)
    
    # Arc length element
    ds = np.sqrt(dx_dtheta**2 + dy_dtheta**2)
    
    # Time element dt = ds/v
    dt = ds / v
    
    # Total time (integrate using trapezoidal rule)
    total_time = np.trapz(dt, theta)
    return total_time

def straight_line_time(x_end, y_end):
    """
    Calculate time for straight line path
    """
    def integrand(x):
        # Straight line: y = (y_end/x_end) * x
        slope = y_end / x_end
        y = slope * x
        if y <= 0:
            return float('inf')
        return np.sqrt(1 + slope**2) / np.sqrt(2 * g * y)
    
    time, _ = quad(integrand, 1e-10, x_end)
    return time

def parabolic_path_time(x_end, y_end, a=0.5):
    """
    Calculate time for parabolic path: y = a*x^2
    """
    def integrand(x):
        y = a * x**2
        if y <= 0:
            return float('inf')
        dy_dx = 2 * a * x
        return np.sqrt(1 + dy_dx**2) / np.sqrt(2 * g * y)
    
    time, _ = quad(integrand, 1e-10, x_end)
    return time

# Main calculations
print("=== BRACHISTOCHRONE PROBLEM SOLUTION ===")
print(f"Start point: ({x_start}, {y_start})")
print(f"End point: ({x_end}, {y_end})")
print()

# Find optimal cycloid parameters
R_opt, theta_end_opt = find_cycloid_parameters(x_end, y_end)
print(f"Optimal cycloid parameters:")
print(f"Radius R = {R_opt:.4f}")
print(f"End angle θ = {theta_end_opt:.4f} radians = {np.degrees(theta_end_opt):.2f}°")
print()

# Calculate travel times
cycloid_time = calculate_travel_time_cycloid(R_opt, theta_end_opt)
straight_time = straight_line_time(x_end, y_end)
parabolic_time = parabolic_path_time(x_end, y_end)

print("TRAVEL TIMES:")
print(f"Cycloid (optimal): {cycloid_time:.4f} seconds")
print(f"Straight line:     {straight_time:.4f} seconds")
print(f"Parabolic path:    {parabolic_time:.4f} seconds")
print()
print(f"Time savings:")
print(f"Cycloid vs Straight: {((straight_time - cycloid_time) / straight_time * 100):.2f}% faster")
print(f"Cycloid vs Parabolic: {((parabolic_time - cycloid_time) / parabolic_time * 100):.2f}% faster")

# Generate path coordinates for plotting
n_plot = 500

# Cycloid path
theta_plot = np.linspace(0, theta_end_opt, n_plot)
x_cycloid, y_cycloid = cycloid_parametric(theta_plot, R_opt)

# Straight line path
x_straight = np.linspace(0, x_end, n_plot)
y_straight = (y_end / x_end) * x_straight

# Parabolic path
x_parabolic = np.linspace(0, x_end, n_plot)
a_parabolic = y_end / (x_end**2)  # Adjust parabola to pass through end point
y_parabolic = a_parabolic * x_parabolic**2

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

# Plot 1: Path comparison
ax1.plot(x_cycloid, y_cycloid, 'r-', linewidth=3, label=f'Cycloid (t={cycloid_time:.3f}s)')
ax1.plot(x_straight, y_straight, 'b--', linewidth=2, label=f'Straight (t={straight_time:.3f}s)')
ax1.plot(x_parabolic, y_parabolic, 'g:', linewidth=2, label=f'Parabolic (t={parabolic_time:.3f}s)')
ax1.plot([x_start, x_end], [y_start, y_end], 'ko', markersize=8)
ax1.set_xlabel('x (m)')
ax1.set_ylabel('y (m)')
ax1.set_title('Path Comparison: Brachistochrone Problem')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.invert_yaxis()  # Invert y-axis to show downward motion

# Plot 2: Speed along cycloid path
theta_speed = np.linspace(1e-6, theta_end_opt, n_plot)
x_speed, y_speed = cycloid_parametric(theta_speed, R_opt)
speed_cycloid = np.sqrt(2 * g * y_speed)

ax2.plot(x_speed, speed_cycloid, 'r-', linewidth=2)
ax2.set_xlabel('x (m)')
ax2.set_ylabel('Speed (m/s)')
ax2.set_title('Speed Along Cycloid Path')
ax2.grid(True, alpha=0.3)

# Plot 3: Cycloid construction visualization
circle_positions = np.linspace(0, x_end, 8)
ax3.plot(x_cycloid, y_cycloid, 'r-', linewidth=3, label='Cycloid')

# Show rolling circle construction
for i, x_pos in enumerate(circle_positions[1:]):
    # Find corresponding theta
    theta_pos = theta_end_opt * x_pos / x_end
    x_center, y_center = x_pos, R_opt
    
    circle = plt.Circle((x_center, y_center), R_opt, fill=False, 
                       color='blue', alpha=0.3, linewidth=1)
    ax3.add_patch(circle)
    
    # Mark the tracing point
    x_trace, y_trace = cycloid_parametric(theta_pos, R_opt)
    ax3.plot(x_trace, y_trace, 'ro', markersize=4, alpha=0.7)

ax3.axhline(y=R_opt, color='k', linestyle='-', alpha=0.5, label='Rolling line')
ax3.set_xlabel('x (m)')
ax3.set_ylabel('y (m)')
ax3.set_title('Cycloid Construction (Rolling Circle)')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_aspect('equal')

# Plot 4: Time comparison bar chart
paths = ['Cycloid\n(Optimal)', 'Straight\nLine', 'Parabolic\nPath']
times = [cycloid_time, straight_time, parabolic_time]
colors = ['red', 'blue', 'green']

bars = ax4.bar(paths, times, color=colors, alpha=0.7)
ax4.set_ylabel('Travel Time (seconds)')
ax4.set_title('Travel Time Comparison')
ax4.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar, time in zip(bars, times):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height + 0.001,
             f'{time:.3f}s', ha='center', va='bottom', fontweight='bold')

plt.tight_layout()
plt.show()

# Additional analysis: Calculus of variations insight
print("\n=== MATHEMATICAL INSIGHT ===")
print("The Euler-Lagrange equation for this problem leads to:")
print("d/dx[∂L/∂y'] - ∂L/∂y = 0")
print("where L = √(1 + (y')²) / √(2gy)")
print()
print("This yields the first integral:")
print("√(1 + (y')²) / √(2gy) = constant = 1/√(2gR)")
print()
print("Solving this differential equation gives the cycloid solution:")
print("x = R(θ - sin θ)")
print("y = R(1 - cos θ)")
print(f"with R = {R_opt:.4f} m for our specific boundary conditions.")

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize_scalar, minimize
from scipy.special import jv  # Bessel function
import seaborn as sns

# Set style for better plots
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

class QuantumClock:
    """
    A class to simulate and optimize quantum clock performance
    """
    
    def __init__(self, N_ions=10, gamma=1e-3, omega_0=1e15):
        """
        Initialize quantum clock parameters
        
        Parameters:
        - N_ions: Number of trapped ions
        - gamma: Decoherence rate (Hz)
        - omega_0: Transition frequency (Hz)
        """
        self.N_ions = N_ions
        self.gamma = gamma  # Decoherence rate
        self.omega_0 = omega_0  # Transition frequency
        self.hbar = 1.055e-34  # Reduced Planck constant
        
    def quantum_fisher_information(self, T, entangled=True):
        """
        Calculate quantum Fisher information for the clock
        
        The QFI determines the ultimate precision limit via Cramér-Rao bound:
        var(ω) ≥ 1/F_Q
        
        For entangled states: F_Q = N²T² (Heisenberg limit)
        For separable states: F_Q = NT² (Standard quantum limit)
        """
        if entangled:
            # Heisenberg-limited scaling with decoherence
            decoherence_factor = np.exp(-self.gamma * T)
            return self.N_ions**2 * T**2 * decoherence_factor
        else:
            # Standard quantum limit
            decoherence_factor = np.exp(-self.gamma * T)
            return self.N_ions * T**2 * decoherence_factor
    
    def frequency_uncertainty(self, T, entangled=True):
        """
        Calculate frequency uncertainty using Cramér-Rao bound
        Δω ≥ 1/√F_Q
        """
        F_Q = self.quantum_fisher_information(T, entangled)
        return 1.0 / np.sqrt(F_Q)
    
    def allan_deviation(self, T, tau_measurement):
        """
        Calculate Allan deviation for clock stability
        
        Allan deviation characterizes frequency stability:
        σ_A(τ) = Δω/ω₀ × √(T/2τ)
        
        where τ is the measurement time
        """
        delta_omega = self.frequency_uncertainty(T, entangled=True)
        return (delta_omega / self.omega_0) * np.sqrt(T / (2 * tau_measurement))
    
    def optimize_interrogation_time(self, T_range=(1e-6, 1e-2)):
        """
        Find optimal interrogation time that minimizes frequency uncertainty
        """
        def objective(T):
            return self.frequency_uncertainty(T, entangled=True)
        
        result = minimize_scalar(objective, bounds=T_range, method='bounded')
        return result.x, result.fun
    
    def sensitivity_analysis(self, param_name, param_range, T_optimal):
        """
        Analyze sensitivity to parameter changes
        """
        original_value = getattr(self, param_name)
        sensitivities = []
        
        for param_value in param_range:
            setattr(self, param_name, param_value)
            uncertainty = self.frequency_uncertainty(T_optimal, entangled=True)
            sensitivities.append(uncertainty)
        
        # Restore original value
        setattr(self, param_name, original_value)
        return np.array(sensitivities)

# Initialize quantum clock system
qclock = QuantumClock(N_ions=50, gamma=1e-4, omega_0=1e15)

# Define time ranges for analysis
T_range = np.logspace(-6, -2, 100)  # 1 μs to 10 ms
tau_range = np.logspace(0, 4, 50)   # 1 s to 10,000 s

print("=== Quantum Clock Optimization Analysis ===")
print(f"System parameters:")
print(f"- Number of ions: {qclock.N_ions}")
print(f"- Decoherence rate: {qclock.gamma:.1e} Hz")
print(f"- Transition frequency: {qclock.omega_0:.1e} Hz")

# Find optimal interrogation time
T_opt, uncertainty_opt = qclock.optimize_interrogation_time()
print(f"\nOptimal interrogation time: {T_opt*1e6:.2f} μs")
print(f"Minimum frequency uncertainty: {uncertainty_opt:.2e} Hz")

# Calculate quantum Fisher information and uncertainties
F_Q_entangled = [qclock.quantum_fisher_information(T, entangled=True) for T in T_range]
F_Q_separable = [qclock.quantum_fisher_information(T, entangled=False) for T in T_range]

uncertainty_entangled = [qclock.frequency_uncertainty(T, entangled=True) for T in T_range]
uncertainty_separable = [qclock.frequency_uncertainty(T, entangled=False) for T in T_range]

# Allan deviation calculation
allan_dev = [qclock.allan_deviation(T_opt, tau) for tau in tau_range]

# Sensitivity analysis
N_ions_range = np.arange(10, 101, 10)
gamma_range = np.logspace(-6, -2, 20)

sensitivity_N = qclock.sensitivity_analysis('N_ions', N_ions_range, T_opt)
sensitivity_gamma = qclock.sensitivity_analysis('gamma', gamma_range, T_opt)

print(f"\nQuantum advantage factor: {uncertainty_separable[50]/uncertainty_entangled[50]:.1f}x")
print("Analysis complete. Generating visualizations...")

# ===============================
# VISUALIZATION SECTION
# ===============================

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

# 1. Quantum Fisher Information comparison
ax1 = plt.subplot(2, 3, 1)
plt.loglog(T_range*1e6, F_Q_entangled, 'b-', linewidth=2, label='Entangled (Heisenberg limit)')
plt.loglog(T_range*1e6, F_Q_separable, 'r--', linewidth=2, label='Separable (SQL)')
plt.xlabel('Interrogation Time (μs)')
plt.ylabel('Quantum Fisher Information')
plt.title('Quantum Fisher Information vs Time')
plt.legend()
plt.grid(True, alpha=0.3)

# 2. Frequency uncertainty comparison
ax2 = plt.subplot(2, 3, 2)
plt.loglog(T_range*1e6, uncertainty_entangled, 'b-', linewidth=2, label='Entangled states')
plt.loglog(T_range*1e6, uncertainty_separable, 'r--', linewidth=2, label='Separable states')
plt.axvline(T_opt*1e6, color='green', linestyle=':', linewidth=2, label=f'Optimal T = {T_opt*1e6:.1f} μs')
plt.xlabel('Interrogation Time (μs)')
plt.ylabel('Frequency Uncertainty (Hz)')
plt.title('Frequency Uncertainty vs Time')
plt.legend()
plt.grid(True, alpha=0.3)

# 3. Allan deviation
ax3 = plt.subplot(2, 3, 3)
plt.loglog(tau_range, allan_dev, 'purple', linewidth=2)
plt.xlabel('Measurement Time τ (s)')
plt.ylabel('Allan Deviation σ_A(τ)')
plt.title('Clock Stability (Allan Deviation)')
plt.grid(True, alpha=0.3)

# 4. Sensitivity to number of ions
ax4 = plt.subplot(2, 3, 4)
plt.plot(N_ions_range, sensitivity_N, 'go-', linewidth=2, markersize=6)
plt.xlabel('Number of Ions')
plt.ylabel('Frequency Uncertainty (Hz)')
plt.title('Sensitivity to Ion Number')
plt.grid(True, alpha=0.3)

# 5. Sensitivity to decoherence rate
ax5 = plt.subplot(2, 3, 5)
plt.semilogx(gamma_range, sensitivity_gamma, 'ro-', linewidth=2, markersize=4)
plt.xlabel('Decoherence Rate γ (Hz)')
plt.ylabel('Frequency Uncertainty (Hz)')
plt.title('Sensitivity to Decoherence')
plt.grid(True, alpha=0.3)

# 6. 3D surface plot: Uncertainty vs N_ions and T
ax6 = plt.subplot(2, 3, 6, projection='3d')
N_mesh = np.arange(10, 51, 5)
T_mesh = np.logspace(-6, -3, 20)
N_grid, T_grid = np.meshgrid(N_mesh, T_mesh)
Z_grid = np.zeros_like(N_grid)

for i, N in enumerate(N_mesh):
    temp_clock = QuantumClock(N_ions=N, gamma=qclock.gamma, omega_0=qclock.omega_0)
    for j, T in enumerate(T_mesh):
        Z_grid[j, i] = temp_clock.frequency_uncertainty(T, entangled=True)

surf = ax6.plot_surface(N_grid, T_grid*1e6, Z_grid, cmap='viridis', alpha=0.8)
ax6.set_xlabel('Number of Ions')
ax6.set_ylabel('Time (μs)')
ax6.set_zlabel('Uncertainty (Hz)')
ax6.set_title('3D Optimization Surface')

plt.tight_layout()
plt.show()

# ===============================
# DETAILED ANALYSIS OUTPUT
# ===============================

print("\n=== DETAILED RESULTS ===")
print(f"Optimal Parameters:")
print(f"- Interrogation time: {T_opt*1e6:.3f} μs")
print(f"- Minimum uncertainty: {uncertainty_opt:.3e} Hz")
print(f"- Relative precision: {uncertainty_opt/qclock.omega_0:.3e}")

print(f"\nQuantum Advantage:")
idx_opt = np.argmin(np.abs(T_range - T_opt))
quantum_advantage = uncertainty_separable[idx_opt] / uncertainty_entangled[idx_opt]
print(f"- Entangled vs separable states: {quantum_advantage:.1f}x improvement")

print(f"\nScaling Analysis:")
print(f"- Heisenberg scaling: Δω ∝ 1/N (N = number of ions)")
print(f"- Standard quantum limit: Δω ∝ 1/√N")
print(f"- Decoherence impact: exp(-γT) suppression")

# Calculate theoretical limits
theoretical_limit = 1.0 / (qclock.N_ions * T_opt)
print(f"\nTheoretical vs Achieved:")
print(f"- Theoretical limit: {theoretical_limit:.3e} Hz")
print(f"- Achieved precision: {uncertainty_opt:.3e} Hz")
print(f"- Efficiency: {theoretical_limit/uncertainty_opt*100:.1f}%")

def quantum_fisher_information(self, T, entangled=True):
    if entangled:
        decoherence_factor = np.exp(-self.gamma * T)
        return self.N_ions**2 * T**2 * decoherence_factor

def objective(T):
    return self.frequency_uncertainty(T, entangled=True)