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

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

class PlanetaryProtectionOptimizer:
    """
    Optimizer for planetary protection protocols in sample return missions.
    
    This class models the contamination risk and optimizes the protection
    strategy considering multiple sterilization methods and containment barriers.
    """
    
    def __init__(self):
        # Mission parameters
        self.P0 = 1e-3  # Initial contamination probability (0.1%)
        self.max_budget = 10e6  # Maximum budget: $10 million
        self.max_time = 72  # Maximum time: 72 hours
        
        # Cost parameters (in millions of dollars)
        self.cost_per_barrier = 1.2  # Cost per containment barrier
        self.cost_sterilization_base = 0.5  # Base sterilization cost
        self.cost_per_temp_unit = 0.01  # Cost per degree Celsius
        self.cost_per_hour = 0.05  # Operating cost per hour
        self.cost_uv_per_hour = 0.08  # UV equipment cost per hour
        
        # Effectiveness parameters
        self.barrier_effectiveness = 0.95  # Each barrier reduces risk by 95%
        self.temp_effectiveness_coef = 0.004  # Temperature effectiveness coefficient
        self.duration_effectiveness_coef = 0.15  # Duration effectiveness coefficient
        self.uv_effectiveness_coef = 0.12  # UV effectiveness coefficient
        
    def contamination_probability(self, x):
        """
        Calculate total contamination probability.
        
        Parameters:
        x: array [n_barriers, temp, duration, uv_time]
        
        Returns:
        Total contamination probability
        """
        n_barriers, temp, duration, uv_time = x
        
        # Barrier effectiveness (multiplicative)
        P_barrier = self.P0 * ((1 - self.barrier_effectiveness) ** n_barriers)
        
        # Temperature sterilization effectiveness (exponential decay)
        # Higher temperature and longer duration reduce contamination
        eta_temp = 1 - np.exp(-self.temp_effectiveness_coef * (temp - 120))
        
        # Duration effectiveness (logarithmic)
        eta_duration = 1 - np.exp(-self.duration_effectiveness_coef * duration)
        
        # UV sterilization effectiveness
        eta_uv = 1 - np.exp(-self.uv_effectiveness_coef * uv_time)
        
        # Combined probability (independent events)
        P_total = P_barrier * (1 - eta_temp) * (1 - eta_duration) * (1 - eta_uv)
        
        return P_total
    
    def total_cost(self, x):
        """Calculate total mission cost in millions of dollars."""
        n_barriers, temp, duration, uv_time = x
        
        cost = (self.cost_per_barrier * n_barriers + 
                self.cost_sterilization_base +
                self.cost_per_temp_unit * temp +
                self.cost_per_hour * duration +
                self.cost_uv_per_hour * uv_time)
        
        return cost
    
    def total_time(self, x):
        """Calculate total mission time in hours."""
        n_barriers, temp, duration, uv_time = x
        
        # Barrier installation time (5 hours per barrier)
        # Plus sterilization duration and UV exposure time
        time = 5 * n_barriers + duration + uv_time
        
        return time
    
    def objective_function(self, x):
        """
        Objective function to minimize (contamination probability).
        Includes penalty terms for constraint violations.
        """
        # Calculate contamination probability
        P = self.contamination_probability(x)
        
        # Add penalty for constraint violations
        penalty = 0
        
        # Budget constraint
        cost = self.total_cost(x)
        if cost > self.max_budget:
            penalty += 1e6 * (cost - self.max_budget)
        
        # Time constraint
        time = self.total_time(x)
        if time > self.max_time:
            penalty += 1e3 * (time - self.max_time)
        
        return P + penalty
    
    def optimize(self):
        """
        Perform optimization using differential evolution algorithm.
        
        Returns:
        Optimal solution and results dictionary
        """
        # Define bounds for each variable
        # [n_barriers, temp, duration, uv_time]
        bounds = [
            (1, 5),      # Number of barriers (integer, but treated as continuous)
            (120, 200),  # Temperature (°C)
            (0.5, 10),   # Sterilization duration (hours)
            (0, 8)       # UV exposure time (hours)
        ]
        
        # Optimize using differential evolution
        result = differential_evolution(
            self.objective_function,
            bounds,
            seed=42,
            maxiter=1000,
            popsize=15,
            atol=1e-10,
            tol=1e-10
        )
        
        # Round number of barriers to nearest integer
        optimal_solution = result.x.copy()
        optimal_solution[0] = round(optimal_solution[0])
        
        # Calculate final metrics
        results = {
            'solution': optimal_solution,
            'contamination_prob': self.contamination_probability(optimal_solution),
            'cost': self.total_cost(optimal_solution),
            'time': self.total_time(optimal_solution),
            'n_barriers': int(optimal_solution[0]),
            'temperature': optimal_solution[1],
            'duration': optimal_solution[2],
            'uv_time': optimal_solution[3]
        }
        
        return results
    
    def sensitivity_analysis(self, optimal_solution):
        """
        Perform sensitivity analysis by varying each parameter.
        
        Returns:
        Dictionary containing sensitivity data for each parameter
        """
        sensitivity = {}
        
        # Temperature sensitivity
        temps = np.linspace(120, 200, 50)
        contamination_temp = []
        for t in temps:
            x = optimal_solution.copy()
            x[1] = t
            contamination_temp.append(self.contamination_probability(x))
        sensitivity['temperature'] = (temps, contamination_temp)
        
        # Duration sensitivity
        durations = np.linspace(0.5, 10, 50)
        contamination_duration = []
        for d in durations:
            x = optimal_solution.copy()
            x[2] = d
            contamination_duration.append(self.contamination_probability(x))
        sensitivity['duration'] = (durations, contamination_duration)
        
        # UV time sensitivity
        uv_times = np.linspace(0, 8, 50)
        contamination_uv = []
        for u in uv_times:
            x = optimal_solution.copy()
            x[3] = u
            contamination_uv.append(self.contamination_probability(x))
        sensitivity['uv_time'] = (uv_times, contamination_uv)
        
        # Number of barriers sensitivity
        barriers = np.arange(1, 6)
        contamination_barriers = []
        for b in barriers:
            x = optimal_solution.copy()
            x[0] = b
            contamination_barriers.append(self.contamination_probability(x))
        sensitivity['barriers'] = (barriers, contamination_barriers)
        
        return sensitivity
    
    def trade_space_analysis(self):
        """
        Analyze the trade space between cost, time, and contamination risk.
        
        Returns:
        Arrays for visualization
        """
        n_samples = 1000
        costs = []
        times = []
        contaminations = []
        
        # Generate random samples within constraints
        for _ in range(n_samples):
            x = np.array([
                np.random.randint(1, 6),
                np.random.uniform(120, 200),
                np.random.uniform(0.5, 10),
                np.random.uniform(0, 8)
            ])
            
            costs.append(self.total_cost(x))
            times.append(self.total_time(x))
            contaminations.append(self.contamination_probability(x))
        
        return np.array(costs), np.array(times), np.array(contaminations)


# Main execution
print("=" * 70)
print("PLANETARY PROTECTION PROTOCOL OPTIMIZATION")
print("Mars Sample Return Mission")
print("=" * 70)

# Create optimizer instance
optimizer = PlanetaryProtectionOptimizer()

# Perform optimization
print("\n[1] Running optimization algorithm...")
results = optimizer.optimize()

# Display results
print("\n" + "=" * 70)
print("OPTIMAL SOLUTION")
print("=" * 70)
print(f"Number of Containment Barriers: {results['n_barriers']}")
print(f"Sterilization Temperature: {results['temperature']:.2f} °C")
print(f"Sterilization Duration: {results['duration']:.2f} hours")
print(f"UV Exposure Time: {results['uv_time']:.2f} hours")
print(f"\nContamination Probability: {results['contamination_prob']:.2e} ({results['contamination_prob']*100:.6f}%)")
print(f"Total Cost: ${results['cost']:.2f} million")
print(f"Total Time: {results['time']:.2f} hours")
print("=" * 70)

# Perform sensitivity analysis
print("\n[2] Performing sensitivity analysis...")
sensitivity = optimizer.sensitivity_analysis(results['solution'])

# Perform trade space analysis
print("\n[3] Analyzing design trade space...")
costs, times, contaminations = optimizer.trade_space_analysis()

# Visualization
print("\n[4] Generating visualizations...")

fig = plt.figure(figsize=(18, 12))

# Plot 1: Sensitivity to Temperature
ax1 = plt.subplot(2, 3, 1)
temps, cont_temp = sensitivity['temperature']
ax1.plot(temps, np.array(cont_temp) * 1e6, 'b-', linewidth=2)
ax1.axvline(results['temperature'], color='r', linestyle='--', linewidth=2, label='Optimal')
ax1.set_xlabel('Temperature (°C)', fontsize=11, fontweight='bold')
ax1.set_ylabel('Contamination Probability (×10⁻⁶)', fontsize=11, fontweight='bold')
ax1.set_title('Sensitivity to Sterilization Temperature', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Plot 2: Sensitivity to Duration
ax2 = plt.subplot(2, 3, 2)
durations, cont_dur = sensitivity['duration']
ax2.plot(durations, np.array(cont_dur) * 1e6, 'g-', linewidth=2)
ax2.axvline(results['duration'], color='r', linestyle='--', linewidth=2, label='Optimal')
ax2.set_xlabel('Sterilization Duration (hours)', fontsize=11, fontweight='bold')
ax2.set_ylabel('Contamination Probability (×10⁻⁶)', fontsize=11, fontweight='bold')
ax2.set_title('Sensitivity to Sterilization Duration', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
ax2.legend()

# Plot 3: Sensitivity to UV Exposure
ax3 = plt.subplot(2, 3, 3)
uv_times, cont_uv = sensitivity['uv_time']
ax3.plot(uv_times, np.array(cont_uv) * 1e6, 'm-', linewidth=2)
ax3.axvline(results['uv_time'], color='r', linestyle='--', linewidth=2, label='Optimal')
ax3.set_xlabel('UV Exposure Time (hours)', fontsize=11, fontweight='bold')
ax3.set_ylabel('Contamination Probability (×10⁻⁶)', fontsize=11, fontweight='bold')
ax3.set_title('Sensitivity to UV Radiation', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.legend()

# Plot 4: Sensitivity to Number of Barriers
ax4 = plt.subplot(2, 3, 4)
barriers, cont_bar = sensitivity['barriers']
ax4.bar(barriers, np.array(cont_bar) * 1e6, color='orange', alpha=0.7, edgecolor='black')
ax4.axvline(results['n_barriers'], color='r', linestyle='--', linewidth=2, label='Optimal')
ax4.set_xlabel('Number of Containment Barriers', fontsize=11, fontweight='bold')
ax4.set_ylabel('Contamination Probability (×10⁻⁶)', fontsize=11, fontweight='bold')
ax4.set_title('Sensitivity to Containment Barriers', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3, axis='y')
ax4.legend()

# Plot 5: Cost vs Contamination Risk Trade-off
ax5 = plt.subplot(2, 3, 5)
scatter = ax5.scatter(costs, contaminations * 1e6, c=times, cmap='viridis', 
                      alpha=0.5, s=30, edgecolors='black', linewidth=0.5)
ax5.scatter(results['cost'], results['contamination_prob'] * 1e6, 
           color='red', s=200, marker='*', edgecolors='black', linewidth=2, 
           label='Optimal', zorder=5)
cbar = plt.colorbar(scatter, ax=ax5)
cbar.set_label('Time (hours)', fontsize=10, fontweight='bold')
ax5.set_xlabel('Cost ($ millions)', fontsize=11, fontweight='bold')
ax5.set_ylabel('Contamination Probability (×10⁻⁶)', fontsize=11, fontweight='bold')
ax5.set_title('Cost vs Risk Trade-off', fontsize=12, fontweight='bold')
ax5.grid(True, alpha=0.3)
ax5.legend()

# Plot 6: 3D Trade Space
ax6 = plt.subplot(2, 3, 6, projection='3d')
scatter3d = ax6.scatter(costs, times, contaminations * 1e6, 
                        c=contaminations * 1e6, cmap='coolwarm', 
                        alpha=0.4, s=20, edgecolors='black', linewidth=0.3)
ax6.scatter(results['cost'], results['time'], results['contamination_prob'] * 1e6,
           color='red', s=300, marker='*', edgecolors='black', linewidth=2, 
           label='Optimal Solution')
ax6.set_xlabel('Cost ($ millions)', fontsize=10, fontweight='bold')
ax6.set_ylabel('Time (hours)', fontsize=10, fontweight='bold')
ax6.set_zlabel('Contamination (×10⁻⁶)', fontsize=10, fontweight='bold')
ax6.set_title('3D Design Trade Space', fontsize=12, fontweight='bold')
ax6.legend()

plt.tight_layout()
plt.savefig('planetary_protection_optimization.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n[5] Optimization complete!")
print("=" * 70)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import differential_evolution
import seaborn as sns

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

class PrebioticChemistry:
    """
    Simulates a prebiotic chemical reaction network with the following species:
    0: Simple organics (H2CO, HCN, etc.)
    1: Amino acid precursors
    2: Nucleotide precursors
    3: Lipid precursors
    4: Simple peptides
    5: Simple oligonucleotides
    6: Complex catalytic molecules
    """
    
    def __init__(self, temperature=300, energy_flux=1.0):
        self.n_species = 7
        self.temperature = temperature  # Kelvin
        self.energy_flux = energy_flux  # Energy availability (normalized)
        
        # Reaction rate constants (influenced by temperature and energy)
        self.k_base = np.array([
            0.1,   # Simple organic formation
            0.05,  # Amino acid precursor formation
            0.05,  # Nucleotide precursor formation
            0.03,  # Lipid precursor formation
            0.02,  # Peptide formation
            0.02,  # Oligonucleotide formation
            0.01,  # Catalytic molecule formation
            0.015, # Autocatalytic enhancement
        ])
        
        # Temperature dependence (Arrhenius-like)
        self.activation_energy = 50  # kJ/mol
        self.R = 8.314e-3  # kJ/(mol·K)
        
    def arrhenius_factor(self):
        """Temperature-dependent rate modification"""
        return np.exp(-self.activation_energy / (self.R * self.temperature))
    
    def reaction_network(self, concentrations, t):
        """
        Defines the chemical reaction network ODEs
        """
        C = concentrations
        dC = np.zeros(self.n_species)
        
        # Effective rate constants
        k = self.k_base * self.arrhenius_factor() * self.energy_flux
        
        # Reaction 1: Simple organics formation (from environment)
        dC[0] = k[0] * (1.0 - C[0]) - k[1] * C[0]
        
        # Reaction 2: Amino acid precursors from simple organics
        dC[1] = k[1] * C[0] - k[4] * C[1]**2 - 0.01 * C[1]
        
        # Reaction 3: Nucleotide precursors from simple organics
        dC[2] = k[2] * C[0] - k[5] * C[2]**2 - 0.01 * C[2]
        
        # Reaction 4: Lipid precursors from simple organics
        dC[3] = k[3] * C[0] - 0.005 * C[3]
        
        # Reaction 5: Peptide formation (polymerization)
        dC[4] = k[4] * C[1]**2 - 0.015 * C[4] + k[7] * C[6] * C[1]
        
        # Reaction 6: Oligonucleotide formation
        dC[5] = k[5] * C[2]**2 - 0.015 * C[5] + k[7] * C[6] * C[2]
        
        # Reaction 7: Catalytic molecule formation (autocatalytic)
        # This is key: complex molecules catalyze their own formation
        autocatalytic_boost = 1.0 + 2.0 * C[6]
        dC[6] = k[6] * C[4] * C[5] * autocatalytic_boost - 0.02 * C[6]
        
        return dC
    
    def simulate(self, initial_conditions, t_span):
        """Run the simulation"""
        solution = odeint(self.reaction_network, initial_conditions, t_span)
        return solution
    
    def calculate_self_organization_score(self, solution):
        """
        Calculate metrics for self-organization:
        - Autocatalytic strength: presence of catalytic molecules
        - Diversity: Shannon entropy of species distribution
        - Complexity: network connectivity and high-level molecules
        """
        final_state = solution[-1, :]
        
        # Autocatalytic strength (weight on catalytic molecules)
        A = final_state[6] * 10.0
        
        # Diversity (Shannon entropy)
        concentrations_norm = final_state / (np.sum(final_state) + 1e-10)
        concentrations_norm = concentrations_norm[concentrations_norm > 1e-10]
        D = -np.sum(concentrations_norm * np.log(concentrations_norm + 1e-10))
        
        # Complexity (weighted sum of complex molecules)
        C = final_state[4] * 2.0 + final_state[5] * 2.0 + final_state[6] * 5.0
        
        # Combined score
        alpha, beta, gamma = 1.0, 0.5, 1.0
        S = alpha * A + beta * D + gamma * C
        
        return S, A, D, C

def objective_function(initial_conditions, chemistry, t_span):
    """
    Objective function for optimization
    Returns negative score (since we minimize)
    """
    # Ensure non-negative concentrations
    initial_conditions = np.abs(initial_conditions)
    
    try:
        solution = chemistry.simulate(initial_conditions, t_span)
        score, _, _, _ = chemistry.calculate_self_organization_score(solution)
        return -score  # Negative because we minimize
    except:
        return 1e10  # Penalty for failed simulations

# Main simulation parameters
t_span = np.linspace(0, 100, 1000)  # Time in arbitrary units
chemistry = PrebioticChemistry(temperature=300, energy_flux=1.0)

# Initial guess (small random concentrations)
initial_guess = np.random.rand(chemistry.n_species) * 0.1

# Optimization bounds (0 to 2.0 for each species)
bounds = [(0, 2.0) for _ in range(chemistry.n_species)]

print("=" * 70)
print("ORIGIN OF LIFE: OPTIMIZING INITIAL CONDITIONS")
print("=" * 70)
print("\nSearching for optimal initial concentrations...")
print("This may take a minute...\n")

# Run optimization
result = differential_evolution(
    objective_function,
    bounds,
    args=(chemistry, t_span),
    maxiter=100,
    popsize=15,
    seed=42,
    atol=1e-4,
    tol=1e-4
)

optimal_initial = result.x
print("Optimization complete!")
print(f"Optimal self-organization score: {-result.fun:.4f}\n")

# Compare optimal vs random initial conditions
random_initial = np.random.rand(chemistry.n_species) * 0.5

# Simulate both scenarios
solution_optimal = chemistry.simulate(optimal_initial, t_span)
solution_random = chemistry.simulate(random_initial, t_span)

# Calculate scores over time
scores_optimal = []
scores_random = []
A_vals_opt, D_vals_opt, C_vals_opt = [], [], []
A_vals_rnd, D_vals_rnd, C_vals_rnd = [], [], []

for i in range(len(t_span)):
    s_opt, a_opt, d_opt, c_opt = chemistry.calculate_self_organization_score(solution_optimal[:i+1])
    s_rnd, a_rnd, d_rnd, c_rnd = chemistry.calculate_self_organization_score(solution_random[:i+1])
    scores_optimal.append(s_opt)
    scores_random.append(s_rnd)
    A_vals_opt.append(a_opt)
    D_vals_opt.append(d_opt)
    C_vals_opt.append(c_opt)
    A_vals_rnd.append(a_rnd)
    D_vals_rnd.append(d_rnd)
    C_vals_rnd.append(c_rnd)

# Print detailed results
print("=" * 70)
print("INITIAL CONDITIONS COMPARISON")
print("=" * 70)
species_names = ['Simple Organics', 'Amino Acid Pre.', 'Nucleotide Pre.',
                 'Lipid Pre.', 'Peptides', 'Oligonucleotides', 'Catalytic Mol.']

print("\nOptimal Initial Conditions:")
for i, name in enumerate(species_names):
    print(f"  {name:20s}: {optimal_initial[i]:.4f}")

print("\nRandom Initial Conditions:")
for i, name in enumerate(species_names):
    print(f"  {name:20s}: {random_initial[i]:.4f}")

print("\n" + "=" * 70)
print("FINAL STATE COMPARISON (t=100)")
print("=" * 70)
print("\nOptimal System:")
for i, name in enumerate(species_names):
    print(f"  {name:20s}: {solution_optimal[-1, i]:.4f}")

print("\nRandom System:")
for i, name in enumerate(species_names):
    print(f"  {name:20s}: {solution_random[-1, i]:.4f}")

# Final scores
score_opt_final, A_opt, D_opt, C_opt = chemistry.calculate_self_organization_score(solution_optimal)
score_rnd_final, A_rnd, D_rnd, C_rnd = chemistry.calculate_self_organization_score(solution_random)

print("\n" + "=" * 70)
print("SELF-ORGANIZATION METRICS")
print("=" * 70)
print(f"\nOptimal System:")
print(f"  Autocatalytic Strength (A): {A_opt:.4f}")
print(f"  Diversity (D):              {D_opt:.4f}")
print(f"  Complexity (C):             {C_opt:.4f}")
print(f"  Total Score (S):            {score_opt_final:.4f}")

print(f"\nRandom System:")
print(f"  Autocatalytic Strength (A): {A_rnd:.4f}")
print(f"  Diversity (D):              {D_rnd:.4f}")
print(f"  Complexity (C):             {C_rnd:.4f}")
print(f"  Total Score (S):            {score_rnd_final:.4f}")

print(f"\nImprovement: {(score_opt_final/score_rnd_final - 1)*100:.1f}%")

# ===== VISUALIZATION =====
plt.style.use('seaborn-v0_8-darkgrid')
fig = plt.figure(figsize=(16, 12))

# Plot 1: Concentration evolution (Optimal)
ax1 = plt.subplot(3, 3, 1)
for i in range(chemistry.n_species):
    ax1.plot(t_span, solution_optimal[:, i], label=species_names[i], linewidth=2)
ax1.set_xlabel('Time (arbitrary units)', fontsize=10)
ax1.set_ylabel('Concentration', fontsize=10)
ax1.set_title('Optimal System: Concentration Evolution', fontsize=12, fontweight='bold')
ax1.legend(fontsize=7, loc='best')
ax1.grid(True, alpha=0.3)

# Plot 2: Concentration evolution (Random)
ax2 = plt.subplot(3, 3, 2)
for i in range(chemistry.n_species):
    ax2.plot(t_span, solution_random[:, i], label=species_names[i], linewidth=2)
ax2.set_xlabel('Time (arbitrary units)', fontsize=10)
ax2.set_ylabel('Concentration', fontsize=10)
ax2.set_title('Random System: Concentration Evolution', fontsize=12, fontweight='bold')
ax2.legend(fontsize=7, loc='best')
ax2.grid(True, alpha=0.3)

# Plot 3: Self-organization score comparison
ax3 = plt.subplot(3, 3, 3)
ax3.plot(t_span, scores_optimal, label='Optimal Initial', linewidth=3, color='green')
ax3.plot(t_span, scores_random, label='Random Initial', linewidth=3, color='red', linestyle='--')
ax3.set_xlabel('Time (arbitrary units)', fontsize=10)
ax3.set_ylabel('Self-Organization Score', fontsize=10)
ax3.set_title('Self-Organization Score Comparison', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# Plot 4: Catalytic molecules (key for life!)
ax4 = plt.subplot(3, 3, 4)
ax4.plot(t_span, solution_optimal[:, 6], label='Optimal', linewidth=3, color='purple')
ax4.plot(t_span, solution_random[:, 6], label='Random', linewidth=3, color='orange', linestyle='--')
ax4.set_xlabel('Time (arbitrary units)', fontsize=10)
ax4.set_ylabel('Concentration', fontsize=10)
ax4.set_title('Catalytic Molecules: Key to Self-Organization', fontsize=12, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

# Plot 5: Initial conditions bar chart
ax5 = plt.subplot(3, 3, 5)
x = np.arange(len(species_names))
width = 0.35
ax5.bar(x - width/2, optimal_initial, width, label='Optimal', color='green', alpha=0.7)
ax5.bar(x + width/2, random_initial, width, label='Random', color='red', alpha=0.7)
ax5.set_xlabel('Species', fontsize=10)
ax5.set_ylabel('Initial Concentration', fontsize=10)
ax5.set_title('Initial Conditions Comparison', fontsize=12, fontweight='bold')
ax5.set_xticks(x)
ax5.set_xticklabels([s[:10] for s in species_names], rotation=45, ha='right', fontsize=8)
ax5.legend(fontsize=10)
ax5.grid(True, alpha=0.3, axis='y')

# Plot 6: Final state bar chart
ax6 = plt.subplot(3, 3, 6)
ax6.bar(x - width/2, solution_optimal[-1, :], width, label='Optimal', color='green', alpha=0.7)
ax6.bar(x + width/2, solution_random[-1, :], width, label='Random', color='red', alpha=0.7)
ax6.set_xlabel('Species', fontsize=10)
ax6.set_ylabel('Final Concentration', fontsize=10)
ax6.set_title('Final State Comparison (t=100)', fontsize=12, fontweight='bold')
ax6.set_xticks(x)
ax6.set_xticklabels([s[:10] for s in species_names], rotation=45, ha='right', fontsize=8)
ax6.legend(fontsize=10)
ax6.grid(True, alpha=0.3, axis='y')

# Plot 7: Metric decomposition (Optimal)
ax7 = plt.subplot(3, 3, 7)
ax7.plot(t_span, A_vals_opt, label='Autocatalytic (A)', linewidth=2, color='blue')
ax7.plot(t_span, D_vals_opt, label='Diversity (D)', linewidth=2, color='orange')
ax7.plot(t_span, C_vals_opt, label='Complexity (C)', linewidth=2, color='green')
ax7.set_xlabel('Time (arbitrary units)', fontsize=10)
ax7.set_ylabel('Metric Value', fontsize=10)
ax7.set_title('Optimal System: Metric Decomposition', fontsize=12, fontweight='bold')
ax7.legend(fontsize=10)
ax7.grid(True, alpha=0.3)

# Plot 8: Metric decomposition (Random)
ax8 = plt.subplot(3, 3, 8)
ax8.plot(t_span, A_vals_rnd, label='Autocatalytic (A)', linewidth=2, color='blue', linestyle='--')
ax8.plot(t_span, D_vals_rnd, label='Diversity (D)', linewidth=2, color='orange', linestyle='--')
ax8.plot(t_span, C_vals_rnd, label='Complexity (C)', linewidth=2, color='green', linestyle='--')
ax8.set_xlabel('Time (arbitrary units)', fontsize=10)
ax8.set_ylabel('Metric Value', fontsize=10)
ax8.set_title('Random System: Metric Decomposition', fontsize=12, fontweight='bold')
ax8.legend(fontsize=10)
ax8.grid(True, alpha=0.3)

# Plot 9: Phase space (Complex molecules vs Catalytic)
ax9 = plt.subplot(3, 3, 9)
complex_opt = solution_optimal[:, 4] + solution_optimal[:, 5]
complex_rnd = solution_random[:, 4] + solution_random[:, 5]
ax9.plot(complex_opt, solution_optimal[:, 6], linewidth=2, color='green', label='Optimal')
ax9.plot(complex_rnd, solution_random[:, 6], linewidth=2, color='red', linestyle='--', label='Random')
ax9.scatter([complex_opt[0]], [solution_optimal[0, 6]], s=100, c='green', marker='o', label='Opt Start', zorder=5)
ax9.scatter([complex_opt[-1]], [solution_optimal[-1, 6]], s=100, c='darkgreen', marker='*', label='Opt End', zorder=5)
ax9.scatter([complex_rnd[0]], [solution_random[0, 6]], s=100, c='red', marker='o', label='Rnd Start', zorder=5)
ax9.scatter([complex_rnd[-1]], [solution_random[-1, 6]], s=100, c='darkred', marker='*', label='Rnd End', zorder=5)
ax9.set_xlabel('Total Complex Molecules (Peptides + Oligonucleotides)', fontsize=10)
ax9.set_ylabel('Catalytic Molecules', fontsize=10)
ax9.set_title('Phase Space: Path to Self-Organization', fontsize=12, fontweight='bold')
ax9.legend(fontsize=8)
ax9.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('origin_of_life_optimization.png', dpi=150, bbox_inches='tight')
print("\n" + "=" * 70)
print("Figure saved as 'origin_of_life_optimization.png'")
print("=" * 70)
plt.show()

print("\n" + "=" * 70)
print("ANALYSIS COMPLETE")
print("=" * 70)

autocatalytic_boost = 1.0 + 2.0 * C[6]
dC[6] = k[6] * C[4] * C[5] * autocatalytic_boost - 0.02 * C[6]

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import seaborn as sns

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

class EvolutionSimulator:
    """
    Simulates evolution with adjustable selection pressure.
    
    Parameters:
    -----------
    pop_size : int
        Population size
    genome_length : int
        Number of genes per organism
    mutation_rate : float
        Probability of mutation per gene
    mutation_sigma : float
        Standard deviation of mutation effect
    selection_sigma : float
        Controls selection pressure strength
    optimal_trait : float
        Environmental optimum for trait value
    """
    
    def __init__(self, pop_size=200, genome_length=10, 
                 mutation_rate=0.1, mutation_sigma=0.5,
                 selection_sigma=2.0, optimal_trait=5.0):
        self.pop_size = pop_size
        self.genome_length = genome_length
        self.mutation_rate = mutation_rate
        self.mutation_sigma = mutation_sigma
        self.selection_sigma = selection_sigma
        self.optimal_trait = optimal_trait
        
        # Initialize population with random genomes
        self.population = np.random.randn(pop_size, genome_length)
        
        # History tracking
        self.mean_trait_history = []
        self.diversity_history = []
        self.fitness_history = []
        
    def calculate_trait(self, genome):
        """Calculate phenotypic trait as mean of genome values"""
        return np.mean(genome)
    
    def fitness_function(self, trait_value):
        """
        Gaussian fitness function centered at optimal_trait.
        
        f(x) = exp(-(x - θ)² / (2σ²))
        """
        return np.exp(-((trait_value - self.optimal_trait) ** 2) / 
                      (2 * self.selection_sigma ** 2))
    
    def calculate_fitness(self):
        """Calculate fitness for all organisms"""
        traits = np.array([self.calculate_trait(genome) 
                          for genome in self.population])
        fitness = np.array([self.fitness_function(trait) 
                           for trait in traits])
        return fitness, traits
    
    def selection(self, fitness):
        """
        Fitness-proportionate selection (roulette wheel selection).
        Returns indices of selected parents.
        """
        # Normalize fitness to probabilities
        total_fitness = np.sum(fitness)
        if total_fitness == 0:
            probabilities = np.ones(len(fitness)) / len(fitness)
        else:
            probabilities = fitness / total_fitness
        
        # Select parents
        selected_indices = np.random.choice(
            len(self.population), 
            size=self.pop_size, 
            replace=True, 
            p=probabilities
        )
        return selected_indices
    
    def mutate(self, genome):
        """Apply mutations to genome"""
        mutation_mask = np.random.random(self.genome_length) < self.mutation_rate
        mutations = np.random.randn(self.genome_length) * self.mutation_sigma
        genome = genome + mutation_mask * mutations
        return genome
    
    def reproduce(self, parent_indices):
        """Create new generation through reproduction and mutation"""
        new_population = []
        for idx in parent_indices:
            # Copy parent genome
            child_genome = self.population[idx].copy()
            # Apply mutations
            child_genome = self.mutate(child_genome)
            new_population.append(child_genome)
        
        self.population = np.array(new_population)
    
    def step(self):
        """Execute one generation of evolution"""
        # Calculate fitness
        fitness, traits = self.calculate_fitness()
        
        # Record statistics
        self.mean_trait_history.append(np.mean(traits))
        self.diversity_history.append(np.var(traits))
        self.fitness_history.append(np.mean(fitness))
        
        # Selection and reproduction
        parent_indices = self.selection(fitness)
        self.reproduce(parent_indices)
    
    def run(self, generations):
        """Run simulation for specified generations"""
        for _ in range(generations):
            self.step()

# Experiment: Compare different selection pressures
def run_experiments():
    """
    Run simulations with different selection pressure values.
    Compare genetic diversity and stability.
    """
    
    # Selection pressure values to test
    # Smaller sigma = stronger selection pressure
    selection_sigmas = [0.5, 1.0, 2.0, 4.0, 8.0]
    generations = 150
    
    results = {}
    
    for sigma in selection_sigmas:
        print(f"Running simulation with selection_sigma = {sigma}")
        sim = EvolutionSimulator(
            pop_size=200,
            genome_length=10,
            mutation_rate=0.1,
            mutation_sigma=0.5,
            selection_sigma=sigma,
            optimal_trait=5.0
        )
        sim.run(generations)
        
        results[sigma] = {
            'mean_trait': np.array(sim.mean_trait_history),
            'diversity': np.array(sim.diversity_history),
            'fitness': np.array(sim.fitness_history),
            'final_diversity': sim.diversity_history[-1],
            'final_fitness': sim.fitness_history[-1],
            'convergence_gen': np.argmax(np.array(sim.fitness_history) > 0.9 * np.max(sim.fitness_history)) 
                              if np.max(sim.fitness_history) > 0.1 else generations
        }
    
    return results, generations, selection_sigmas

# Run experiments
results, generations, selection_sigmas = run_experiments()

# Visualization
fig = plt.figure(figsize=(16, 12))
gs = GridSpec(3, 3, figure=fig, hspace=0.3, wspace=0.3)

# Color palette
colors = sns.color_palette("husl", len(selection_sigmas))

# Plot 1: Mean Trait Evolution
ax1 = fig.add_subplot(gs[0, :2])
for i, sigma in enumerate(selection_sigmas):
    ax1.plot(results[sigma]['mean_trait'], 
             label=f'σ = {sigma}', color=colors[i], linewidth=2)
ax1.axhline(y=5.0, color='red', linestyle='--', linewidth=2, label='Optimal Trait')
ax1.set_xlabel('Generation', fontsize=12)
ax1.set_ylabel('Mean Trait Value', fontsize=12)
ax1.set_title('Trait Evolution Under Different Selection Pressures', fontsize=14, fontweight='bold')
ax1.legend(loc='best')
ax1.grid(True, alpha=0.3)

# Plot 2: Genetic Diversity Over Time
ax2 = fig.add_subplot(gs[1, :2])
for i, sigma in enumerate(selection_sigmas):
    ax2.plot(results[sigma]['diversity'], 
             label=f'σ = {sigma}', color=colors[i], linewidth=2)
ax2.set_xlabel('Generation', fontsize=12)
ax2.set_ylabel('Trait Variance (Diversity)', fontsize=12)
ax2.set_title('Genetic Diversity Dynamics', fontsize=14, fontweight='bold')
ax2.legend(loc='best')
ax2.grid(True, alpha=0.3)

# Plot 3: Mean Fitness Over Time
ax3 = fig.add_subplot(gs[2, :2])
for i, sigma in enumerate(selection_sigmas):
    ax3.plot(results[sigma]['fitness'], 
             label=f'σ = {sigma}', color=colors[i], linewidth=2)
ax3.set_xlabel('Generation', fontsize=12)
ax3.set_ylabel('Mean Fitness', fontsize=12)
ax3.set_title('Population Fitness Evolution', fontsize=14, fontweight='bold')
ax3.legend(loc='best')
ax3.grid(True, alpha=0.3)

# Plot 4: Final Diversity vs Selection Pressure
ax4 = fig.add_subplot(gs[0, 2])
final_diversities = [results[sigma]['final_diversity'] for sigma in selection_sigmas]
ax4.bar(range(len(selection_sigmas)), final_diversities, color=colors, alpha=0.7)
ax4.set_xticks(range(len(selection_sigmas)))
ax4.set_xticklabels([f'{s}' for s in selection_sigmas])
ax4.set_xlabel('Selection Sigma (σ)', fontsize=11)
ax4.set_ylabel('Final Diversity', fontsize=11)
ax4.set_title('Final Genetic Diversity', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3, axis='y')

# Plot 5: Final Fitness vs Selection Pressure
ax5 = fig.add_subplot(gs[1, 2])
final_fitness = [results[sigma]['final_fitness'] for sigma in selection_sigmas]
ax5.bar(range(len(selection_sigmas)), final_fitness, color=colors, alpha=0.7)
ax5.set_xticks(range(len(selection_sigmas)))
ax5.set_xticklabels([f'{s}' for s in selection_sigmas])
ax5.set_xlabel('Selection Sigma (σ)', fontsize=11)
ax5.set_ylabel('Final Mean Fitness', fontsize=11)
ax5.set_title('Final Population Fitness', fontsize=12, fontweight='bold')
ax5.grid(True, alpha=0.3, axis='y')

# Plot 6: Convergence Speed
ax6 = fig.add_subplot(gs[2, 2])
convergence_gens = [results[sigma]['convergence_gen'] for sigma in selection_sigmas]
ax6.bar(range(len(selection_sigmas)), convergence_gens, color=colors, alpha=0.7)
ax6.set_xticks(range(len(selection_sigmas)))
ax6.set_xticklabels([f'{s}' for s in selection_sigmas])
ax6.set_xlabel('Selection Sigma (σ)', fontsize=11)
ax6.set_ylabel('Generations to 90% Max Fitness', fontsize=11)
ax6.set_title('Convergence Speed', fontsize=12, fontweight='bold')
ax6.grid(True, alpha=0.3, axis='y')

plt.suptitle('Evolution Simulation: Selection Pressure Optimization', 
             fontsize=16, fontweight='bold', y=0.995)

plt.tight_layout()
plt.show()

# Summary Statistics
print("\n" + "="*80)
print("SUMMARY STATISTICS")
print("="*80)
print(f"{'Selection σ':<15} {'Final Diversity':<20} {'Final Fitness':<20} {'Convergence Gen':<20}")
print("-"*80)
for sigma in selection_sigmas:
    print(f"{sigma:<15.1f} {results[sigma]['final_diversity']:<20.4f} "
          f"{results[sigma]['final_fitness']:<20.4f} {results[sigma]['convergence_gen']:<20}")
print("="*80)

# Calculate and display optimal parameter
diversity_fitness_product = [(sigma, 
                              results[sigma]['final_diversity'] * results[sigma]['final_fitness']) 
                             for sigma in selection_sigmas]
optimal_sigma = max(diversity_fitness_product, key=lambda x: x[1])[0]
print(f"\nOptimal Selection Pressure (maximizing diversity × fitness): σ = {optimal_sigma}")
print("="*80)

def fitness_function(self, trait_value):
    return np.exp(-((trait_value - self.optimal_trait) ** 2) / 
                  (2 * self.selection_sigma ** 2))

def mutate(self, genome):
    mutation_mask = np.random.random(self.genome_length) < self.mutation_rate
    mutations = np.random.randn(self.genome_length) * self.mutation_sigma
    genome = genome + mutation_mask * mutations
    return genome

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

# Physical constants
SIGMA = 5.67e-8  # Stefan-Boltzmann constant (W m^-2 K^-4)
L_SUN = 3.828e26  # Solar luminosity (W)
AU = 1.496e11  # Astronomical Unit (m)

class HabitableZoneCalculator:
    """
    Calculate habitable zone boundaries considering:
    - Stellar radiation characteristics
    - Planetary albedo
    - Greenhouse effect
    - Atmospheric composition
    """
    
    def __init__(self, stellar_luminosity, stellar_temp, star_name="Star"):
        """
        Initialize with stellar properties
        
        Parameters:
        -----------
        stellar_luminosity : float
            Luminosity in solar luminosities
        stellar_temp : float
            Effective temperature in Kelvin
        star_name : str
            Name of the star for labeling
        """
        self.L_star = stellar_luminosity * L_SUN  # Convert to Watts
        self.T_star = stellar_temp
        self.name = star_name
        
    def calculate_distance(self, T_planet, albedo, greenhouse_factor=1.0):
        """
        Calculate orbital distance for a given planet temperature
        
        The energy balance equation:
        T_eff = [(L_star * (1-A) * greenhouse_factor) / (16 * π * σ * d²)]^(1/4)
        
        Solving for d:
        d = sqrt[(L_star * (1-A) * greenhouse_factor) / (16 * π * σ * T_eff^4)]
        
        Parameters:
        -----------
        T_planet : float
            Target planet effective temperature (K)
        albedo : float
            Planetary albedo (0-1)
        greenhouse_factor : float
            Multiplicative factor for greenhouse warming (default=1.0)
            
        Returns:
        --------
        float : Distance in AU
        """
        numerator = self.L_star * (1 - albedo) * greenhouse_factor
        denominator = 16 * np.pi * SIGMA * (T_planet ** 4)
        distance_m = np.sqrt(numerator / denominator)
        return distance_m / AU  # Convert to AU
    
    def calculate_hz_boundaries(self, scenarios):
        """
        Calculate HZ boundaries for different scenarios
        
        Parameters:
        -----------
        scenarios : dict
            Dictionary with scenario names as keys and parameters as values
            Each scenario should have: T_inner, T_outer, albedo, greenhouse
            
        Returns:
        --------
        dict : Boundaries for each scenario
        """
        results = {}
        
        for scenario_name, params in scenarios.items():
            inner = self.calculate_distance(
                params['T_inner'], 
                params['albedo'],
                params.get('greenhouse_inner', 1.0)
            )
            outer = self.calculate_distance(
                params['T_outer'],
                params['albedo'],
                params.get('greenhouse_outer', 1.0)
            )
            
            results[scenario_name] = {
                'inner': inner,
                'outer': outer,
                'width': outer - inner,
                'params': params
            }
            
        return results

# Define different scenarios for HZ calculation
# These represent different assumptions about planetary characteristics

scenarios = {
    'Conservative': {
        'T_inner': 373,  # Water boiling point (runaway greenhouse)
        'T_outer': 273,  # Water freezing point (no greenhouse)
        'albedo': 0.3,   # Earth-like albedo
        'greenhouse_inner': 1.0,
        'greenhouse_outer': 1.0,
        'color': 'green',
        'description': 'Basic liquid water range'
    },
    'Optimistic': {
        'T_inner': 373,
        'T_outer': 200,  # With strong CO2 greenhouse
        'albedo': 0.2,   # Lower albedo (darker surface)
        'greenhouse_inner': 1.0,
        'greenhouse_outer': 2.5,  # Strong greenhouse at outer edge
        'color': 'blue',
        'description': 'With greenhouse warming'
    },
    'Recent Venus': {
        'T_inner': 335,  # Recent Venus threshold (Kopparapu et al. 2013)
        'T_outer': 273,
        'albedo': 0.3,
        'greenhouse_inner': 1.1,
        'greenhouse_outer': 1.0,
        'color': 'orange',
        'description': 'Moist greenhouse limit'
    },
    'Maximum Greenhouse': {
        'T_inner': 373,
        'T_outer': 188,  # Maximum CO2 greenhouse (early Mars)
        'albedo': 0.25,
        'greenhouse_inner': 1.0,
        'greenhouse_outer': 3.0,  # Very strong greenhouse
        'color': 'cyan',
        'description': 'Maximum CO2 warming'
    }
}

# Define different stellar types to analyze
stellar_types = {
    'M-dwarf (Cool)': {'luminosity': 0.04, 'temp': 3200, 'color': 'red'},
    'K-dwarf': {'luminosity': 0.3, 'temp': 4500, 'color': 'orange'},
    'Sun (G-dwarf)': {'luminosity': 1.0, 'temp': 5778, 'color': 'yellow'},
    'F-dwarf (Hot)': {'luminosity': 2.5, 'temp': 6500, 'color': 'lightblue'},
}

# Calculate HZ for each stellar type
all_results = {}

for star_name, star_props in stellar_types.items():
    hz_calc = HabitableZoneCalculator(
        star_props['luminosity'],
        star_props['temp'],
        star_name
    )
    all_results[star_name] = hz_calc.calculate_hz_boundaries(scenarios)

# Print detailed results
print("=" * 80)
print("HABITABLE ZONE ANALYSIS: Optimal Boundary Definitions")
print("=" * 80)
print()

for star_name, star_results in all_results.items():
    print(f"\n{'=' * 80}")
    print(f"STAR TYPE: {star_name}")
    print(f"Luminosity: {stellar_types[star_name]['luminosity']:.2f} L_sun")
    print(f"Temperature: {stellar_types[star_name]['temp']} K")
    print(f"{'=' * 80}\n")
    
    for scenario, bounds in star_results.items():
        print(f"  {scenario} Scenario:")
        print(f"    Description: {scenarios[scenario]['description']}")
        print(f"    Inner boundary: {bounds['inner']:.4f} AU")
        print(f"    Outer boundary: {bounds['outer']:.4f} AU")
        print(f"    Zone width: {bounds['width']:.4f} AU")
        print(f"    Parameters:")
        print(f"      - Inner temp: {bounds['params']['T_inner']} K")
        print(f"      - Outer temp: {bounds['params']['T_outer']} K")
        print(f"      - Albedo: {bounds['params']['albedo']}")
        print()

# Visualization 1: HZ boundaries for different stellar types
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
axes = axes.flatten()

for idx, (star_name, star_results) in enumerate(all_results.items()):
    ax = axes[idx]
    
    y_pos = 0
    for scenario, bounds in star_results.items():
        # Draw HZ rectangle
        width = bounds['outer'] - bounds['inner']
        rect = Rectangle(
            (bounds['inner'], y_pos),
            width, 0.8,
            facecolor=scenarios[scenario]['color'],
            alpha=0.5,
            edgecolor='black',
            linewidth=2
        )
        ax.add_patch(rect)
        
        # Add labels
        mid_point = bounds['inner'] + width / 2
        ax.text(mid_point, y_pos + 0.4, scenario,
                ha='center', va='center', fontsize=9, fontweight='bold')
        ax.text(bounds['inner'], y_pos - 0.15, f"{bounds['inner']:.3f}",
                ha='center', fontsize=8)
        ax.text(bounds['outer'], y_pos - 0.15, f"{bounds['outer']:.3f}",
                ha='center', fontsize=8)
        
        y_pos += 1
    
    # Mark Earth's position for Sun
    if star_name == 'Sun (G-dwarf)':
        ax.axvline(x=1.0, color='green', linestyle='--', linewidth=2, alpha=0.7)
        ax.text(1.0, y_pos - 0.5, 'Earth (1.0 AU)', rotation=90,
                va='bottom', ha='right', fontsize=10, color='green', fontweight='bold')
    
    ax.set_xlim(0, max([b['outer'] for b in star_results.values()]) * 1.1)
    ax.set_ylim(-0.5, y_pos)
    ax.set_xlabel('Distance from Star (AU)', fontsize=11, fontweight='bold')
    ax.set_ylabel('Scenario', fontsize=11, fontweight='bold')
    ax.set_title(f'{star_name}\nL = {stellar_types[star_name]["luminosity"]} L☉, T = {stellar_types[star_name]["temp"]} K',
                 fontsize=12, fontweight='bold')
    ax.grid(True, alpha=0.3)
    ax.set_yticks([])

plt.tight_layout()
plt.savefig('hz_boundaries_by_star.png', dpi=300, bbox_inches='tight')
plt.show()

# Visualization 2: Comparison of HZ widths
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

star_names = list(all_results.keys())
scenario_names = list(scenarios.keys())
x = np.arange(len(star_names))
width = 0.2

for idx, scenario in enumerate(scenario_names):
    inner_bounds = [all_results[star][scenario]['inner'] for star in star_names]
    outer_bounds = [all_results[star][scenario]['outer'] for star in star_names]
    zone_widths = [all_results[star][scenario]['width'] for star in star_names]
    
    # Plot inner and outer boundaries
    ax1.bar(x + idx*width - 1.5*width, inner_bounds, width,
            label=f'{scenario} (inner)', alpha=0.7,
            color=scenarios[scenario]['color'], edgecolor='black')
    ax1.bar(x + idx*width - 1.5*width, outer_bounds, width,
            label=f'{scenario} (outer)', alpha=0.4,
            color=scenarios[scenario]['color'], edgecolor='black')
    
    # Plot zone widths
    ax2.bar(x + idx*width - 1.5*width, zone_widths, width,
            label=scenario, alpha=0.7,
            color=scenarios[scenario]['color'], edgecolor='black')

ax1.set_xlabel('Stellar Type', fontsize=12, fontweight='bold')
ax1.set_ylabel('Distance (AU)', fontsize=12, fontweight='bold')
ax1.set_title('Habitable Zone Boundaries by Star Type', fontsize=14, fontweight='bold')
ax1.set_xticks(x)
ax1.set_xticklabels(star_names, rotation=15, ha='right')
ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=8)
ax1.grid(True, alpha=0.3)

ax2.set_xlabel('Stellar Type', fontsize=12, fontweight='bold')
ax2.set_ylabel('Zone Width (AU)', fontsize=12, fontweight='bold')
ax2.set_title('Habitable Zone Width Comparison', fontsize=14, fontweight='bold')
ax2.set_xticks(x)
ax2.set_xticklabels(star_names, rotation=15, ha='right')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('hz_comparison.png', dpi=300, bbox_inches='tight')
plt.show()

# Visualization 3: Temperature-Distance relationship
fig, ax = plt.subplots(figsize=(14, 8))

distances = np.linspace(0.1, 3.0, 1000)

for star_name, star_props in stellar_types.items():
    hz_calc = HabitableZoneCalculator(
        star_props['luminosity'],
        star_props['temp'],
        star_name
    )
    
    # Calculate effective temperature at each distance (with Earth-like albedo)
    temps = []
    for d in distances:
        # Reverse calculation: given distance, find temperature
        T = (hz_calc.L_star * (1 - 0.3) / (16 * np.pi * SIGMA * (d * AU)**2))**(1/4)
        temps.append(T)
    
    ax.plot(distances, temps, linewidth=2.5, label=star_name,
            color=star_props['color'], alpha=0.8)

# Mark key temperature boundaries
ax.axhline(y=373, color='red', linestyle='--', linewidth=2, alpha=0.7, label='Water boiling (373 K)')
ax.axhline(y=273, color='blue', linestyle='--', linewidth=2, alpha=0.7, label='Water freezing (273 K)')
ax.axhline(y=288, color='green', linestyle=':', linewidth=2, alpha=0.7, label='Earth average (288 K)')

# Shade habitable temperature range
ax.axhspan(273, 373, alpha=0.1, color='green', label='Liquid water range')

ax.set_xlabel('Distance from Star (AU)', fontsize=12, fontweight='bold')
ax.set_ylabel('Effective Planet Temperature (K)', fontsize=12, fontweight='bold')
ax.set_title('Planet Temperature vs Orbital Distance\n(Albedo = 0.3, No Greenhouse Effect)',
             fontsize=14, fontweight='bold')
ax.set_xlim(0.1, 3.0)
ax.set_ylim(150, 500)
ax.legend(loc='upper right', fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('temp_distance_relationship.png', dpi=300, bbox_inches='tight')
plt.show()

# Visualization 4: Effect of albedo on HZ boundaries
fig, ax = plt.subplots(figsize=(14, 8))

albedos = np.linspace(0.0, 0.7, 8)
colors = plt.cm.viridis(np.linspace(0, 1, len(albedos)))

# Use Sun as reference
hz_calc = HabitableZoneCalculator(1.0, 5778, "Sun")

for idx, albedo in enumerate(albedos):
    inner_bounds = []
    outer_bounds = []
    
    for scenario in ['Conservative', 'Optimistic']:
        params = scenarios[scenario].copy()
        params['albedo'] = albedo
        
        inner = hz_calc.calculate_distance(
            params['T_inner'], albedo, params.get('greenhouse_inner', 1.0)
        )
        outer = hz_calc.calculate_distance(
            params['T_outer'], albedo, params.get('greenhouse_outer', 1.0)
        )
        
        inner_bounds.append(inner)
        outer_bounds.append(outer)
    
    # Plot as error bar style
    x_pos = idx
    ax.plot([x_pos, x_pos], [inner_bounds[0], outer_bounds[0]],
            linewidth=8, color=colors[idx], alpha=0.6, label=f'A = {albedo:.2f}')
    ax.plot([x_pos, x_pos], [inner_bounds[1], outer_bounds[1]],
            linewidth=8, color=colors[idx], alpha=0.3, linestyle='--')

ax.axhline(y=1.0, color='green', linestyle=':', linewidth=2, alpha=0.5, label='Earth (1.0 AU)')
ax.set_xlabel('Albedo Value (0 = dark, 1 = reflective)', fontsize=12, fontweight='bold')
ax.set_ylabel('Distance from Sun (AU)', fontsize=12, fontweight='bold')
ax.set_title('Effect of Planetary Albedo on Habitable Zone\nSolid: Conservative | Dashed: Optimistic',
             fontsize=14, fontweight='bold')
ax.set_xticks(range(len(albedos)))
ax.set_xticklabels([f'{a:.2f}' for a in albedos])
ax.legend(loc='upper left', fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('albedo_effect.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "=" * 80)
print("ANALYSIS COMPLETE")
print("=" * 80)
print("\nKey Findings:")
print("1. Cooler stars (M-dwarfs) have HZ much closer to the star")
print("2. Hotter stars (F-dwarfs) have HZ farther from the star")
print("3. Greenhouse effect significantly expands the outer boundary")
print("4. Albedo strongly affects HZ position (darker = farther out)")
print("5. HZ width varies significantly with stellar type")
print("\nGraphs saved: hz_boundaries_by_star.png, hz_comparison.png,")
print("              temp_distance_relationship.png, albedo_effect.png")
print("=" * 80)