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

# Constants
g = 1.62              # Lunar gravity (m/s^2)
m = 1000              # mass of lander (kg)
T = 10                # total time (s)
N = 100               # number of time steps
dt = T / N
h0, v0 = 1000, 0      # initial height and velocity
u_max = 15000         # max thrust (N)

# Initial guess: uniform thrust
u0 = np.ones(N) * (m * g)

# Bounds for thrust: 0 <= u <= u_max
bounds = [(0, u_max) for _ in range(N)]

# Dynamics integration
def simulate(u):
    h = np.zeros(N+1)
    v = np.zeros(N+1)
    h[0] = h0
    v[0] = v0
    for i in range(N):
        a = u[i] / m - g
        v[i+1] = v[i] + a * dt
        h[i+1] = h[i] - v[i] * dt
    return h, v

# Objective: minimize total thrust (fuel surrogate)
def objective(u):
    return np.sum(u) * dt

# Constraints for final state: height = 0, velocity = 0
def constraint_final_state(u):
    h, v = simulate(u)
    return [h[-1], v[-1]]

# Define constraint dict
constraints = {'type': 'eq', 'fun': constraint_final_state}

# Solve optimization
result = minimize(objective, u0, bounds=bounds, constraints=constraints)

# Extract result
u_opt = result.x
h_opt, v_opt = simulate(u_opt)
t = np.linspace(0, T, N+1)

# Plotting
plt.figure(figsize=(12, 6))

plt.subplot(3, 1, 1)
plt.plot(t, h_opt, label='Height (m)')
plt.ylabel('Height (m)')
plt.grid(True)
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(t, v_opt, label='Velocity (m/s)', color='orange')
plt.ylabel('Velocity (m/s)')
plt.grid(True)
plt.legend()

plt.subplot(3, 1, 3)
plt.plot(t[:-1], u_opt, label='Thrust (N)', color='green')
plt.ylabel('Thrust (N)')
plt.xlabel('Time (s)')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.integrate import quad
import seaborn as sns

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

# Physical constants
c = 299792.458  # Speed of light in km/s

class DarkEnergyModel:
    """
    Lambda-CDM dark energy model for cosmological parameter optimization
    """
    
    def __init__(self):
        self.name = "Lambda-CDM"
    
    def hubble_parameter(self, z, h0, omega_m, omega_lambda):
        """
        Calculate the Hubble parameter H(z) at redshift z
        
        Parameters:
        z: redshift
        h0: Hubble constant (km/s/Mpc)
        omega_m: matter density parameter
        omega_lambda: dark energy density parameter
        """
        return h0 * np.sqrt(omega_m * (1 + z)**3 + omega_lambda)
    
    def E_function(self, z, omega_m, omega_lambda):
        """
        Dimensionless Hubble parameter E(z) = H(z)/H0
        """
        return np.sqrt(omega_m * (1 + z)**3 + omega_lambda)
    
    def luminosity_distance(self, z, h0, omega_m, omega_lambda):
        """
        Calculate luminosity distance using numerical integration
        """
        def integrand(zp):
            return 1.0 / self.E_function(zp, omega_m, omega_lambda)
        
        # Handle single value or array
        if np.isscalar(z):
            if z == 0:
                return 0.0
            integral, _ = quad(integrand, 0, z)
            return (c * (1 + z) / h0) * integral
        else:
            distances = []
            for zi in z:
                if zi == 0:
                    distances.append(0.0)
                else:
                    integral, _ = quad(integrand, 0, zi)
                    distances.append((c * (1 + zi) / h0) * integral)
            return np.array(distances)
    
    def distance_modulus(self, z, h0, omega_m, omega_lambda):
        """
        Calculate distance modulus μ(z) = 5*log10(dL/10pc)
        """
        dl = self.luminosity_distance(z, h0, omega_m, omega_lambda)
        # Convert from Mpc to pc and calculate distance modulus
        dl_pc = dl * 1e6  # Convert Mpc to pc
        return 5 * np.log10(dl_pc / 10.0)

class SupernovaDataGenerator:
    """
    Generate synthetic Type Ia supernova data for testing
    """
    
    def __init__(self, model):
        self.model = model
    
    def generate_data(self, n_points=50, z_max=2.0, noise_level=0.1):
        """
        Generate synthetic supernova data with realistic noise
        
        Parameters:
        n_points: number of data points
        z_max: maximum redshift
        noise_level: standard deviation of observational errors
        """
        # True cosmological parameters (Planck 2018 values)
        true_h0 = 67.4
        true_omega_m = 0.315
        true_omega_lambda = 0.685
        
        # Generate redshift points (more dense at low z, like real surveys)
        z_data = np.random.exponential(0.3, n_points)
        z_data = z_data[z_data <= z_max]
        z_data = np.sort(z_data)
        
        # Calculate theoretical distance moduli
        mu_theory = self.model.distance_modulus(z_data, true_h0, true_omega_m, true_omega_lambda)
        
        # Add realistic noise
        mu_errors = np.random.normal(0, noise_level, len(z_data))
        mu_observed = mu_theory + mu_errors
        
        # Observational uncertainties (typical for SN Ia)
        mu_uncertainties = np.random.uniform(0.05, 0.2, len(z_data))
        
        return {
            'redshift': z_data,
            'distance_modulus': mu_observed,
            'errors': mu_uncertainties,
            'true_params': {'h0': true_h0, 'omega_m': true_omega_m, 'omega_lambda': true_omega_lambda}
        }

class CosmologyOptimizer:
    """
    Optimize cosmological parameters using chi-squared minimization
    """
    
    def __init__(self, model, data):
        self.model = model
        self.data = data
    
    def chi_squared(self, params):
        """
        Calculate chi-squared statistic for given parameters
        """
        h0, omega_m, omega_lambda = params
        
        # Physical constraints
        if h0 <= 0 or omega_m <= 0 or omega_lambda <= 0:
            return 1e10
        if omega_m + omega_lambda < 0.5 or omega_m + omega_lambda > 1.5:
            return 1e10
        
        # Calculate theoretical predictions
        mu_theory = self.model.distance_modulus(
            self.data['redshift'], h0, omega_m, omega_lambda
        )
        
        # Calculate chi-squared
        residuals = (self.data['distance_modulus'] - mu_theory) / self.data['errors']
        chi2 = np.sum(residuals**2)
        
        return chi2
    
    def optimize(self, initial_guess=None):
        """
        Find best-fit parameters using optimization
        """
        if initial_guess is None:
            initial_guess = [70.0, 0.3, 0.7]  # Reasonable starting values
        
        # Set parameter bounds
        bounds = [(50, 100), (0.1, 0.5), (0.5, 1.0)]
        
        # Perform optimization
        result = minimize(
            self.chi_squared,
            initial_guess,
            method='L-BFGS-B',
            bounds=bounds,
            options={'maxiter': 1000}
        )
        
        return {
            'best_fit_params': result.x,
            'chi_squared': result.fun,
            'success': result.success,
            'message': result.message
        }
    
    def parameter_confidence_intervals(self, best_fit_params, n_samples=1000):
        """
        Estimate parameter uncertainties using bootstrap-like sampling
        """
        chi2_min = self.chi_squared(best_fit_params)
        
        # Sample around best-fit parameters
        h0_samples = np.random.normal(best_fit_params[0], 2.0, n_samples)
        omega_m_samples = np.random.normal(best_fit_params[1], 0.05, n_samples)
        omega_lambda_samples = np.random.normal(best_fit_params[2], 0.05, n_samples)
        
        valid_samples = []
        chi2_values = []
        
        for i in range(n_samples):
            params = [h0_samples[i], omega_m_samples[i], omega_lambda_samples[i]]
            chi2 = self.chi_squared(params)
            
            # Accept samples within reasonable chi-squared range
            if chi2 < chi2_min + 10:  # Delta chi-squared = 10 for rough 3-sigma
                valid_samples.append(params)
                chi2_values.append(chi2)
        
        valid_samples = np.array(valid_samples)
        
        if len(valid_samples) > 10:
            uncertainties = np.std(valid_samples, axis=0)
            return uncertainties
        else:
            return np.array([1.0, 0.01, 0.01])  # Default uncertainties

# Initialize the dark energy model
model = DarkEnergyModel()

# Generate synthetic supernova data
data_generator = SupernovaDataGenerator(model)
sn_data = data_generator.generate_data(n_points=100, z_max=1.5, noise_level=0.15)

print("Generated Supernova Dataset:")
print(f"Number of supernovae: {len(sn_data['redshift'])}")
print(f"Redshift range: {sn_data['redshift'].min():.3f} - {sn_data['redshift'].max():.3f}")
print(f"True parameters:")
print(f"  H₀ = {sn_data['true_params']['h0']:.1f} km/s/Mpc")
print(f"  Ωₘ = {sn_data['true_params']['omega_m']:.3f}")
print(f"  ΩΛ = {sn_data['true_params']['omega_lambda']:.3f}")

# Initialize optimizer
optimizer = CosmologyOptimizer(model, sn_data)

# Perform optimization
print("\nOptimizing cosmological parameters...")
result = optimizer.optimize(initial_guess=[68.0, 0.32, 0.68])

if result['success']:
    print("Optimization successful!")
    best_params = result['best_fit_params']
    print(f"Best-fit parameters:")
    print(f"  H₀ = {best_params[0]:.2f} km/s/Mpc")
    print(f"  Ωₘ = {best_params[1]:.4f}")
    print(f"  ΩΛ = {best_params[2]:.4f}")
    print(f"  χ² = {result['chi_squared']:.2f}")
    
    # Estimate uncertainties
    uncertainties = optimizer.parameter_confidence_intervals(best_params)
    print(f"\nParameter uncertainties (1σ estimates):")
    print(f"  σ(H₀) = ±{uncertainties[0]:.2f} km/s/Mpc")
    print(f"  σ(Ωₘ) = ±{uncertainties[1]:.4f}")
    print(f"  σ(ΩΛ) = ±{uncertainties[2]:.4f}")
else:
    print("Optimization failed:", result['message'])

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

# Plot 1: Hubble diagram (distance modulus vs redshift)
z_theory = np.linspace(0.01, 1.5, 100)
mu_true = model.distance_modulus(z_theory, 
                                sn_data['true_params']['h0'],
                                sn_data['true_params']['omega_m'],
                                sn_data['true_params']['omega_lambda'])
mu_fit = model.distance_modulus(z_theory, *best_params)

ax1.errorbar(sn_data['redshift'], sn_data['distance_modulus'], 
             yerr=sn_data['errors'], fmt='o', alpha=0.7, color='blue', 
             markersize=4, label='Simulated SN Ia data')
ax1.plot(z_theory, mu_true, 'r-', linewidth=2, label='True model')
ax1.plot(z_theory, mu_fit, 'g--', linewidth=2, label='Best fit')
ax1.set_xlabel('Redshift (z)')
ax1.set_ylabel('Distance Modulus (μ)')
ax1.set_title('Hubble Diagram: Type Ia Supernovae')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Residuals
mu_data_fit = model.distance_modulus(sn_data['redshift'], *best_params)
residuals = sn_data['distance_modulus'] - mu_data_fit
standardized_residuals = residuals / sn_data['errors']

ax2.errorbar(sn_data['redshift'], standardized_residuals, 
             yerr=np.ones_like(sn_data['errors']), fmt='o', alpha=0.7, color='red', markersize=4)
ax2.axhline(y=0, color='black', linestyle='-', alpha=0.8)
ax2.axhline(y=2, color='gray', linestyle='--', alpha=0.5)
ax2.axhline(y=-2, color='gray', linestyle='--', alpha=0.5)
ax2.set_xlabel('Redshift (z)')
ax2.set_ylabel('Standardized Residuals (σ)')
ax2.set_title('Fit Residuals')
ax2.grid(True, alpha=0.3)
ax2.set_ylim(-4, 4)

# Plot 3: Hubble parameter evolution
H_true = [model.hubble_parameter(z, sn_data['true_params']['h0'],
                                sn_data['true_params']['omega_m'],
                                sn_data['true_params']['omega_lambda']) for z in z_theory]
H_fit = [model.hubble_parameter(z, *best_params) for z in z_theory]

ax3.plot(z_theory, H_true, 'r-', linewidth=2, label='True H(z)')
ax3.plot(z_theory, H_fit, 'g--', linewidth=2, label='Best-fit H(z)')
ax3.set_xlabel('Redshift (z)')
ax3.set_ylabel('Hubble Parameter H(z) [km/s/Mpc]')
ax3.set_title('Hubble Parameter Evolution')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Parameter comparison
param_names = ['H₀\n[km/s/Mpc]', 'Ωₘ', 'ΩΛ']
true_values = [sn_data['true_params']['h0'], 
               sn_data['true_params']['omega_m'], 
               sn_data['true_params']['omega_lambda']]
fit_values = best_params
fit_errors = uncertainties

x_pos = np.arange(len(param_names))
ax4.bar(x_pos - 0.2, true_values, 0.4, label='True values', color='red', alpha=0.7)
ax4.bar(x_pos + 0.2, fit_values, 0.4, yerr=fit_errors, 
        label='Best fit ± 1σ', color='green', alpha=0.7, capsize=5)
ax4.set_xlabel('Parameters')
ax4.set_ylabel('Parameter Values')
ax4.set_title('Parameter Recovery')
ax4.set_xticks(x_pos)
ax4.set_xticklabels(param_names)
ax4.legend()
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Additional analysis: Chi-squared surface
print("\nGenerating chi-squared contour analysis...")

# Create parameter grids for contour plot
h0_range = np.linspace(best_params[0] - 5, best_params[0] + 5, 30)
omega_m_range = np.linspace(max(0.1, best_params[1] - 0.1), 
                           min(0.5, best_params[1] + 0.1), 30)

H0_grid, Om_grid = np.meshgrid(h0_range, omega_m_range)
chi2_grid = np.zeros_like(H0_grid)

# Fix omega_lambda at best-fit value for 2D slice
omega_lambda_fixed = best_params[2]

for i in range(len(omega_m_range)):
    for j in range(len(h0_range)):
        params = [H0_grid[i,j], Om_grid[i,j], omega_lambda_fixed]
        chi2_grid[i,j] = optimizer.chi_squared(params)

# Create contour plot
fig, ax = plt.subplots(1, 1, figsize=(10, 8))

# Calculate confidence level contours
chi2_min = result['chi_squared']
levels = [chi2_min + 2.3, chi2_min + 6.17, chi2_min + 11.8]  # 1σ, 2σ, 3σ for 2 parameters

contour = ax.contour(H0_grid, Om_grid, chi2_grid, levels=levels, 
                    colors=['green', 'orange', 'red'], linestyles=['-', '--', ':'])
ax.clabel(contour, inline=True, fontsize=10, fmt='%.1f')

# Mark true and best-fit values
ax.plot(sn_data['true_params']['h0'], sn_data['true_params']['omega_m'], 
        'r*', markersize=15, label='True values')
ax.plot(best_params[0], best_params[1], 'go', markersize=10, label='Best fit')

ax.set_xlabel('H₀ [km/s/Mpc]')
ax.set_ylabel('Ωₘ')
ax.set_title('χ² Confidence Contours (ΩΛ fixed)')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Summary statistics
print(f"\nFinal Analysis Summary:")
print(f"Data points used: {len(sn_data['redshift'])}")
print(f"Degrees of freedom: {len(sn_data['redshift']) - 3}")
print(f"Reduced χ²: {result['chi_squared']/(len(sn_data['redshift']) - 3):.3f}")
print(f"Parameter recovery accuracy:")
print(f"  H₀: {abs(best_params[0] - sn_data['true_params']['h0'])/sn_data['true_params']['h0']*100:.2f}% difference")
print(f"  Ωₘ: {abs(best_params[1] - sn_data['true_params']['omega_m'])/sn_data['true_params']['omega_m']*100:.2f}% difference")
print(f"  ΩΛ: {abs(best_params[2] - sn_data['true_params']['omega_lambda'])/sn_data['true_params']['omega_lambda']*100:.2f}% difference")

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

# Physical constants
STEFAN_BOLTZMANN = 5.67e-8  # W/m²K⁴
SOLAR_CONSTANT = 1361       # W/m² (solar flux at Earth orbit)

class SpacecraftThermalOptimizer:
    """
    Spacecraft thermal design optimization considering:
    - Multi-layer insulation (MLI)
    - Radiative heat transfer
    - Weight constraints
    - Cost optimization
    """
    
    def __init__(self):
        # Spacecraft parameters
        self.surface_area = 10.0        # m² (total surface area)
        self.electronics_power = 500    # W (internal heat generation)
        
        # Environmental conditions
        self.T_space = 4               # K (deep space temperature)
        self.T_sun_side = 393          # K (120°C, sun-facing side)
        self.T_shade_side = 173        # K (-100°C, shaded side)
        
        # Component temperature limits
        self.T_min_electronics = 253   # K (-20°C)
        self.T_max_electronics = 343   # K (70°C)
        self.T_min_battery = 273       # K (0°C)
        self.T_max_battery = 323       # K (50°C)
        
        # Material properties
        self.mli_properties = {
            'thermal_conductivity': 0.002,  # W/m·K (very low for MLI)
            'density': 50,                  # kg/m³
            'cost_per_kg': 10000,          # $/kg
            'emissivity_low': 0.03,        # Low-emissivity side
            'emissivity_high': 0.8         # High-emissivity side
        }
        
        # Mission costs
        self.launch_cost_per_kg = 5000  # $/kg
        self.thermal_control_base_cost = 50000  # $
        self.penalty_cost_per_degree = 1000     # $/K (temperature violation)
    
    def calculate_heat_transfer(self, thickness, emissivity):
        """
        Calculate heat transfer through MLI and radiation
        
        Args:
            thickness: MLI thickness in meters
            emissivity: Effective emissivity of outer surface
        
        Returns:
            dict: Heat transfer rates and temperatures
        """
        # Conductive heat transfer through MLI
        k = self.mli_properties['thermal_conductivity']
        q_conduction = (k * self.surface_area * 
                       (self.T_sun_side - self.T_shade_side)) / thickness
        
        # Radiative heat transfer (simplified model)
        # Assuming spacecraft acts as a radiator to space
        T_avg_exterior = np.sqrt((self.T_sun_side**2 + self.T_shade_side**2) / 2)
        
        q_radiation_out = (STEFAN_BOLTZMANN * self.surface_area * emissivity * 
                          (T_avg_exterior**4 - self.T_space**4))
        
        # Solar heating (absorbed)
        absorptivity = 0.3  # Typical for spacecraft surfaces
        q_solar = SOLAR_CONSTANT * self.surface_area * absorptivity * 0.5  # 50% sun exposure
        
        # Heat balance for internal temperature
        q_net = q_solar + self.electronics_power - q_radiation_out - q_conduction
        
        # Estimate internal temperature (simplified)
        # This is a simplified model - real spacecraft use detailed thermal networks
        thermal_mass = 1000  # J/K (thermal capacitance)
        T_internal = self.T_space + (q_net / (STEFAN_BOLTZMANN * self.surface_area * 0.1))
        
        return {
            'q_conduction': q_conduction,
            'q_radiation_out': q_radiation_out,
            'q_solar': q_solar,
            'q_net': q_net,
            'T_internal': T_internal,
            'T_electronics': T_internal + 10,  # Electronics run hotter
            'T_battery': T_internal + 5        # Battery temperature
        }
    
    def calculate_mass_and_cost(self, thickness):
        """Calculate MLI mass and associated costs"""
        volume = self.surface_area * thickness
        mass = volume * self.mli_properties['density']
        
        launch_cost = mass * self.launch_cost_per_kg
        material_cost = mass * self.mli_properties['cost_per_kg']
        
        return mass, launch_cost + material_cost
    
    def temperature_penalty(self, temperatures):
        """Calculate penalty for temperature violations"""
        penalty = 0
        
        T_electronics = temperatures['T_electronics']
        T_battery = temperatures['T_battery']
        
        # Electronics temperature penalties
        if T_electronics < self.T_min_electronics:
            penalty += (self.T_min_electronics - T_electronics) * self.penalty_cost_per_degree
        elif T_electronics > self.T_max_electronics:
            penalty += (T_electronics - self.T_max_electronics) * self.penalty_cost_per_degree
        
        # Battery temperature penalties
        if T_battery < self.T_min_battery:
            penalty += (self.T_min_battery - T_battery) * self.penalty_cost_per_degree
        elif T_battery > self.T_max_battery:
            penalty += (T_battery - self.T_max_battery) * self.penalty_cost_per_degree
        
        return penalty
    
    def objective_function(self, design_vars):
        """
        Objective function to minimize total mission cost
        
        Args:
            design_vars: [thickness (m), emissivity]
        
        Returns:
            float: Total cost ($)
        """
        thickness, emissivity = design_vars
        
        # Constraints
        if thickness <= 0.001 or thickness > 0.1:  # 1mm to 10cm
            return 1e10
        if emissivity <= 0.01 or emissivity > 1.0:
            return 1e10
        
        # Calculate heat transfer
        heat_results = self.calculate_heat_transfer(thickness, emissivity)
        
        # Calculate mass and cost
        mass, material_launch_cost = self.calculate_mass_and_cost(thickness)
        
        # Temperature penalty
        temp_penalty = self.temperature_penalty(heat_results)
        
        # Total cost
        total_cost = (material_launch_cost + 
                     self.thermal_control_base_cost + 
                     temp_penalty)
        
        return total_cost
    
    def optimize_design(self):
        """Perform optimization to find optimal MLI thickness and emissivity"""
        # Initial guess: moderate thickness and emissivity
        x0 = [0.02, 0.5]  # 2cm thickness, 0.5 emissivity
        
        # Bounds
        bounds = [(0.001, 0.1), (0.01, 1.0)]
        
        # Optimization
        result = minimize(self.objective_function, x0, 
                         method='L-BFGS-B', bounds=bounds)
        
        return result
    
    def analyze_design_space(self):
        """Analyze the design space to understand trade-offs"""
        # Create parameter grids
        thickness_range = np.linspace(0.005, 0.08, 50)
        emissivity_range = np.linspace(0.1, 0.9, 40)
        
        # Initialize result arrays
        cost_matrix = np.zeros((len(thickness_range), len(emissivity_range)))
        temp_matrix = np.zeros((len(thickness_range), len(emissivity_range)))
        mass_matrix = np.zeros((len(thickness_range), len(emissivity_range)))
        
        # Calculate for each combination
        for i, thickness in enumerate(thickness_range):
            for j, emissivity in enumerate(emissivity_range):
                cost = self.objective_function([thickness, emissivity])
                heat_results = self.calculate_heat_transfer(thickness, emissivity)
                mass, _ = self.calculate_mass_and_cost(thickness)
                
                cost_matrix[i, j] = cost if cost < 1e9 else np.nan
                temp_matrix[i, j] = heat_results['T_electronics']
                mass_matrix[i, j] = mass
        
        return thickness_range, emissivity_range, cost_matrix, temp_matrix, mass_matrix

# Initialize and run optimization
optimizer = SpacecraftThermalOptimizer()

print("🚀 SPACECRAFT THERMAL DESIGN OPTIMIZATION")
print("=" * 50)

# Perform optimization
print("\n📊 Running optimization...")
result = optimizer.optimize_design()

optimal_thickness, optimal_emissivity = result.x
optimal_cost = result.fun

print(f"\n✅ OPTIMAL DESIGN FOUND:")
print(f"   MLI Thickness: {optimal_thickness*1000:.1f} mm")
print(f"   Surface Emissivity: {optimal_emissivity:.3f}")
print(f"   Total Mission Cost: ${optimal_cost:,.0f}")

# Analyze optimal design performance
heat_results = optimizer.calculate_heat_transfer(optimal_thickness, optimal_emissivity)
mass, material_cost = optimizer.calculate_mass_and_cost(optimal_thickness)

print(f"\n🌡️  THERMAL PERFORMANCE:")
print(f"   Electronics Temperature: {heat_results['T_electronics']:.1f} K ({heat_results['T_electronics']-273:.1f}°C)")
print(f"   Battery Temperature: {heat_results['T_battery']:.1f} K ({heat_results['T_battery']-273:.1f}°C)")
print(f"   Heat Conduction: {heat_results['q_conduction']:.1f} W")
print(f"   Heat Radiation: {heat_results['q_radiation_out']:.1f} W")

print(f"\n⚖️  MASS AND COST BREAKDOWN:")
print(f"   MLI Mass: {mass:.2f} kg")
print(f"   Launch + Material Cost: ${material_cost:,.0f}")
print(f"   Temperature Penalty: ${optimizer.temperature_penalty(heat_results):,.0f}")

print(f"\n🎯 DESIGN VERIFICATION:")
temp_ok = (optimizer.T_min_electronics <= heat_results['T_electronics'] <= optimizer.T_max_electronics and
           optimizer.T_min_battery <= heat_results['T_battery'] <= optimizer.T_max_battery)
print(f"   Temperature Constraints: {'✅ SATISFIED' if temp_ok else '❌ VIOLATED'}")

# Analyze design space
print(f"\n🔍 Analyzing design space...")
thickness_range, emissivity_range, cost_matrix, temp_matrix, mass_matrix = optimizer.analyze_design_space()

# Create comprehensive visualizations
fig = plt.figure(figsize=(20, 15))

# 1. Cost contour plot
ax1 = plt.subplot(2, 3, 1)
T, E = np.meshgrid(thickness_range * 1000, emissivity_range)  # Convert to mm
valid_costs = np.where(np.isfinite(cost_matrix.T), cost_matrix.T, np.nan)
contour = plt.contourf(T, E, valid_costs, levels=20, cmap='viridis')
plt.colorbar(contour, label='Total Cost ($)')
plt.contour(T, E, valid_costs, levels=10, colors='white', alpha=0.6, linewidths=0.5)
plt.plot(optimal_thickness*1000, optimal_emissivity, 'r*', markersize=15, label='Optimal Design')
plt.xlabel('MLI Thickness (mm)')
plt.ylabel('Surface Emissivity')
plt.title('Total Mission Cost Optimization')
plt.legend()
plt.grid(True, alpha=0.3)

# 2. Temperature contour plot
ax2 = plt.subplot(2, 3, 2)
temp_contour = plt.contourf(T, E, temp_matrix.T, levels=20, cmap='coolwarm')
plt.colorbar(temp_contour, label='Electronics Temperature (K)')
# Add temperature constraint lines
plt.contour(T, E, temp_matrix.T, levels=[optimizer.T_min_electronics, optimizer.T_max_electronics], 
           colors=['blue', 'red'], linewidths=2, linestyles=['--', '--'])
plt.plot(optimal_thickness*1000, optimal_emissivity, 'r*', markersize=15, label='Optimal Design')
plt.xlabel('MLI Thickness (mm)')
plt.ylabel('Surface Emissivity')
plt.title('Electronics Temperature Distribution')
plt.legend()
plt.grid(True, alpha=0.3)

# 3. Mass contour plot
ax3 = plt.subplot(2, 3, 3)
mass_contour = plt.contourf(T, E, mass_matrix.T, levels=20, cmap='plasma')
plt.colorbar(mass_contour, label='MLI Mass (kg)')
plt.contour(T, E, mass_matrix.T, levels=10, colors='white', alpha=0.6, linewidths=0.5)
plt.plot(optimal_thickness*1000, optimal_emissivity, 'r*', markersize=15, label='Optimal Design')
plt.xlabel('MLI Thickness (mm)')
plt.ylabel('Surface Emissivity')
plt.title('MLI Mass Distribution')
plt.legend()
plt.grid(True, alpha=0.3)

# 4. Trade-off analysis: Cost vs Temperature
ax4 = plt.subplot(2, 3, 4)
thickness_test = np.linspace(0.005, 0.08, 100)
costs = []
temps = []
masses = []

for t in thickness_test:
    cost = optimizer.objective_function([t, optimal_emissivity])
    if cost < 1e9:
        heat_res = optimizer.calculate_heat_transfer(t, optimal_emissivity)
        mass, _ = optimizer.calculate_mass_and_cost(t)
        costs.append(cost)
        temps.append(heat_res['T_electronics'])
        masses.append(mass)
    else:
        costs.append(np.nan)
        temps.append(np.nan)
        masses.append(np.nan)

plt.plot([t*1000 for t in thickness_test], costs, 'b-', linewidth=2, label='Total Cost')
plt.axvline(optimal_thickness*1000, color='red', linestyle='--', alpha=0.7, label='Optimal Thickness')
plt.xlabel('MLI Thickness (mm)')
plt.ylabel('Total Cost ($)', color='blue')
plt.title('Cost vs Thickness Trade-off')
plt.grid(True, alpha=0.3)

ax4_twin = ax4.twinx()
ax4_twin.plot([t*1000 for t in thickness_test], [t-273 for t in temps], 'g-', linewidth=2, label='Electronics Temp')
ax4_twin.axhline(optimizer.T_min_electronics-273, color='orange', linestyle=':', alpha=0.7, label='Min Temp Limit')
ax4_twin.axhline(optimizer.T_max_electronics-273, color='orange', linestyle=':', alpha=0.7, label='Max Temp Limit')
ax4_twin.set_ylabel('Temperature (°C)', color='green')
ax4_twin.legend(loc='upper right')
ax4.legend(loc='upper left')

# 5. Heat transfer breakdown
ax5 = plt.subplot(2, 3, 5)
heat_components = ['Solar Input', 'Electronics', 'Conduction Loss', 'Radiation Loss']
heat_values = [heat_results['q_solar'], optimizer.electronics_power, 
               -heat_results['q_conduction'], -heat_results['q_radiation_out']]
colors = ['gold', 'orange', 'lightblue', 'lightcoral']

bars = plt.bar(heat_components, heat_values, color=colors, alpha=0.8)
plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)
plt.ylabel('Heat Flow (W)')
plt.title('Heat Transfer Component Analysis')
plt.xticks(rotation=45)
plt.grid(True, alpha=0.3, axis='y')

# Add value labels on bars
for bar, value in zip(bars, heat_values):
    plt.text(bar.get_x() + bar.get_width()/2, value + (5 if value > 0 else -15), 
             f'{value:.0f}W', ha='center', va='bottom' if value > 0 else 'top')

# 6. Design sensitivity analysis
ax6 = plt.subplot(2, 3, 6)
# Vary thickness around optimal
thickness_sensitivity = np.linspace(optimal_thickness*0.5, optimal_thickness*2, 50)
cost_sensitivity = []
temp_sensitivity = []

for t in thickness_sensitivity:
    cost = optimizer.objective_function([t, optimal_emissivity])
    heat_res = optimizer.calculate_heat_transfer(t, optimal_emissivity)
    cost_sensitivity.append(cost if cost < 1e9 else np.nan)
    temp_sensitivity.append(heat_res['T_electronics'])

plt.plot([t*1000 for t in thickness_sensitivity], cost_sensitivity, 'b-', linewidth=2, label='Cost')
plt.axvline(optimal_thickness*1000, color='red', linestyle='--', alpha=0.7, label='Optimal')
plt.xlabel('MLI Thickness (mm)')
plt.ylabel('Total Cost ($)', color='blue')
plt.title('Design Sensitivity Analysis')
plt.grid(True, alpha=0.3)

ax6_twin = ax6.twinx()
ax6_twin.plot([t*1000 for t in thickness_sensitivity], [t-273 for t in temp_sensitivity], 
              'g-', linewidth=2, label='Temperature')
ax6_twin.set_ylabel('Temperature (°C)', color='green')
ax6_twin.axhline(optimizer.T_min_electronics-273, color='orange', linestyle=':', alpha=0.5)
ax6_twin.axhline(optimizer.T_max_electronics-273, color='orange', linestyle=':', alpha=0.5)

plt.tight_layout()
plt.suptitle('SPACECRAFT THERMAL DESIGN OPTIMIZATION ANALYSIS', 
             fontsize=16, fontweight='bold', y=0.98)
plt.show()

# Summary statistics
print(f"\n📈 DESIGN SPACE ANALYSIS:")
print(f"   Design Space Explored: {len(thickness_range)} × {len(emissivity_range)} = {len(thickness_range)*len(emissivity_range)} combinations")
print(f"   Feasible Designs: {np.sum(np.isfinite(cost_matrix))}")
print(f"   Cost Range: ${np.nanmin(cost_matrix):,.0f} - ${np.nanmax(cost_matrix):,.0f}")
print(f"   Temperature Range: {np.nanmin(temp_matrix):.1f}K - {np.nanmax(temp_matrix):.1f}K")
print(f"   Mass Range: {np.nanmin(mass_matrix):.2f}kg - {np.nanmax(mass_matrix):.2f}kg")

print(f"\n🎯 OPTIMIZATION CONVERGENCE:")
print(f"   Optimization Success: {'✅ YES' if result.success else '❌ NO'}")
print(f"   Function Evaluations: {result.nfev}")
print(f"   Final Message: {result.message}")

# Conductive heat transfer through MLI
q_conduction = (k * self.surface_area * 
               (self.T_sun_side - self.T_shade_side)) / thickness

q_radiation_out = (STEFAN_BOLTZMANN * self.surface_area * emissivity * 
                  (T_avg_exterior**4 - self.T_space**4))

import numpy as np
import matplotlib.pyplot as plt

# Constants
R_e = 6371  # Earth radius in km
h = 550     # Satellite altitude in km (LEO)
R = R_e + h

# Angular radius (coverage half-angle)
theta = np.arccos(R_e / R)

# Coverage area per satellite (spherical cap)
A_sat = 2 * np.pi * R_e**2 * (1 - np.cos(theta))
A_earth = 4 * np.pi * R_e**2

# Brute-force search: estimate required number of satellites
min_satellites = int(np.ceil(A_earth / A_sat))
print(f"Minimum satellites needed for full coverage: {min_satellites}")

# Visualizing coverage (simplified 2D projection)
def plot_constellation(n_sats, title="Satellite Footprints on Earth"):
    fig, ax = plt.subplots(figsize=(8, 4))
    ax.set_title(title)
    ax.set_xlim(-180, 180)
    ax.set_ylim(-90, 90)
    ax.set_xlabel("Longitude")
    ax.set_ylabel("Latitude")

    # Distribute satellites evenly in latitude bands
    lats = np.linspace(-60, 60, int(np.sqrt(n_sats)))
    lons = np.linspace(-180, 180, int(np.sqrt(n_sats)))
    for lat in lats:
        for lon in lons:
            circle = plt.Circle((lon, lat), np.degrees(theta), color='skyblue', alpha=0.4)
            ax.add_patch(circle)
    ax.grid(True)
    plt.show()

plot_constellation(min_satellites)

import numpy as np
import matplotlib.pyplot as plt

# Gravitational constant in kpc (km/s)^2 / Msun
G = 4.302e-6

# Create radial distance values in kiloparsecs (kpc)
radii = np.linspace(0.1, 30, 100)  # Avoid r = 0 to prevent division by zero

# Simulate a realistic rotation curve (in km/s)
# Rising quickly then flattening, e.g., Milky Way-like
def rotation_curve(r):
    v_max = 220  # km/s
    r_turn = 5   # kpc
    return v_max * r / np.sqrt(r**2 + r_turn**2)

velocities = rotation_curve(radii)

# Calculate enclosed mass M(r) from v(r)
mass_enclosed = velocities**2 * radii / G  # In solar masses

# Plotting the rotation curve
plt.figure(figsize=(14, 5))

plt.subplot(1, 2, 1)
plt.plot(radii, velocities, label='Rotation Curve', color='blue')
plt.xlabel('Radius [kpc]')
plt.ylabel('Velocity [km/s]')
plt.title('Synthetic Galactic Rotation Curve')
plt.grid(True)
plt.legend()

# Plotting the mass distribution
plt.subplot(1, 2, 2)
plt.plot(radii, mass_enclosed, label='Enclosed Mass', color='darkgreen')
plt.xlabel('Radius [kpc]')
plt.ylabel(r'Mass $M(<r)$ [$M_\odot$]')
plt.title('Enclosed Mass vs Radius')
plt.grid(True)
plt.legend()

plt.tight_layout()
plt.show()

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

# Physical constants (in geometric units where G = c = 1)
class BlackHoleAccretion:
    def __init__(self, mass_bh=1.0):
        """
        Initialize black hole accretion system
        mass_bh: Black hole mass in solar masses (normalized to 1)
        """
        self.M = mass_bh  # Black hole mass
        self.rs = 2 * self.M  # Schwarzschild radius
        self.r_isco = 6 * self.M  # Innermost stable circular orbit
        
    def efficiency(self, r):
        """
        Calculate energy extraction efficiency at radius r
        η = 1 - √(1 - 2M/r)
        """
        if r <= self.rs:
            return 0  # Inside event horizon
        return 1 - np.sqrt(1 - 2*self.M/r)
    
    def orbital_velocity(self, r):
        """
        Keplerian orbital velocity: v = √(M/r)
        """
        if r <= self.rs:
            return 0
        return np.sqrt(self.M/r)
    
    def accretion_rate_density(self, r, mdot_0=1.0, alpha=1.5):
        """
        Accretion rate density profile: ρ(r) ∝ r^(-α)
        Normalized so that total accretion rate is mdot_0
        """
        if r <= self.rs:
            return 0
        # Cutoff at ISCO to avoid unphysical solutions
        r_min = max(r, self.r_isco)
        return mdot_0 * (r_min/self.r_isco)**(-alpha)
    
    def luminosity_density(self, r, mdot_0=1.0, alpha=1.5):
        """
        Luminosity density: L(r) = η(r) × ṁ(r) × c²
        """
        return self.efficiency(r) * self.accretion_rate_density(r, mdot_0, alpha)
    
    def total_luminosity(self, r_min=None, r_max=100, mdot_0=1.0, alpha=1.5):
        """
        Integrate luminosity from r_min to r_max
        """
        if r_min is None:
            r_min = self.r_isco
        
        def integrand(r):
            return self.luminosity_density(r, mdot_0, alpha) * r  # Include r for proper integration
        
        result, _ = quad(integrand, r_min, r_max)
        return result
    
    def optimize_accretion_radius(self, r_range=(3, 50)):
        """
        Find radius that maximizes luminosity density
        """
        def neg_luminosity(r):
            return -self.luminosity_density(r)
        
        result = minimize_scalar(neg_luminosity, bounds=r_range, method='bounded')
        return result.x, -result.fun
    
    def optimize_disk_parameters(self, alpha_range=(0.5, 3.0), mdot_range=(0.1, 10.0)):
        """
        Optimize both disk parameters (α, ṁ₀) for maximum total luminosity
        """
        def neg_total_lum(params):
            alpha, mdot_0 = params
            return -self.total_luminosity(mdot_0=mdot_0, alpha=alpha)
        
        # Initial guess
        x0 = [1.5, 1.0]
        bounds = [alpha_range, mdot_range]
        
        result = minimize(neg_total_lum, x0, bounds=bounds, method='L-BFGS-B')
        return result.x, -result.fun

# Create black hole system
bh = BlackHoleAccretion(mass_bh=1.0)

print("Black Hole Accretion Optimization Analysis")
print("=" * 50)
print(f"Black hole mass: {bh.M} M☉")
print(f"Schwarzschild radius: {bh.rs:.2f} GM/c²")
print(f"ISCO radius: {bh.r_isco:.2f} GM/c²")
print()

# Find optimal radius for maximum luminosity density
r_opt, lum_max = bh.optimize_accretion_radius()
print(f"Optimal accretion radius: {r_opt:.2f} GM/c²")
print(f"Maximum luminosity density: {lum_max:.4f}")
print(f"Efficiency at optimal radius: {bh.efficiency(r_opt):.4f}")
print()

# Optimize disk parameters
params_opt, total_lum_max = bh.optimize_disk_parameters()
alpha_opt, mdot_opt = params_opt
print(f"Optimal disk parameters:")
print(f"  α (density exponent): {alpha_opt:.3f}")
print(f"  ṁ₀ (base accretion rate): {mdot_opt:.3f}")
print(f"  Maximum total luminosity: {total_lum_max:.4f}")

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

# Plot 1: Efficiency vs Radius
r_range = np.linspace(2.1, 20, 1000)
efficiency = [bh.efficiency(r) for r in r_range]

ax1.plot(r_range, efficiency, 'b-', linewidth=2, label='Efficiency η(r)')
ax1.axvline(bh.r_isco, color='r', linestyle='--', label=f'ISCO (r = {bh.r_isco})')
ax1.axvline(r_opt, color='g', linestyle=':', label=f'Optimal r = {r_opt:.2f}')
ax1.set_xlabel('Radius r [GM/c²]')
ax1.set_ylabel('Efficiency η')
ax1.set_title('Energy Extraction Efficiency vs Radius')
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.text(0.05, 0.95, r'$\eta = 1 - \sqrt{1 - \frac{2M}{r}}$', 
         transform=ax1.transAxes, fontsize=12, verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

# Plot 2: Luminosity Density vs Radius
lum_density = [bh.luminosity_density(r) for r in r_range]

ax2.plot(r_range, lum_density, 'purple', linewidth=2)
ax2.axvline(bh.r_isco, color='r', linestyle='--', label='ISCO')
ax2.axvline(r_opt, color='g', linestyle=':', label=f'Maximum at r = {r_opt:.2f}')
ax2.scatter([r_opt], [lum_max], color='red', s=100, zorder=5, label=f'Max = {lum_max:.4f}')
ax2.set_xlabel('Radius r [GM/c²]')
ax2.set_ylabel('Luminosity Density')
ax2.set_title('Luminosity Density Distribution')
ax2.grid(True, alpha=0.3)
ax2.legend()

# Plot 3: Parameter Space Optimization
alpha_vals = np.linspace(0.5, 3.0, 50)
mdot_vals = np.linspace(0.1, 5.0, 50)
A, M = np.meshgrid(alpha_vals, mdot_vals)
L = np.zeros_like(A)

for i, alpha in enumerate(alpha_vals):
    for j, mdot in enumerate(mdot_vals):
        L[j, i] = bh.total_luminosity(mdot_0=mdot, alpha=alpha)

contour = ax3.contourf(A, M, L, levels=20, cmap='viridis')
ax3.scatter([alpha_opt], [mdot_opt], color='red', s=100, marker='*', 
           label=f'Optimum: α={alpha_opt:.2f}, ṁ₀={mdot_opt:.2f}')
ax3.set_xlabel('Density Exponent α')
ax3.set_ylabel('Base Accretion Rate ṁ₀')
ax3.set_title('Total Luminosity in Parameter Space')
plt.colorbar(contour, ax=ax3, label='Total Luminosity')
ax3.legend()

# Plot 4: Comparison of Different Scenarios
r_comparison = np.linspace(3, 30, 200)
scenarios = [
    (1.0, 1.0, 'Standard', 'blue'),
    (1.5, 1.0, 'Steeper Profile', 'green'),
    (1.0, 2.0, 'Higher Rate', 'orange'),
    (alpha_opt, mdot_opt, 'Optimized', 'red')
]

for alpha, mdot, label, color in scenarios:
    lum_profile = [bh.luminosity_density(r, mdot_0=mdot, alpha=alpha) for r in r_comparison]
    ax4.plot(r_comparison, lum_profile, color=color, linewidth=2, label=label)

ax4.axvline(bh.r_isco, color='black', linestyle='--', alpha=0.5, label='ISCO')
ax4.set_xlabel('Radius r [GM/c²]')
ax4.set_ylabel('Luminosity Density')
ax4.set_title('Comparison of Accretion Scenarios')
ax4.grid(True, alpha=0.3)
ax4.legend()
ax4.set_xlim(3, 25)

plt.tight_layout()
plt.show()

# Additional analysis: Energy considerations
print("\nDetailed Energy Analysis:")
print("-" * 30)
r_test_points = [bh.r_isco, r_opt, 10, 20, 50]
for r in r_test_points:
    eff = bh.efficiency(r)
    v_orb = bh.orbital_velocity(r)
    lum_dens = bh.luminosity_density(r)
    print(f"r = {r:5.1f}: η = {eff:.4f}, v_orb = {v_orb:.4f}c, L_density = {lum_dens:.4f}")

# Calculate binding energy and compare
print(f"\nBinding energy at ISCO: {bh.efficiency(bh.r_isco):.4f} mc²")
print(f"This means {bh.efficiency(bh.r_isco)*100:.1f}% of rest mass energy is extractable")

print(f"\nOptimal configuration summary:")
print(f"- Maximum efficiency occurs at r = {r_opt:.2f} GM/c²")
print(f"- This is {r_opt/bh.r_isco:.2f} times the ISCO radius")
print(f"- Energy extraction efficiency: {bh.efficiency(r_opt)*100:.1f}%")
print(f"- Orbital velocity: {bh.orbital_velocity(r_opt)*100:.1f}% of light speed")