2026_mcm_b/p4/environmental_impact.py

"""
Task 4: Environmental Impact Analysis for Moon Colony Construction
使用液氧甲烷(LOX/CH4)燃料的环境影响评估

分析不同建设方案对地球环境的影响，并给出最小化环境影响的模型调整策略。
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib import rcParams
from dataclasses import dataclass
from typing import List, Dict, Tuple
import pandas as pd

# 设置字体
rcParams['font.sans-serif'] = ['Arial Unicode MS', 'DejaVu Sans', 'SimHei']
rcParams['axes.unicode_minus'] = False

# ============== 物理常数 ==============
G0 = 9.81  # m/s²
OMEGA_EARTH = 7.27e-5  # rad/s
R_EARTH = 6.371e6  # m

# ============== 任务参数 (与任务一一致) ==============
TOTAL_PAYLOAD = 100e6  # 100 million metric tons

# ============== 太空电梯参数 ==============
NUM_ELEVATORS = 3
ELEVATOR_CAPACITY_PER_YEAR = 179000  # metric tons per elevator per year
TOTAL_ELEVATOR_CAPACITY = NUM_ELEVATORS * ELEVATOR_CAPACITY_PER_YEAR  # 537,000 tons/year
ELEVATOR_SPECIFIC_ENERGY = 157.2e9  # J per metric ton (157.2 MJ/kg)

# 电梯能量分解
ELEVATOR_LIFT_ENERGY = 57.0e9  # J/ton (电力提升能)
ELEVATOR_TOP_ROCKET_ENERGY = 100.2e9  # J/ton (顶端火箭推进)

# ============== 液氧甲烷燃料参数 (LOX/CH4) ==============
@dataclass
class MethaneFuelParams:
    """液氧甲烷燃料参数"""
    name: str = "LOX/CH4"
    isp: float = 360  # 比冲 (秒), Raptor发动机级别
    specific_energy: float = 12.9e6  # 燃料比能量 (J/kg)
    ox_fuel_ratio: float = 3.5  # 氧燃比 O2:CH4
    
    @property
    def exhaust_velocity(self) -> float:
        """排气速度 (m/s)"""
        return self.isp * G0
    
    @property
    def ch4_mass_fraction(self) -> float:
        """甲烷在燃料中的质量分数"""
        return 1 / (1 + self.ox_fuel_ratio)  # ≈ 0.222
    
    @property
    def lox_mass_fraction(self) -> float:
        """液氧在燃料中的质量分数"""
        return self.ox_fuel_ratio / (1 + self.ox_fuel_ratio)  # ≈ 0.778


# ============== 排放因子 (基于化学计量) ==============
@dataclass
class EmissionFactors:
    """排放因子 (kg排放物/kg燃料)"""
    # 燃烧排放: CH4 + 2O2 → CO2 + 2H2O
    # 每kg CH4产生: 44/16 = 2.75 kg CO2, 36/16 = 2.25 kg H2O
    
    fuel_params: MethaneFuelParams = None
    
    def __post_init__(self):
        if self.fuel_params is None:
            self.fuel_params = MethaneFuelParams()
    
    @property
    def co2_per_kg_fuel(self) -> float:
        """每kg燃料的CO2排放 (kg)"""
        return self.fuel_params.ch4_mass_fraction * 2.75  # ≈ 0.611 kg CO2/kg燃料
    
    @property
    def h2o_per_kg_fuel(self) -> float:
        """每kg燃料的H2O排放 (kg)"""
        return self.fuel_params.ch4_mass_fraction * 2.25  # ≈ 0.500 kg H2O/kg燃料
    
    # 燃料生产排放
    ch4_production_co2: float = 2.5  # kg CO2/kg CH4 (天然气制甲烷，含开采和加工)
    lox_production_kwh: float = 0.5  # kWh/kg LOX (空气分离)
    
    # 电网碳强度
    grid_carbon_intensity: float = 0.475  # kg CO2/kWh (全球平均, IEA 2023)
    grid_carbon_intensity_renewable: float = 0.05  # kg CO2/kWh (可再生能源)
    
    # 平流层排放比例
    stratosphere_emission_fraction: float = 0.40  # 40%的排放进入平流层
    stratosphere_h2o_residence_time: float = 3.0  # 年 (平流层水蒸气滞留时间)


# ============== 发射场定义 ==============
@dataclass
class LaunchSite:
    name: str
    short_name: str
    latitude: float
    max_launches_per_day: int = 1
    
    @property
    def abs_latitude(self) -> float:
        return abs(self.latitude)
    
    @property
    def rotation_velocity(self) -> float:
        return OMEGA_EARTH * R_EARTH * np.cos(np.radians(self.abs_latitude))
    
    @property
    def delta_v_loss(self) -> float:
        v_equator = OMEGA_EARTH * R_EARTH
        return v_equator - self.rotation_velocity


# 10个发射场 (按纬度排序)
LAUNCH_SITES = sorted([
    LaunchSite("Kourou (French Guiana)", "Kourou", 5.2),
    LaunchSite("Satish Dhawan (India)", "SDSC", 13.7),
    LaunchSite("Boca Chica (Texas)", "Texas", 26.0),
    LaunchSite("Cape Canaveral (Florida)", "Florida", 28.5),
    LaunchSite("Vandenberg (California)", "California", 34.7),
    LaunchSite("Wallops (Virginia)", "Virginia", 37.8),
    LaunchSite("Taiyuan (China)", "Taiyuan", 38.8),
    LaunchSite("Mahia (New Zealand)", "Mahia", 39.3),
    LaunchSite("Baikonur (Kazakhstan)", "Baikonur", 45.6),
    LaunchSite("Kodiak (Alaska)", "Alaska", 57.4),
], key=lambda x: x.abs_latitude)


# ============== 火箭参数 ==============
PAYLOAD_PER_LAUNCH = 125  # metric tons per launch
ALPHA = 0.10  # 结构系数
NUM_STAGES = 3
DELTA_V_BASE = 13300  # m/s (赤道发射到月球)

FUEL_PARAMS = MethaneFuelParams()
EMISSION_FACTORS = EmissionFactors(fuel_params=FUEL_PARAMS)


# ============== 核心计算函数 ==============

def fuel_ratio_multistage(delta_v: float, fuel_params: MethaneFuelParams = FUEL_PARAMS) -> float:
    """多级火箭燃料/载荷比"""
    ve = fuel_params.exhaust_velocity
    delta_v_per_stage = delta_v / NUM_STAGES
    R_stage = np.exp(delta_v_per_stage / ve)
    
    denominator = 1 - ALPHA * (R_stage - 1)
    if denominator <= 0:
        return np.inf
    
    k_stage = (R_stage - 1) / denominator
    
    total_fuel_ratio = 0
    remaining_ratio = 1.0
    
    for _ in range(NUM_STAGES):
        fuel_this_stage = remaining_ratio * k_stage
        total_fuel_ratio += fuel_this_stage
        remaining_ratio *= (1 + k_stage * (1 + ALPHA))
    
    return total_fuel_ratio


def rocket_fuel_per_ton(site: LaunchSite) -> float:
    """每吨载荷的燃料消耗 (kg燃料/ton载荷)"""
    delta_v = DELTA_V_BASE + site.delta_v_loss
    k = fuel_ratio_multistage(delta_v)
    return k * 1000  # kg fuel per metric ton payload


def rocket_energy_per_ton(site: LaunchSite) -> float:
    """火箭发射每吨载荷的能量消耗 (J/ton)"""
    return rocket_fuel_per_ton(site) * FUEL_PARAMS.specific_energy


# ============== 环境影响计算 ==============

@dataclass
class EnvironmentalImpact:
    """环境影响结果"""
    scenario_name: str
    completion_years: float
    total_payload_tons: float
    
    # 能量消耗
    total_energy_PJ: float
    rocket_energy_PJ: float
    elevator_energy_PJ: float
    
    # 燃料消耗
    total_fuel_Mt: float  # 百万吨
    rocket_fuel_Mt: float
    elevator_top_rocket_fuel_Mt: float
    
    # 排放量
    co2_combustion_Mt: float  # 燃烧直接排放
    co2_fuel_production_Mt: float  # 燃料生产排放
    co2_electricity_Mt: float  # 电力消耗排放
    h2o_total_Mt: float
    h2o_stratosphere_Mt: float  # 平流层水蒸气
    
    # 发射统计
    rocket_launches: int
    
    @property
    def co2_total_Mt(self) -> float:
        """总CO2排放 (百万吨)"""
        return self.co2_combustion_Mt + self.co2_fuel_production_Mt + self.co2_electricity_Mt
    
    @property
    def co2_per_ton_payload(self) -> float:
        """每吨载荷的CO2排放 (kg)"""
        return self.co2_total_Mt * 1e9 / self.total_payload_tons
    
    @property
    def annual_co2_Mt(self) -> float:
        """年均CO2排放 (百万吨/年)"""
        return self.co2_total_Mt / self.completion_years if self.completion_years > 0 else 0


def calculate_scenario_impact(
    scenario_name: str,
    completion_years: float,
    elevator_fraction: float,
    rocket_sites_used: List[LaunchSite],
    emission_factors: EmissionFactors = EMISSION_FACTORS,
    renewable_energy_fraction: float = 0.0
) -> EnvironmentalImpact:
    """
    计算单个情景的环境影响
    
    参数:
        scenario_name: 情景名称
        completion_years: 完成年限
        elevator_fraction: 电梯运输比例 (0-1)
        rocket_sites_used: 使用的火箭发射场列表
        emission_factors: 排放因子
        renewable_energy_fraction: 可再生能源占比 (0-1)
    """
    # 载荷分配
    elevator_payload = TOTAL_PAYLOAD * elevator_fraction
    rocket_payload = TOTAL_PAYLOAD * (1 - elevator_fraction)
    
    # ===== 电梯部分 =====
    elevator_energy = elevator_payload * ELEVATOR_SPECIFIC_ENERGY  # J
    elevator_lift_energy = elevator_payload * ELEVATOR_LIFT_ENERGY  # J (电力)
    elevator_top_rocket_energy = elevator_payload * ELEVATOR_TOP_ROCKET_ENERGY  # J (化学能)
    
    # 电梯顶端火箭燃料 (使用较小的燃料比，因为已在高轨道)
    # 顶端ΔV约1800 m/s (LOI + 轨道调整)
    elevator_top_delta_v = 1800
    elevator_top_fuel_ratio = fuel_ratio_multistage(elevator_top_delta_v)
    elevator_top_rocket_fuel = elevator_payload * 1000 * elevator_top_fuel_ratio  # kg
    
    # ===== 火箭部分 =====
    if rocket_payload > 0 and len(rocket_sites_used) > 0:
        launches_per_site = int(np.ceil(rocket_payload / PAYLOAD_PER_LAUNCH / len(rocket_sites_used)))
        total_launches = launches_per_site * len(rocket_sites_used)
        
        # 按发射场计算燃料消耗 (简化：均匀分配)
        rocket_fuel_total = 0
        rocket_energy_total = 0
        for site in rocket_sites_used:
            site_payload = rocket_payload / len(rocket_sites_used)
            site_fuel = site_payload * 1000 * fuel_ratio_multistage(DELTA_V_BASE + site.delta_v_loss)
            rocket_fuel_total += site_fuel
            rocket_energy_total += site_fuel * FUEL_PARAMS.specific_energy
    else:
        total_launches = 0
        rocket_fuel_total = 0
        rocket_energy_total = 0
    
    # ===== 总燃料消耗 =====
    total_fuel = rocket_fuel_total + elevator_top_rocket_fuel  # kg
    
    # ===== 排放计算 =====
    # 1. 燃烧直接排放
    co2_combustion = total_fuel * emission_factors.co2_per_kg_fuel  # kg
    h2o_combustion = total_fuel * emission_factors.h2o_per_kg_fuel  # kg
    
    # 2. 燃料生产排放
    ch4_mass = total_fuel * FUEL_PARAMS.ch4_mass_fraction
    lox_mass = total_fuel * FUEL_PARAMS.lox_mass_fraction
    
    co2_ch4_production = ch4_mass * emission_factors.ch4_production_co2
    co2_lox_production = lox_mass * emission_factors.lox_production_kwh * emission_factors.grid_carbon_intensity
    co2_fuel_production = co2_ch4_production + co2_lox_production
    
    # 3. 电梯电力消耗排放
    elevator_electricity_kwh = elevator_lift_energy / 3.6e6  # J -> kWh
    grid_intensity = (emission_factors.grid_carbon_intensity * (1 - renewable_energy_fraction) + 
                     emission_factors.grid_carbon_intensity_renewable * renewable_energy_fraction)
    co2_electricity = elevator_electricity_kwh * grid_intensity
    
    # 4. 平流层水蒸气 (仅火箭发射部分)
    # 电梯顶端释放点已在GEO以上，不经过平流层
    h2o_rocket = rocket_fuel_total * emission_factors.h2o_per_kg_fuel
    h2o_stratosphere = h2o_rocket * emission_factors.stratosphere_emission_fraction
    
    return EnvironmentalImpact(
        scenario_name=scenario_name,
        completion_years=completion_years,
        total_payload_tons=TOTAL_PAYLOAD,
        total_energy_PJ=(rocket_energy_total + elevator_energy) / 1e15,
        rocket_energy_PJ=rocket_energy_total / 1e15,
        elevator_energy_PJ=elevator_energy / 1e15,
        total_fuel_Mt=total_fuel / 1e9,
        rocket_fuel_Mt=rocket_fuel_total / 1e9,
        elevator_top_rocket_fuel_Mt=elevator_top_rocket_fuel / 1e9,
        co2_combustion_Mt=co2_combustion / 1e9,
        co2_fuel_production_Mt=co2_fuel_production / 1e9,
        co2_electricity_Mt=co2_electricity / 1e9,
        h2o_total_Mt=h2o_combustion / 1e9,
        h2o_stratosphere_Mt=h2o_stratosphere / 1e9,
        rocket_launches=total_launches
    )


# ============== 情景定义 ==============

def define_scenarios() -> List[Dict]:
    """
    定义四个基准情景 (基于任务一结果)
    """
    scenarios = []
    
    # 情景1: 纯火箭 (全部10个发射场)
    rocket_only_years = TOTAL_PAYLOAD / (len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH)
    scenarios.append({
        'name': 'Rocket Only',
        'completion_years': rocket_only_years,  # ~219年
        'elevator_fraction': 0.0,
        'rocket_sites': LAUNCH_SITES.copy(),
        'description': '纯火箭方案：10个发射场满负荷'
    })
    
    # 情景2: 纯电梯
    elevator_only_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY
    scenarios.append({
        'name': 'Elevator Only',
        'completion_years': elevator_only_years,  # ~186年
        'elevator_fraction': 1.0,
        'rocket_sites': [],
        'description': '纯电梯方案：3个电梯满负荷'
    })
    
    # 情景3: 混合方案 (最短时间)
    combined_capacity = TOTAL_ELEVATOR_CAPACITY + len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH
    combined_min_years = TOTAL_PAYLOAD / combined_capacity
    elevator_payload_min = TOTAL_ELEVATOR_CAPACITY * combined_min_years
    elevator_fraction_min = elevator_payload_min / TOTAL_PAYLOAD
    scenarios.append({
        'name': 'Combined (Min Time)',
        'completion_years': combined_min_years,  # ~101年
        'elevator_fraction': elevator_fraction_min,  # ~54%
        'rocket_sites': LAUNCH_SITES.copy(),
        'description': '混合方案(最短时间)：电梯+全部火箭站点'
    })
    
    # 情景4: 混合方案 (膝点优化)
    knee_years = 139  # 从任务一Pareto分析
    elevator_payload_knee = min(TOTAL_ELEVATOR_CAPACITY * knee_years, TOTAL_PAYLOAD)
    elevator_fraction_knee = elevator_payload_knee / TOTAL_PAYLOAD
    # 膝点时只需要部分低纬度发射场
    scenarios.append({
        'name': 'Combined (Knee Point)',
        'completion_years': knee_years,
        'elevator_fraction': elevator_fraction_knee,  # ~74.6%
        'rocket_sites': LAUNCH_SITES[:5],  # 只使用5个低纬度站点
        'description': '混合方案(膝点)：电梯优先+低纬度火箭'
    })
    
    return scenarios


def analyze_all_scenarios(renewable_fraction: float = 0.0) -> pd.DataFrame:
    """分析所有情景"""
    scenarios = define_scenarios()
    results = []
    
    for scenario in scenarios:
        impact = calculate_scenario_impact(
            scenario_name=scenario['name'],
            completion_years=scenario['completion_years'],
            elevator_fraction=scenario['elevator_fraction'],
            rocket_sites_used=scenario['rocket_sites'],
            renewable_energy_fraction=renewable_fraction
        )
        
        results.append({
            'Scenario': impact.scenario_name,
            'Years': impact.completion_years,
            'Elevator %': impact.elevator_energy_PJ / impact.total_energy_PJ * 100 if impact.total_energy_PJ > 0 else 0,
            'Energy (PJ)': impact.total_energy_PJ,
            'Fuel (Mt)': impact.total_fuel_Mt,
            'CO2 Total (Mt)': impact.co2_total_Mt,
            'CO2 Combustion (Mt)': impact.co2_combustion_Mt,
            'CO2 Production (Mt)': impact.co2_fuel_production_Mt,
            'CO2 Electricity (Mt)': impact.co2_electricity_Mt,
            'H2O Total (Mt)': impact.h2o_total_Mt,
            'H2O Stratosphere (Mt)': impact.h2o_stratosphere_Mt,
            'Annual CO2 (Mt/yr)': impact.annual_co2_Mt,
            'CO2/ton (kg)': impact.co2_per_ton_payload,
            'Launches': impact.rocket_launches,
        })
    
    return pd.DataFrame(results)


# ============== 模型优化：环境约束优化 ==============

def optimize_for_environment(
    max_annual_co2_Mt: float = 50.0,
    min_completion_years: float = 100.0,
    renewable_fraction: float = 0.5
) -> Dict:
    """
    带环境约束的优化模型
    
    目标: 最小化 总CO2排放
    约束:
        - 完成时间 >= min_completion_years
        - 年均CO2 <= max_annual_co2_Mt
    
    决策变量:
        - 电梯运输比例 (0-1)
        - 使用的火箭站点数量 (0-10)
    """
    best_result = None
    best_co2 = float('inf')
    
    # 网格搜索
    for elevator_frac in np.linspace(0, 1, 21):
        for n_sites in range(0, 11):
            # 计算完成时间
            elevator_capacity = TOTAL_ELEVATOR_CAPACITY
            rocket_capacity = n_sites * 365 * PAYLOAD_PER_LAUNCH if n_sites > 0 else 0
            
            elevator_payload = TOTAL_PAYLOAD * elevator_frac
            rocket_payload = TOTAL_PAYLOAD * (1 - elevator_frac)
            
            if elevator_frac == 1.0:
                years = elevator_payload / elevator_capacity if elevator_capacity > 0 else float('inf')
            elif elevator_frac == 0.0:
                years = rocket_payload / rocket_capacity if rocket_capacity > 0 else float('inf')
            else:
                # 电梯和火箭并行
                elevator_years = elevator_payload / elevator_capacity
                rocket_years = rocket_payload / rocket_capacity if rocket_capacity > 0 else float('inf')
                years = max(elevator_years, rocket_years)
            
            if years < min_completion_years:
                continue
            
            # 计算环境影响
            sites_used = LAUNCH_SITES[:n_sites] if n_sites > 0 else []
            impact = calculate_scenario_impact(
                scenario_name='Optimization',
                completion_years=years,
                elevator_fraction=elevator_frac,
                rocket_sites_used=sites_used,
                renewable_energy_fraction=renewable_fraction
            )
            
            # 检查约束
            if impact.annual_co2_Mt > max_annual_co2_Mt:
                continue
            
            # 更新最优解
            if impact.co2_total_Mt < best_co2:
                best_co2 = impact.co2_total_Mt
                best_result = {
                    'elevator_fraction': elevator_frac,
                    'n_sites': n_sites,
                    'years': years,
                    'impact': impact
                }
    
    return best_result


# ============== 可视化 ==============

def plot_environmental_comparison(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """绘制环境影响对比图（2张子图版本，马卡龙色系论文风格）"""
    df = analyze_all_scenarios()
    
    # 马卡龙色系
    MACARON = {
        'pink': '#F8B4C0',       # 樱花粉
        'mint': '#A8E6CF',       # 薄荷绿
        'lavender': '#C5A3FF',   # 薰衣草紫
        'peach': '#FFCBA4',      # 蜜桃橙
        'sky': '#A8D8EA',        # 天空蓝
        'cream': '#FFEAA7',      # 奶油黄
        'rose': '#DDA0DD',       # 玫瑰粉
        'sage': '#B5EAD7',       # 鼠尾草绿
    }
    
    # 论文风格设置
    plt.rcParams.update({
        'font.size': 11,
        'axes.labelsize': 12,
        'axes.titlesize': 13,
        'xtick.labelsize': 10,
        'ytick.labelsize': 10,
        'legend.fontsize': 9,
        'axes.linewidth': 0.8,
        'axes.edgecolor': '#666666',
    })
    
    fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))
    fig.patch.set_facecolor('white')
    
    scenarios = df['Scenario'].tolist()
    x = np.arange(len(scenarios))
    width = 0.35
    
    # ========== 左图: CO2排放分解 + H2O排放 ==========
    ax1 = axes[0]
    ax1.set_facecolor('white')
    
    # CO2数据（左侧柱子）
    co2_combustion = df['CO2 Combustion (Mt)'].values
    co2_production = df['CO2 Production (Mt)'].values
    co2_electricity = df['CO2 Electricity (Mt)'].values
    
    bars1 = ax1.bar(x - width/2, co2_combustion, width, label='Combustion', 
                    color=MACARON['pink'], edgecolor='#666666', linewidth=0.5)
    bars2 = ax1.bar(x - width/2, co2_production, width, bottom=co2_combustion, 
                    label='Fuel Prod.', color=MACARON['mint'], edgecolor='#666666', linewidth=0.5)
    bars3 = ax1.bar(x - width/2, co2_electricity, width, bottom=co2_combustion+co2_production, 
                    label='Electricity', color=MACARON['sky'], edgecolor='#666666', linewidth=0.5)
    
    ax1.set_ylabel('CO₂ Emissions (Mt)', color='#333333')
    ax1.tick_params(axis='y', labelcolor='#333333')
    
    # 添加CO2总量标签
    totals = df['CO2 Total (Mt)'].values
    for i, total in enumerate(totals):
        ax1.text(i - width/2, total + 300, f'{total:.0f}', ha='center', fontsize=9, fontweight='bold')
    
    # H2O数据（右侧柱子）
    ax1_twin = ax1.twinx()
    
    h2o_surface = df['H2O Total (Mt)'].values - df['H2O Stratosphere (Mt)'].values
    h2o_strat = df['H2O Stratosphere (Mt)'].values
    
    bars4 = ax1_twin.bar(x + width/2, h2o_surface, width, label='Troposphere', 
                         color=MACARON['sage'], edgecolor='#666666', linewidth=0.5)
    bars5 = ax1_twin.bar(x + width/2, h2o_strat, width, bottom=h2o_surface, 
                         label='Stratosphere', color=MACARON['rose'], edgecolor='#666666', linewidth=0.5)
    
    ax1_twin.set_ylabel('H₂O Emissions (Mt)', color='#7B68A0')
    ax1_twin.tick_params(axis='y', labelcolor='#7B68A0')
    
    ax1.set_title('(a) CO₂ & H₂O Emissions', fontweight='bold', pad=10)
    ax1.set_xticks(x)
    ax1.set_xticklabels(scenarios, rotation=20, ha='right')
    ax1.grid(True, alpha=0.2, axis='y', linestyle='--')
    ax1.set_axisbelow(True)
    
    # 合并图例
    lines1, labels1 = ax1.get_legend_handles_labels()
    lines2, labels2 = ax1_twin.get_legend_handles_labels()
    ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right', 
               framealpha=0.95, edgecolor='#cccccc', ncol=2)
    
    # ========== 右图: 年均CO2排放 + 碳强度 ==========
    ax2 = axes[1]
    ax2.set_facecolor('white')
    
    annual_co2 = df['Annual CO2 (Mt/yr)'].values
    co2_per_ton = df['CO2/ton (kg)'].values
    
    # 年均CO2（左侧柱子）- 使用渐变马卡龙色
    colors_annual = [MACARON['pink'], MACARON['mint'], MACARON['sky'], MACARON['sage']]
    bars_annual = ax2.bar(x - width/2, annual_co2, width, color=colors_annual, 
                          edgecolor='#666666', linewidth=0.5)
    
    ax2.set_ylabel('Annual CO₂ (Mt/yr)', color='#333333')
    ax2.tick_params(axis='y', labelcolor='#333333')
    
    # 添加年均CO2标签
    global_annual_co2 = 37000
    for i, val in enumerate(annual_co2):
        pct = val / global_annual_co2 * 100
        ax2.text(i - width/2, val + 2, f'{val:.1f}\n({pct:.2f}%)', ha='center', fontsize=8)
    
    # 碳强度（右侧柱子）
    ax2_twin = ax2.twinx()
    
    bars_intensity = ax2_twin.bar(x + width/2, co2_per_ton, width, color=MACARON['cream'], 
                                   edgecolor='#666666', linewidth=0.5)
    
    ax2_twin.set_ylabel('CO₂ per ton payload (kg)', color='#B8860B')
    ax2_twin.tick_params(axis='y', labelcolor='#B8860B')
    
    # 添加碳强度标签
    for i, val in enumerate(co2_per_ton):
        ax2_twin.text(i + width/2, val + 3000, f'{val:.0f}', ha='center', fontsize=9, fontweight='bold')
    
    ax2.set_title('(b) Annual Emissions & Carbon Intensity', fontweight='bold', pad=10)
    ax2.set_xticks(x)
    ax2.set_xticklabels(scenarios, rotation=20, ha='right')
    ax2.grid(True, alpha=0.2, axis='y', linestyle='--')
    ax2.set_axisbelow(True)
    
    # 图例说明
    from matplotlib.patches import Patch
    legend_elements = [Patch(facecolor=MACARON['pink'], edgecolor='#666666', label='Annual CO₂'),
                       Patch(facecolor=MACARON['cream'], edgecolor='#666666', label='Carbon Intensity')]
    ax2.legend(handles=legend_elements, loc='upper right', framealpha=0.95, edgecolor='#cccccc')
    
    plt.tight_layout()
    plt.savefig(f'{save_dir}/environmental_comparison.png', dpi=300, bbox_inches='tight', 
                facecolor='white', edgecolor='none')
    print(f"环境影响对比图已保存至: {save_dir}/environmental_comparison.png")
    
    return df


def plot_radar_chart(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """单独绘制雷达图（马卡龙色系，无标题，论文风格）"""
    df = analyze_all_scenarios()
    
    # 马卡龙色系
    MACARON = {
        'pink': '#F8B4C0',       # 樱花粉
        'mint': '#A8E6CF',       # 薄荷绿
        'lavender': '#C5A3FF',   # 薰衣草紫
        'sky': '#A8D8EA',        # 天空蓝
    }
    colors = [MACARON['pink'], MACARON['mint'], MACARON['sky'], MACARON['lavender']]
    
    # 论文风格设置
    plt.rcParams.update({
        'font.size': 11,
        'axes.labelsize': 12,
        'xtick.labelsize': 11,
        'ytick.labelsize': 10,
        'legend.fontsize': 10,
    })
    
    fig, ax = plt.subplots(figsize=(5, 5), subplot_kw=dict(projection='polar'))
    fig.patch.set_facecolor('white')
    
    categories = ['Time\n(inv)', 'Energy\n(inv)', 'CO₂\n(inv)', 'Launches\n(inv)', 'Strat. H₂O\n(inv)']
    N = len(categories)
    
    # 归一化数据 (取逆以使"更好"在外围)
    max_vals = df[['Years', 'Energy (PJ)', 'CO2 Total (Mt)', 'Launches', 'H2O Stratosphere (Mt)']].max()
    
    angles = np.linspace(0, 2*np.pi, N, endpoint=False).tolist()
    angles += angles[:1]
    
    for i, (idx, row) in enumerate(df.iterrows()):
        values = [
            1 - row['Years'] / max_vals['Years'],
            1 - row['Energy (PJ)'] / max_vals['Energy (PJ)'],
            1 - row['CO2 Total (Mt)'] / max_vals['CO2 Total (Mt)'],
            1 - row['Launches'] / max_vals['Launches'] if max_vals['Launches'] > 0 else 1,
            1 - row['H2O Stratosphere (Mt)'] / max_vals['H2O Stratosphere (Mt)'] if max_vals['H2O Stratosphere (Mt)'] > 0 else 1
        ]
        values += values[:1]
        ax.plot(angles, values, 'o-', linewidth=2, label=row['Scenario'], 
                color=colors[i], markersize=5)
        ax.fill(angles, values, alpha=0.2, color=colors[i])
    
    ax.set_xticks(angles[:-1])
    ax.set_xticklabels(categories, fontsize=10)
    ax.set_ylim(0, 1)
    
    # 设置径向刻度
    ax.set_yticks([0.2, 0.4, 0.6, 0.8, 1.0])
    ax.set_yticklabels(['0.2', '0.4', '0.6', '0.8', '1.0'], fontsize=8, color='#666666')
    
    # 图例放在右侧
    ax.legend(loc='upper left', bbox_to_anchor=(1.05, 1.0), framealpha=0.95, 
              edgecolor='#cccccc')
    
    # 网格样式
    ax.grid(True, alpha=0.3, linestyle='--')
    ax.spines['polar'].set_color('#666666')
    ax.spines['polar'].set_linewidth(0.8)
    
    plt.tight_layout()
    plt.savefig(f'{save_dir}/radar_chart.png', dpi=300, bbox_inches='tight',
                facecolor='white', edgecolor='none')
    print(f"雷达图已保存至: {save_dir}/radar_chart.png")
    
    return df


def plot_marginal_benefit(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """单独绘制边际收益分析图（马卡龙色系，无标题，论文风格）"""
    
    # 马卡龙色系
    MACARON = {
        'sky': '#A8D8EA',        # 天空蓝
        'pink': '#F8B4C0',       # 樱花粉
        'mint': '#A8E6CF',       # 薄荷绿
    }
    
    # 论文风格设置
    plt.rcParams.update({
        'font.size': 11,
        'axes.labelsize': 12,
        'xtick.labelsize': 10,
        'ytick.labelsize': 10,
        'legend.fontsize': 10,
        'axes.linewidth': 0.8,
    })
    
    # 生成Pareto数据
    years_range = np.linspace(101, 200, 100)
    
    results = []
    for years in years_range:
        # 计算电梯运输比例
        elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * years, TOTAL_PAYLOAD)
        elevator_frac = elevator_payload / TOTAL_PAYLOAD
        
        # 确定使用的发射场
        rocket_payload = TOTAL_PAYLOAD - elevator_payload
        if rocket_payload > 0:
            rocket_capacity_per_site = 365 * years * PAYLOAD_PER_LAUNCH
            n_sites_needed = int(np.ceil(rocket_payload / rocket_capacity_per_site))
            n_sites_needed = min(n_sites_needed, 10)
            sites = LAUNCH_SITES[:n_sites_needed]
        else:
            sites = []
        
        impact = calculate_scenario_impact(
            scenario_name=f'{years:.0f}yr',
            completion_years=years,
            elevator_fraction=elevator_frac,
            rocket_sites_used=sites
        )
        
        results.append({
            'Years': years,
            'Energy (PJ)': impact.total_energy_PJ,
            'CO2 (Mt)': impact.co2_total_Mt,
        })
    
    df = pd.DataFrame(results)
    
    # 计算边际减少率
    T = df['Years'].values
    E = df['Energy (PJ)'].values
    C = df['CO2 (Mt)'].values
    
    dE_dT = np.gradient(E, T)
    dC_dT = np.gradient(C, T)
    
    energy_reduction_rate = -dE_dT  # 每多1年节省的能量
    co2_reduction_rate = -dC_dT     # 每多1年减少的CO2
    
    # 绘图（黄金比例）
    fig, ax = plt.subplots(figsize=(8, 8*0.618))
    fig.patch.set_facecolor('white')
    ax.set_facecolor('white')
    
    # 能量减少率
    ax.plot(T, energy_reduction_rate, color=MACARON['sky'], linewidth=2, 
            label='Energy reduction rate')
    ax.fill_between(T, energy_reduction_rate, alpha=0.2, color=MACARON['sky'])
    
    # CO2减少率（使用第二Y轴）
    ax2 = ax.twinx()
    ax2.plot(T, co2_reduction_rate, color=MACARON['pink'], linewidth=2, 
             label='CO₂ reduction rate')
    ax2.fill_between(T, co2_reduction_rate, alpha=0.2, color=MACARON['pink'])
    
    # 推荐点
    recommended_year = 186
    ax.axvline(x=recommended_year, color=MACARON['mint'], linestyle='--', linewidth=2,
               label=f'Recommended: {recommended_year}yr')
    
    ax.set_xlabel('Completion Time (years)')
    ax.set_ylabel('Energy Reduction (PJ/yr)', color='#5B9BD5')
    ax.tick_params(axis='y', labelcolor='#5B9BD5')
    
    ax2.set_ylabel('CO₂ Reduction (Mt/yr)', color='#E57373')
    ax2.tick_params(axis='y', labelcolor='#E57373')
    
    ax.grid(True, alpha=0.2, linestyle='--')
    ax.set_axisbelow(True)
    
    # 合并图例
    lines1, labels1 = ax.get_legend_handles_labels()
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax.legend(lines1 + lines2, labels1 + labels2, loc='upper right', 
              framealpha=0.95, edgecolor='#cccccc')
    
    plt.tight_layout()
    plt.savefig(f'{save_dir}/marginal_benefit.png', dpi=300, bbox_inches='tight',
                facecolor='white', edgecolor='none')
    print(f"边际收益图已保存至: {save_dir}/marginal_benefit.png")
    
    return df


def plot_pareto_with_environment(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """绘制时间-能量-环境三目标Pareto图"""
    
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # 生成Pareto数据
    years_range = np.linspace(101, 220, 50)
    
    results = []
    for years in years_range:
        # 计算电梯运输比例
        elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * years, TOTAL_PAYLOAD)
        elevator_frac = elevator_payload / TOTAL_PAYLOAD
        
        # 确定使用的发射场
        rocket_payload = TOTAL_PAYLOAD - elevator_payload
        if rocket_payload > 0:
            rocket_capacity_per_site = 365 * years * PAYLOAD_PER_LAUNCH
            n_sites_needed = int(np.ceil(rocket_payload / rocket_capacity_per_site))
            n_sites_needed = min(n_sites_needed, 10)
            sites = LAUNCH_SITES[:n_sites_needed]
        else:
            sites = []
        
        impact = calculate_scenario_impact(
            scenario_name=f'{years:.0f}yr',
            completion_years=years,
            elevator_fraction=elevator_frac,
            rocket_sites_used=sites
        )
        
        results.append({
            'Years': years,
            'Energy (PJ)': impact.total_energy_PJ,
            'CO2 (Mt)': impact.co2_total_Mt,
            'Elevator %': elevator_frac * 100
        })
    
    df = pd.DataFrame(results)
    
    # ========== 图1: 时间 vs CO2 ==========
    ax1 = axes[0]
    
    scatter = ax1.scatter(df['Years'], df['CO2 (Mt)'], c=df['Elevator %'], 
                         cmap='RdYlGn', s=50, alpha=0.8)
    plt.colorbar(scatter, ax=ax1, label='Elevator %')
    
    ax1.set_xlabel('Completion Time (years)', fontsize=12)
    ax1.set_ylabel('Total CO₂ Emissions (Mt)', fontsize=12)
    ax1.set_title('Pareto Front: Time vs CO₂ Emissions\n(Color: Elevator Usage)', fontsize=14)
    ax1.grid(True, alpha=0.3)
    
    # 标记关键点
    ax1.axvline(x=139, color='red', linestyle='--', alpha=0.7, label='Knee Point (139yr)')
    ax1.axvline(x=186, color='green', linestyle='--', alpha=0.7, label='Elevator Only (186yr)')
    ax1.legend()
    
    # ========== 图2: 能量 vs CO2 ==========
    ax2 = axes[1]
    
    scatter = ax2.scatter(df['Energy (PJ)'], df['CO2 (Mt)'], c=df['Years'], 
                         cmap='viridis', s=50, alpha=0.8)
    plt.colorbar(scatter, ax=ax2, label='Years')
    
    ax2.set_xlabel('Total Energy (PJ)', fontsize=12)
    ax2.set_ylabel('Total CO₂ Emissions (Mt)', fontsize=12)
    ax2.set_title('Pareto Front: Energy vs CO₂ Emissions\n(Color: Completion Time)', fontsize=14)
    ax2.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(f'{save_dir}/pareto_environmental.png', dpi=150, bbox_inches='tight')
    print(f"Pareto环境分析图已保存至: {save_dir}/pareto_environmental.png")
    
    return df


def plot_mitigation_strategies(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """绘制减排策略对比图"""
    
    # 基准情景 (混合膝点)
    baseline = calculate_scenario_impact(
        scenario_name='Baseline',
        completion_years=139,
        elevator_fraction=0.746,
        rocket_sites_used=LAUNCH_SITES[:5],
        renewable_energy_fraction=0.0
    )
    
    # 策略1: 50%可再生能源
    strategy1 = calculate_scenario_impact(
        scenario_name='50% Renewable',
        completion_years=139,
        elevator_fraction=0.746,
        rocket_sites_used=LAUNCH_SITES[:5],
        renewable_energy_fraction=0.5
    )
    
    # 策略2: 100%可再生能源
    strategy2 = calculate_scenario_impact(
        scenario_name='100% Renewable',
        completion_years=139,
        elevator_fraction=0.746,
        rocket_sites_used=LAUNCH_SITES[:5],
        renewable_energy_fraction=1.0
    )
    
    # 策略3: 延长工期到纯电梯
    strategy3 = calculate_scenario_impact(
        scenario_name='Elevator Only',
        completion_years=186,
        elevator_fraction=1.0,
        rocket_sites_used=[],
        renewable_energy_fraction=0.0
    )
    
    # 策略4: 纯电梯 + 100%可再生
    strategy4 = calculate_scenario_impact(
        scenario_name='Elevator + 100% Renewable',
        completion_years=186,
        elevator_fraction=1.0,
        rocket_sites_used=[],
        renewable_energy_fraction=1.0
    )
    
    strategies = [baseline, strategy1, strategy2, strategy3, strategy4]
    
    fig, ax = plt.subplots(figsize=(12, 7))
    
    names = [s.scenario_name for s in strategies]
    co2_vals = [s.co2_total_Mt for s in strategies]
    reductions = [(baseline.co2_total_Mt - s.co2_total_Mt) / baseline.co2_total_Mt * 100 
                  for s in strategies]
    
    x = np.arange(len(names))
    colors = ['#FF6B6B', '#FFE66D', '#4ECDC4', '#45B7D1', '#96CEB4']
    
    bars = ax.bar(x, co2_vals, color=colors)
    
    ax.set_ylabel('Total CO₂ Emissions (Mt)', fontsize=12)
    ax.set_title('CO₂ Reduction Strategies Comparison\n(Knee Point baseline: 139 years)', fontsize=14)
    ax.set_xticks(x)
    ax.set_xticklabels(names, rotation=15, ha='right')
    ax.grid(True, alpha=0.3, axis='y')
    
    # 添加数值和减排比例
    for i, (bar, val, red) in enumerate(zip(bars, co2_vals, reductions)):
        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 20, 
               f'{val:.0f} Mt\n({red:+.1f}%)', ha='center', fontsize=10)
    
    # 添加基准线
    ax.axhline(y=baseline.co2_total_Mt, color='red', linestyle='--', alpha=0.5, 
               label=f'Baseline: {baseline.co2_total_Mt:.0f} Mt')
    ax.legend()
    
    plt.tight_layout()
    plt.savefig(f'{save_dir}/mitigation_strategies.png', dpi=150, bbox_inches='tight')
    print(f"减排策略对比图已保存至: {save_dir}/mitigation_strategies.png")


def plot_comprehensive_summary(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """绘制综合汇总图"""
    
    df = analyze_all_scenarios()
    
    fig = plt.figure(figsize=(16, 12))
    
    # 创建4x4网格
    gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)
    
    scenarios = df['Scenario'].tolist()
    
    # ========== 左上: 雷达图 ==========
    ax1 = fig.add_subplot(gs[0, 0], projection='polar')
    
    categories = ['Time\n(inv)', 'Energy\n(inv)', 'CO2\n(inv)', 'Launches\n(inv)', 'Strat H2O\n(inv)']
    N = len(categories)
    
    # 归一化数据 (取逆以使"更好"在外围)
    max_vals = df[['Years', 'Energy (PJ)', 'CO2 Total (Mt)', 'Launches', 'H2O Stratosphere (Mt)']].max()
    
    angles = np.linspace(0, 2*np.pi, N, endpoint=False).tolist()
    angles += angles[:1]
    
    colors = ['#FF6B6B', '#4ECDC4', '#45B7D1', '#96CEB4']
    
    for i, (idx, row) in enumerate(df.iterrows()):
        values = [
            1 - row['Years'] / max_vals['Years'],
            1 - row['Energy (PJ)'] / max_vals['Energy (PJ)'],
            1 - row['CO2 Total (Mt)'] / max_vals['CO2 Total (Mt)'],
            1 - row['Launches'] / max_vals['Launches'] if max_vals['Launches'] > 0 else 1,
            1 - row['H2O Stratosphere (Mt)'] / max_vals['H2O Stratosphere (Mt)'] if max_vals['H2O Stratosphere (Mt)'] > 0 else 1
        ]
        values += values[:1]
        ax1.plot(angles, values, 'o-', linewidth=2, label=row['Scenario'], color=colors[i])
        ax1.fill(angles, values, alpha=0.15, color=colors[i])
    
    ax1.set_xticks(angles[:-1])
    ax1.set_xticklabels(categories, fontsize=9)
    ax1.set_title('Multi-Criteria Comparison\n(outer = better)', fontsize=12, pad=20)
    ax1.legend(loc='upper right', bbox_to_anchor=(1.3, 1.0), fontsize=8)
    
    # ========== 右上: 汇总表格 ==========
    ax2 = fig.add_subplot(gs[0, 1:])
    ax2.axis('off')
    
    table_data = []
    for idx, row in df.iterrows():
        table_data.append([
            row['Scenario'],
            f"{row['Years']:.0f}",
            f"{row['Energy (PJ)']:.0f}",
            f"{row['CO2 Total (Mt)']:.0f}",
            f"{row['H2O Stratosphere (Mt)']:.1f}",
            f"{row['Annual CO2 (Mt/yr)']:.1f}",
        ])
    
    table = ax2.table(
        cellText=table_data,
        colLabels=['Scenario', 'Years', 'Energy(PJ)', 'CO2(Mt)', 'Strat.H2O(Mt)', 'Annual CO2'],
        loc='center',
        cellLoc='center'
    )
    table.auto_set_font_size(False)
    table.set_fontsize(10)
    table.scale(1.2, 1.8)
    
    # 高亮最优值
    for i in range(len(df)):
        for j in range(1, 6):
            cell = table[(i+1, j)]
            if j == 1 and df.iloc[i]['Years'] == df['Years'].min():
                cell.set_facecolor('#90EE90')
            elif j == 2 and df.iloc[i]['Energy (PJ)'] == df['Energy (PJ)'].min():
                cell.set_facecolor('#90EE90')
            elif j == 3 and df.iloc[i]['CO2 Total (Mt)'] == df['CO2 Total (Mt)'].min():
                cell.set_facecolor('#90EE90')
    
    ax2.set_title('Scenario Summary Table\n(Green = Best in Category)', fontsize=12, pad=10)
    
    # ========== 中间行: CO2分解 + 时间线 ==========
    ax3 = fig.add_subplot(gs[1, :2])
    
    x = np.arange(len(scenarios))
    width = 0.25
    
    ax3.bar(x - width, df['CO2 Combustion (Mt)'], width, label='Combustion', color='#FF6B6B')
    ax3.bar(x, df['CO2 Production (Mt)'], width, label='Fuel Production', color='#4ECDC4')
    ax3.bar(x + width, df['CO2 Electricity (Mt)'], width, label='Electricity', color='#45B7D1')
    
    ax3.set_ylabel('CO₂ (Mt)', fontsize=11)
    ax3.set_title('CO₂ Emissions by Source', fontsize=12)
    ax3.set_xticks(x)
    ax3.set_xticklabels(scenarios, rotation=15, ha='right')
    ax3.legend()
    ax3.grid(True, alpha=0.3, axis='y')
    
    # ========== 右侧: 时间-CO2权衡 ==========
    ax4 = fig.add_subplot(gs[1, 2])
    
    ax4.scatter(df['Years'], df['CO2 Total (Mt)'], s=200, c=colors, edgecolors='black', linewidth=2)
    
    for i, row in df.iterrows():
        ax4.annotate(row['Scenario'].replace(' ', '\n'), 
                    (row['Years'], row['CO2 Total (Mt)']),
                    textcoords="offset points", xytext=(0, 15), ha='center', fontsize=8)
    
    ax4.set_xlabel('Completion Time (years)', fontsize=11)
    ax4.set_ylabel('Total CO₂ (Mt)', fontsize=11)
    ax4.set_title('Time-CO₂ Trade-off', fontsize=12)
    ax4.grid(True, alpha=0.3)
    
    # ========== 底行: 结论文字 ==========
    ax5 = fig.add_subplot(gs[2, :])
    ax5.axis('off')
    
    conclusions = """
    KEY FINDINGS (LOX/CH4 Fuel):
    
    1. Rocket Only: Highest CO₂ emissions ({:.0f} Mt) due to large fuel consumption ({:.1f} Mt fuel).
       - Annual emissions: {:.1f} Mt/yr (≈{:.3f}% of global annual emissions)
       - Stratospheric H₂O injection: {:.1f} Mt (potential climate impact)
    
    2. Elevator Only: Lowest CO₂ emissions ({:.0f} Mt), but longest timeline ({:.0f} years).
       - CO₂ reduction vs Rocket: {:.0f}% (primary savings from avoiding rocket fuel production)
       - No stratospheric H₂O injection from main transport
    
    3. Combined (Knee Point): Best trade-off at 139 years with {:.0f} Mt CO₂.
       - 74.6% payload via elevator, 25.4% via low-latitude rockets
       - Reasonable timeline with significant environmental benefit
    
    4. RECOMMENDATIONS for Minimizing Environmental Impact:
       a) Prioritize space elevator transport (saves ~{:.0f}% CO₂ vs pure rocket)
       b) Use renewable energy for elevator operation (reduces electricity emissions)
       c) If rockets needed, prefer low-latitude sites (lower fuel per kg)
       d) Consider carbon capture for CH4 production phase
    """.format(
        df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0],
        df.loc[df['Scenario']=='Rocket Only', 'Fuel (Mt)'].values[0],
        df.loc[df['Scenario']=='Rocket Only', 'Annual CO2 (Mt/yr)'].values[0],
        df.loc[df['Scenario']=='Rocket Only', 'Annual CO2 (Mt/yr)'].values[0] / 37000 * 100,
        df.loc[df['Scenario']=='Rocket Only', 'H2O Stratosphere (Mt)'].values[0],
        df.loc[df['Scenario']=='Elevator Only', 'CO2 Total (Mt)'].values[0],
        df.loc[df['Scenario']=='Elevator Only', 'Years'].values[0],
        (1 - df.loc[df['Scenario']=='Elevator Only', 'CO2 Total (Mt)'].values[0] / 
         df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]) * 100,
        df.loc[df['Scenario']=='Combined (Knee Point)', 'CO2 Total (Mt)'].values[0],
        (1 - df.loc[df['Scenario']=='Elevator Only', 'CO2 Total (Mt)'].values[0] / 
         df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]) * 100
    )
    
    ax5.text(0.02, 0.95, conclusions, transform=ax5.transAxes, fontsize=10,
            verticalalignment='top', fontfamily='monospace',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    plt.suptitle('Moon Colony Environmental Impact Analysis\n'
                 '(100 Million Tons Payload, LOX/CH4 Fuel)', fontsize=14, y=0.98)
    
    plt.savefig(f'{save_dir}/comprehensive_summary.png', dpi=150, bbox_inches='tight')
    print(f"综合汇总图已保存至: {save_dir}/comprehensive_summary.png")


# ============== 报告生成 ==============

def generate_report(save_dir: str = '/Volumes/Files/code/mm/20260130_b/p4'):
    """生成分析报告"""
    df = analyze_all_scenarios()
    
    # 优化结果
    opt_result = optimize_for_environment(
        max_annual_co2_Mt=50.0,
        min_completion_years=120.0,
        renewable_fraction=0.5
    )
    
    report = []
    report.append("=" * 80)
    report.append("TASK 4: ENVIRONMENTAL IMPACT ANALYSIS")
    report.append("Moon Colony Construction - 100 Million Metric Tons")
    report.append("Fuel: LOX/CH4 (Liquid Oxygen / Methane)")
    report.append("=" * 80)
    
    report.append("\n" + "=" * 80)
    report.append("1. FUEL PARAMETERS (LOX/CH4)")
    report.append("=" * 80)
    report.append(f"  Specific Impulse (Isp): {FUEL_PARAMS.isp} s")
    report.append(f"  Exhaust Velocity: {FUEL_PARAMS.exhaust_velocity:.0f} m/s")
    report.append(f"  Specific Energy: {FUEL_PARAMS.specific_energy/1e6:.1f} MJ/kg")
    report.append(f"  O2/CH4 Mass Ratio: {FUEL_PARAMS.ox_fuel_ratio}:1")
    report.append(f"  CH4 Mass Fraction: {FUEL_PARAMS.ch4_mass_fraction*100:.1f}%")
    report.append(f"\n  Combustion Reaction: CH4 + 2O2 → CO2 + 2H2O")
    report.append(f"  CO2 per kg fuel: {EMISSION_FACTORS.co2_per_kg_fuel:.3f} kg")
    report.append(f"  H2O per kg fuel: {EMISSION_FACTORS.h2o_per_kg_fuel:.3f} kg")
    
    report.append("\n" + "=" * 80)
    report.append("2. SCENARIO COMPARISON")
    report.append("=" * 80)
    
    for idx, row in df.iterrows():
        report.append(f"\n  {row['Scenario']}:")
        report.append(f"    Completion Time: {row['Years']:.0f} years")
        report.append(f"    Total Energy: {row['Energy (PJ)']:.0f} PJ")
        report.append(f"    Total Fuel: {row['Fuel (Mt)']:.2f} Mt")
        report.append(f"    CO2 Emissions: {row['CO2 Total (Mt)']:.0f} Mt")
        report.append(f"      - Combustion: {row['CO2 Combustion (Mt)']:.0f} Mt")
        report.append(f"      - Fuel Production: {row['CO2 Production (Mt)']:.0f} Mt")
        report.append(f"      - Electricity: {row['CO2 Electricity (Mt)']:.1f} Mt")
        report.append(f"    H2O Emissions: {row['H2O Total (Mt)']:.2f} Mt")
        report.append(f"      - Stratospheric: {row['H2O Stratosphere (Mt)']:.2f} Mt")
        report.append(f"    Annual CO2: {row['Annual CO2 (Mt/yr)']:.2f} Mt/yr ({row['Annual CO2 (Mt/yr)']/370:.3f}% of global)")
        report.append(f"    CO2 per ton payload: {row['CO2/ton (kg)']:.0f} kg/ton")
    
    report.append("\n" + "=" * 80)
    report.append("3. ENVIRONMENTAL IMPACT RANKING")
    report.append("=" * 80)
    
    df_sorted = df.sort_values('CO2 Total (Mt)')
    report.append("\n  By Total CO2 (lowest first):")
    for i, (idx, row) in enumerate(df_sorted.iterrows(), 1):
        report.append(f"    {i}. {row['Scenario']}: {row['CO2 Total (Mt)']:.0f} Mt")
    
    rocket_co2 = df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]
    report.append("\n  CO2 Reduction vs Rocket Only:")
    for idx, row in df.iterrows():
        if row['Scenario'] != 'Rocket Only':
            reduction = (rocket_co2 - row['CO2 Total (Mt)']) / rocket_co2 * 100
            report.append(f"    {row['Scenario']}: -{reduction:.1f}%")
    
    report.append("\n" + "=" * 80)
    report.append("4. OPTIMIZATION FOR MINIMUM ENVIRONMENTAL IMPACT")
    report.append("=" * 80)
    
    if opt_result:
        report.append(f"\n  Optimization Constraints:")
        report.append(f"    - Minimum completion time: 120 years")
        report.append(f"    - Maximum annual CO2: 50 Mt/yr")
        report.append(f"    - Renewable energy fraction: 50%")
        report.append(f"\n  Optimal Solution:")
        report.append(f"    - Elevator fraction: {opt_result['elevator_fraction']*100:.1f}%")
        report.append(f"    - Rocket sites used: {opt_result['n_sites']}")
        report.append(f"    - Completion time: {opt_result['years']:.0f} years")
        report.append(f"    - Total CO2: {opt_result['impact'].co2_total_Mt:.0f} Mt")
        report.append(f"    - Annual CO2: {opt_result['impact'].annual_co2_Mt:.2f} Mt/yr")
    
    report.append("\n" + "=" * 80)
    report.append("5. MODEL ADJUSTMENT RECOMMENDATIONS")
    report.append("=" * 80)
    report.append("""
  To minimize environmental impact, adjust the model as follows:

  A. OBJECTIVE FUNCTION MODIFICATION:
     Original: min(α×Time + β×Energy)
     Modified: min(α×Time + β×Energy + γ×CO2_total + δ×H2O_stratosphere)
     
     Where γ and δ are weights for environmental costs.

  B. RECOMMENDED STRATEGIES (by effectiveness):
     
     1. Prioritize Elevator Transport (CO2 reduction: ~70-80%)
        - Elevator has no direct combustion emissions
        - Top rocket stage emissions much smaller than full rocket launch
        
     2. Use Renewable Energy for Elevator (CO2 reduction: 5-10%)
        - Replace grid electricity with solar/wind
        - Reduces electricity-related emissions
        
     3. Low-Latitude Launch Sites Only (CO2 reduction: 3-5%)
        - Lower ΔV requirements = less fuel
        - Kourou (5.2°) vs Alaska (57.4°): ~10% fuel difference
        
     4. Green Methane Production (CO2 reduction: 15-20%)
        - Use bio-methane or synthetic methane from CO2 capture
        - Reduces fuel production emissions
        
     5. Extend Timeline (CO2 reduction: proportional)
        - More time → more elevator use → less rocket use
        - Trade-off: slower colony establishment

  C. QUANTIFIED IMPACT:
     - Baseline (Rocket Only): {:.0f} Mt CO2
     - Best Case (Elevator + 100% Renewable): {:.0f} Mt CO2
     - Maximum Reduction: {:.1f}%
""".format(
        df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0],
        calculate_scenario_impact('Best', 186, 1.0, [], renewable_energy_fraction=1.0).co2_total_Mt,
        (1 - calculate_scenario_impact('Best', 186, 1.0, [], renewable_energy_fraction=1.0).co2_total_Mt / 
         df.loc[df['Scenario']=='Rocket Only', 'CO2 Total (Mt)'].values[0]) * 100
    ))
    
    report.append("\n" + "=" * 80)
    report.append("6. STRATOSPHERIC IMPACT ANALYSIS")
    report.append("=" * 80)
    report.append(f"""
  Stratospheric Water Vapor Injection (Rocket Only scenario):
  
  - Total H2O emitted: {df.loc[df['Scenario']=='Rocket Only', 'H2O Total (Mt)'].values[0]:.2f} Mt
  - Stratospheric fraction: {EMISSION_FACTORS.stratosphere_emission_fraction*100:.0f}%
  - Stratospheric H2O: {df.loc[df['Scenario']=='Rocket Only', 'H2O Stratosphere (Mt)'].values[0]:.2f} Mt
  - Residence time: {EMISSION_FACTORS.stratosphere_h2o_residence_time:.0f} years
  
  Potential Climate Effects:
  - Stratospheric H2O has greenhouse effect (radiative forcing)
  - May affect ozone chemistry
  - Elevator-only scenario eliminates this impact entirely
""")
    
    report.append("=" * 80)
    
    # 保存报告
    report_text = '\n'.join(report)
    with open(f'{save_dir}/environmental_report.txt', 'w', encoding='utf-8') as f:
        f.write(report_text)
    print(f"分析报告已保存至: {save_dir}/environmental_report.txt")
    
    # 保存CSV
    df.to_csv(f'{save_dir}/scenario_comparison.csv', index=False)
    print(f"数据表已保存至: {save_dir}/scenario_comparison.csv")
    
    return report_text


# ============== 主程序 ==============

if __name__ == "__main__":
    print("=" * 80)
    print("TASK 4: ENVIRONMENTAL IMPACT ANALYSIS")
    print("Moon Colony Construction with LOX/CH4 Fuel")
    print("=" * 80)
    
    save_dir = '/Volumes/Files/code/mm/20260130_b/p4'
    
    # 1. 基础情景对比
    print("\n1. Analyzing baseline scenarios...")
    df = plot_environmental_comparison(save_dir)
    print("\nScenario Comparison:")
    print(df.to_string(index=False))
    
    # 2. Pareto分析
    print("\n2. Generating Pareto analysis...")
    plot_pareto_with_environment(save_dir)
    
    # 3. 减排策略对比
    print("\n3. Analyzing mitigation strategies...")
    plot_mitigation_strategies(save_dir)
    
    # 4. 综合汇总
    print("\n4. Generating comprehensive summary...")
    plot_comprehensive_summary(save_dir)
    
    # 5. 生成报告
    print("\n5. Generating report...")
    generate_report(save_dir)
    
    print("\n" + "=" * 80)
    print("Analysis complete! All outputs saved to:", save_dir)
    print("=" * 80)