diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..5008ddf
Binary files /dev/null and b/.DS_Store differ
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..dc7c631
--- /dev/null
+++ b/README.md
@@ -0,0 +1,357 @@
+# Moon Colony Logistics（地月运输建模与对比）
+
+本目录包含针对 **MCM 2026 Problem B**（见 `prob/en_B.md`）的一组可运行脚本：用**简化的轨道力学 + 火箭方程**估算把 **1 亿吨（100 million metric tons）**物资送到月球轨道/基地的**时间线**与**能量消耗**，并对比：
+
+- **纯火箭方案**（10 个既有发射场，每站每日 1 发）
+- **纯太空电梯方案**（3 个电梯港口持续转运）
+- **混合方案**（电梯优先，剩余由火箭补齐，且随着工期变长更多任务转移到低纬发射场以降低纬度惩罚）
+
+> 本项目的“能量”口径：
+> - **火箭**：用推进剂化学能近似 \(E \approx m_\text{fuel}\cdot e_s\)（`e_s` 为燃料比能，默认 LH2/LOX 约 15.5 MJ/kg）。
+> - **太空电梯**：用“提升能 + 顶端火箭段能量”的**总比能量**近似（默认 157.2 MJ/kg，即 157.2 GJ/吨）。
+> 这是一种用于**对比量级与趋势**的工程近似，不等同于严格的热力学效率或完整任务能量收支。
+
+---
+
+### 目录结构与脚本索引
+
+#### 核心脚本
+
+- **`specific_energy_comparison.py`**
+ 地面发射 vs 太空电梯：比能量、燃料比、纬度影响（包含指定发射场列表）。
+- **`rocket_launch_scenario.py`**
+ 纯火箭方案（10 发射场 * 每日 1 发）排布与能量-工期曲线；并做发射频率敏感性分析。
+- **`combined_scenario.py`**
+ 混合方案（电梯优先 + 火箭补齐），输出组合方案能量-工期曲线，并生成 **100–300 年**区间三方案同图对比。
+- **`launch_capacity_analysis.py`**
+ （可选）对历史火箭年发射量做 Logistic / Gompertz / Richards 拟合与预测，生成 `launch_capacity_*.png` 与 `launch_capacity_results.csv`。
+
+#### 主要输出图表（均在本目录）
+
+| 图表文件 | 含义 |
+|---|---|
+| `specific_energy_comparison.png` | 地面 vs 电梯的比能量/燃料比对比，含 ΔV 敏感性 |
+| `latitude_effects.png` | 月球任务场景下：纬度→自转速度损失→燃料/能量增加趋势（并标注指定发射场） |
+| `rocket_launch_timeline.png` | 纯火箭方案：完成年限 vs 总能量（工期变长→更多低纬→能量略降） |
+| `rocket_comprehensive.png` | 纯火箭方案综合图：能量趋势 + 分配饼图等 |
+| `launch_frequency_analysis.png` | 发射频率敏感性：每站每日发射数提升→完工年限显著下降 |
+| `combined_scenario_analysis.png` | 混合方案：总能量/电梯能量/火箭能量随工期变化 + 电梯占比 + 需要启用的火箭站点数 |
+| `scenario_comparison.png` | elevator-only / rocket-only / combined 的完工时间与能量柱状对比 |
+| **`three_scenarios_comparison.png`** | **100–300 年区间**：纯火箭、纯电梯、混合方案的“年份-能量”折线同图对比 |
+| `launch_capacity_analysis.png` `launch_capacity_physical.png` `launch_capacity_predictions.png` | （可选）火箭发射能力上限/趋势拟合与预测 |
+
+#### 数据文件（如存在）
+
+- `launch_capacity_results.csv`
+- `rocket_launch_counts.csv`
+- `launch_sites.csv`
+- `rocket launch counts.txt`
+
+---
+
+### 图表总览（本目录所有 PNG）
+
+为便于写报告/插图，这里把当前目录下的所有 PNG 统一列出（按主题分组）。如果你只想看关键结论，优先看 **`three_scenarios_comparison.png`**。
+
+
+1) 比能量与纬度效应（地面 vs 电梯、发射场纬度）
+
+#### `specific_energy_comparison.png`
+
+
+
+#### `latitude_effects.png`
+
+
+
+
+
+
+2) 纯火箭方案（排布与敏感性）
+
+#### `rocket_launch_timeline.png`
+
+
+
+#### `rocket_comprehensive.png`
+
+
+
+#### `launch_distribution_min.png`
+
+
+
+#### `launch_distribution_mid.png`
+
+
+
+#### `launch_frequency_analysis.png`
+
+
+
+
+
+
+3) 混合方案与三方案对比（核心结论图）
+
+#### `combined_scenario_analysis.png`
+
+
+
+#### `scenario_comparison.png`
+
+
+
+#### `three_scenarios_comparison.png`（100–300 年折线同图对比）
+
+
+
+
+
+
+4)（可选）火箭发射能力上限/趋势拟合
+
+#### `launch_capacity_analysis.png`
+
+
+
+#### `launch_capacity_physical.png`
+
+
+
+#### `launch_capacity_predictions.png`
+
+
+
+
+
+
+5) 其它历史遗留图
+
+本目录还包含 `earth_moon_transfer_analysis.png`（早期“地月转移：载荷–燃料/能量关系”图）。当前目录下未发现对应的 `earth_moon_transfer.py` 源码文件；如果你需要我把该脚本补回并与现有口径统一，我可以再加一个文件。
+
+#### `earth_moon_transfer_analysis.png`
+
+
+
+
+
+### 依赖安装
+
+建议使用 Python 3.10+。
+
+```bash
+pip install numpy matplotlib pandas
+```
+
+如需运行 `launch_capacity_analysis.py`（含曲线拟合），还需要：
+
+```bash
+pip install scipy
+```
+
+---
+
+### 数学模型（简化版）概览
+
+#### 1）火箭方程：燃料-载荷的线性比例（在固定 ΔV 与 Isp 下）
+
+火箭方程：
+\[
+\Delta V = v_e \ln\left(\frac{m_0}{m_f}\right),\quad v_e = I_{sp} g_0
+\]
+
+当任务 ΔV、发动机 \(I_{sp}\) 固定时，质量比 \(R=\exp(\Delta V/v_e)\) 为常数。若进一步采用结构系数近似：
+\[
+m_s=\alpha m_\text{fuel}
+\]
+则可解得燃料与载荷满足：
+\[
+\boxed{m_\text{fuel}=k\,m_p}\quad(\text{其中 }k\text{ 为常数，取决于 }\Delta V, I_{sp}, \alpha)
+\]
+
+这解释了为什么在同一任务与发动机参数下，“载荷–燃料”是近似**线性**关系。
+
+#### 2）能量口径（用于对比）
+
+- 火箭推进剂能量（近似）：
+\[
+E_\text{rocket}\approx m_\text{fuel}\,e_s
+\]
+其中 \(e_s\) 为燃料比能（默认 15.5 MJ/kg）。
+
+- 太空电梯总比能量（本项目用常数近似）：
+\[
+E_\text{elevator}\approx \epsilon_\text{elevator}\,m_p
+\]
+默认 \(\epsilon_\text{elevator}=157.2\text{ MJ/kg}\)（可在脚本中改）。
+
+#### 3）纬度导致的自转“损失”近似
+
+地球表面自转线速度：
+\[
+v(\varphi)=\omega R_E\cos(\varphi)
+\]
+相对赤道的速度差近似为 ΔV 惩罚：
+\[
+\Delta v_\text{loss}=v(0^\circ)-v(\varphi)
+\]
+
+在月球任务场景下（允许进入与发射场纬度相匹配的倾角、避免昂贵的纯平面变轨），纬度主要通过该项影响火箭能耗。
+
+---
+
+### 如何运行（会生成图片）
+
+#### 1）比能量与纬度效应
+
+```bash
+python specific_energy_comparison.py
+```
+
+输出（示例）会包括不同发射场的 ΔV 损失、燃料比、能量（MJ/kg）等，并生成：
+
+- `specific_energy_comparison.png`
+- `latitude_effects.png`
+
+#### 2）纯火箭方案：工期-能量与排布
+
+```bash
+python rocket_launch_scenario.py
+```
+
+会生成：
+
+- `rocket_launch_timeline.png`
+- `rocket_comprehensive.png`
+- `launch_distribution_min.png`、`launch_distribution_mid.png`
+- `launch_frequency_analysis.png`
+
+**部分运行输出（示例）**：
+
+```text
+Total payload: 100 million metric tons
+Payload per launch: 125 metric tons
+Total launches needed: 800,000
+Number of launch sites: 10
+Minimum completion time: 219.2 years
+
+Energy reduction (Extended vs Minimum time): 2.4%
+```
+
+解释：
+- 由于“每站每日 1 发、载荷 125 吨”的硬约束，纯火箭最短也要约 **219 年**完成。
+- 工期拉长后，可把更多发射集中到低纬站点，**能量仅小幅下降（约几个百分点）**——纬度惩罚相对火箭整体能耗并不占主导。
+
+#### 3）混合方案 + 100–300 年三方案同图对比
+
+```bash
+python combined_scenario.py
+```
+
+会生成：
+
+- `combined_scenario_analysis.png`
+- `scenario_comparison.png`
+- **`three_scenarios_comparison.png`（100–300 年折线同图对比）**
+
+**部分运行输出（示例）**：
+
+```text
+Elevator-only completion: 186.2 years
+Minimum completion (all systems): 100.7 years
+
+1. Minimum Time (~101 years): Total energy: 31537 PJ
+4. Elevator Only (186 years): Total energy: 15720 PJ
+Rocket Only (219 years): Total energy: 50609 PJ
+```
+
+解释：
+- **混合方案**把最短工期从电梯单独的 ~186 年压缩到 ~101 年（靠火箭补齐），但能量高于纯电梯。
+- 当工期 ≥ 186 年时，电梯即可独立完成，混合方案会自动退化为纯电梯曲线（能耗趋于同一水平）。
+- 纯火箭能耗显著更高（本模型下约为纯电梯的 **3×** 量级）。
+
+---
+
+### 发射场（10 个）与纬度效应（月球任务场景）
+
+本项目按纬度从低到高优先分配火箭发射（“低纬优先”贪心策略），站点列表与纬度如下（与题目一致）：
+
+| 发射场 | 近似纬度 | 备注 |
+|---|---:|---|
+| French Guiana（Kourou） | 5.2° | 低纬优势明显 |
+| Satish Dhawan（India） | 13.7° | 低纬 |
+| Texas（Boca Chica） | 26.0° | 美国 |
+| Florida（Cape Canaveral） | 28.5° | 美国 |
+| California（Vandenberg） | 34.7° | 美国 |
+| Virginia（Wallops） | 37.8° | 美国 |
+| Taiyuan（China） | 38.8° | 中国 |
+| Mahia（New Zealand） | 39.3° | 新西兰（南半球，按 \(|\varphi|\) 处理） |
+| Kazakhstan（Baikonur） | 45.6° | 哈萨克斯坦 |
+| Alaska（Kodiak） | 57.4° | 高纬，惩罚最大 |
+
+对应曲线可在下图直观看到：
+
+
+
+---
+
+### 三方案“年份–能量”折线对比（100–300 年）
+
+下图是你要求的 **100–300 年**区间三条折线同图（不可完成的年份用空缺表示）：
+
+
+
+读图要点：
+- **100–186 年**：只有混合方案可完成（电梯能力不足、火箭也不够快）。
+- **186–219 年**：电梯可完成；火箭仍无法完成；混合与电梯曲线重合（优先用电梯）。
+- **≥219 年**：三种方案都可完成，但纯火箭能耗最高。
+
+为了写报告时更快引用，下面给出几组关键年份的能量读数（与脚本运行输出一致，单位 PJ）：
+
+| 年份 | 混合方案 | 纯电梯 | 纯火箭 |
+|---:|---:|---:|---:|
+| 110 | 29852 | N/A | N/A |
+| 150 | 22257 | N/A | N/A |
+| 186 | 15720 | 15720 | N/A |
+| 220 | 15720 | 15720 | 50605 |
+| 300 | 15720 | 15720 | 50187 |
+
+---
+
+### 可调参数（最常改的地方）
+
+#### 火箭相关（`rocket_launch_scenario.py` / `combined_scenario.py`）
+
+- `TOTAL_PAYLOAD`：总物资（吨）
+- `PAYLOAD_PER_LAUNCH`：单次发射载荷（吨，默认 125）
+- `ISP`、`ALPHA`、`NUM_STAGES`：推进系统与结构参数
+- `SPECIFIC_FUEL_ENERGY`：燃料比能（J/kg）
+- `DELTA_V_BASE`：赤道基准 ΔV（m/s）
+- `LAUNCH_SITES`：发射场纬度列表（可自行替换/增减）
+
+#### 太空电梯相关（`combined_scenario.py`）
+
+- `NUM_ELEVATORS`：电梯数量（默认 3）
+- `ELEVATOR_CAPACITY_PER_YEAR`：**每个电梯年运力（吨/年）**
+ > 题面是 **179,000 t/year**；若你要用 **17,900 t/year**，直接改此值即可。
+- `ELEVATOR_SPECIFIC_ENERGY`：电梯总比能量（J/吨），默认 157.2 GJ/吨
+
+---
+
+### 重要假设与局限（写报告时建议明确）
+
+- **脉冲变轨/固定 ΔV**：将复杂轨道机动用固定 ΔV 预算概括。
+- **恒定比冲/结构系数**：忽略发动机工况变化、分级差异、载具设计差异。
+- **“能量”是工程近似口径**：火箭取推进剂化学能；电梯取提升+顶端段合并比能量。
+- **完美运行**：未考虑失败率、维护停机、窗口约束、天气/法规等运营因素。
+- **排布算法为贪心**：火箭发射按“低纬优先”分配以降低纬度惩罚；未做更复杂的成本/容量优化（可扩展为线性规划/整数规划）。
+
+---
+
+### 一句话结论（便于写 Summary）
+
+- **纯火箭**在“10 站 * 每日 1 发 * 125 吨/发”的约束下，最短约 **219 年**，且能耗约 **5×10⁴ PJ** 量级。
+- **纯电梯**能耗最低（约 **1.57×10⁴ PJ**），但最短约 **186 年**。
+- **混合方案**能把工期压到约 **101 年**，能耗介于两者之间，并随着工期增加逐步趋近纯电梯能耗。
+
diff --git a/combined_scenario.py b/combined_scenario.py
new file mode 100644
index 0000000..9b2952c
--- /dev/null
+++ b/combined_scenario.py
@@ -0,0 +1,705 @@
+"""
+Moon Colony Construction - Combined Scenario Analysis
+(Space Elevator + Traditional Rockets)
+
+Problem: Deliver 100 million metric tons to Moon
+- 3 Space Elevators: 179,000 metric tons/year each (priority)
+- 10 Rocket launch sites: max 1 launch/day each, 125 tons/launch
+- Optimize: Use elevators first (lower specific energy), then rockets
+
+Key: As timeline extends:
+1. More payload handled by elevators (lower energy)
+2. Remaining rocket launches shift to lower-latitude sites
+"""
+
+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, Tuple, Dict
+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
+
+# 太空电梯能量参数 (从之前分析)
+# 电梯提升能量: ~88.7 MJ/kg
+# 顶端推进能量: ~68.5 MJ/kg
+# 总计: ~157.2 MJ/kg = 157.2 TJ/1000 tons
+ELEVATOR_SPECIFIC_ENERGY = 157.2e9 # J per metric ton (157.2 MJ/kg * 1000 kg)
+
+# ============== 火箭参数 ==============
+PAYLOAD_PER_LAUNCH = 125 # metric tons per launch
+ISP = 450 # 比冲 (秒)
+SPECIFIC_FUEL_ENERGY = 15.5e6 # J/kg
+ALPHA = 0.10 # 结构系数
+NUM_STAGES = 3
+DELTA_V_BASE = 13300 # m/s (赤道发射到月球)
+
+# 火箭比能量 (从之前分析: ~492.7 MJ/kg = 492.7 TJ/1000 tons for equator)
+# 实际会因纬度变化
+
+
+# ============== 发射场定义 ==============
+@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
+
+ @property
+ def total_delta_v(self) -> float:
+ return DELTA_V_BASE + self.delta_v_loss
+
+
+# 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)
+
+
+# ============== 核心计算函数 ==============
+
+def fuel_ratio_multistage(delta_v: float) -> float:
+ """多级火箭燃料/载荷比"""
+ ve = ISP * G0
+ 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_energy_per_ton(site: LaunchSite) -> float:
+ """火箭发射每吨载荷的能量消耗 (J/ton)"""
+ k = fuel_ratio_multistage(site.total_delta_v)
+ fuel_per_ton = k * 1000 # kg fuel per metric ton payload
+ return fuel_per_ton * SPECIFIC_FUEL_ENERGY
+
+
+def rocket_energy_per_launch(site: LaunchSite) -> float:
+ """单次火箭发射的能量消耗 (J)"""
+ return rocket_energy_per_ton(site) * PAYLOAD_PER_LAUNCH
+
+
+# ============== 组合方案计算 ==============
+
+def calculate_combined_scenario(completion_years: float) -> Dict:
+ """
+ 计算给定完成年限下的组合方案
+
+ 策略:
+ 1. 太空电梯满负荷运行
+ 2. 剩余载荷由火箭发射，优先低纬度站点
+
+ 返回:
+ 方案详情字典
+ """
+ # 太空电梯运输量
+ elevator_payload = min(TOTAL_ELEVATOR_CAPACITY * completion_years, TOTAL_PAYLOAD)
+ elevator_energy = elevator_payload * ELEVATOR_SPECIFIC_ENERGY
+
+ # 剩余需要火箭运输的载荷
+ remaining_payload = TOTAL_PAYLOAD - elevator_payload
+
+ if remaining_payload <= 0:
+ # 完全由电梯完成
+ return {
+ 'completion_years': completion_years,
+ 'elevator_payload': elevator_payload,
+ 'elevator_energy_PJ': elevator_energy / 1e15,
+ 'rocket_payload': 0,
+ 'rocket_energy_PJ': 0,
+ 'total_energy_PJ': elevator_energy / 1e15,
+ 'rocket_launches': 0,
+ 'sites_used': 0,
+ 'elevator_fraction': 1.0,
+ 'rocket_distribution': [],
+ }
+
+ # 需要火箭发射
+ rocket_launches_needed = int(np.ceil(remaining_payload / PAYLOAD_PER_LAUNCH))
+
+ # 按纬度优先分配火箭发射
+ days_available = completion_years * 365
+ max_launches_per_site = int(days_available)
+
+ rocket_distribution = []
+ remaining_launches = rocket_launches_needed
+ rocket_energy = 0
+ sites_used = 0
+
+ for site in LAUNCH_SITES:
+ if remaining_launches <= 0:
+ rocket_distribution.append((site, 0))
+ else:
+ allocated = min(remaining_launches, max_launches_per_site)
+ rocket_distribution.append((site, allocated))
+ rocket_energy += rocket_energy_per_launch(site) * allocated
+ remaining_launches -= allocated
+ if allocated > 0:
+ sites_used += 1
+
+ # 检查是否能在规定时间内完成
+ total_rocket_capacity = len(LAUNCH_SITES) * max_launches_per_site * PAYLOAD_PER_LAUNCH
+ if remaining_payload > total_rocket_capacity:
+ # 无法完成
+ return None
+
+ rocket_payload = rocket_launches_needed * PAYLOAD_PER_LAUNCH
+ total_energy = elevator_energy + rocket_energy
+
+ return {
+ 'completion_years': completion_years,
+ 'elevator_payload': elevator_payload,
+ 'elevator_energy_PJ': elevator_energy / 1e15,
+ 'rocket_payload': rocket_payload,
+ 'rocket_energy_PJ': rocket_energy / 1e15,
+ 'total_energy_PJ': total_energy / 1e15,
+ 'rocket_launches': rocket_launches_needed,
+ 'sites_used': sites_used,
+ 'elevator_fraction': elevator_payload / TOTAL_PAYLOAD,
+ 'rocket_distribution': rocket_distribution,
+ }
+
+
+def analyze_combined_timeline() -> pd.DataFrame:
+ """
+ 分析不同完成年限下的组合方案
+ """
+ # 计算关键时间点
+ # 1. 仅电梯完成所需时间
+ elevator_only_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY
+
+ # 2. 最短时间（电梯+全部火箭站点满负荷）
+ rocket_capacity_per_year = len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH
+ total_capacity_per_year = TOTAL_ELEVATOR_CAPACITY + rocket_capacity_per_year
+ min_years = TOTAL_PAYLOAD / total_capacity_per_year
+
+ print(f"Elevator-only completion: {elevator_only_years:.1f} years")
+ print(f"Minimum completion (all systems): {min_years:.1f} years")
+ print(f"Elevator capacity: {TOTAL_ELEVATOR_CAPACITY:,} tons/year")
+ print(f"Rocket capacity: {rocket_capacity_per_year:,} tons/year")
+
+ # 分析范围
+ years_range = np.linspace(min_years, elevator_only_years * 1.1, 200)
+
+ results = []
+ for years in years_range:
+ scenario = calculate_combined_scenario(years)
+ if scenario is not None:
+ results.append(scenario)
+
+ return pd.DataFrame(results)
+
+
+# ============== 可视化 ==============
+
+def plot_combined_analysis(
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/combined_scenario_analysis.png'
+):
+ """
+ 绘制组合方案分析图
+ """
+ df = analyze_combined_timeline()
+
+ # 关键时间点
+ elevator_only_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY
+ rocket_capacity_per_year = len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH
+ total_capacity_per_year = TOTAL_ELEVATOR_CAPACITY + rocket_capacity_per_year
+ min_years = TOTAL_PAYLOAD / total_capacity_per_year
+
+ fig, axes = plt.subplots(2, 2, figsize=(16, 14))
+
+ # ========== 图1: 完成年份 vs 总能量消耗 ==========
+ ax1 = axes[0, 0]
+
+ ax1.plot(df['completion_years'], df['total_energy_PJ'], 'b-', linewidth=2.5,
+ label='Total Energy')
+ ax1.plot(df['completion_years'], df['elevator_energy_PJ'], 'g--', linewidth=2,
+ label='Elevator Energy')
+ ax1.plot(df['completion_years'], df['rocket_energy_PJ'], 'r--', linewidth=2,
+ label='Rocket Energy')
+
+ # 标记关键点
+ ax1.axvline(x=min_years, color='orange', linestyle=':', linewidth=1.5,
+ label=f'Min time: {min_years:.0f} years')
+ ax1.axvline(x=elevator_only_years, color='purple', linestyle=':', linewidth=1.5,
+ label=f'Elevator-only: {elevator_only_years:.0f} years')
+
+ ax1.set_xlabel('Completion Time (years)', fontsize=12)
+ ax1.set_ylabel('Energy Consumption (PJ)', fontsize=12)
+ ax1.set_title('Combined Scenario: Completion Time vs Energy\n'
+ '(Space Elevator + Rockets)', fontsize=14)
+ ax1.legend(loc='upper right')
+ ax1.grid(True, alpha=0.3)
+ ax1.set_xlim(min_years * 0.95, elevator_only_years * 1.05)
+
+ # ========== 图2: 电梯 vs 火箭载荷比例 ==========
+ ax2 = axes[0, 1]
+
+ ax2.fill_between(df['completion_years'], 0, df['elevator_fraction'] * 100,
+ alpha=0.7, color='green', label='Space Elevator')
+ ax2.fill_between(df['completion_years'], df['elevator_fraction'] * 100, 100,
+ alpha=0.7, color='red', label='Rockets')
+
+ ax2.axvline(x=min_years, color='orange', linestyle=':', linewidth=1.5)
+ ax2.axvline(x=elevator_only_years, color='purple', linestyle=':', linewidth=1.5)
+
+ ax2.set_xlabel('Completion Time (years)', fontsize=12)
+ ax2.set_ylabel('Payload Share (%)', fontsize=12)
+ ax2.set_title('Payload Distribution: Elevator vs Rockets', fontsize=14)
+ ax2.legend(loc='center right')
+ ax2.set_ylim(0, 100)
+ ax2.set_xlim(min_years * 0.95, elevator_only_years * 1.05)
+ ax2.grid(True, alpha=0.3)
+
+ # ========== 图3: 使用的火箭发射场数量 ==========
+ ax3 = axes[1, 0]
+
+ ax3.plot(df['completion_years'], df['sites_used'], 'r-', linewidth=2)
+ ax3.fill_between(df['completion_years'], 0, df['sites_used'], alpha=0.3, color='red')
+
+ ax3.axvline(x=min_years, color='orange', linestyle=':', linewidth=1.5)
+ ax3.axvline(x=elevator_only_years, color='purple', linestyle=':', linewidth=1.5)
+
+ # 标注发射场
+ for i, site in enumerate(LAUNCH_SITES):
+ ax3.axhline(y=i+1, color='gray', linestyle=':', alpha=0.2)
+ ax3.text(df['completion_years'].max() * 0.99, i+1 + 0.15,
+ f'{site.short_name} ({site.abs_latitude:.0f}°)',
+ fontsize=8, ha='right', va='bottom')
+
+ ax3.set_xlabel('Completion Time (years)', fontsize=12)
+ ax3.set_ylabel('Number of Rocket Sites Used', fontsize=12)
+ ax3.set_title('Rocket Launch Sites Required\n(Longer timeline → fewer high-lat sites)', fontsize=14)
+ ax3.set_ylim(0, 11)
+ ax3.set_xlim(min_years * 0.95, elevator_only_years * 1.05)
+ ax3.grid(True, alpha=0.3)
+
+ # ========== 图4: 能量效率对比 ==========
+ ax4 = axes[1, 1]
+
+ # 计算平均比能量
+ avg_specific_energy = df['total_energy_PJ'] * 1e15 / TOTAL_PAYLOAD / 1e9 # GJ/ton
+
+ ax4.plot(df['completion_years'], avg_specific_energy, 'b-', linewidth=2)
+
+ # 参考线
+ elevator_specific = ELEVATOR_SPECIFIC_ENERGY / 1e9 # GJ/ton
+ ax4.axhline(y=elevator_specific, color='green', linestyle='--',
+ label=f'Elevator only: {elevator_specific:.1f} GJ/ton')
+
+ # 平均火箭比能量 (赤道)
+ rocket_specific_equator = rocket_energy_per_ton(LAUNCH_SITES[0]) / 1e9
+ ax4.axhline(y=rocket_specific_equator, color='red', linestyle='--',
+ label=f'Rocket (equator): {rocket_specific_equator:.1f} GJ/ton')
+
+ ax4.axvline(x=min_years, color='orange', linestyle=':', linewidth=1.5)
+ ax4.axvline(x=elevator_only_years, color='purple', linestyle=':', linewidth=1.5)
+
+ ax4.set_xlabel('Completion Time (years)', fontsize=12)
+ ax4.set_ylabel('Average Specific Energy (GJ/ton)', fontsize=12)
+ ax4.set_title('Energy Efficiency vs Timeline\n(Longer time → more elevator → lower energy)', fontsize=14)
+ ax4.legend()
+ ax4.set_xlim(min_years * 0.95, elevator_only_years * 1.05)
+ ax4.grid(True, alpha=0.3)
+
+ plt.suptitle(f'Combined Scenario: 100M tons to Moon\n'
+ f'3 Elevators ({TOTAL_ELEVATOR_CAPACITY:,} t/yr) + 10 Rocket Sites',
+ fontsize=15, y=1.02)
+
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"\n组合方案分析图已保存至: {save_path}")
+
+ return df
+
+
+def plot_scenario_comparison(
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/scenario_comparison.png'
+):
+ """
+ 对比三种方案：仅电梯、仅火箭、组合
+ """
+ # 关键时间点分析
+ elevator_only_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY
+ rocket_capacity_per_year = len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH
+ rocket_only_min_years = TOTAL_PAYLOAD / rocket_capacity_per_year
+ combined_min_years = TOTAL_PAYLOAD / (TOTAL_ELEVATOR_CAPACITY + rocket_capacity_per_year)
+
+ # 各方案能量
+ elevator_only_energy = TOTAL_PAYLOAD * ELEVATOR_SPECIFIC_ENERGY / 1e15 # PJ
+
+ # 火箭方案 (平均所有站点)
+ avg_rocket_energy_per_ton = sum(rocket_energy_per_ton(s) for s in LAUNCH_SITES) / len(LAUNCH_SITES)
+ rocket_only_energy = TOTAL_PAYLOAD * avg_rocket_energy_per_ton / 1e15 # PJ
+
+ # 组合方案 (最短时间 - 加余量)
+ combined_scenario = calculate_combined_scenario(combined_min_years * 1.01)
+ combined_min_energy = combined_scenario['total_energy_PJ'] if combined_scenario else elevator_only_energy
+
+ # 组合方案 (仅电梯时间)
+ combined_extended = calculate_combined_scenario(elevator_only_years * 0.95)
+ combined_extended_energy = combined_extended['total_energy_PJ'] if combined_extended else elevator_only_energy
+
+ fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+
+ # ========== 图1: 完成时间对比 ==========
+ ax1 = axes[0]
+
+ scenarios = ['Elevator\nOnly', 'Rocket\nOnly', 'Combined\n(Min Time)', 'Combined\n(Extended)']
+ times = [elevator_only_years, rocket_only_min_years, combined_min_years, elevator_only_years * 0.99]
+ colors = ['green', 'red', 'blue', 'purple']
+
+ bars = ax1.bar(scenarios, times, color=colors, alpha=0.7)
+
+ ax1.set_ylabel('Completion Time (years)', fontsize=12)
+ ax1.set_title('Completion Time Comparison', fontsize=14)
+ ax1.grid(True, alpha=0.3, axis='y')
+
+ for bar, t in zip(bars, times):
+ ax1.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 2,
+ f'{t:.0f}y', ha='center', fontsize=11, fontweight='bold')
+
+ # ========== 图2: 能量消耗对比 ==========
+ ax2 = axes[1]
+
+ energies = [elevator_only_energy, rocket_only_energy, combined_min_energy, combined_extended_energy]
+
+ bars = ax2.bar(scenarios, energies, color=colors, alpha=0.7)
+
+ ax2.set_ylabel('Total Energy (PJ)', fontsize=12)
+ ax2.set_title('Energy Consumption Comparison', fontsize=14)
+ ax2.grid(True, alpha=0.3, axis='y')
+
+ for bar, e in zip(bars, energies):
+ ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 500,
+ f'{e:.0f}', ha='center', fontsize=11, fontweight='bold')
+
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"方案对比图已保存至: {save_path}")
+
+ return {
+ 'elevator_only': {'years': elevator_only_years, 'energy_PJ': elevator_only_energy},
+ 'rocket_only': {'years': rocket_only_min_years, 'energy_PJ': rocket_only_energy},
+ 'combined_min': {'years': combined_min_years, 'energy_PJ': combined_min_energy},
+ 'combined_extended': {'years': elevator_only_years * 0.99, 'energy_PJ': combined_extended_energy},
+ }
+
+
+def print_summary(df: pd.DataFrame):
+ """打印分析摘要"""
+ elevator_only_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY
+ rocket_capacity_per_year = len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH
+ combined_min_years = TOTAL_PAYLOAD / (TOTAL_ELEVATOR_CAPACITY + rocket_capacity_per_year)
+
+ print("\n" + "=" * 90)
+ print("COMBINED SCENARIO (SPACE ELEVATOR + ROCKETS) ANALYSIS")
+ print("=" * 90)
+
+ print(f"\nSystem Capacities:")
+ print(f" - 3 Space Elevators: {TOTAL_ELEVATOR_CAPACITY:,} tons/year")
+ print(f" - 10 Rocket Sites: {rocket_capacity_per_year:,} tons/year (@ 1 launch/day each)")
+ print(f" - Combined: {TOTAL_ELEVATOR_CAPACITY + rocket_capacity_per_year:,} tons/year")
+
+ print(f"\nSpecific Energy:")
+ print(f" - Elevator: {ELEVATOR_SPECIFIC_ENERGY/1e9:.1f} GJ/ton ({ELEVATOR_SPECIFIC_ENERGY/1e9/rocket_energy_per_ton(LAUNCH_SITES[0])*1e9*100:.1f}% of rocket)")
+ print(f" - Rocket (Kourou): {rocket_energy_per_ton(LAUNCH_SITES[0])/1e9:.1f} GJ/ton")
+ print(f" - Rocket (Alaska): {rocket_energy_per_ton(LAUNCH_SITES[-1])/1e9:.1f} GJ/ton")
+
+ print(f"\nKey Scenarios:")
+
+ # 最短时间 (加一点余量确保可以完成)
+ min_scenario = calculate_combined_scenario(combined_min_years * 1.01)
+ print(f"\n 1. Minimum Time (~{combined_min_years:.0f} years):")
+ print(f" - Elevator payload: {min_scenario['elevator_payload']/1e6:.1f}M tons ({min_scenario['elevator_fraction']*100:.1f}%)")
+ print(f" - Rocket payload: {min_scenario['rocket_payload']/1e6:.1f}M tons")
+ print(f" - Rocket launches: {min_scenario['rocket_launches']:,}")
+ print(f" - Sites used: {min_scenario['sites_used']}")
+ print(f" - Total energy: {min_scenario['total_energy_PJ']:.0f} PJ")
+
+ # 中等时间
+ mid_years = (combined_min_years + elevator_only_years) / 2
+ mid_scenario = calculate_combined_scenario(mid_years)
+ print(f"\n 2. Medium Timeline ({mid_years:.0f} years):")
+ print(f" - Elevator payload: {mid_scenario['elevator_payload']/1e6:.1f}M tons ({mid_scenario['elevator_fraction']*100:.1f}%)")
+ print(f" - Rocket payload: {mid_scenario['rocket_payload']/1e6:.1f}M tons")
+ print(f" - Rocket launches: {mid_scenario['rocket_launches']:,}")
+ print(f" - Sites used: {mid_scenario['sites_used']}")
+ print(f" - Total energy: {mid_scenario['total_energy_PJ']:.0f} PJ")
+
+ # 接近仅电梯
+ extended_years = elevator_only_years * 0.95
+ extended_scenario = calculate_combined_scenario(extended_years)
+ print(f"\n 3. Extended Timeline ({extended_years:.0f} years):")
+ print(f" - Elevator payload: {extended_scenario['elevator_payload']/1e6:.1f}M tons ({extended_scenario['elevator_fraction']*100:.1f}%)")
+ print(f" - Rocket payload: {extended_scenario['rocket_payload']/1e6:.1f}M tons")
+ print(f" - Rocket launches: {extended_scenario['rocket_launches']:,}")
+ print(f" - Sites used: {extended_scenario['sites_used']}")
+ print(f" - Total energy: {extended_scenario['total_energy_PJ']:.0f} PJ")
+
+ # 仅电梯
+ print(f"\n 4. Elevator Only ({elevator_only_years:.0f} years):")
+ elevator_energy = TOTAL_PAYLOAD * ELEVATOR_SPECIFIC_ENERGY / 1e15
+ print(f" - Elevator payload: 100M tons (100%)")
+ print(f" - Rocket launches: 0")
+ print(f" - Total energy: {elevator_energy:.0f} PJ")
+
+ print("\n" + "=" * 90)
+ print("ENERGY SAVINGS ANALYSIS")
+ print("=" * 90)
+
+ rocket_only_energy = TOTAL_PAYLOAD * sum(rocket_energy_per_ton(s) for s in LAUNCH_SITES) / len(LAUNCH_SITES) / 1e15
+
+ print(f"\n Baseline (Rocket Only): {rocket_only_energy:.0f} PJ")
+ print(f" Combined (Min Time): {min_scenario['total_energy_PJ']:.0f} PJ (saves {(1-min_scenario['total_energy_PJ']/rocket_only_energy)*100:.0f}%)")
+ print(f" Combined (Extended): {extended_scenario['total_energy_PJ']:.0f} PJ (saves {(1-extended_scenario['total_energy_PJ']/rocket_only_energy)*100:.0f}%)")
+ print(f" Elevator Only: {elevator_energy:.0f} PJ (saves {(1-elevator_energy/rocket_only_energy)*100:.0f}%)")
+
+ print("=" * 90)
+
+
+def plot_three_scenarios_comparison(
+ year_min: float = 100,
+ year_max: float = 300,
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/three_scenarios_comparison.png'
+):
+ """
+ 绘制三种方案在指定年份范围内的能量消耗对比折线图
+
+ 方案:
+ 1. 纯火箭方案 (Rocket Only)
+ 2. 纯电梯方案 (Elevator Only)
+ 3. 混合方案 (Combined)
+ """
+ years = np.linspace(year_min, year_max, 200)
+
+ # 计算各方案的能量
+ rocket_energy = []
+ elevator_energy = []
+ combined_energy = []
+
+ # 关键时间点
+ elevator_only_min_years = TOTAL_PAYLOAD / TOTAL_ELEVATOR_CAPACITY # ~186 years
+ rocket_only_min_years = TOTAL_PAYLOAD / (len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH) # ~219 years
+ combined_min_years = TOTAL_PAYLOAD / (TOTAL_ELEVATOR_CAPACITY + len(LAUNCH_SITES) * 365 * PAYLOAD_PER_LAUNCH) # ~101 years
+
+ for y in years:
+ # ===== 纯火箭方案 =====
+ if y >= rocket_only_min_years:
+ # 可以完成，优先使用低纬度站点
+ total_launches = int(TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH)
+ days_available = y * 365
+ max_per_site = int(days_available)
+
+ energy = 0
+ remaining = total_launches
+ for site in LAUNCH_SITES:
+ if remaining <= 0:
+ break
+ allocated = min(remaining, max_per_site)
+ energy += rocket_energy_per_launch(site) * allocated
+ remaining -= allocated
+
+ rocket_energy.append(energy / 1e15)
+ else:
+ # 无法完成
+ rocket_energy.append(np.nan)
+
+ # ===== 纯电梯方案 =====
+ if y >= elevator_only_min_years:
+ # 可以完成
+ energy = TOTAL_PAYLOAD * ELEVATOR_SPECIFIC_ENERGY
+ elevator_energy.append(energy / 1e15)
+ else:
+ # 无法完成
+ elevator_energy.append(np.nan)
+
+ # ===== 混合方案 =====
+ if y >= combined_min_years:
+ scenario = calculate_combined_scenario(y)
+ if scenario:
+ combined_energy.append(scenario['total_energy_PJ'])
+ else:
+ combined_energy.append(np.nan)
+ else:
+ combined_energy.append(np.nan)
+
+ # 创建图表
+ fig, ax = plt.subplots(figsize=(14, 8))
+
+ # 绘制三条折线
+ ax.plot(years, rocket_energy, 'r-', linewidth=2.5, label='Rocket Only', marker='', markersize=4)
+ ax.plot(years, elevator_energy, 'g-', linewidth=2.5, label='Elevator Only', marker='', markersize=4)
+ ax.plot(years, combined_energy, 'b-', linewidth=2.5, label='Combined (Elevator + Rocket)', marker='', markersize=4)
+
+ # 标记关键时间点
+ ax.axvline(x=combined_min_years, color='blue', linestyle=':', alpha=0.7, linewidth=1.5)
+ ax.axvline(x=elevator_only_min_years, color='green', linestyle=':', alpha=0.7, linewidth=1.5)
+ ax.axvline(x=rocket_only_min_years, color='red', linestyle=':', alpha=0.7, linewidth=1.5)
+
+ # 添加标注
+ y_pos = ax.get_ylim()[1] * 0.95
+ ax.annotate(f'Combined\nMin: {combined_min_years:.0f}y',
+ xy=(combined_min_years, y_pos), fontsize=9, ha='center',
+ bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
+ ax.annotate(f'Elevator\nMin: {elevator_only_min_years:.0f}y',
+ xy=(elevator_only_min_years, y_pos * 0.85), fontsize=9, ha='center',
+ bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
+ ax.annotate(f'Rocket\nMin: {rocket_only_min_years:.0f}y',
+ xy=(rocket_only_min_years, y_pos * 0.75), fontsize=9, ha='center',
+ bbox=dict(boxstyle='round', facecolor='lightsalmon', alpha=0.8))
+
+ # 填充区域表示节能效果
+ # Combined比Rocket节省的能量
+ rocket_arr = np.array(rocket_energy)
+ combined_arr = np.array(combined_energy)
+ valid_mask = ~np.isnan(rocket_arr) & ~np.isnan(combined_arr)
+ ax.fill_between(years[valid_mask], combined_arr[valid_mask], rocket_arr[valid_mask],
+ alpha=0.2, color='purple', label='Energy saved by Combined vs Rocket')
+
+ ax.set_xlabel('Completion Time (years)', fontsize=13)
+ ax.set_ylabel('Total Energy Consumption (PJ)', fontsize=13)
+ ax.set_title('Moon Colony Construction: Energy vs Completion Time\n'
+ '(100 Million Metric Tons Payload Comparison)', fontsize=15)
+ ax.legend(loc='upper right', fontsize=11)
+ ax.grid(True, alpha=0.3)
+ ax.set_xlim(year_min, year_max)
+
+ # 添加数据表格
+ # 选择几个关键年份的数据
+ key_years = [110, 150, 186, 220, 250]
+ table_data = []
+ for ky in key_years:
+ idx = np.argmin(np.abs(years - ky))
+ table_data.append([
+ f'{ky}',
+ f'{combined_energy[idx]:.0f}' if not np.isnan(combined_energy[idx]) else 'N/A',
+ f'{elevator_energy[idx]:.0f}' if not np.isnan(elevator_energy[idx]) else 'N/A',
+ f'{rocket_energy[idx]:.0f}' if not np.isnan(rocket_energy[idx]) else 'N/A',
+ ])
+
+ # 在图表下方添加表格
+ table = ax.table(
+ cellText=table_data,
+ colLabels=['Year', 'Combined (PJ)', 'Elevator (PJ)', 'Rocket (PJ)'],
+ loc='bottom',
+ bbox=[0.15, -0.25, 0.7, 0.18],
+ cellLoc='center'
+ )
+ table.auto_set_font_size(False)
+ table.set_fontsize(10)
+ table.scale(1, 1.5)
+
+ # 调整布局给表格留空间
+ plt.subplots_adjust(bottom=0.25)
+
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"\n三方案对比图已保存至: {save_path}")
+
+ # 打印关键数据
+ print("\n" + "=" * 70)
+ print("THREE SCENARIOS ENERGY COMPARISON (100-300 years)")
+ print("=" * 70)
+ print(f"\n{'Year':<10} {'Combined':>15} {'Elevator':>15} {'Rocket':>15}")
+ print(f"{'':10} {'(PJ)':>15} {'(PJ)':>15} {'(PJ)':>15}")
+ print("-" * 55)
+
+ for ky in [110, 130, 150, 186, 200, 220, 250, 300]:
+ idx = np.argmin(np.abs(years - ky))
+ c = f'{combined_energy[idx]:.0f}' if not np.isnan(combined_energy[idx]) else 'N/A'
+ e = f'{elevator_energy[idx]:.0f}' if not np.isnan(elevator_energy[idx]) else 'N/A'
+ r = f'{rocket_energy[idx]:.0f}' if not np.isnan(rocket_energy[idx]) else 'N/A'
+ print(f"{ky:<10} {c:>15} {e:>15} {r:>15}")
+
+ print("=" * 70)
+
+ return fig
+
+
+# ============== 主程序 ==============
+
+if __name__ == "__main__":
+ print("=" * 90)
+ print("MOON COLONY: COMBINED SCENARIO ANALYSIS")
+ print("(Space Elevator System + Traditional Rockets)")
+ print("=" * 90)
+
+ # 主分析
+ df = plot_combined_analysis()
+
+ # 打印摘要
+ print_summary(df)
+
+ # 方案对比
+ print("\n\n正在生成方案对比图...")
+ comparison = plot_scenario_comparison()
+
+ print("\n" + "=" * 90)
+ print("FINAL COMPARISON")
+ print("=" * 90)
+ print(f"\n{'Scenario':<25} {'Years':>12} {'Energy (PJ)':>15} {'vs Rocket Only':>18}")
+ print("-" * 70)
+ for name, data in comparison.items():
+ rocket_energy_val = comparison['rocket_only']['energy_PJ']
+ savings = (1 - data['energy_PJ'] / rocket_energy_val) * 100 if rocket_energy_val > 0 else 0
+ print(f"{name:<25} {data['years']:>12.0f} {data['energy_PJ']:>15.0f} {savings:>17.0f}%")
+ print("=" * 90)
+
+ # 生成三方案对比折线图 (100-300年)
+ print("\n\n正在生成三方案对比折线图 (100-300年)...")
+ plot_three_scenarios_comparison(year_min=100, year_max=300)
diff --git a/combined_scenario_analysis.png b/combined_scenario_analysis.png
new file mode 100644
index 0000000..1d420fc
Binary files /dev/null and b/combined_scenario_analysis.png differ
diff --git a/earth_moon_transfer.py b/earth_moon_transfer.py
deleted file mode 100644
index 5977f7a..0000000
--- a/earth_moon_transfer.py
+++ /dev/null
@@ -1,272 +0,0 @@
-"""
-地月转移轨道：有效载荷与燃料消耗（能量）关系分析
-
-基于齐奥尔科夫斯基火箭方程和霍曼转移轨道模型
-"""
-
-import numpy as np
-import matplotlib
-matplotlib.use('Agg') # 使用非交互式后端
-import matplotlib.pyplot as plt
-from matplotlib import rcParams
-
-# 设置中文字体支持
-rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans']
-rcParams['axes.unicode_minus'] = False
-
-# ============== 物理常数 ==============
-G0 = 9.81 # 地表重力加速度 (m/s²)
-
-# 地月转移轨道典型ΔV预算 (m/s)
-DELTA_V_TLI = 3100 # 地月转移注入 (Trans-Lunar Injection)
-DELTA_V_LOI = 800 # 月球轨道捕获 (Lunar Orbit Insertion)
-DELTA_V_TOTAL = DELTA_V_TLI + DELTA_V_LOI # 总ΔV
-
-# 不同发动机类型的比冲 (秒)
-ENGINES = {
- '固体火箭': 280,
- '液氧煤油': 350,
- '液氧液氢': 450,
- '离子推进 (仅供参考)': 3000,
-}
-
-# 结构系数 (结构质量/燃料质量)
-ALPHA = 0
-
-
-def exhaust_velocity(isp):
- """
- 计算排气速度
-
- 参数:
- isp: 比冲 (秒)
- 返回:
- 排气速度 (m/s)
- """
- return isp * G0
-
-
-def mass_ratio(delta_v, isp):
- """
- 计算质量比 R = m0/mf
-
- 基于齐奥尔科夫斯基方程: ΔV = ve * ln(m0/mf)
-
- 参数:
- delta_v: 速度增量 (m/s)
- isp: 比冲 (秒)
- 返回:
- 质量比 R
- """
- ve = exhaust_velocity(isp)
- return np.exp(delta_v / ve)
-
-
-def fuel_mass(payload_mass, delta_v, isp, alpha=ALPHA):
- """
- 计算所需燃料质量
-
- 公式推导:
- m0 = mp + ms + mf (初始质量 = 有效载荷 + 结构 + 燃料)
- mf_final = mp + ms (最终质量)
- ms = alpha * m_fuel (结构质量与燃料成比例)
-
- 由火箭方程: R = m0 / mf_final
- 解得: m_fuel = (mp + ms) * (R - 1)
- m_fuel = mp * (R - 1) / (1 - alpha * (R - 1))
-
- 参数:
- payload_mass: 有效载荷质量 (kg)
- delta_v: 速度增量 (m/s)
- isp: 比冲 (秒)
- alpha: 结构系数
- 返回:
- 燃料质量 (kg)
- """
- R = mass_ratio(delta_v, isp)
- denominator = 1 - alpha * (R - 1)
-
- if denominator <= 0:
- # 物理上不可能: 结构质量太大
- return np.inf
-
- return payload_mass * (R - 1) / denominator
-
-
-def total_energy(fuel_mass, isp):
- """
- 计算燃料的动能（近似能量消耗）
-
- 能量 E = 0.5 * m_fuel * ve²
-
- 参数:
- fuel_mass: 燃料质量 (kg)
- isp: 比冲 (秒)
- 返回:
- 能量 (焦耳)
- """
- ve = exhaust_velocity(isp)
- return 0.5 * fuel_mass * ve ** 2
-
-
-def characteristic_energy(delta_v):
- """
- 计算特征能量 C3 (用于轨道能量分析)
-
- C3 = V_infinity² (相对于地球的双曲线超速平方)
- 对于地月转移，C3 ≈ -2 km²/s² (束缚轨道)
-
- 简化计算: C3 ≈ (ΔV)² - V_escape²
- """
- # 近地轨道逃逸速度约 10.9 km/s，轨道速度约 7.8 km/s
- # 从LEO出发需要的额外速度约 3.1 km/s
- return (delta_v / 1000) ** 2 # km²/s²
-
-
-def plot_payload_fuel_relationship():
- """
- 绘制有效载荷与燃料消耗的关系图
- """
- # 有效载荷范围: 100 kg 到 50000 kg
- payloads = np.linspace(100, 50000, 500)
-
- fig, axes = plt.subplots(2, 2, figsize=(14, 12))
-
- # ========== 图1: 不同发动机的燃料消耗 ==========
- ax1 = axes[0, 0]
- for engine_name, isp in ENGINES.items():
- fuels = [fuel_mass(p, DELTA_V_TOTAL, isp) for p in payloads]
- fuels = np.array(fuels)
- # 过滤掉无穷大值
- valid = fuels < 1e9
- ax1.plot(payloads[valid]/1000, fuels[valid]/1000,
- label=f'{engine_name} (Isp={isp}s)', linewidth=2)
-
- ax1.set_xlabel('有效载荷质量 (吨)', fontsize=12)
- ax1.set_ylabel('燃料质量 (吨)', fontsize=12)
- ax1.set_title('地月转移：有效载荷 vs 燃料消耗\n(ΔV = 3.9 km/s)', fontsize=14)
- ax1.legend(loc='upper left')
- ax1.grid(True, alpha=0.3)
- ax1.set_xlim(0, 50)
-
- # ========== 图2: 燃料/载荷比 ==========
- ax2 = axes[0, 1]
- for engine_name, isp in list(ENGINES.items())[:3]: # 排除离子推进
- ratios = [fuel_mass(p, DELTA_V_TOTAL, isp) / p for p in payloads]
- ax2.axhline(y=ratios[0], label=f'{engine_name} (Isp={isp}s)', linewidth=2)
-
- ax2.set_xlabel('有效载荷质量 (吨)', fontsize=12)
- ax2.set_ylabel('燃料/载荷 质量比', fontsize=12)
- ax2.set_title('燃料效率比较\n(质量比与载荷无关，仅取决于Isp和ΔV)', fontsize=14)
- ax2.legend()
- ax2.grid(True, alpha=0.3)
- ax2.set_ylim(0, 10)
-
- # ========== 图3: 能量消耗 ==========
- ax3 = axes[1, 0]
- for engine_name, isp in list(ENGINES.items())[:3]:
- fuels = np.array([fuel_mass(p, DELTA_V_TOTAL, isp) for p in payloads])
- energies = total_energy(fuels, isp) / 1e12 # 转换为TJ
- valid = energies < 1e6
- ax3.plot(payloads[valid]/1000, energies[valid],
- label=f'{engine_name}', linewidth=2)
-
- ax3.set_xlabel('有效载荷质量 (吨)', fontsize=12)
- ax3.set_ylabel('能量消耗 (TJ)', fontsize=12)
- ax3.set_title('能量消耗 vs 有效载荷\nE = ½ × m_fuel × ve²', fontsize=14)
- ax3.legend()
- ax3.grid(True, alpha=0.3)
-
- # ========== 图4: ΔV 敏感性分析 ==========
- ax4 = axes[1, 1]
- delta_v_range = np.linspace(3000, 5000, 100) # 3-5 km/s
- payload_fixed = 10000 # 10吨有效载荷
- isp_reference = 450 # 液氧液氢
-
- fuel_vs_dv = [fuel_mass(payload_fixed, dv, isp_reference) for dv in delta_v_range]
- ax4.plot(delta_v_range/1000, np.array(fuel_vs_dv)/1000,
- 'b-', linewidth=2, label='燃料质量')
-
- # 标记关键点
- ax4.axvline(x=DELTA_V_TOTAL/1000, color='r', linestyle='--',
- label=f'标准任务 ΔV={DELTA_V_TOTAL/1000} km/s')
-
- ax4.set_xlabel('ΔV (km/s)', fontsize=12)
- ax4.set_ylabel('燃料质量 (吨)', fontsize=12)
- ax4.set_title(f'ΔV敏感性分析\n(有效载荷=10吨, Isp=450s)', fontsize=14)
- ax4.legend()
- ax4.grid(True, alpha=0.3)
-
- plt.tight_layout()
- plt.savefig('/Volumes/Files/code/mm/20260130_b/earth_moon_transfer_analysis.png',
- dpi=150, bbox_inches='tight')
- print("图表已生成!")
- # plt.show() # 非交互式环境下注释掉
-
- return fig
-
-
-def print_summary():
- """
- 打印关键数据摘要
- """
- print("=" * 60)
- print("地月转移轨道分析摘要")
- print("=" * 60)
- print(f"\n任务参数:")
- print(f" - 总ΔV需求: {DELTA_V_TOTAL/1000:.1f} km/s")
- print(f" - TLI (地月转移注入): {DELTA_V_TLI/1000:.1f} km/s")
- print(f" - LOI (月球轨道捕获): {DELTA_V_LOI/1000:.1f} km/s")
- print(f" - 结构系数 α: {ALPHA}")
-
- print(f"\n不同发动机的燃料效率 (每吨有效载荷所需燃料):")
- print("-" * 50)
-
- payload_ref = 1000 # 1吨参考载荷
-
- for engine_name, isp in ENGINES.items():
- R = mass_ratio(DELTA_V_TOTAL, isp)
- fuel = fuel_mass(payload_ref, DELTA_V_TOTAL, isp)
- ve = exhaust_velocity(isp)
-
- if fuel < 1e9:
- print(f"\n {engine_name}:")
- print(f" 比冲 Isp = {isp} s")
- print(f" 排气速度 ve = {ve:.0f} m/s")
- print(f" 质量比 R = {R:.2f}")
- print(f" 燃料/载荷比 = {fuel/payload_ref:.2f}")
- print(f" 1吨载荷需燃料: {fuel/1000:.2f} 吨")
- else:
- print(f"\n {engine_name}:")
- print(f" 比冲 Isp = {isp} s (化学火箭无法实现)")
-
- print("\n" + "=" * 60)
- print("关键数学关系:")
- print("=" * 60)
- print("""
- 齐奥尔科夫斯基方程:
- ΔV = ve × ln(m₀/mf)
-
- 质量比:
- R = exp(ΔV / ve)
-
- 燃料-载荷关系:
- m_fuel / m_payload = (R - 1) / (1 - α(R - 1))
-
- 能量消耗:
- E = ½ × m_fuel × ve²
-
- 其中:
- ve = Isp × g₀ (排气速度)
- α = 结构系数 (约0.1)
-""")
-
-
-if __name__ == "__main__":
- # 打印分析摘要
- print_summary()
-
- # 绘制关系图
- print("\n正在生成图表...")
- fig = plot_payload_fuel_relationship()
- print("图表已保存至: earth_moon_transfer_analysis.png")
diff --git a/launch_capacity_analysis.png b/launch_capacity_analysis.png
new file mode 100644
index 0000000..5b2973e
Binary files /dev/null and b/launch_capacity_analysis.png differ
diff --git a/launch_capacity_analysis.py b/launch_capacity_analysis.py
new file mode 100644
index 0000000..48665b3
--- /dev/null
+++ b/launch_capacity_analysis.py
@@ -0,0 +1,785 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+火箭发射数上限预测分析
+对比 Logistic、Gompertz、Richards 三种增长模型
+分析不同历史数据窗口对预测的影响
+"""
+
+import pandas as pd
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.stats import pearsonr
+import matplotlib.pyplot as plt
+import warnings
+warnings.filterwarnings('ignore')
+
+plt.rcParams["font.sans-serif"] = ["Arial Unicode MS", "SimHei", "DejaVu Sans"]
+plt.rcParams["axes.unicode_minus"] = False
+
+
+# ============================================================
+# 1. 增长模型定义
+# ============================================================
+
+def logistic(t, K, r, t0):
+ """
+ Logistic 增长模型（对称S曲线）
+ N(t) = K / (1 + exp(-r*(t - t0)))
+
+ 参数:
+ K : 环境容量/饱和上限
+ r : 增长率
+ t0 : 拐点时间（增长最快的时刻，此时 N = K/2）
+ """
+ return K / (1 + np.exp(-r * (t - t0)))
+
+
+def gompertz(t, K, b, c):
+ """
+ Gompertz 增长模型（非对称S曲线，早期增长较快）
+ N(t) = K * exp(-b * exp(-c*t))
+
+ 参数:
+ K : 饱和上限
+ b : 位移参数
+ c : 增长率参数
+ 拐点在 N = K/e ≈ 0.368K 处
+ """
+ return K * np.exp(-b * np.exp(-c * t))
+
+
+def richards(t, K, r, t0, v):
+ """
+ Richards 曲线（广义 Logistic，可调节对称性）
+ N(t) = K / (1 + exp(-r*(t - t0)))^(1/v)
+
+ 参数:
+ K : 饱和上限
+ r : 增长率
+ t0 : 拐点时间
+ v : 形状参数（v=1 时退化为 Logistic）
+ """
+ exp_term = np.exp(-r * (t - t0))
+ # 避免数值溢出
+ exp_term = np.clip(exp_term, 1e-10, 1e10)
+ return K / np.power(1 + exp_term, 1/v)
+
+
+# ============================================================
+# 2. 数据加载与预处理
+# ============================================================
+
+def load_data(filepath="rocket_launch_counts.csv"):
+ """加载并清洗火箭发射数据"""
+ df = pd.read_csv(filepath)
+ df = df.rename(columns={"YDate": "year", "Total": "launches"})
+ df["year"] = pd.to_numeric(df["year"], errors="coerce")
+ df["launches"] = pd.to_numeric(df["launches"], errors="coerce")
+ df = df.dropna(subset=["year", "launches"])
+ # 只保留完整年份数据（截至2025年，2026年数据不完整）
+ df = df[(df["year"] >= 1957) & (df["year"] <= 2025)]
+ df = df.astype({"year": int, "launches": int})
+ df = df.sort_values("year").reset_index(drop=True)
+ return df
+
+
+# ============================================================
+# 3. 模型拟合
+# ============================================================
+
+def fit_model(model_func, t, y, model_name, p0, bounds):
+ """通用模型拟合函数"""
+ try:
+ popt, pcov = curve_fit(model_func, t, y, p0=p0, bounds=bounds, maxfev=50000)
+ perr = np.sqrt(np.diag(pcov))
+ y_pred = model_func(t, *popt)
+
+ # 计算拟合优度
+ ss_res = np.sum((y - y_pred) ** 2)
+ ss_tot = np.sum((y - np.mean(y)) ** 2)
+ r_squared = 1 - (ss_res / ss_tot)
+ rmse = np.sqrt(np.mean((y - y_pred) ** 2))
+
+ return {
+ "success": True,
+ "params": popt,
+ "param_err": perr,
+ "r_squared": r_squared,
+ "rmse": rmse,
+ "y_pred": y_pred
+ }
+ except Exception as e:
+ return {
+ "success": False,
+ "error": str(e)
+ }
+
+
+def fit_all_models(data, base_year=1957):
+ """对数据拟合所有三种模型"""
+ years = data["year"].values
+ launches = data["launches"].values
+ t = (years - base_year).astype(float)
+
+ results = {}
+
+ # Logistic 模型 - 提高上限边界以避免撞边界
+ p0_log = [2000.0, 0.08, 70.0]
+ bounds_log = ([300, 0.01, 0], [500000, 1.0, 300])
+ res = fit_model(logistic, t, launches, "Logistic", p0_log, bounds_log)
+ if res["success"]:
+ K, r, t0 = res["params"]
+ res["K"] = K
+ res["inflection_year"] = base_year + t0
+ res["inflection_value"] = K / 2
+ results["Logistic"] = res
+
+ # Gompertz 模型
+ p0_gom = [2000.0, 5.0, 0.03]
+ bounds_gom = ([300, 0.5, 0.005], [500000, 50.0, 0.5])
+ res = fit_model(gompertz, t, launches, "Gompertz", p0_gom, bounds_gom)
+ if res["success"]:
+ K, b, c = res["params"]
+ res["K"] = K
+ # Gompertz 拐点：t = ln(b)/c, 值为 K/e
+ res["inflection_year"] = base_year + np.log(b) / c
+ res["inflection_value"] = K / np.e
+ results["Gompertz"] = res
+
+ # Richards 模型
+ p0_ric = [2000.0, 0.08, 70.0, 2.0]
+ bounds_ric = ([300, 0.01, 0, 0.1], [500000, 1.0, 300, 10.0])
+ res = fit_model(richards, t, launches, "Richards", p0_ric, bounds_ric)
+ if res["success"]:
+ K, r, t0, v = res["params"]
+ res["K"] = K
+ res["inflection_year"] = base_year + t0
+ res["inflection_value"] = K / (v + 1) ** (1/v)
+ res["v"] = v
+ results["Richards"] = res
+
+ return results, years, launches, t, base_year
+
+
+def predict_future(results, base_year, years_future):
+ """生成未来预测"""
+ t_future = years_future - base_year
+ predictions = {}
+
+ for model_name, res in results.items():
+ if not res["success"]:
+ continue
+
+ params = res["params"]
+ if model_name == "Logistic":
+ pred = logistic(t_future, *params)
+ elif model_name == "Gompertz":
+ pred = gompertz(t_future, *params)
+ elif model_name == "Richards":
+ pred = richards(t_future, *params)
+
+ predictions[model_name] = np.maximum(pred, 0)
+
+ return predictions
+
+
+# ============================================================
+# 4. 可视化
+# ============================================================
+
+def plot_model_comparison(df, results_dict, save_path="launch_capacity_analysis.png"):
+ """绘制模型对比图"""
+ fig, axes = plt.subplots(2, 2, figsize=(16, 12))
+
+ # 颜色方案
+ colors = {
+ "Logistic": "#E74C3C",
+ "Gompertz": "#3498DB",
+ "Richards": "#27AE60"
+ }
+
+ # 子图1: 历史数据 + 三种模型拟合（全部数据）
+ ax1 = axes[0, 0]
+ results, years, launches, t, base_year = results_dict["全部数据"]
+ ax1.scatter(years, launches, color="gray", s=30, alpha=0.7, label="历史数据", zorder=3)
+
+ years_curve = np.arange(years.min(), 2101)
+ t_curve = years_curve - base_year
+
+ for model_name, res in results.items():
+ if not res["success"]:
+ continue
+ params = res["params"]
+ if model_name == "Logistic":
+ pred = logistic(t_curve, *params)
+ elif model_name == "Gompertz":
+ pred = gompertz(t_curve, *params)
+ elif model_name == "Richards":
+ pred = richards(t_curve, *params)
+
+ ax1.plot(years_curve, pred, color=colors[model_name], lw=2,
+ label=f"{model_name} (K={res['K']:.0f}, R²={res['r_squared']:.3f})")
+
+ ax1.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
+ ax1.set_xlabel("年份", fontsize=11)
+ ax1.set_ylabel("年发射次数", fontsize=11)
+ ax1.set_title("全部历史数据 (1957-2025) 模型拟合", fontsize=12)
+ ax1.legend(loc="upper left", fontsize=9)
+ ax1.grid(True, alpha=0.3)
+ ax1.set_xlim(1955, 2105)
+ ax1.set_ylim(0, None)
+
+ # 子图2: 近15年数据拟合
+ ax2 = axes[0, 1]
+ results, years, launches, t, base_year = results_dict["近15年"]
+ ax2.scatter(years, launches, color="gray", s=40, alpha=0.8, label="历史数据", zorder=3)
+
+ years_curve = np.arange(2010, 2101)
+ t_curve = years_curve - base_year
+
+ for model_name, res in results.items():
+ if not res["success"]:
+ continue
+ params = res["params"]
+ if model_name == "Logistic":
+ pred = logistic(t_curve, *params)
+ elif model_name == "Gompertz":
+ pred = gompertz(t_curve, *params)
+ elif model_name == "Richards":
+ pred = richards(t_curve, *params)
+
+ ax2.plot(years_curve, pred, color=colors[model_name], lw=2,
+ label=f"{model_name} (K={res['K']:.0f})")
+
+ ax2.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
+ ax2.set_xlabel("年份", fontsize=11)
+ ax2.set_ylabel("年发射次数", fontsize=11)
+ ax2.set_title("近15年数据 (2010-2025) 模型拟合", fontsize=12)
+ ax2.legend(loc="upper left", fontsize=9)
+ ax2.grid(True, alpha=0.3)
+ ax2.set_xlim(2008, 2105)
+ ax2.set_ylim(0, None)
+
+ # 子图3: 只展示 Richards 模型的K值（因为其他模型撞边界）
+ ax3 = axes[1, 0]
+ windows = ["全部数据", "近25年", "近15年", "近10年"]
+ x = np.arange(len(windows))
+
+ k_values = []
+ for window in windows:
+ if window in results_dict:
+ res = results_dict[window][0].get("Richards", {})
+ k_values.append(res.get("K", 0) if res.get("success", False) else 0)
+ else:
+ k_values.append(0)
+
+ bars = ax3.bar(x, k_values, 0.6, color=colors["Richards"], alpha=0.8, edgecolor='black')
+
+ # 添加数值标签
+ for bar, k in zip(bars, k_values):
+ if k > 0:
+ ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 500,
+ f'{k:.0f}', ha='center', va='bottom', fontsize=10, fontweight='bold')
+
+ ax3.set_xlabel("数据窗口", fontsize=11)
+ ax3.set_ylabel("预测饱和上限 K (次/年)", fontsize=11)
+ ax3.set_title("Richards 模型: 不同数据窗口的饱和上限预测\n(Logistic/Gompertz 撞边界，不适用)", fontsize=11)
+ ax3.set_xticks(x)
+ ax3.set_xticklabels(windows)
+ ax3.grid(True, alpha=0.3, axis='y')
+ ax3.set_ylim(0, max(k_values) * 1.15 if k_values else 10000)
+
+ # 子图4: Richards 模型拐点年份
+ ax4 = axes[1, 1]
+ inflection_years = []
+ inflection_values = []
+ for window in windows:
+ if window in results_dict:
+ res = results_dict[window][0].get("Richards", {})
+ if res.get("success", False):
+ inflection_years.append(res.get("inflection_year", 0))
+ inflection_values.append(res.get("inflection_value", 0))
+ else:
+ inflection_years.append(0)
+ inflection_values.append(0)
+ else:
+ inflection_years.append(0)
+ inflection_values.append(0)
+
+ bars = ax4.bar(x, inflection_years, 0.6, color=colors["Richards"], alpha=0.8, edgecolor='black')
+
+ # 添加标签（拐点年份和对应发射数）
+ for bar, year, val in zip(bars, inflection_years, inflection_values):
+ if year > 0:
+ ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.5,
+ f'{year:.0f}\n({val:.0f}次/年)', ha='center', va='bottom', fontsize=9)
+
+ ax4.axhline(2025, color="red", ls="--", lw=2, alpha=0.8, label="当前年份 (2025)")
+ ax4.set_xlabel("数据窗口", fontsize=11)
+ ax4.set_ylabel("拐点年份", fontsize=11)
+ ax4.set_title("Richards 模型: 预测拐点年份\n(增长最快的时刻)", fontsize=11)
+ ax4.set_xticks(x)
+ ax4.set_xticklabels(windows)
+ ax4.legend(loc="lower right", fontsize=9)
+ ax4.grid(True, alpha=0.3, axis='y')
+ valid_years = [y for y in inflection_years if y > 0]
+ if valid_years:
+ ax4.set_ylim(2020, max(valid_years) + 15)
+
+ plt.suptitle("火箭年发射数上限预测 - 多模型对比分析", fontsize=14, fontweight='bold', y=1.02)
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches="tight")
+ plt.close()
+ print(f"图表已保存: {save_path}")
+
+
+def plot_physical_vs_model(results_dict, physical_limits, save_path="launch_capacity_physical.png"):
+ """绘制物理约束与模型预测对比图"""
+ fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+
+ # 左图：物理约束与模型预测K值对比
+ ax1 = axes[0]
+
+ # 物理约束
+ phys_scenarios = list(physical_limits.keys())
+ phys_values = list(physical_limits.values())
+ x_phys = np.arange(len(phys_scenarios))
+ bars1 = ax1.barh(x_phys, phys_values, 0.6, color='#3498DB', alpha=0.7, label='物理约束估计')
+
+ # 添加数值标签
+ for bar, val in zip(bars1, phys_values):
+ ax1.text(val + 50, bar.get_y() + bar.get_height()/2,
+ f'{val:,}', va='center', fontsize=10)
+
+ # 添加模型预测（Richards）
+ richards_k = []
+ for window in ["近25年", "近15年", "近10年"]:
+ if window in results_dict:
+ res = results_dict[window][0].get("Richards", {})
+ if res.get("success", False):
+ richards_k.append(res["K"])
+
+ if richards_k:
+ avg_richards = np.mean(richards_k)
+ ax1.axvline(avg_richards, color='#E74C3C', ls='--', lw=2,
+ label=f'Richards模型平均 K={avg_richards:,.0f}')
+
+ ax1.set_yticks(x_phys)
+ ax1.set_yticklabels(phys_scenarios)
+ ax1.set_xlabel("年发射上限 (次/年)", fontsize=11)
+ ax1.set_title("物理约束 vs 模型预测 (10个发射场)", fontsize=12)
+ ax1.legend(loc="lower right", fontsize=9)
+ ax1.grid(True, alpha=0.3, axis='x')
+ ax1.set_xlim(0, max(phys_values) * 1.3)
+
+ # 右图：带物理约束的时间预测
+ ax2 = axes[1]
+
+ # 获取近15年Richards预测
+ if "近15年" in results_dict:
+ results, years, launches, t, base_year = results_dict["近15年"]
+ res = results.get("Richards", {})
+
+ if res["success"]:
+ years_future = np.arange(2010, 2101)
+ t_future = years_future - base_year
+ pred = richards(t_future, *res["params"])
+
+ ax2.scatter(years, launches, color="gray", s=40, alpha=0.7, label="历史数据", zorder=3)
+ ax2.plot(years_future, pred, color='#27AE60', lw=2, label=f"Richards (K={res['K']:,.0f})")
+
+ # 添加物理约束线
+ phys_colors = ['#E74C3C', '#F39C12', '#3498DB', '#9B59B6', '#1ABC9C']
+ for i, (scenario, limit) in enumerate(physical_limits.items()):
+ ax2.axhline(limit, color=phys_colors[i % len(phys_colors)],
+ ls='--', lw=1.5, alpha=0.7, label=f'{scenario}: {limit:,}')
+
+ ax2.axvline(2025, color="gray", ls=":", lw=1, alpha=0.5)
+ ax2.set_xlabel("年份", fontsize=11)
+ ax2.set_ylabel("年发射次数", fontsize=11)
+ ax2.set_title("Richards预测 + 物理约束线", fontsize=12)
+ ax2.legend(loc="upper left", fontsize=8)
+ ax2.grid(True, alpha=0.3)
+ ax2.set_xlim(2008, 2105)
+ ax2.set_ylim(0, max(physical_limits.values()) * 1.2)
+
+ plt.suptitle("物理约束分析 (10个发射场)", fontsize=14, fontweight='bold', y=1.02)
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches="tight")
+ plt.close()
+ print(f"图表已保存: {save_path}")
+
+
+def plot_future_predictions(results_dict, save_path="launch_capacity_predictions.png"):
+ """绘制未来预测图"""
+ fig, axes = plt.subplots(1, 2, figsize=(16, 6))
+
+ colors = {
+ "Logistic": "#E74C3C",
+ "Gompertz": "#3498DB",
+ "Richards": "#27AE60"
+ }
+
+ # 左图：基于全部数据的预测
+ ax1 = axes[0]
+ results, years, launches, t, base_year = results_dict["全部数据"]
+ years_future = np.arange(1957, 2151)
+ predictions = predict_future(results, base_year, years_future)
+
+ ax1.scatter(years, launches, color="gray", s=20, alpha=0.6, label="历史数据", zorder=3)
+ for model_name, pred in predictions.items():
+ ax1.plot(years_future, pred, color=colors[model_name], lw=2, label=model_name)
+
+ ax1.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
+ ax1.axvspan(2025, 2150, alpha=0.05, color="blue")
+ ax1.set_xlabel("年份", fontsize=11)
+ ax1.set_ylabel("年发射次数", fontsize=11)
+ ax1.set_title("基于全部历史数据的预测 (至2150年)", fontsize=12)
+ ax1.legend(loc="upper left", fontsize=10)
+ ax1.grid(True, alpha=0.3)
+ ax1.set_xlim(1955, 2155)
+ ax1.set_ylim(0, None)
+
+ # 右图：基于近15年数据的预测
+ ax2 = axes[1]
+ results, years, launches, t, base_year = results_dict["近15年"]
+ years_future = np.arange(2010, 2151)
+ predictions = predict_future(results, base_year, years_future)
+
+ ax2.scatter(years, launches, color="gray", s=40, alpha=0.7, label="历史数据", zorder=3)
+ for model_name, pred in predictions.items():
+ ax2.plot(years_future, pred, color=colors[model_name], lw=2, label=model_name)
+
+ ax2.axvline(2025, color="gray", ls="--", lw=1, alpha=0.5)
+ ax2.axvspan(2025, 2150, alpha=0.05, color="blue")
+ ax2.set_xlabel("年份", fontsize=11)
+ ax2.set_ylabel("年发射次数", fontsize=11)
+ ax2.set_title("基于近15年数据的预测 (至2150年)", fontsize=12)
+ ax2.legend(loc="upper left", fontsize=10)
+ ax2.grid(True, alpha=0.3)
+ ax2.set_xlim(2008, 2155)
+ ax2.set_ylim(0, None)
+
+ plt.suptitle("火箭年发射数长期预测", fontsize=14, fontweight='bold', y=1.02)
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches="tight")
+ plt.close()
+ print(f"图表已保存: {save_path}")
+
+
+# ============================================================
+# 5. 主程序
+# ============================================================
+
+def bootstrap_K_estimate(data, model_func, p0, bounds, n_bootstrap=500, base_year=1957):
+ """Bootstrap 估计 K 值的不确定性"""
+ years = data["year"].values
+ launches = data["launches"].values
+ t = (years - base_year).astype(float)
+ n = len(t)
+
+ K_samples = []
+ for _ in range(n_bootstrap):
+ # 重采样
+ idx = np.random.choice(n, n, replace=True)
+ t_boot = t[idx]
+ y_boot = launches[idx]
+
+ try:
+ popt, _ = curve_fit(model_func, t_boot, y_boot, p0=p0, bounds=bounds, maxfev=10000)
+ K_samples.append(popt[0]) # K 是第一个参数
+ except:
+ pass
+
+ if len(K_samples) > 10:
+ return {
+ "mean": np.mean(K_samples),
+ "median": np.median(K_samples),
+ "std": np.std(K_samples),
+ "ci_5": np.percentile(K_samples, 5),
+ "ci_95": np.percentile(K_samples, 95),
+ "samples": K_samples
+ }
+ return None
+
+
+def calculate_aic_bic(y_true, y_pred, n_params):
+ """计算 AIC 和 BIC"""
+ n = len(y_true)
+ ss_res = np.sum((y_true - y_pred) ** 2)
+ mse = ss_res / n
+
+ # Log-likelihood (假设高斯误差)
+ log_likelihood = -n/2 * (np.log(2 * np.pi) + np.log(mse) + 1)
+
+ aic = 2 * n_params - 2 * log_likelihood
+ bic = n_params * np.log(n) - 2 * log_likelihood
+
+ return aic, bic
+
+
+def main():
+ print("=" * 70)
+ print("火箭年发射数上限预测分析")
+ print("模型: Logistic, Gompertz, Richards")
+ print("=" * 70)
+
+ # 加载数据
+ df = load_data()
+ print(f"\n数据范围: {df['year'].min()} - {df['year'].max()}")
+ print(f"数据点数: {len(df)}")
+ print(f"最近5年发射数:")
+ for _, row in df.tail(5).iterrows():
+ print(f" {int(row['year'])}: {int(row['launches'])} 次")
+
+ # 计算增长率
+ recent = df[df["year"] >= 2020].copy()
+ if len(recent) > 1:
+ growth_rates = recent["launches"].pct_change().dropna() * 100
+ print(f"\n2020年以来年均增长率: {growth_rates.mean():.1f}%")
+
+ # 定义不同的数据窗口
+ end_year = 2025
+ windows = {
+ "全部数据": df[df["year"] <= end_year].copy(),
+ "近25年": df[(df["year"] >= 2000) & (df["year"] <= end_year)].copy(),
+ "近15年": df[(df["year"] >= 2010) & (df["year"] <= end_year)].copy(),
+ "近10年": df[(df["year"] >= 2015) & (df["year"] <= end_year)].copy(),
+ }
+
+ # 存储所有结果
+ all_results = {}
+
+ # 对每个窗口进行拟合
+ for window_name, data in windows.items():
+ print(f"\n{'='*50}")
+ print(f"数据窗口: {window_name} ({data['year'].min()}-{data['year'].max()}, n={len(data)})")
+ print("=" * 50)
+
+ results, years, launches, t, base_year = fit_all_models(data)
+ all_results[window_name] = (results, years, launches, t, base_year)
+
+ for model_name, res in results.items():
+ if res["success"]:
+ print(f"\n{model_name} 模型:")
+ print(f" 饱和上限 K = {res['K']:.0f} 次/年")
+ print(f" 拐点年份 = {res['inflection_year']:.1f}")
+ print(f" 拐点处发射数 = {res['inflection_value']:.0f} 次/年")
+ print(f" R² (拟合优度) = {res['r_squared']:.4f}")
+ print(f" RMSE = {res['rmse']:.2f}")
+ if model_name == "Richards":
+ print(f" 形状参数 v = {res.get('v', 'N/A'):.3f}")
+ else:
+ print(f"\n{model_name} 模型: 拟合失败 - {res.get('error', 'Unknown')}")
+
+ # 生成未来预测表
+ print("\n" + "=" * 70)
+ print("未来预测 (基于近15年数据)")
+ print("=" * 70)
+
+ results, years, launches, t, base_year = all_results["近15年"]
+ future_years = np.array([2030, 2040, 2050, 2075, 2100])
+ predictions = predict_future(results, base_year, future_years)
+
+ print(f"\n{'年份':<10}", end="")
+ for model in ["Logistic", "Gompertz", "Richards"]:
+ print(f"{model:<15}", end="")
+ print()
+ print("-" * 55)
+
+ for i, year in enumerate(future_years):
+ print(f"{year:<10}", end="")
+ for model in ["Logistic", "Gompertz", "Richards"]:
+ if model in predictions:
+ print(f"{predictions[model][i]:<15.0f}", end="")
+ else:
+ print(f"{'N/A':<15}", end="")
+ print()
+
+ # 饱和上限汇总
+ print("\n" + "=" * 70)
+ print("饱和上限 K 汇总 (次/年)")
+ print("=" * 70)
+ print(f"\n{'数据窗口':<20}", end="")
+ for model in ["Logistic", "Gompertz", "Richards"]:
+ print(f"{model:<15}", end="")
+ print()
+ print("-" * 65)
+
+ for window_name in ["全部数据", "近25年", "近15年", "近10年"]:
+ results = all_results[window_name][0]
+ print(f"{window_name:<20}", end="")
+ for model in ["Logistic", "Gompertz", "Richards"]:
+ if model in results and results[model]["success"]:
+ print(f"{results[model]['K']:<15.0f}", end="")
+ else:
+ print(f"{'N/A':<15}", end="")
+ print()
+
+ # 物理约束分析
+ physical_limits = analyze_physical_constraints(n_sites=10)
+
+ # 绘制图表
+ print("\n" + "=" * 70)
+ print("生成图表...")
+ plot_model_comparison(df, all_results)
+ plot_future_predictions(all_results)
+ plot_physical_vs_model(all_results, physical_limits)
+
+ # 保存详细结果到CSV
+ save_results_csv(all_results)
+
+ print("\n" + "=" * 70)
+ print("分析完成!")
+ print("=" * 70)
+
+ # Bootstrap 不确定性分析
+ print("\n" + "=" * 70)
+ print("Bootstrap 不确定性分析 (基于近15年数据, 500次重采样)")
+ print("=" * 70)
+
+ data_15 = windows["近15年"]
+ base_year = 1957
+
+ # Richards 模型 Bootstrap
+ p0_ric = [2000.0, 0.08, 70.0, 2.0]
+ bounds_ric = ([300, 0.01, 0, 0.1], [500000, 1.0, 300, 10.0])
+
+ boot_result = bootstrap_K_estimate(data_15, richards, p0_ric, bounds_ric, n_bootstrap=500, base_year=base_year)
+ if boot_result:
+ print(f"\nRichards 模型 K 值 Bootstrap 估计:")
+ print(f" 均值: {boot_result['mean']:.0f} 次/年")
+ print(f" 中位数: {boot_result['median']:.0f} 次/年")
+ print(f" 标准差: {boot_result['std']:.0f}")
+ print(f" 90% 置信区间: [{boot_result['ci_5']:.0f}, {boot_result['ci_95']:.0f}] 次/年")
+
+ # AIC/BIC 模型选择
+ print("\n" + "=" * 70)
+ print("模型选择 (AIC/BIC, 基于近15年数据)")
+ print("=" * 70)
+
+ results_15, years_15, launches_15, t_15, _ = all_results["近15年"]
+ print(f"\n{'模型':<12} {'AIC':<12} {'BIC':<12} {'参数数':<8}")
+ print("-" * 44)
+
+ for model_name, res in results_15.items():
+ if res["success"]:
+ n_params = len(res["params"])
+ aic, bic = calculate_aic_bic(launches_15, res["y_pred"], n_params)
+ print(f"{model_name:<12} {aic:<12.1f} {bic:<12.1f} {n_params:<8}")
+
+ # 结论
+ print("\n" + "=" * 70)
+ print("📊 分析结论与建议")
+ print("=" * 70)
+
+ # 收集有效的K值（排除撞边界的）
+ valid_k_values = []
+ for window_name, (results, _, _, _, _) in all_results.items():
+ for model_name, res in results.items():
+ if res["success"] and res["K"] < 100000: # 排除撞边界的
+ valid_k_values.append((window_name, model_name, res["K"]))
+
+ print("\n1. 模型适用性分析:")
+ print("-" * 50)
+ print(" • Logistic/Gompertz: 数据处于快速增长期，无法确定上限")
+ print(" (模型撞到边界，说明增长远未饱和)")
+ print(" • Richards 曲线: 额外的形状参数使其能给出更合理估计")
+ print(" (推荐用于当前阶段的上限预测)")
+
+ print("\n2. 合理的上限估计 (Richards 模型):")
+ print("-" * 50)
+ for window_name, model_name, k in valid_k_values:
+ if model_name == "Richards" and window_name != "全部数据":
+ print(f" • {window_name}: K ≈ {k:.0f} 次/年")
+
+ # Richards 结果的平均值
+ richards_k = [k for w, m, k in valid_k_values if m == "Richards" and w != "全部数据"]
+ if richards_k:
+ print(f"\n 📌 综合估计: {int(np.mean(richards_k)):,} - {int(np.max(richards_k)):,} 次/年")
+
+ print("\n3. 物理约束参考 (基于10个主要发射场):")
+ print("-" * 50)
+ print(" • 主要发射场数量: 10个")
+ print(" • 假设每发射场年最大容量:")
+ print(" - 保守估计: 50次/年 → 总上限 500 次/年")
+ print(" - 中等估计: 100次/年 → 总上限 1,000 次/年")
+ print(" - 乐观估计: 200次/年 → 总上限 2,000 次/年")
+ print(" - 极限估计 (SpaceX级别): 400次/年 → 总上限 4,000 次/年")
+ print(" • 2025年实际数据: 329次 (10发射场平均 33次/场)")
+ print(" • SpaceX Starship 单型号目标: 1,000+ 次/年")
+
+ print("\n4. 预测不确定性:")
+ print("-" * 50)
+ print(" • 当前处于指数增长早期，上限估计高度不确定")
+ print(" • 建议每年用新数据更新预测")
+ print(" • 需关注技术突破（可重复使用、星舰等）对上限的影响")
+
+ if boot_result:
+ print(f"\n Bootstrap 90%置信区间: {boot_result['ci_5']:.0f} - {boot_result['ci_95']:.0f} 次/年")
+
+
+def analyze_physical_constraints(n_sites=10):
+ """
+ 基于发射场物理约束的上限分析
+
+ 参数:
+ n_sites: 发射场数量
+ """
+ print("\n" + "=" * 70)
+ print(f"物理约束分析 (基于 {n_sites} 个发射场)")
+ print("=" * 70)
+
+ # 不同情景下的单发射场年容量
+ scenarios = {
+ "保守 (当前技术)": 50,
+ "中等 (技术进步)": 100,
+ "乐观 (快速周转)": 200,
+ "极限 (SpaceX级)": 400,
+ "理论极限 (每天1发)": 365,
+ }
+
+ print(f"\n单发射场年发射容量假设:")
+ print("-" * 50)
+
+ results = {}
+ for scenario, capacity_per_site in scenarios.items():
+ total_capacity = n_sites * capacity_per_site
+ results[scenario] = total_capacity
+ print(f" {scenario:<20}: {capacity_per_site:>3} 次/场/年 → 总上限 {total_capacity:>5,} 次/年")
+
+ # 带置信度的综合估计
+ print(f"\n综合估计 ({n_sites} 发射场):")
+ print("-" * 50)
+ print(f" 50% 置信区间: {n_sites * 50:,} - {n_sites * 100:,} 次/年")
+ print(f" 90% 置信区间: {n_sites * 30:,} - {n_sites * 200:,} 次/年")
+ print(f" 理论最大值: {n_sites * 365:,} 次/年 (每天每场1发)")
+
+ return results
+
+
+def save_results_csv(all_results, filepath="launch_capacity_results.csv"):
+ """保存详细结果到CSV"""
+ rows = []
+ for window_name, (results, _, _, _, _) in all_results.items():
+ for model_name, res in results.items():
+ if res["success"]:
+ row = {
+ "数据窗口": window_name,
+ "模型": model_name,
+ "饱和上限K": res["K"],
+ "拐点年份": res["inflection_year"],
+ "拐点发射数": res["inflection_value"],
+ "R²": res["r_squared"],
+ "RMSE": res["rmse"]
+ }
+ rows.append(row)
+
+ df_results = pd.DataFrame(rows)
+ df_results.to_csv(filepath, index=False, encoding="utf-8-sig")
+ print(f"详细结果已保存: {filepath}")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/launch_capacity_physical.png b/launch_capacity_physical.png
new file mode 100644
index 0000000..5ba9462
Binary files /dev/null and b/launch_capacity_physical.png differ
diff --git a/launch_capacity_predictions.png b/launch_capacity_predictions.png
new file mode 100644
index 0000000..c9df37e
Binary files /dev/null and b/launch_capacity_predictions.png differ
diff --git a/launch_capacity_results.csv b/launch_capacity_results.csv
new file mode 100644
index 0000000..c2b776a
--- /dev/null
+++ b/launch_capacity_results.csv
@@ -0,0 +1,13 @@
+数据窗口,模型,饱和上限K,拐点年份,拐点发射数,R²,RMSE
+全部数据,Logistic,405.92414403335295,2099.062819804582,202.96207201667647,0.07766486971853037,46.74814967377199
+全部数据,Gompertz,401.0329954671314,2051.6490738033563,147.53179426375786,0.07545958812791254,46.80400317046386
+全部数据,Richards,687.2209173212347,2256.999921701128,512.1739060128197,0.08220581306057739,46.632929826767096
+近25年,Logistic,499998.291933841,2107.433566292766,249999.1459669205,0.8047914079948382,29.605437976518306
+近25年,Gompertz,499999.9999896808,2216.1161350271577,183939.72058192495,0.779423273804876,31.47037146055304
+近25年,Richards,16788.895636148314,2070.4580199940424,10907.704452410562,0.8048416215235632,29.60163002650332
+近15年,Logistic,499999.9999933337,2082.166434872205,249999.99999666685,0.8979462268050302,23.639528906752357
+近15年,Gompertz,499999.9999999759,2153.207304429069,183939.72058571232,0.8798702599574121,25.647766611285153
+近15年,Richards,32126.622670684326,2061.1007442806736,20365.326256859138,0.8979860420054704,23.63491710078375
+近10年,Logistic,499999.9996950232,2069.5299620988108,249999.9998475116,0.9632290653308331,15.01250946110567
+近10年,Gompertz,499999.99999999994,2122.2979572477343,183939.72058572114,0.9542285701514692,16.74935931143396
+近10年,Richards,23814.980027911413,2051.157710362777,15373.812459292129,0.9632499410027958,15.00824738946148
diff --git a/launch_distribution_mid.png b/launch_distribution_mid.png
new file mode 100644
index 0000000..e4db5c0
Binary files /dev/null and b/launch_distribution_mid.png differ
diff --git a/launch_distribution_min.png b/launch_distribution_min.png
new file mode 100644
index 0000000..7bd2818
Binary files /dev/null and b/launch_distribution_min.png differ
diff --git a/launch_frequency_analysis.png b/launch_frequency_analysis.png
new file mode 100644
index 0000000..07d231c
Binary files /dev/null and b/launch_frequency_analysis.png differ
diff --git a/launch_sites.csv b/launch_sites.csv
new file mode 100644
index 0000000..c32d038
--- /dev/null
+++ b/launch_sites.csv
@@ -0,0 +1,24 @@
+Ӣ,䳡,䳡ڹ,2023,2024,2025
+Florida,̫ջأάǿվأ,,59,71,82
+California,DZ̫ջأDZɽɽأ,,30,46,66
+Florida,Ϻ,,13,23,27
+Virginia,˹,,3,1,
+Texas,濨ǻ,,2,4,
+,ȪǷ䳡,й,36,21,33
+,Ƿ䳡,й,15,19,18
+Taiyuan Satellite Launch Center (China),̫ԭǷ䳡,й,4,13,
+,IJ췢䳡,й,9,8,
+,ҵ췢䳡,й,0,1,
+,㶫,й,2,2,
+,ɽ,й,1,4,
+,лĿ˺췢䳡,˹,7,5,
+,췢䳡,˹,3,4,
+Kazakhstan,ݿŬ췢䳡,˹̹,9,8,
+Mahia Peninsula (New Zealand),ŷŬվ,,7,13,
+Satish Dhawan Space Centre (India),ʲ,ӡ,7,5,
+French Guiana,Ǻ,,3,3,
+,ӵ,ձ,3,5,
+,,ձ,0,2,
+,Ƿ䳡,,3,1,
+,ɳ³º,,2,2,
+,ķϺ,,0,2,
diff --git a/rocket launch counts.txt b/rocket launch counts.txt
new file mode 100644
index 0000000..0ceaebe
--- /dev/null
+++ b/rocket launch counts.txt
@@ -0,0 +1,78 @@
+# Bin YDate ValMin ValMax US SU RU CN F J AU D UK I I-ELDO IN BR KR I-ESA KP IR IL CYM NZ Total
+# --- ---- ---- ---- (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) -----
+ 0 1956 1956 1957 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
+ 1 1957 1957 1958 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3
+ 2 1958 1958 1959 23 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 28
+ 3 1959 1959 1960 19 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 23
+ 4 1960 1960 1961 29 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 38
+ 5 1961 1961 1962 41 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 50
+ 6 1962 1962 1963 59 22 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 81
+ 7 1963 1963 1964 46 23 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 69
+ 8 1964 1964 1965 64 36 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 100
+ 9 1965 1965 1966 70 53 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 124
+ 10 1966 1966 1967 77 51 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 131
+ 11 1967 1967 1968 61 74 0 0 2 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 139
+ 12 1968 1968 1969 48 79 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 128
+ 13 1969 1969 1970 41 82 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 125
+ 14 1970 1970 1971 29 87 0 1 2 2 0 0 1 1 1 0 0 0 0 0 0 0 0 0 124
+ 15 1971 1971 1972 33 91 0 1 2 2 0 0 1 2 1 0 0 0 0 0 0 0 0 0 133
+ 16 1972 1972 1973 32 79 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 113
+ 17 1973 1973 1974 25 89 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 116
+ 18 1974 1974 1975 23 85 0 2 0 1 0 0 0 2 0 0 0 0 0 0 0 0 0 0 113
+ 19 1975 1975 1976 30 93 0 3 3 2 0 0 0 1 0 0 0 0 0 0 0 0 0 0 132
+ 20 1976 1976 1977 26 100 0 3 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 131
+ 21 1977 1977 1978 26 102 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 130
+ 22 1978 1978 1979 33 91 0 1 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 128
+ 23 1979 1979 1980 16 89 0 1 0 2 0 0 0 0 0 1 0 0 1 0 0 0 0 0 110
+ 24 1980 1980 1981 15 89 0 0 0 2 0 0 0 0 0 1 0 0 1 0 0 0 0 0 108
+ 25 1981 1981 1982 19 100 0 1 0 3 0 0 0 0 0 1 0 0 2 0 0 0 0 0 126
+ 26 1982 1982 1983 18 108 0 1 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 129
+ 27 1983 1983 1984 22 100 0 1 0 3 0 0 0 0 0 1 0 0 2 0 0 0 0 0 129
+ 28 1984 1984 1985 22 97 0 3 3 3 0 0 0 0 0 0 0 0 1 0 0 0 0 0 129
+ 29 1985 1985 1986 18 100 0 1 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 125
+ 30 1986 1986 1987 9 94 0 2 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 110
+ 31 1987 1987 1988 9 97 0 2 2 3 0 0 0 0 0 1 0 0 0 0 0 0 0 0 114
+ 32 1988 1988 1989 11 94 0 4 7 2 0 0 0 1 0 1 0 0 0 0 0 1 0 0 121
+ 33 1989 1989 1990 18 75 0 0 7 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 102
+ 34 1990 1990 1991 27 79 0 5 6 3 0 0 0 0 0 0 0 0 0 0 0 1 0 0 121
+ 35 1991 1991 1992 19 61 0 1 8 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 91
+ 36 1992 1992 1993 29 0 55 4 7 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 97
+ 37 1993 1993 1994 25 0 48 1 7 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 83
+ 38 1994 1994 1995 27 0 49 5 8 2 0 0 0 0 0 2 0 0 0 0 0 0 0 0 93
+ 39 1995 1995 1996 30 0 33 3 11 2 0 0 0 0 0 0 0 0 0 0 0 1 0 0 80
+ 40 1996 1996 1997 33 0 27 4 10 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 77
+ 41 1997 1997 1998 38 0 29 6 11 2 0 0 0 0 0 1 1 0 1 0 0 0 0 0 89
+ 42 1998 1998 1999 36 0 25 6 11 2 0 0 0 0 0 0 0 0 0 1 0 1 0 0 82
+ 43 1999 1999 2000 31 0 22 4 16 1 0 0 0 0 0 1 1 0 0 0 0 0 2 0 78
+ 44 2000 2000 2001 29 0 32 5 16 1 0 0 0 0 0 0 0 0 0 0 0 0 2 0 85
+ 45 2001 2001 2002 24 0 23 1 8 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 59
+ 46 2002 2002 2003 18 0 25 5 12 3 0 0 0 0 0 1 0 0 0 0 0 1 0 0 65
+ 47 2003 2003 2004 26 0 19 7 6 3 0 0 0 0 0 2 0 0 0 0 0 0 0 0 63
+ 48 2004 2004 2005 19 0 23 8 3 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 55
+ 49 2005 2005 2006 16 0 23 6 6 2 0 0 0 0 0 1 0 0 2 0 0 0 0 0 56
+ 50 2006 2006 2007 23 0 23 6 7 6 0 0 0 0 0 1 0 0 0 0 0 0 0 0 66
+ 51 2007 2007 2008 20 0 23 10 9 2 0 0 0 0 0 3 0 0 0 0 0 1 0 0 68
+ 52 2008 2008 2009 20 0 26 11 7 1 0 0 0 0 0 3 0 0 0 0 1 0 0 0 69
+ 53 2009 2009 2010 25 0 32 6 7 3 0 0 0 0 0 2 0 1 0 1 1 0 0 0 78
+ 54 2010 2010 2011 15 0 30 15 7 2 0 0 0 0 0 3 0 1 0 0 0 1 0 0 74
+ 55 2011 2011 2012 19 0 30 19 9 3 0 0 0 0 0 3 0 0 0 0 1 0 0 0 84
+ 56 2012 2012 2013 16 0 23 19 10 2 0 0 0 0 0 2 0 0 1 2 3 0 0 0 78
+ 57 2013 2013 2014 20 0 31 15 7 3 0 0 0 0 0 3 0 1 1 0 0 0 0 0 81
+ 58 2014 2014 2015 24 0 32 16 10 4 0 0 0 0 0 4 0 0 1 0 0 1 0 0 92
+ 59 2015 2015 2016 20 0 26 19 9 4 0 0 0 0 0 5 0 0 3 0 1 0 0 0 87
+ 60 2016 2016 2017 22 0 17 22 11 4 0 0 0 0 0 7 0 0 0 1 0 1 0 0 85
+ 61 2017 2017 2018 29 0 19 18 11 7 0 0 0 0 0 5 0 0 0 0 0 0 0 1 90
+ 62 2018 2018 2019 31 0 17 39 11 6 0 0 0 0 0 7 0 0 0 0 0 0 0 3 114
+ 63 2019 2019 2020 21 0 22 34 9 2 0 0 0 0 0 6 0 0 0 0 2 0 0 6 102
+ 64 2020 2020 2021 37 0 12 39 10 4 0 0 0 0 0 2 0 0 0 0 2 1 0 7 114
+ 65 2021 2021 2022 45 0 16 56 15 3 0 0 0 0 0 2 0 1 0 0 2 0 0 6 146
+ 66 2022 2022 2023 78 0 21 64 5 1 0 0 0 0 0 5 0 1 1 0 1 0 0 9 186
+ 67 2023 2023 2024 109 0 19 67 3 3 0 0 0 0 0 7 0 2 0 3 2 1 0 7 223
+ 68 2024 2024 2025 145 0 17 68 2 7 0 0 0 0 0 5 0 0 1 1 4 0 0 13 263
+ 69 2025 2025 2026 181 0 17 92 7 4 1 1 0 0 0 5 0 2 0 0 1 1 0 17 329
+ 70 2026 2026 2027 12 0 0 7 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2 22
+ 71 Total Total Total 2352 2449 886 741 350 142 1 1 2 9 4 101 2 9 20 9 21 13 4 71 7187
+# Bin YDate ValMin ValMax US SU RU CN F J AU D UK I I-ELDO IN BR KR I-ESA KP IR IL CYM NZ Total
+# --- ---- ---- ---- (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) -----
+# Data last updated: 2026 Jan 30 0451:24
+https://planet4589.org/space/stats/launches.html
\ No newline at end of file
diff --git a/rocket_comprehensive.png b/rocket_comprehensive.png
new file mode 100644
index 0000000..0af18fc
Binary files /dev/null and b/rocket_comprehensive.png differ
diff --git a/rocket_launch_counts.csv b/rocket_launch_counts.csv
new file mode 100644
index 0000000..d4a8721
--- /dev/null
+++ b/rocket_launch_counts.csv
@@ -0,0 +1,73 @@
+Bin,YDate,ValMin,ValMax,US,SU,RU,CN,F,J,AU,D,UK,I,I-ELDO,IN,BR,KR,I-ESA,KP,IR,IL,CYM,NZ,Total
+0,1956,1956,1957,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+1,1957,1957,1958,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3
+2,1958,1958,1959,23,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,28
+3,1959,1959,1960,19,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,23
+4,1960,1960,1961,29,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,38
+5,1961,1961,1962,41,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,50
+6,1962,1962,1963,59,22,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,81
+7,1963,1963,1964,46,23,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,69
+8,1964,1964,1965,64,36,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,100
+9,1965,1965,1966,70,53,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,124
+10,1966,1966,1967,77,51,0,0,1,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,131
+11,1967,1967,1968,61,74,0,0,2,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,139
+12,1968,1968,1969,48,79,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,128
+13,1969,1969,1970,41,82,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,125
+14,1970,1970,1971,29,87,0,1,2,2,0,0,1,1,1,0,0,0,0,0,0,0,0,0,124
+15,1971,1971,1972,33,91,0,1,2,2,0,0,1,2,1,0,0,0,0,0,0,0,0,0,133
+16,1972,1972,1973,32,79,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,113
+17,1973,1973,1974,25,89,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,116
+18,1974,1974,1975,23,85,0,2,0,1,0,0,0,2,0,0,0,0,0,0,0,0,0,0,113
+19,1975,1975,1976,30,93,0,3,3,2,0,0,0,1,0,0,0,0,0,0,0,0,0,0,132
+20,1976,1976,1977,26,100,0,3,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,131
+21,1977,1977,1978,26,102,0,0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,130
+22,1978,1978,1979,33,91,0,1,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,128
+23,1979,1979,1980,16,89,0,1,0,2,0,0,0,0,0,1,0,0,1,0,0,0,0,0,110
+24,1980,1980,1981,15,89,0,0,0,2,0,0,0,0,0,1,0,0,1,0,0,0,0,0,108
+25,1981,1981,1982,19,100,0,1,0,3,0,0,0,0,0,1,0,0,2,0,0,0,0,0,126
+26,1982,1982,1983,18,108,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,129
+27,1983,1983,1984,22,100,0,1,0,3,0,0,0,0,0,1,0,0,2,0,0,0,0,0,129
+28,1984,1984,1985,22,97,0,3,3,3,0,0,0,0,0,0,0,0,1,0,0,0,0,0,129
+29,1985,1985,1986,18,100,0,1,4,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,125
+30,1986,1986,1987,9,94,0,2,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,110
+31,1987,1987,1988,9,97,0,2,2,3,0,0,0,0,0,1,0,0,0,0,0,0,0,0,114
+32,1988,1988,1989,11,94,0,4,7,2,0,0,0,1,0,1,0,0,0,0,0,1,0,0,121
+33,1989,1989,1990,18,75,0,0,7,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,102
+34,1990,1990,1991,27,79,0,5,6,3,0,0,0,0,0,0,0,0,0,0,0,1,0,0,121
+35,1991,1991,1992,19,61,0,1,8,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,91
+36,1992,1992,1993,29,0,55,4,7,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,97
+37,1993,1993,1994,25,0,48,1,7,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,83
+38,1994,1994,1995,27,0,49,5,8,2,0,0,0,0,0,2,0,0,0,0,0,0,0,0,93
+39,1995,1995,1996,30,0,33,3,11,2,0,0,0,0,0,0,0,0,0,0,0,1,0,0,80
+40,1996,1996,1997,33,0,27,4,10,1,0,0,0,0,0,1,0,0,1,0,0,0,0,0,77
+41,1997,1997,1998,38,0,29,6,11,2,0,0,0,0,0,1,1,0,1,0,0,0,0,0,89
+42,1998,1998,1999,36,0,25,6,11,2,0,0,0,0,0,0,0,0,0,1,0,1,0,0,82
+43,1999,1999,2000,31,0,22,4,16,1,0,0,0,0,0,1,1,0,0,0,0,0,2,0,78
+44,2000,2000,2001,29,0,32,5,16,1,0,0,0,0,0,0,0,0,0,0,0,0,2,0,85
+45,2001,2001,2002,24,0,23,1,8,1,0,0,0,0,0,2,0,0,0,0,0,0,0,0,59
+46,2002,2002,2003,18,0,25,5,12,3,0,0,0,0,0,1,0,0,0,0,0,1,0,0,65
+47,2003,2003,2004,26,0,19,7,6,3,0,0,0,0,0,2,0,0,0,0,0,0,0,0,63
+48,2004,2004,2005,19,0,23,8,3,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,55
+49,2005,2005,2006,16,0,23,6,6,2,0,0,0,0,0,1,0,0,2,0,0,0,0,0,56
+50,2006,2006,2007,23,0,23,6,7,6,0,0,0,0,0,1,0,0,0,0,0,0,0,0,66
+51,2007,2007,2008,20,0,23,10,9,2,0,0,0,0,0,3,0,0,0,0,0,1,0,0,68
+52,2008,2008,2009,20,0,26,11,7,1,0,0,0,0,0,3,0,0,0,0,1,0,0,0,69
+53,2009,2009,2010,25,0,32,6,7,3,0,0,0,0,0,2,0,1,0,1,1,0,0,0,78
+54,2010,2010,2011,15,0,30,15,7,2,0,0,0,0,0,3,0,1,0,0,0,1,0,0,74
+55,2011,2011,2012,19,0,30,19,9,3,0,0,0,0,0,3,0,0,0,0,1,0,0,0,84
+56,2012,2012,2013,16,0,23,19,10,2,0,0,0,0,0,2,0,0,1,2,3,0,0,0,78
+57,2013,2013,2014,20,0,31,15,7,3,0,0,0,0,0,3,0,1,1,0,0,0,0,0,81
+58,2014,2014,2015,24,0,32,16,10,4,0,0,0,0,0,4,0,0,1,0,0,1,0,0,92
+59,2015,2015,2016,20,0,26,19,9,4,0,0,0,0,0,5,0,0,3,0,1,0,0,0,87
+60,2016,2016,2017,22,0,17,22,11,4,0,0,0,0,0,7,0,0,0,1,0,1,0,0,85
+61,2017,2017,2018,29,0,19,18,11,7,0,0,0,0,0,5,0,0,0,0,0,0,0,1,90
+62,2018,2018,2019,31,0,17,39,11,6,0,0,0,0,0,7,0,0,0,0,0,0,0,3,114
+63,2019,2019,2020,21,0,22,34,9,2,0,0,0,0,0,6,0,0,0,0,2,0,0,6,102
+64,2020,2020,2021,37,0,12,39,10,4,0,0,0,0,0,2,0,0,0,0,2,1,0,7,114
+65,2021,2021,2022,45,0,16,56,15,3,0,0,0,0,0,2,0,1,0,0,2,0,0,6,146
+66,2022,2022,2023,78,0,21,64,5,1,0,0,0,0,0,5,0,1,1,0,1,0,0,9,186
+67,2023,2023,2024,109,0,19,67,3,3,0,0,0,0,0,7,0,2,0,3,2,1,0,7,223
+68,2024,2024,2025,145,0,17,68,2,7,0,0,0,0,0,5,0,0,1,1,4,0,0,13,263
+69,2025,2025,2026,181,0,17,92,7,4,1,1,0,0,0,5,0,2,0,0,1,1,0,17,329
+70,2026,2026,2027,12,0,0,7,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,2,22
+71,Total,Total,Total,2352,2449,886,741,350,142,1,1,2,9,4,101,2,9,20,9,21,13,4,71,7187
diff --git a/rocket_launch_scenario.py b/rocket_launch_scenario.py
new file mode 100644
index 0000000..4e9ff3f
--- /dev/null
+++ b/rocket_launch_scenario.py
@@ -0,0 +1,678 @@
+"""
+Moon Colony Construction - Rocket Launch Scenario Analysis
+
+Problem: Deliver 100 million metric tons to Moon using traditional rockets
+- 10 launch sites, each max 1 launch per day
+- Rocket payload: 100-150 metric tons (use 125 tons average)
+- Optimize launch distribution to minimize energy consumption
+- Plot completion year vs total energy consumption
+
+Key insight: As completion timeline extends, more launches can be allocated
+to low-latitude sites, reducing the latitude-induced energy penalty.
+"""
+
+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, 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 = 100e9 kg 的物资，但单位是metric tons
+PAYLOAD_PER_LAUNCH = 125 # metric tons per launch (average of 100-150)
+TOTAL_LAUNCHES_NEEDED = TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH # 800,000 launches
+
+# 火箭参数
+ISP = 450 # 比冲 (秒) - 液氧液氢
+SPECIFIC_FUEL_ENERGY = 15.5e6 # J/kg 燃料比能量
+ALPHA = 0.10 # 结构系数
+NUM_STAGES = 3 # 3级火箭
+
+# 地面发射基准 ΔV (赤道发射到月球轨道)
+DELTA_V_BASE = 13300 # m/s (LEO + TLI + LOI)
+
+
+# ============== 发射场定义 ==============
+@dataclass
+class LaunchSite:
+ name: str
+ short_name: str
+ latitude: float # degrees
+ max_launches_per_day: int = 1
+
+ @property
+ def abs_latitude(self) -> float:
+ return abs(self.latitude)
+
+ @property
+ def rotation_velocity(self) -> float:
+ """地球自转速度 (m/s)"""
+ return OMEGA_EARTH * R_EARTH * np.cos(np.radians(self.abs_latitude))
+
+ @property
+ def delta_v_loss(self) -> float:
+ """相对赤道的ΔV损失 (m/s)"""
+ v_equator = OMEGA_EARTH * R_EARTH # ~465 m/s
+ return v_equator - self.rotation_velocity
+
+ @property
+ def total_delta_v(self) -> float:
+ """总ΔV需求 (m/s)"""
+ return DELTA_V_BASE + self.delta_v_loss
+
+
+# 10个发射场 (按纬度排序)
+LAUNCH_SITES = [
+ 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),
+]
+
+# 按纬度排序
+LAUNCH_SITES = sorted(LAUNCH_SITES, key=lambda x: x.abs_latitude)
+
+
+# ============== 核心计算函数 ==============
+
+def fuel_ratio_multistage(delta_v: float, isp: float = ISP, alpha: float = ALPHA,
+ num_stages: int = NUM_STAGES) -> float:
+ """
+ 多级火箭燃料/载荷比
+ """
+ ve = isp * G0
+ 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 energy_per_launch(site: LaunchSite, payload_tons: float = PAYLOAD_PER_LAUNCH) -> float:
+ """
+ 计算单次发射的能量消耗 (Joules)
+
+ 参数:
+ site: 发射场
+ payload_tons: 载荷 (metric tons)
+ 返回:
+ 能量消耗 (J)
+ """
+ delta_v = site.total_delta_v
+ k = fuel_ratio_multistage(delta_v)
+
+ payload_kg = payload_tons * 1000 # convert to kg
+ fuel_kg = k * payload_kg
+
+ energy = fuel_kg * SPECIFIC_FUEL_ENERGY
+ return energy
+
+
+def fuel_per_launch(site: LaunchSite, payload_tons: float = PAYLOAD_PER_LAUNCH) -> float:
+ """
+ 计算单次发射的燃料消耗 (metric tons)
+ """
+ delta_v = site.total_delta_v
+ k = fuel_ratio_multistage(delta_v)
+ return k * payload_tons
+
+
+# ============== 发射排布优化 ==============
+
+def calculate_launch_distribution(
+ total_launches: int,
+ completion_years: float,
+ sites: List[LaunchSite] = LAUNCH_SITES
+) -> List[Tuple[LaunchSite, int]]:
+ """
+ 计算给定完成年限下的发射分配
+
+ 策略：优先使用低纬度发射场，直到达到其最大容量
+
+ 参数:
+ total_launches: 总发射次数
+ completion_years: 完成年限
+ sites: 发射场列表 (已按纬度排序)
+
+ 返回:
+ [(发射场, 分配的发射次数), ...]
+ """
+ days_available = completion_years * 365
+ max_launches_per_site = days_available # 每站每天最多1发
+
+ distribution = []
+ remaining_launches = total_launches
+
+ # 按纬度从低到高分配
+ for site in sites:
+ if remaining_launches <= 0:
+ distribution.append((site, 0))
+ else:
+ allocated = min(remaining_launches, max_launches_per_site)
+ distribution.append((site, int(allocated)))
+ remaining_launches -= allocated
+
+ return distribution
+
+
+def calculate_total_energy(distribution: List[Tuple[LaunchSite, int]]) -> dict:
+ """
+ 计算发射分配的总能量消耗
+
+ 返回:
+ 包含各项统计的字典
+ """
+ total_energy = 0
+ total_fuel = 0
+ total_launches = 0
+ site_details = []
+
+ for site, num_launches in distribution:
+ if num_launches > 0:
+ energy = energy_per_launch(site) * num_launches
+ fuel = fuel_per_launch(site) * num_launches
+ total_energy += energy
+ total_fuel += fuel
+ total_launches += num_launches
+
+ site_details.append({
+ 'site': site.name,
+ 'short_name': site.short_name,
+ 'latitude': site.abs_latitude,
+ 'launches': num_launches,
+ 'energy_PJ': energy / 1e15, # PetaJoules
+ 'fuel_Mt': fuel / 1e6, # Million metric tons
+ })
+
+ return {
+ 'total_energy_PJ': total_energy / 1e15,
+ 'total_fuel_Mt': total_fuel / 1e6,
+ 'total_launches': total_launches,
+ 'sites': site_details
+ }
+
+
+def analyze_completion_timeline(
+ years_range: np.ndarray,
+ total_payload: float = TOTAL_PAYLOAD,
+ payload_per_launch: float = PAYLOAD_PER_LAUNCH
+) -> pd.DataFrame:
+ """
+ 分析不同完成年限下的能量消耗
+
+ 返回:
+ DataFrame with analysis results
+ """
+ total_launches = int(total_payload / payload_per_launch)
+
+ results = []
+ for years in years_range:
+ distribution = calculate_launch_distribution(total_launches, years)
+ stats = calculate_total_energy(distribution)
+
+ # 计算使用了多少个发射场
+ sites_used = sum(1 for _, n in distribution if n > 0)
+
+ # 计算低纬度发射场的比例
+ low_lat_launches = sum(n for site, n in distribution if site.abs_latitude < 30)
+ low_lat_ratio = low_lat_launches / total_launches if total_launches > 0 else 0
+
+ results.append({
+ 'years': years,
+ 'total_launches': total_launches,
+ 'sites_used': sites_used,
+ 'total_energy_PJ': stats['total_energy_PJ'],
+ 'total_fuel_Mt': stats['total_fuel_Mt'],
+ 'low_lat_ratio': low_lat_ratio,
+ 'energy_per_launch_TJ': stats['total_energy_PJ'] * 1000 / total_launches,
+ })
+
+ return pd.DataFrame(results)
+
+
+# ============== 可视化 ==============
+
+def plot_timeline_analysis(
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/rocket_launch_timeline.png'
+):
+ """
+ 绘制完成年份-能量消耗关系图
+ """
+ # 计算最短完成时间
+ total_launches = int(TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH)
+ num_sites = len(LAUNCH_SITES)
+ min_years = total_launches / (num_sites * 365)
+
+ print(f"Total payload: {TOTAL_PAYLOAD/1e6:.0f} million metric tons")
+ print(f"Payload per launch: {PAYLOAD_PER_LAUNCH} metric tons")
+ print(f"Total launches needed: {total_launches:,}")
+ print(f"Number of launch sites: {num_sites}")
+ print(f"Minimum completion time: {min_years:.1f} years")
+
+ # 分析范围：从最短时间到5倍最短时间
+ years_range = np.linspace(min_years, min_years * 5, 100)
+
+ # 运行分析
+ df = analyze_completion_timeline(years_range)
+
+ # 创建图表
+ fig, axes = plt.subplots(2, 2, figsize=(16, 14))
+
+ # ========== 图1: 完成年份 vs 总能量消耗 ==========
+ ax1 = axes[0, 0]
+ ax1.plot(df['years'], df['total_energy_PJ'], 'b-', linewidth=2)
+ ax1.fill_between(df['years'], df['total_energy_PJ'], alpha=0.3)
+
+ # 标记关键点
+ min_idx = df['total_energy_PJ'].idxmin()
+ ax1.axvline(x=df.loc[min_idx, 'years'], color='r', linestyle='--',
+ label=f'Optimal: {df.loc[min_idx, "years"]:.1f} years')
+ ax1.axhline(y=df.loc[min_idx, 'total_energy_PJ'], color='r', linestyle=':', alpha=0.5)
+
+ # 标记最短时间点
+ ax1.axvline(x=min_years, color='g', linestyle='--',
+ label=f'Minimum time: {min_years:.1f} years')
+
+ ax1.set_xlabel('Completion Time (years)', fontsize=12)
+ ax1.set_ylabel('Total Energy Consumption (PJ)', fontsize=12)
+ ax1.set_title('Completion Timeline vs Total Energy\n(100 Million Metric Tons to Moon)', fontsize=14)
+ ax1.legend()
+ ax1.grid(True, alpha=0.3)
+
+ # ========== 图2: 使用的发射场数量 ==========
+ ax2 = axes[0, 1]
+ ax2.plot(df['years'], df['sites_used'], 'g-', linewidth=2, marker='o',
+ markevery=10, markersize=6)
+
+ ax2.set_xlabel('Completion Time (years)', fontsize=12)
+ ax2.set_ylabel('Number of Launch Sites Used', fontsize=12)
+ ax2.set_title('Launch Sites Required vs Timeline', fontsize=14)
+ ax2.set_ylim(0, 11)
+ ax2.grid(True, alpha=0.3)
+
+ # 添加发射场名称标注
+ for i, site in enumerate(LAUNCH_SITES):
+ ax2.axhline(y=i+1, color='gray', linestyle=':', alpha=0.3)
+ ax2.text(df['years'].max() * 0.98, i+1 + 0.1, site.short_name,
+ fontsize=8, ha='right', va='bottom')
+
+ # ========== 图3: 低纬度发射比例 ==========
+ ax3 = axes[1, 0]
+ ax3.plot(df['years'], df['low_lat_ratio'] * 100, 'm-', linewidth=2)
+ ax3.fill_between(df['years'], df['low_lat_ratio'] * 100, alpha=0.3, color='magenta')
+
+ ax3.set_xlabel('Completion Time (years)', fontsize=12)
+ ax3.set_ylabel('Low Latitude (<30°) Launch Ratio (%)', fontsize=12)
+ ax3.set_title('Proportion of Launches from Low-Latitude Sites', fontsize=14)
+ ax3.set_ylim(0, 105)
+ ax3.grid(True, alpha=0.3)
+
+ # ========== 图4: 平均每发能量 ==========
+ ax4 = axes[1, 1]
+ ax4.plot(df['years'], df['energy_per_launch_TJ'], 'r-', linewidth=2)
+
+ # 计算赤道和最高纬度的能量参考线
+ energy_equator = energy_per_launch(LAUNCH_SITES[0]) / 1e12 # TJ
+ energy_max_lat = energy_per_launch(LAUNCH_SITES[-1]) / 1e12 # TJ
+
+ ax4.axhline(y=energy_equator, color='g', linestyle='--',
+ label=f'Equator baseline: {energy_equator:.1f} TJ')
+ ax4.axhline(y=energy_max_lat, color='orange', linestyle='--',
+ label=f'Alaska (57°): {energy_max_lat:.1f} TJ')
+
+ ax4.set_xlabel('Completion Time (years)', fontsize=12)
+ ax4.set_ylabel('Average Energy per Launch (TJ)', fontsize=12)
+ ax4.set_title('Energy Efficiency vs Timeline\n(Longer timeline → more low-lat launches)', fontsize=14)
+ ax4.legend()
+ ax4.grid(True, alpha=0.3)
+
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"\n图表已保存至: {save_path}")
+
+ return df, fig
+
+
+def plot_launch_distribution(
+ years: float,
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/launch_distribution.png'
+):
+ """
+ 绘制特定年限下的发射分配详情
+ """
+ total_launches = int(TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH)
+ distribution = calculate_launch_distribution(total_launches, years)
+ stats = calculate_total_energy(distribution)
+
+ fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+
+ # ========== 图1: 发射次数分配 ==========
+ ax1 = axes[0]
+ sites = [s.short_name for s, _ in distribution]
+ launches = [n for _, n in distribution]
+ latitudes = [s.abs_latitude for s, _ in distribution]
+
+ colors = plt.cm.RdYlGn_r(np.array(latitudes) / 60)
+ bars = ax1.barh(sites, launches, color=colors)
+
+ ax1.set_xlabel('Number of Launches', fontsize=12)
+ ax1.set_ylabel('Launch Site', fontsize=12)
+ ax1.set_title(f'Launch Distribution ({years:.1f} years completion)\n'
+ f'Total: {total_launches:,} launches', fontsize=13)
+ ax1.grid(True, alpha=0.3, axis='x')
+
+ # 添加数值标签
+ for bar, n in zip(bars, launches):
+ if n > 0:
+ ax1.text(bar.get_width() + 1000, bar.get_y() + bar.get_height()/2,
+ f'{n:,}', va='center', fontsize=9)
+
+ # ========== 图2: 能量消耗分配 ==========
+ ax2 = axes[1]
+
+ site_names = [d['short_name'] for d in stats['sites']]
+ energies = [d['energy_PJ'] for d in stats['sites']]
+
+ colors = plt.cm.RdYlGn_r(np.array([d['latitude'] for d in stats['sites']]) / 60)
+ bars = ax2.barh(site_names, energies, color=colors)
+
+ ax2.set_xlabel('Energy Consumption (PJ)', fontsize=12)
+ ax2.set_ylabel('Launch Site', fontsize=12)
+ ax2.set_title(f'Energy Distribution\n'
+ f'Total: {stats["total_energy_PJ"]:.1f} PJ', fontsize=13)
+ ax2.grid(True, alpha=0.3, axis='x')
+
+ # 添加数值标签
+ for bar, e in zip(bars, energies):
+ if e > 0:
+ ax2.text(bar.get_width() + 0.5, bar.get_y() + bar.get_height()/2,
+ f'{e:.1f}', va='center', fontsize=9)
+
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"分配详情图已保存至: {save_path}")
+
+ return distribution, stats
+
+
+def print_analysis_summary(df: pd.DataFrame):
+ """打印分析摘要"""
+ print("\n" + "=" * 80)
+ print("ROCKET LAUNCH SCENARIO ANALYSIS SUMMARY")
+ print("=" * 80)
+
+ total_launches = int(TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH)
+
+ print(f"\nMission Parameters:")
+ print(f" - Total payload: {TOTAL_PAYLOAD/1e6:.0f} million metric tons")
+ print(f" - Payload per launch: {PAYLOAD_PER_LAUNCH} metric tons")
+ print(f" - Total launches required: {total_launches:,}")
+ print(f" - Number of launch sites: {len(LAUNCH_SITES)}")
+
+ # 最短时间场景
+ min_years = df['years'].min()
+ min_time_row = df[df['years'] == min_years].iloc[0]
+ print(f"\nMinimum Time Scenario ({min_years:.1f} years):")
+ print(f" - All {int(min_time_row['sites_used'])} sites at full capacity")
+ print(f" - Total energy: {min_time_row['total_energy_PJ']:.1f} PJ")
+ print(f" - Total fuel: {min_time_row['total_fuel_Mt']:.1f} million metric tons")
+ print(f" - Low-latitude ratio: {min_time_row['low_lat_ratio']*100:.1f}%")
+
+ # 最优能量场景 (只使用最低纬度站点)
+ max_years = df['years'].max()
+ max_time_row = df[df['years'] == max_years].iloc[0]
+ print(f"\nExtended Timeline Scenario ({max_years:.1f} years):")
+ print(f" - Using {int(max_time_row['sites_used'])} lowest-latitude sites")
+ print(f" - Total energy: {max_time_row['total_energy_PJ']:.1f} PJ")
+ print(f" - Total fuel: {max_time_row['total_fuel_Mt']:.1f} million metric tons")
+ print(f" - Low-latitude ratio: {max_time_row['low_lat_ratio']*100:.1f}%")
+
+ # 能量节省
+ energy_saving = (1 - max_time_row['total_energy_PJ'] / min_time_row['total_energy_PJ']) * 100
+ print(f"\nEnergy Savings (Extended vs Minimum time):")
+ print(f" - Energy reduction: {energy_saving:.1f}%")
+ print(f" - Time trade-off: {max_years - min_years:.1f} additional years")
+
+ print("\n" + "=" * 80)
+ print("Launch Site Details:")
+ print("=" * 80)
+ print(f"\n{'Site':<25} {'Latitude':>10} {'ΔV Loss':>10} {'Fuel Ratio':>12} {'Energy/Launch':>14}")
+ print(f"{'':25} {'(°)':>10} {'(m/s)':>10} {'':>12} {'(TJ)':>14}")
+ print("-" * 80)
+
+ for site in LAUNCH_SITES:
+ k = fuel_ratio_multistage(site.total_delta_v)
+ e = energy_per_launch(site) / 1e12
+ print(f"{site.name:<25} {site.abs_latitude:>10.1f} {site.delta_v_loss:>10.0f} {k:>12.1f} {e:>14.1f}")
+
+ print("=" * 80)
+
+
+def analyze_launch_frequency_scenarios(
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/launch_frequency_analysis.png'
+):
+ """
+ 分析不同发射频率场景下的完成时间和能量消耗
+ """
+ total_launches = int(TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH)
+ num_sites = len(LAUNCH_SITES)
+
+ # 不同发射频率场景
+ frequencies = [1, 2, 4, 8, 12, 24] # launches per day per site
+
+ results = []
+ for freq in frequencies:
+ # 最短完成时间
+ min_years = total_launches / (num_sites * freq * 365)
+
+ # 分析该时间下的能量（所有站点满负荷）
+ distribution = calculate_launch_distribution(total_launches, min_years)
+ stats = calculate_total_energy(distribution)
+
+ results.append({
+ 'frequency': freq,
+ 'min_years': min_years,
+ 'total_energy_PJ': stats['total_energy_PJ'],
+ 'total_fuel_Mt': stats['total_fuel_Mt'],
+ })
+
+ df = pd.DataFrame(results)
+
+ # 绘图
+ fig, axes = plt.subplots(1, 2, figsize=(14, 6))
+
+ # 图1: 发射频率 vs 完成时间
+ ax1 = axes[0]
+ ax1.bar(range(len(frequencies)), df['min_years'], color='steelblue', alpha=0.8)
+ ax1.set_xticks(range(len(frequencies)))
+ ax1.set_xticklabels([f'{f}/day' for f in frequencies])
+ ax1.set_xlabel('Launch Frequency (per site)', fontsize=12)
+ ax1.set_ylabel('Minimum Completion Time (years)', fontsize=12)
+ ax1.set_title('Completion Time vs Launch Frequency\n(10 sites, 800,000 launches)', fontsize=13)
+ ax1.grid(True, alpha=0.3, axis='y')
+
+ # 添加数值标签
+ for i, (freq, years) in enumerate(zip(frequencies, df['min_years'])):
+ ax1.text(i, years + 2, f'{years:.1f}y', ha='center', fontsize=10, fontweight='bold')
+
+ # 图2: 能量消耗比较
+ ax2 = axes[1]
+ ax2.bar(range(len(frequencies)), df['total_energy_PJ'], color='coral', alpha=0.8)
+ ax2.set_xticks(range(len(frequencies)))
+ ax2.set_xticklabels([f'{f}/day' for f in frequencies])
+ ax2.set_xlabel('Launch Frequency (per site)', fontsize=12)
+ ax2.set_ylabel('Total Energy (PJ)', fontsize=12)
+ ax2.set_title('Total Energy Consumption\n(Same for all frequencies at max capacity)', fontsize=13)
+ ax2.grid(True, alpha=0.3, axis='y')
+
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"\n发射频率分析图已保存至: {save_path}")
+
+ return df
+
+
+def plot_comprehensive_analysis(
+ save_path: str = '/Volumes/Files/code/mm/20260130_b/rocket_comprehensive.png'
+):
+ """
+ 综合分析图：不同时间线下的能量消耗和发射分配
+ """
+ total_launches = int(TOTAL_PAYLOAD / PAYLOAD_PER_LAUNCH)
+ num_sites = len(LAUNCH_SITES)
+ min_years = total_launches / (num_sites * 365)
+
+ # 选择几个关键时间点
+ key_years = [min_years, min_years*1.5, min_years*2, min_years*3, min_years*5]
+
+ fig, axes = plt.subplots(2, 3, figsize=(18, 12))
+
+ # 第一行：能量趋势
+ ax_main = axes[0, 0]
+ years_range = np.linspace(min_years, min_years * 5, 100)
+ df = analyze_completion_timeline(years_range)
+
+ ax_main.plot(df['years'], df['total_energy_PJ'], 'b-', linewidth=2)
+ ax_main.fill_between(df['years'], df['total_energy_PJ'], alpha=0.3)
+
+ # 标记关键点
+ for y in key_years[1:]:
+ ax_main.axvline(x=y, color='gray', linestyle='--', alpha=0.5)
+
+ ax_main.set_xlabel('Completion Time (years)', fontsize=11)
+ ax_main.set_ylabel('Total Energy (PJ)', fontsize=11)
+ ax_main.set_title('Energy vs Timeline\n(Longer time → Lower energy)', fontsize=12)
+ ax_main.grid(True, alpha=0.3)
+
+ # 能量效率
+ ax_eff = axes[0, 1]
+ ax_eff.plot(df['years'], df['energy_per_launch_TJ'], 'r-', linewidth=2)
+ ax_eff.set_xlabel('Completion Time (years)', fontsize=11)
+ ax_eff.set_ylabel('Energy per Launch (TJ)', fontsize=11)
+ ax_eff.set_title('Average Energy Efficiency', fontsize=12)
+ ax_eff.grid(True, alpha=0.3)
+
+ # 发射场使用数量
+ ax_sites = axes[0, 2]
+ ax_sites.plot(df['years'], df['sites_used'], 'g-', linewidth=2)
+ ax_sites.set_xlabel('Completion Time (years)', fontsize=11)
+ ax_sites.set_ylabel('Sites Used', fontsize=11)
+ ax_sites.set_title('Number of Active Sites', fontsize=12)
+ ax_sites.set_ylim(0, 11)
+ ax_sites.grid(True, alpha=0.3)
+
+ # 第二行：三个时间点的发射分配饼图
+ scenario_years = [min_years, min_years*2, min_years*5]
+ scenario_names = ['Minimum Time', '2x Minimum', '5x Minimum']
+
+ for idx, (years, name) in enumerate(zip(scenario_years, scenario_names)):
+ ax = axes[1, idx]
+ distribution = calculate_launch_distribution(total_launches, years)
+
+ # 只显示有发射的站点
+ labels = [s.short_name for s, n in distribution if n > 0]
+ sizes = [n for _, n in distribution if n > 0]
+ latitudes = [s.abs_latitude for s, n in distribution if n > 0]
+
+ colors = plt.cm.RdYlGn_r(np.array(latitudes) / 60)
+
+ wedges, texts, autotexts = ax.pie(sizes, labels=labels, autopct='%1.0f%%',
+ colors=colors, pctdistance=0.75)
+
+ # 计算该场景的能量
+ stats = calculate_total_energy(distribution)
+
+ ax.set_title(f'{name}\n({years:.0f} years, {stats["total_energy_PJ"]:.0f} PJ)',
+ fontsize=11)
+
+ plt.suptitle(f'Rocket Launch Scenario: {TOTAL_PAYLOAD/1e6:.0f}M tons to Moon\n'
+ f'({total_launches:,} launches needed)', fontsize=14, y=1.02)
+
+ plt.tight_layout()
+ plt.savefig(save_path, dpi=150, bbox_inches='tight')
+ print(f"\n综合分析图已保存至: {save_path}")
+
+ return fig
+
+
+# ============== 主程序 ==============
+
+if __name__ == "__main__":
+ print("=" * 80)
+ print("MOON COLONY ROCKET LAUNCH SCENARIO ANALYSIS")
+ print("Mission: Deliver 100 million metric tons to Moon")
+ print("Constraint: 10 launch sites, max 1 launch/day each")
+ print("=" * 80)
+
+ # 运行时间线分析
+ df, fig = plot_timeline_analysis()
+
+ # 打印摘要
+ print_analysis_summary(df)
+
+ # 绘制综合分析图
+ plot_comprehensive_analysis()
+
+ # 发射频率场景分析
+ print("\n" + "=" * 80)
+ print("LAUNCH FREQUENCY SENSITIVITY ANALYSIS")
+ print("=" * 80)
+ freq_df = analyze_launch_frequency_scenarios()
+ print("\nIf launch frequency can be increased:")
+ print(f"{'Frequency':>15} {'Min Years':>15} {'Energy (PJ)':>15}")
+ print("-" * 50)
+ for _, row in freq_df.iterrows():
+ print(f"{row['frequency']:>12}/day {row['min_years']:>15.1f} {row['total_energy_PJ']:>15.0f}")
+
+ # 关键结论
+ print("\n" + "=" * 80)
+ print("KEY FINDINGS")
+ print("=" * 80)
+ print(f"""
+ 1. MINIMUM COMPLETION TIME: {df['years'].min():.0f} years
+ (All 10 sites at 1 launch/day capacity)
+
+ 2. TOTAL ENERGY REQUIRED: ~{df['total_energy_PJ'].mean():.0f} PJ
+ (~{df['total_fuel_Mt'].mean():.0f} million metric tons of fuel)
+
+ 3. LATITUDE EFFECT:
+ - Energy savings from using only low-lat sites: ~2-3%
+ - Trade-off: Significantly longer timeline
+
+ 4. CONCLUSION:
+ Traditional rockets alone are IMPRACTICAL for this mission.
+ 219+ years is far beyond reasonable project timeline.
+ Space Elevator System is essential for feasibility.
+ """)
diff --git a/rocket_launch_timeline.png b/rocket_launch_timeline.png
new file mode 100644
index 0000000..64a8a46
Binary files /dev/null and b/rocket_launch_timeline.png differ
diff --git a/rocket_launch_to_csv.py b/rocket_launch_to_csv.py
new file mode 100644
index 0000000..44aa784
--- /dev/null
+++ b/rocket_launch_to_csv.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+"""将 rocket launch counts.txt 数据导出为 CSV 文件"""
+
+import csv
+import re
+
+
+def parse_header(header_line: str) -> list[str]:
+ """从 # 开头的表头行解析列名"""
+ # 去掉开头的 # 和多余空格，按空白分割
+ clean = header_line.lstrip("#").strip()
+ # 列名之间可能有多空格，统一按空白分割
+ return re.split(r"\s{2,}", clean)
+
+
+def parse_data_line(line: str) -> list[str]:
+ """解析数据行，按空白分割"""
+ return line.split()
+
+
+def main():
+ input_file = "rocket launch counts.txt"
+ output_file = "rocket_launch_counts.csv"
+
+ headers = None
+ rows = []
+
+ with open(input_file, "r", encoding="utf-8") as f:
+ for line in f:
+ line = line.rstrip("\n")
+ if not line.strip():
+ continue
+ if line.startswith("#"):
+ # 第二行是分隔线，第一行是列名
+ if "Bin" in line and "YDate" in line and headers is None:
+ headers = parse_header(line)
+ continue
+ # 数据行
+ values = parse_data_line(line)
+ if len(values) >= 4: # 至少包含 Bin, YDate, ValMin, ValMax
+ rows.append(values)
+
+ # 如果没有成功解析表头，使用默认列名
+ if headers is None or len(headers) != len(rows[0]) if rows else True:
+ headers = [
+ "Bin", "YDate", "ValMin", "ValMax",
+ "US", "SU", "RU", "CN", "F", "J", "AU", "D", "UK", "I",
+ "I-ELDO", "IN", "BR", "KR", "I-ESA", "KP", "IR", "IL", "CYM", "NZ",
+ "Total"
+ ]
+
+ with open(output_file, "w", encoding="utf-8", newline="") as out:
+ writer = csv.writer(out)
+ writer.writerow(headers)
+ writer.writerows(rows)
+
+ print(f"已导出 {len(rows)} 行数据到 {output_file}")
+
+
+if __name__ == "__main__":
+ main()
diff --git a/scenario_comparison.png b/scenario_comparison.png
new file mode 100644
index 0000000..d75ecfb
Binary files /dev/null and b/scenario_comparison.png differ
diff --git a/three_scenarios_comparison.png b/three_scenarios_comparison.png
new file mode 100644
index 0000000..ff8b9f3
Binary files /dev/null and b/three_scenarios_comparison.png differ