commit 53521059227568f5a90adaee8f3297c7e023bde6 Author: ntnt Date: Sat Jan 31 00:26:34 2026 +0800 first commit diff --git a/earth_moon_transfer.py b/earth_moon_transfer.py new file mode 100644 index 0000000..5977f7a --- /dev/null +++ b/earth_moon_transfer.py @@ -0,0 +1,272 @@ +""" +地月转移轨道：有效载荷与燃料消耗（能量）关系分析 + +基于齐奥尔科夫斯基火箭方程和霍曼转移轨道模型 +""" + +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/earth_moon_transfer_analysis.png b/earth_moon_transfer_analysis.png new file mode 100644 index 0000000..7ed8202 Binary files /dev/null and b/earth_moon_transfer_analysis.png differ diff --git a/latitude_effects.png b/latitude_effects.png new file mode 100644 index 0000000..3e6456d Binary files /dev/null and b/latitude_effects.png differ diff --git a/prob/en_B.md b/prob/en_B.md new file mode 100644 index 0000000..f897388 --- /dev/null +++ b/prob/en_B.md @@ -0,0 +1,59 @@ +# 2026 MCM + +# Problem B: Creating a Moon Colony Using a Space Elevator System + +![](images/d56c8f8a71b712570aacdb699e2e19450e8b1a708748d394f92cbee59747b814.jpg) + +Imagine a future where it's possible for anyone to visit space by taking a leisurely and scenic ride from the Equator to Earth's orbit and then catching a routine, safe, and inexpensive rocket flight to the Moon, Mars, or beyond. In this future, we could build lush, green, and beautiful space habitats with artificial gravity, where people would vacation, work, or even live. These habitats would alleviate pressure on Earth's delicate, overworked, and fragile ecosystems. The technology to enable these events would provide humankind with limitless, safe, routine, environmentally friendly, efficient, and global access to space. To achieve these goals, some people envision a Space Elevator System, powered by electricity, offering a scalable infrastructure for interplanetary logistics, commerce, and exploration. + +At its final operating configuration, the Space Elevator System would comprise three Galactic Harbours, ideally separated by 120 degrees around the equator. Each Galactic Harbour would include a single Earth port with two $100,000\mathrm{km}$ -long tethers connected to two apex anchors, with multiple space elevators operating together, each capable of lifting massive payloads daily from Earth to geosynchronous orbit (GEO) and beyond to the apex anchor where they can be loaded on a rocket and delivered anywhere using much less fuel. + +The Moon Colony Management (MCM) Agency is preparing to build a Moon Colony with an estimated 100,000 people beginning in the year 2050, after completion of the Space Elevator System. It is estimated that the Moon Colony will need about 100 million metric tons of materials. Additionally, water and supplies will routinely need to be sent to sustain the Moon's population once the colony is complete. To get to the Moon, the Galactic Harbour must send material in two steps: first, from the Earth port to the apex anchor via a space elevator, and second, from the apex anchor to the Moon Colony via a rocket. The MCM Agency anticipates that the Galactic Harbour will provide an advanced lift system capable of moving 179,000 metric tons every year, while generating no atmospheric pollution. + +The agency is also considering using traditional rockets to supply material for construction and supplies to the Moon Colony. The Earth current has ten rocket launch sites: Alaska, California, Texas, Florida, and Virginia (United States), Kazakhstan, French Guiana, Satish Dhawan Space Centre (India), Taiyuan Satellite Launch Center (China), and Mahia Peninsula (New Zealand). + +A rocket would require a single step from the rocket launch site on Earth to the Moon Colony. By 2050 it is estimated that rockets will be able to carry 100-150 metric tons of payload to the Moon using advanced Falcon Heavy launches. You may assume perfect conditions for both the Galactic Harbour system (e.g., no swaying of the tether) and rocket launches (e.g., no failed launches). You should consider the cost and timeline to deliver the materials from the surface of the Earth to the Moon Colony site for the different scenarios. + +# Your Task: + +Your task is to utilize a mathematical model to determine the cost and associated timeline in order to transport material to build a 100,000 person Moon Colony starting in 2050. You will need to compare the Modern-Day Space Elevator System's three Galactic Harbours to traditional rockets launched from selected rocket bases. + +# Your model should include: + +1. Consideration of three different scenarios for how the 100 million metric tons of materials will be delivered to build the 100,000-person Moon Colony; + +a. using the Space Elevator System's three Galactic Harbor's alone, +b. traditional rocket launches from existing bases alone (you may choose which facilities to use), or, +c. some combination of the two methods. + +2. To what extent does your solution(s) change if the transportation systems are not in perfect working order (e.g, swaying of the tether, rockets fail, elevators break, etc.).? +3. Investigate the water needs for a one-year period once the 100,000-person Moon Colony is fully operational. Use your delivery model to understand the additional cost and timeline needed to ensure the colony has sufficient water for one full year after the Moon Colony is inhabited. +4. Discuss the impact on the Earth's environment for achieving the 100,000-person Moon Colony under the different scenarios. How would you adjust your model to minimize the environmental impact? +5. Write a one-page letter recommending a course of action to the fictional MCM Agency to build and sustain a 100,000-person Moon Colony. + +Your PDF solution of no more than 25 total pages should include: + +One-page Summary Sheet. +Table of Contents. +- Your complete solution. +One-page letter to MCM Agency +References list. +AI Use Report (If used does not count toward the 25-page limit.) + +Note: There is no specific required minimum page length for a complete MCM submission. You may use up to 25 total pages for all your solution work and any additional information you want to include (for example: drawings, diagrams, calculations, tables). Partial solutions are accepted. We permit the careful use of AI such as ChatGPT, although it is not necessary to create a solution to this problem. If you choose to utilize a generative AI, you must follow the COMAP AI use policy. This will result in an additional AI use report that you must add to the end of your PDF solution file and does not count toward the 25 total page limit for your solution. + +# Glossary + +Space Elevator System is comprised of three Galactic Harbours plus additional support facilities. + +Galactic Harbour is comprised of two apex anchors each connected by two tethers to a single Earth Port. + +Earth Port is the location on the Earth that provides surface support for the Galactic Harbour. + +Tethers are $100,000\mathrm{km}$ long graphene material that links the Earth port and apex anchors in the Space Elevator System. + +Apex Anchor is the counterweight in space at the end of the $100,000\mathrm{km}$ tether. + +Geosynchronous orbit (GEO) is approximately $35,786\mathrm{km}$ above the surface of the Earth where the orbital period to circle Earth is 24 hours, matching Earth's rotation so it stays over the same longitude each day. + +Moon Colony is a habitat on the moon with the capacity to support 100,000-persons. \ No newline at end of file diff --git a/specific_energy_comparison.png b/specific_energy_comparison.png new file mode 100644 index 0000000..eb71a39 Binary files /dev/null and b/specific_energy_comparison.png differ diff --git a/specific_energy_comparison.py b/specific_energy_comparison.py new file mode 100644 index 0000000..43e8823 --- /dev/null +++ b/specific_energy_comparison.py @@ -0,0 +1,1008 @@ +""" +地面发射 vs 太空电梯发射：比能量计算与对比 + +可配置参数进行不同场景的能量分析 +""" + +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 Optional + +# 设置中文字体支持 +rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] +rcParams['axes.unicode_minus'] = False + + +# ============== 物理常数 ============== +@dataclass +class PhysicalConstants: + """物理常数""" + G: float = 6.674e-11 # 万有引力常数 (m³/kg/s²) + M_earth: float = 5.972e24 # 地球质量 (kg) + M_moon: float = 7.342e22 # 月球质量 (kg) + R_earth: float = 6.371e6 # 地球半径 (m) + R_moon: float = 1.737e6 # 月球半径 (m) + g0: float = 9.81 # 地表重力加速度 (m/s²) + d_earth_moon: float = 3.844e8 # 地月平均距离 (m) + omega_earth: float = 7.27e-5 # 地球自转角速度 (rad/s) + + @property + def GM_earth(self) -> float: + """地球引力参数""" + return self.G * self.M_earth + + @property + def GM_moon(self) -> float: + """月球引力参数""" + return self.G * self.M_moon + + +# ============== 发动机参数 ============== +@dataclass +class EngineParams: + """发动机/推进系统参数""" + name: str = "液氧液氢" + isp: float = 450 # 比冲 (秒) + specific_energy: float = 15.5e6 # 燃料比能量 (J/kg) + + @property + def exhaust_velocity(self) -> float: + """排气速度 (m/s)""" + return self.isp * 9.81 + + +# ============== 发射场参数 ============== +@dataclass +class LaunchSite: + """发射场参数""" + name: str = "赤道发射场" + latitude: float = 0.0 # 纬度 (度) + + # 常见发射场纬度参考: + # 赤道: 0° + # 文昌 (中国): 19.6° + # 卡纳维拉尔角 (美国): 28.5° + # 拜科努尔 (哈萨克斯坦): 45.6° + # 酒泉 (中国): 40.9° + # 普列谢茨克 (俄罗斯): 62.9° + + @property + def latitude_rad(self) -> float: + """纬度 (弧度)""" + return np.radians(self.latitude) + + +# 预定义的发射场 +LAUNCH_SITES = { + # 赤道参考 + '赤道': LaunchSite("Equator (参考)", 0.0), + + # 法属圭亚那 (接近赤道) + 'Kourou': LaunchSite("Kourou (French Guiana)", 5.2), + + # 印度 + 'SDSC': LaunchSite("Satish Dhawan (India)", 13.7), + + # 美国 + 'Texas': LaunchSite("Boca Chica (Texas)", 26.0), + 'Florida': LaunchSite("Cape Canaveral (Florida)", 28.5), + 'California': LaunchSite("Vandenberg (California)", 34.7), + 'Virginia': LaunchSite("Wallops (Virginia)", 37.8), + 'Alaska': LaunchSite("Kodiak (Alaska)", 57.4), + + # 中国 + 'Taiyuan': LaunchSite("Taiyuan (China)", 38.8), + + # 新西兰 + 'Mahia': LaunchSite("Mahia (New Zealand)", -39.3), # 南半球 + + # 哈萨克斯坦 + 'Baikonur': LaunchSite("Baikonur (Kazakhstan)", 45.6), +} + + +# ============== 任务参数 ============== +@dataclass +class MissionParams: + """任务参数""" + # 结构系数 + alpha: float = 0.10 + + # 目标轨道 (环月轨道高度) + lunar_orbit_altitude: float = 100e3 # 100 km + + # 地面发射参数 (赤道发射基准值) + leo_altitude: float = 400e3 # LEO高度 400 km + delta_v_ground_to_leo_base: float = 9400 # 赤道发射到LEO的基准ΔV (m/s) + delta_v_tli: float = 3100 # 地月转移注入 (m/s) + delta_v_loi: float = 800 # 月球轨道捕获 (m/s) + + # 目标轨道倾角 (度), 0 = 赤道轨道 + target_inclination: float = 0.0 + + # 太空电梯参数 + elevator_length: float = 1.0e8 # 电梯长度 (m), 默认10万km + + @property + def total_delta_v_ground(self) -> float: + """地面发射总ΔV (赤道发射基准)""" + return self.delta_v_ground_to_leo_base + self.delta_v_tli + self.delta_v_loi + + +# ============== 纬度相关计算函数 ============== + +def earth_rotation_velocity(latitude: float, constants: PhysicalConstants = PhysicalConstants()) -> float: + """ + 计算地球表面某纬度的自转线速度 + + 参数: + latitude: 纬度 (度) + constants: 物理常数 + 返回: + 自转线速度 (m/s) + """ + lat_rad = np.radians(latitude) + return constants.omega_earth * constants.R_earth * np.cos(lat_rad) + + +def delta_v_rotation_loss(latitude: float, constants: PhysicalConstants = PhysicalConstants()) -> float: + """ + 计算相对赤道发射损失的自转速度贡献 + + 参数: + latitude: 发射场纬度 (度) + constants: 物理常数 + 返回: + 损失的ΔV (m/s) + """ + v_equator = earth_rotation_velocity(0, constants) # 赤道速度 ~465 m/s + v_latitude = earth_rotation_velocity(latitude, constants) + return v_equator - v_latitude + + +def delta_v_inclination_change( + launch_latitude: float, + target_inclination: float, + orbital_velocity: float = 7800 # LEO轨道速度 (m/s) +) -> float: + """ + 计算轨道倾角改变所需的ΔV + + 从发射场纬度的轨道转到目标倾角轨道 + 注意: 最小轨道倾角 = 发射场纬度 + + 参数: + launch_latitude: 发射场纬度 (度) + target_inclination: 目标轨道倾角 (度) + orbital_velocity: 轨道速度 (m/s) + 返回: + 倾角改变ΔV (m/s) + """ + # 发射场纬度决定最小可达倾角 + min_inclination = abs(launch_latitude) + + if target_inclination < min_inclination: + # 需要进行纯倾角改变机动 (非常昂贵!) + # ΔV = 2 * v * sin(Δi/2) + delta_i = np.radians(min_inclination - target_inclination) + return 2 * orbital_velocity * np.sin(delta_i / 2) + else: + # 可以通过发射方位角直接进入目标倾角，无额外ΔV + return 0 + + +def total_delta_v_with_latitude( + launch_site: LaunchSite, + mission: MissionParams, + constants: PhysicalConstants = PhysicalConstants() +) -> dict: + """ + 计算考虑纬度因素的总ΔV + + 参数: + launch_site: 发射场 + mission: 任务参数 + constants: 物理常数 + 返回: + 包含各项ΔV的字典 + """ + # 1. 自转速度损失 + dv_rotation_loss = delta_v_rotation_loss(launch_site.latitude, constants) + + # 2. 倾角改变损失 (如果目标是低倾角轨道) + dv_inclination = delta_v_inclination_change( + launch_site.latitude, + mission.target_inclination, + orbital_velocity=7800 + ) + + # 3. 总地面到LEO的ΔV + dv_to_leo = mission.delta_v_ground_to_leo_base + dv_rotation_loss + + # 4. 总ΔV + dv_total = dv_to_leo + dv_inclination + mission.delta_v_tli + mission.delta_v_loi + + return { + 'launch_site': launch_site.name, + 'latitude': launch_site.latitude, + 'rotation_velocity': earth_rotation_velocity(launch_site.latitude, constants), + 'dv_rotation_loss': dv_rotation_loss, + 'dv_inclination_change': dv_inclination, + 'dv_to_leo': dv_to_leo, + 'dv_tli': mission.delta_v_tli, + 'dv_loi': mission.delta_v_loi, + 'dv_total': dv_total + } + + +# ============== 核心计算函数 ============== + +def mass_ratio(delta_v: float, engine: EngineParams) -> float: + """ + 计算质量比 R = m0/mf + + 参数: + delta_v: 速度增量 (m/s) + engine: 发动机参数 + 返回: + 质量比 R + """ + return np.exp(delta_v / engine.exhaust_velocity) + + +def fuel_to_payload_ratio(delta_v: float, engine: EngineParams, alpha: float = 0.1) -> float: + """ + 计算燃料与载荷的质量比（单级火箭） + + 参数: + delta_v: 速度增量 (m/s) + engine: 发动机参数 + alpha: 结构系数 + 返回: + m_fuel / m_payload + """ + R = mass_ratio(delta_v, engine) + denominator = 1 - alpha * (R - 1) + + if denominator <= 0: + return np.inf + + return (R - 1) / denominator + + +def fuel_to_payload_ratio_multistage( + delta_v: float, + engine: EngineParams, + alpha: float = 0.1, + num_stages: int = 3 +) -> float: + """ + 计算燃料与载荷的质量比（多级火箭） + + 假设各级均分ΔV，使用相同发动机和结构系数 + + 参数: + delta_v: 总速度增量 (m/s) + engine: 发动机参数 + alpha: 结构系数 + num_stages: 级数 + 返回: + m_fuel_total / m_payload + """ + # 每级分配的ΔV + delta_v_per_stage = delta_v / num_stages + + # 每级的质量比 + R_stage = mass_ratio(delta_v_per_stage, engine) + + # 每级的燃料比 (相对于该级有效载荷) + denominator = 1 - alpha * (R_stage - 1) + if denominator <= 0: + return np.inf + + k_stage = (R_stage - 1) / denominator + + # 结构质量比 (相对于燃料) + # 每级末质量 = 载荷 + 结构 = 载荷 + alpha * 燃料 + # 每级初质量 = 载荷 + 结构 + 燃料 = 载荷 * (1 + k_stage * (1 + alpha)) + + # 对于多级火箭，总质量比是各级质量比的乘积 + # 第i级的"载荷"是第i+1级及以上所有质量 + + # 简化计算：假设各级结构系数相同 + # 级比 = (1 + k_stage * (1 + alpha)) + stage_ratio = 1 + k_stage * (1 + alpha) + + # 总初始质量/最终载荷 = stage_ratio ^ num_stages + total_ratio = stage_ratio ** num_stages + + # 总燃料 = 总初始质量 - 载荷 - 总结构 + # 近似：总燃料/载荷 ≈ total_ratio - 1 (忽略结构细节) + + # 更精确的计算 + 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 ground_launch_specific_energy( + engine: EngineParams, + mission: MissionParams, + constants: PhysicalConstants = PhysicalConstants(), + num_stages: int = 3 +) -> dict: + """ + 计算地面发射的比能量（使用多级火箭） + + 参数: + engine: 发动机参数 + mission: 任务参数 + constants: 物理常数 + num_stages: 火箭级数 (默认3级) + + 返回: + 包含各项比能量的字典 (单位: J/kg 载荷) + """ + # 总ΔV + delta_v_total = mission.total_delta_v_ground + + # 燃料/载荷比 (使用多级火箭模型) + k = fuel_to_payload_ratio_multistage(delta_v_total, engine, mission.alpha, num_stages) + + # 推进剂化学能 (J/kg 载荷) + propellant_energy = k * engine.specific_energy + + # 载荷轨道能量增益 + r_moon_orbit = constants.R_moon + mission.lunar_orbit_altitude + + # 从地面静止到月球轨道的比能量变化 + E_initial = -constants.GM_earth / constants.R_earth + E_final = -constants.GM_moon / (2 * r_moon_orbit) - constants.GM_earth / constants.d_earth_moon + payload_orbital_energy = E_final - E_initial + + # 效率 + efficiency = payload_orbital_energy / propellant_energy if propellant_energy > 0 and not np.isinf(propellant_energy) else 0 + + return { + 'delta_v': delta_v_total, + 'fuel_ratio': k, + 'num_stages': num_stages, + 'propellant_energy': propellant_energy, + 'payload_orbital_energy': payload_orbital_energy, + 'total_energy': propellant_energy, # 地面发射全靠推进剂 + 'efficiency': efficiency, + 'method': f'地面发射({num_stages}级)' + } + + +def ground_launch_specific_energy_with_latitude( + engine: EngineParams, + mission: MissionParams, + launch_site: LaunchSite, + constants: PhysicalConstants = PhysicalConstants(), + num_stages: int = 3 +) -> dict: + """ + 计算考虑发射场纬度的地面发射比能量 + + 参数: + engine: 发动机参数 + mission: 任务参数 + launch_site: 发射场 + constants: 物理常数 + num_stages: 火箭级数 (默认3级) + + 返回: + 包含各项比能量的字典 (单位: J/kg 载荷) + """ + # 计算纬度相关的ΔV + dv_info = total_delta_v_with_latitude(launch_site, mission, constants) + delta_v_total = dv_info['dv_total'] + + # 燃料/载荷比 (使用多级火箭模型) + k = fuel_to_payload_ratio_multistage(delta_v_total, engine, mission.alpha, num_stages) + + # 推进剂化学能 (J/kg 载荷) + propellant_energy = k * engine.specific_energy + + # 载荷轨道能量增益 + r_moon_orbit = constants.R_moon + mission.lunar_orbit_altitude + E_initial = -constants.GM_earth / constants.R_earth + E_final = -constants.GM_moon / (2 * r_moon_orbit) - constants.GM_earth / constants.d_earth_moon + payload_orbital_energy = E_final - E_initial + + # 效率 + efficiency = payload_orbital_energy / propellant_energy if propellant_energy > 0 and not np.isinf(propellant_energy) else 0 + + return { + 'launch_site': launch_site.name, + 'latitude': launch_site.latitude, + 'rotation_velocity': dv_info['rotation_velocity'], + 'dv_rotation_loss': dv_info['dv_rotation_loss'], + 'dv_inclination_change': dv_info['dv_inclination_change'], + 'delta_v': delta_v_total, + 'fuel_ratio': k, + 'num_stages': num_stages, + 'propellant_energy': propellant_energy, + 'payload_orbital_energy': payload_orbital_energy, + 'total_energy': propellant_energy, + 'efficiency': efficiency, + 'method': f'地面发射({launch_site.name}, {num_stages}级)' + } + + +def compare_launch_sites( + engine: EngineParams = EngineParams(), + mission: MissionParams = MissionParams(), + constants: PhysicalConstants = PhysicalConstants(), + num_stages: int = 3 +) -> list: + """ + 比较不同发射场的燃料和能量需求 + + 返回: + 各发射场的分析结果列表 (按纬度绝对值排序) + """ + results = [] + for name, site in LAUNCH_SITES.items(): + # 使用纬度绝对值进行计算（南北半球对称） + abs_lat = abs(site.latitude) + site_normalized = LaunchSite(site.name, abs_lat) + result = ground_launch_specific_energy_with_latitude( + engine, mission, site_normalized, constants, num_stages + ) + result['original_latitude'] = site.latitude + result['abs_latitude'] = abs_lat + results.append(result) + + return sorted(results, key=lambda x: x['abs_latitude']) + + +def print_latitude_comparison(results: list): + """打印不同纬度发射场的对比""" + print("=" * 110) + print("Launch Site Latitude Effects on Fuel and Energy") + print("=" * 110) + + print(f"\n{'Launch Site':<30} {'Lat':>6} {'V_rot':>8} {'ΔV_loss':>8} {'ΔV_total':>9} {'Fuel/PL':>9} {'Energy':>10}") + print(f"{'':30} {'(°)':>6} {'(m/s)':>8} {'(m/s)':>8} {'(km/s)':>9} {'Ratio':>9} {'(MJ/kg)':>10}") + print("-" * 110) + + for r in results: + lat = r.get('abs_latitude', abs(r['latitude'])) + print(f"{r['launch_site']:<30} {lat:>6.1f} {r['rotation_velocity']:>8.0f} {r['dv_rotation_loss']:>8.0f} {r['delta_v']/1000:>9.2f} {r['fuel_ratio']:>9.1f} {r['total_energy']/1e6:>10.1f}") + + print("-" * 110) + + # 计算相对赤道的损失 + equator_result = next((r for r in results if abs(r['latitude']) < 0.1), results[0]) + print(f"\nExtra consumption relative to Equator launch:") + for r in results: + lat = r.get('abs_latitude', abs(r['latitude'])) + if lat > 0.1: + fuel_increase = (r['fuel_ratio'] / equator_result['fuel_ratio'] - 1) * 100 + energy_increase = (r['total_energy'] / equator_result['total_energy'] - 1) * 100 + print(f" {r['launch_site']:<30}: Fuel +{fuel_increase:>5.1f}%, Energy +{energy_increase:>5.1f}%") + + print("=" * 110) + + +def plot_latitude_effects( + engine: EngineParams = EngineParams(), + mission: MissionParams = MissionParams(), + constants: PhysicalConstants = PhysicalConstants(), + save_path: str = '/Volumes/Files/code/mm/20260130_b/latitude_effects.png', + lunar_mission: bool = True +): + """ + 绘制纬度影响图表 + + 参数: + lunar_mission: 如果True，使用月球任务场景（不需要轨道平面改变） + """ + + # 月球任务场景：设置高目标倾角，避免轨道平面改变惩罚 + if lunar_mission: + mission_plot = MissionParams( + alpha=mission.alpha, + lunar_orbit_altitude=mission.lunar_orbit_altitude, + delta_v_ground_to_leo_base=mission.delta_v_ground_to_leo_base, + delta_v_tli=mission.delta_v_tli, + delta_v_loi=mission.delta_v_loi, + target_inclination=90.0, # 允许任意倾角 + elevator_length=mission.elevator_length + ) + title_suffix = "\n(Lunar Mission - 无轨道平面改变)" + else: + mission_plot = mission + title_suffix = "\n(目标: 赤道轨道)" + + fig, axes = plt.subplots(2, 2, figsize=(16, 14)) + + # 连续纬度范围 + latitudes = np.linspace(0, 65, 100) + + # ========== 图1: 自转速度 vs 纬度 ========== + ax1 = axes[0, 0] + rotation_velocities = [earth_rotation_velocity(lat, constants) for lat in latitudes] + ax1.plot(latitudes, rotation_velocities, 'b-', linewidth=2) + + # 标记发射场 (使用绝对值纬度) + for name, site in LAUNCH_SITES.items(): + abs_lat = abs(site.latitude) + v = earth_rotation_velocity(abs_lat, constants) + ax1.plot(abs_lat, v, 'ro', markersize=8) + # 简化标签显示 + label = site.name.split('(')[0].strip() + ax1.annotate(label, (abs_lat, v), textcoords="offset points", + xytext=(3, 3), fontsize=8, rotation=15) + + ax1.set_xlabel('Latitude |φ| (°)', fontsize=12) + ax1.set_ylabel('Rotation Velocity (m/s)', fontsize=12) + ax1.set_title('Earth Surface Rotation Velocity vs Latitude\nv = ω×R×cos(φ)', fontsize=13) + ax1.grid(True, alpha=0.3) + ax1.set_xlim(0, 65) + + # ========== 图2: ΔV损失 vs 纬度 ========== + ax2 = axes[0, 1] + dv_losses = [delta_v_rotation_loss(lat, constants) for lat in latitudes] + ax2.plot(latitudes, dv_losses, 'r-', linewidth=2) + + for name, site in LAUNCH_SITES.items(): + abs_lat = abs(site.latitude) + dv = delta_v_rotation_loss(abs_lat, constants) + ax2.plot(abs_lat, dv, 'bo', markersize=8) + label = site.name.split('(')[0].strip() + ax2.annotate(label, (abs_lat, dv), textcoords="offset points", + xytext=(3, 3), fontsize=8, rotation=15) + + ax2.set_xlabel('Latitude |φ| (°)', fontsize=12) + ax2.set_ylabel('ΔV Loss (m/s)', fontsize=12) + ax2.set_title('Rotation Velocity Loss vs Latitude\n(Relative to Equator)', fontsize=13) + ax2.grid(True, alpha=0.3) + ax2.set_xlim(0, 65) + + # ========== 图3: 燃料比 vs 纬度 ========== + ax3 = axes[1, 0] + + fuel_ratios = [] + for lat in latitudes: + site = LaunchSite("temp", lat) + result = ground_launch_specific_energy_with_latitude(engine, mission_plot, site, constants) + fuel_ratios.append(result['fuel_ratio']) + + ax3.plot(latitudes, fuel_ratios, 'g-', linewidth=2, label='Fuel Ratio vs Latitude') + + for name, site in LAUNCH_SITES.items(): + abs_lat = abs(site.latitude) + site_temp = LaunchSite(site.name, abs_lat) + result = ground_launch_specific_energy_with_latitude(engine, mission_plot, site_temp, constants) + ax3.plot(abs_lat, result['fuel_ratio'], 'mo', markersize=8) + label = site.name.split('(')[0].strip() + ax3.annotate(label, (abs_lat, result['fuel_ratio']), + textcoords="offset points", xytext=(3, 3), fontsize=8, rotation=15) + + ax3.set_xlabel('Latitude |φ| (°)', fontsize=12) + ax3.set_ylabel('Fuel / Payload Mass Ratio', fontsize=12) + ax3.set_title(f'Fuel Requirement vs Launch Latitude{title_suffix}', fontsize=13) + ax3.grid(True, alpha=0.3) + ax3.set_xlim(0, 65) + + # ========== 图4: 能量对比柱状图 ========== + ax4 = axes[1, 1] + + # 按纬度绝对值排序 + results = [] + for name, site in LAUNCH_SITES.items(): + abs_lat = abs(site.latitude) + site_temp = LaunchSite(site.name, abs_lat) + result = ground_launch_specific_energy_with_latitude(engine, mission_plot, site_temp, constants) + result['abs_latitude'] = abs_lat + results.append(result) + + results = sorted(results, key=lambda x: x['abs_latitude']) + + # 简化标签 + sites = [r['launch_site'].split('(')[0].strip() for r in results] + fuel_ratios_bar = [r['fuel_ratio'] for r in results] + abs_lats = [r['abs_latitude'] for r in results] + + # 颜色映射：纬度越高颜色越深 + colors = plt.cm.RdYlGn_r(np.array(abs_lats) / 65) + bars = ax4.bar(range(len(sites)), fuel_ratios_bar, color=colors) + + ax4.set_ylabel('Fuel / Payload Mass Ratio', fontsize=12) + ax4.set_title(f'Fuel Requirement by Launch Site{title_suffix}', fontsize=13) + ax4.set_xticks(range(len(sites))) + ax4.set_xticklabels(sites, rotation=45, ha='right', fontsize=9) + ax4.grid(True, alpha=0.3, axis='y') + + # 添加数值标签和纬度 + for i, (bar, ratio, lat) in enumerate(zip(bars, fuel_ratios_bar, abs_lats)): + ax4.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.3, + f'{ratio:.1f}', ha='center', fontsize=9, fontweight='bold') + ax4.text(bar.get_x() + bar.get_width()/2, 1, + f'{lat:.0f}°', ha='center', fontsize=8, color='white') + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"纬度影响图表已保存至: {save_path}") + + return fig + + +# ============== 电梯发射比能量计算 ============== + +def elevator_launch_specific_energy( + engine: EngineParams, + mission: MissionParams, + constants: PhysicalConstants = PhysicalConstants(), + release_height: Optional[float] = None +) -> dict: + """ + 计算太空电梯发射的比能量 + + 参数: + engine: 发动机参数 + mission: 任务参数 + constants: 物理常数 + release_height: 释放高度 (m)，None则使用电梯顶端 + + 返回: + 包含各项比能量的字典 (单位: J/kg 载荷) + """ + # 释放点参数 + if release_height is None: + release_height = mission.elevator_length + + r_release = constants.R_earth + release_height + + # 释放点速度 (随地球自转) + v_release = constants.omega_earth * r_release + + # 释放点的逃逸速度 + v_escape = np.sqrt(2 * constants.GM_earth / r_release) + + # 释放点的圆轨道速度 + v_circular = np.sqrt(constants.GM_earth / r_release) + + # 计算C3 (特征能量) + C3 = v_release**2 - v_escape**2 + + # 电梯提升能量 (从地面到释放点) + # 包括势能变化和动能获得 + lift_energy = ( + constants.GM_earth / constants.R_earth + - constants.GM_earth / r_release + + 0.5 * v_release**2 + ) + + # 计算所需ΔV + if C3 > 0: + # 已超过逃逸速度，需要减速 + # 目标：进入地月转移轨道 (C3 ≈ -2 km²/s²) + C3_target = -2e6 # -2 km²/s² 转换为 m²/s² + v_needed = np.sqrt(max(0, v_escape**2 + C3_target)) + delta_v_adjust = abs(v_release - v_needed) + maneuver_type = "减速" + else: + # 未达逃逸速度，需要加速 + # 计算到达月球轨道所需的额外速度 + v_needed = np.sqrt(v_escape**2 - 2e6) # 地月转移 + delta_v_adjust = v_needed - v_release + maneuver_type = "加速" + + # 加上月球轨道捕获 + delta_v_total = abs(delta_v_adjust) + mission.delta_v_loi + + # 燃料/载荷比 + k = fuel_to_payload_ratio(delta_v_total, engine, mission.alpha) + + # 推进剂化学能 + propellant_energy = k * engine.specific_energy + + # 载荷轨道能量增益 (与地面发射相同) + r_moon_orbit = constants.R_moon + mission.lunar_orbit_altitude + E_initial = -constants.GM_earth / constants.R_earth + E_final = -constants.GM_moon / (2 * r_moon_orbit) - constants.GM_earth / constants.d_earth_moon + payload_orbital_energy = E_final - E_initial + + # 总能量 = 电梯提升 + 推进剂 + total_energy = lift_energy + propellant_energy + + # 效率 + efficiency = payload_orbital_energy / total_energy if total_energy > 0 else 0 + + return { + 'release_height': release_height, + 'release_velocity': v_release, + 'escape_velocity': v_escape, + 'C3': C3, + 'maneuver_type': maneuver_type, + 'delta_v': delta_v_total, + 'fuel_ratio': k, + 'lift_energy': lift_energy, + 'propellant_energy': propellant_energy, + 'payload_orbital_energy': payload_orbital_energy, + 'total_energy': total_energy, + 'efficiency': efficiency, + 'method': '电梯发射' + } + + +# ============== 对比分析函数 ============== + +def compare_launch_methods( + engine: EngineParams = EngineParams(), + mission: MissionParams = MissionParams(), + constants: PhysicalConstants = PhysicalConstants() +) -> dict: + """ + 对比两种发射方式 + + 返回: + 对比结果字典 + """ + ground = ground_launch_specific_energy(engine, mission, constants) + elevator = elevator_launch_specific_energy(engine, mission, constants) + + return { + 'ground': ground, + 'elevator': elevator, + 'fuel_saving': 1 - elevator['fuel_ratio'] / ground['fuel_ratio'], + 'energy_saving': 1 - elevator['total_energy'] / ground['total_energy'], + 'delta_v_saving': 1 - elevator['delta_v'] / ground['delta_v'] + } + + +def print_comparison(comparison: dict): + """打印对比结果""" + ground = comparison['ground'] + elevator = comparison['elevator'] + + ground_label = ground['method'] + elevator_label = elevator['method'] + + print("=" * 70) + print("地面发射 vs 太空电梯发射：比能量对比") + print("=" * 70) + + print(f"\n{'参数':<25} {ground_label:>18} {elevator_label:>18}") + print("-" * 70) + + print(f"{'ΔV (km/s)':<25} {ground['delta_v']/1000:>18.2f} {elevator['delta_v']/1000:>18.2f}") + print(f"{'燃料/载荷比':<25} {ground['fuel_ratio']:>18.2f} {elevator['fuel_ratio']:>18.2f}") + print(f"{'电梯提升能量 (MJ/kg)':<25} {'0':>18} {elevator['lift_energy']/1e6:>18.1f}") + print(f"{'推进剂能量 (MJ/kg)':<25} {ground['propellant_energy']/1e6:>18.1f} {elevator['propellant_energy']/1e6:>18.1f}") + print(f"{'总能量 (MJ/kg)':<25} {ground['total_energy']/1e6:>18.1f} {elevator['total_energy']/1e6:>18.1f}") + print(f"{'载荷能量增益 (MJ/kg)':<25} {ground['payload_orbital_energy']/1e6:>18.1f} {elevator['payload_orbital_energy']/1e6:>18.1f}") + print(f"{'系统效率':<25} {ground['efficiency']*100:>17.1f}% {elevator['efficiency']*100:>17.1f}%") + + print("\n" + "-" * 70) + print(f"{'ΔV 节省':<25} {comparison['delta_v_saving']*100:>18.1f}%") + print(f"{'燃料节省':<25} {comparison['fuel_saving']*100:>18.1f}%") + print(f"{'总能量节省':<25} {comparison['energy_saving']*100:>18.1f}%") + print("=" * 70) + + +def plot_comparison( + engine: EngineParams = EngineParams(), + mission: MissionParams = MissionParams(), + constants: PhysicalConstants = PhysicalConstants(), + save_path: str = '/Volumes/Files/code/mm/20260130_b/specific_energy_comparison.png' +): + """绘制对比图表""" + + fig, axes = plt.subplots(2, 2, figsize=(14, 12)) + + # ========== 图1: 不同电梯高度的ΔV需求 ========== + ax1 = axes[0, 0] + heights = np.linspace(3.6e7, 1.5e8, 100) # 36,000 km 到 150,000 km + + delta_vs = [] + for h in heights: + mission_temp = MissionParams(elevator_length=h) + result = elevator_launch_specific_energy(engine, mission_temp, constants, h) + delta_vs.append(result['delta_v'] / 1000) + + ax1.plot(heights/1e6, delta_vs, 'b-', linewidth=2, label='电梯发射') + ax1.axhline(y=mission.total_delta_v_ground/1000, color='r', linestyle='--', + linewidth=2, label='地面发射') + ax1.axvline(x=35.786, color='g', linestyle=':', label='GEO高度') + + ax1.set_xlabel('电梯释放高度 (千km)', fontsize=12) + ax1.set_ylabel('ΔV 需求 (km/s)', fontsize=12) + ax1.set_title('释放高度 vs ΔV需求', fontsize=14) + ax1.legend() + ax1.grid(True, alpha=0.3) + + # ========== 图2: 燃料/载荷比对比 ========== + ax2 = axes[0, 1] + + fuel_ratios_elevator = [] + for h in heights: + mission_temp = MissionParams(elevator_length=h) + result = elevator_launch_specific_energy(engine, mission_temp, constants, h) + fuel_ratios_elevator.append(result['fuel_ratio']) + + ground_result = ground_launch_specific_energy(engine, mission, constants, num_stages=3) + + ax2.plot(heights/1e6, fuel_ratios_elevator, 'b-', linewidth=2, label='电梯发射') + ax2.axhline(y=ground_result['fuel_ratio'], color='r', linestyle='--', + linewidth=2, label='地面发射(3级)') + + ax2.set_xlabel('电梯释放高度 (千km)', fontsize=12) + ax2.set_ylabel('燃料/载荷 质量比', fontsize=12) + ax2.set_title('释放高度 vs 燃料效率', fontsize=14) + ax2.legend() + ax2.grid(True, alpha=0.3) + ax2.set_ylim(0, 30) + + # ========== 图3: 能量构成对比 (柱状图) ========== + ax3 = axes[1, 0] + + comparison = compare_launch_methods(engine, mission, constants) + + categories = ['地面发射', '电梯发射'] + lift_energy = [0, comparison['elevator']['lift_energy']/1e6] + prop_energy = [comparison['ground']['propellant_energy']/1e6, + comparison['elevator']['propellant_energy']/1e6] + + x = np.arange(len(categories)) + width = 0.35 + + bars1 = ax3.bar(x, lift_energy, width, label='电梯提升能量', color='green') + bars2 = ax3.bar(x, prop_energy, width, bottom=lift_energy, label='推进剂能量', color='orange') + + ax3.set_ylabel('比能量 (MJ/kg)', fontsize=12) + ax3.set_title('能量构成对比', fontsize=14) + ax3.set_xticks(x) + ax3.set_xticklabels(categories) + ax3.legend() + ax3.grid(True, alpha=0.3, axis='y') + + # 添加数值标签 + for i, (l, p) in enumerate(zip(lift_energy, prop_energy)): + total = l + p + ax3.text(i, total + 5, f'{total:.0f}', ha='center', fontsize=11, fontweight='bold') + + # ========== 图4: 不同发动机对比 ========== + ax4 = axes[1, 1] + + engines = [ + EngineParams(name="固体火箭", isp=280, specific_energy=5e6), + EngineParams(name="液氧煤油", isp=350, specific_energy=10.3e6), + EngineParams(name="液氧液氢", isp=450, specific_energy=15.5e6), + ] + + x_pos = np.arange(len(engines)) + ground_ratios = [] + elevator_ratios = [] + + for eng in engines: + g = ground_launch_specific_energy(eng, mission, constants, num_stages=3) + e = elevator_launch_specific_energy(eng, mission, constants) + ground_ratios.append(min(g['fuel_ratio'], 100)) # 限制显示范围 + elevator_ratios.append(e['fuel_ratio']) + + width = 0.35 + ax4.bar(x_pos - width/2, ground_ratios, width, label='地面发射(3级)', color='red', alpha=0.7) + ax4.bar(x_pos + width/2, elevator_ratios, width, label='电梯发射', color='blue', alpha=0.7) + + ax4.set_ylabel('燃料/载荷 质量比', fontsize=12) + ax4.set_title('不同发动机的燃料效率对比', fontsize=14) + ax4.set_xticks(x_pos) + ax4.set_xticklabels([e.name for e in engines]) + ax4.legend() + ax4.grid(True, alpha=0.3, axis='y') + + plt.tight_layout() + plt.savefig(save_path, dpi=150, bbox_inches='tight') + print(f"图表已保存至: {save_path}") + + return fig + + +# ============== 主程序 ============== + +if __name__ == "__main__": + # 可以在这里修改参数 + + # 发动机参数 + engine = EngineParams( + name="液氧液氢", + isp=450, # 比冲 (秒) + specific_energy=15.5e6 # 燃料比能量 (J/kg) + ) + + # 任务参数 + mission = MissionParams( + alpha=0.10, # 结构系数 + lunar_orbit_altitude=100e3, # 环月轨道高度 (m) + delta_v_ground_to_leo_base=9400, # 赤道发射到LEO的基准ΔV (m/s) + delta_v_tli=3100, # 地月转移注入 (m/s) + delta_v_loi=800, # 月球轨道捕获 (m/s) + target_inclination=0.0, # 目标轨道倾角 (度) + elevator_length=1.0e8 # 电梯长度 10万km (m) + ) + + # 物理常数 (通常不需要修改) + constants = PhysicalConstants() + + # ========== 1. 地面发射 vs 电梯发射对比 ========== + comparison = compare_launch_methods(engine, mission, constants) + print_comparison(comparison) + + # 绘制对比图表 + print("\n正在生成对比图表...") + plot_comparison(engine, mission, constants) + + # 电梯发射详细信息 + print("\n" + "=" * 70) + print("电梯发射详细参数:") + print("=" * 70) + elevator = comparison['elevator'] + print(f" 释放高度: {elevator['release_height']/1e6:.1f} 千km") + print(f" 释放点速度: {elevator['release_velocity']/1000:.2f} km/s") + print(f" 当地逃逸速度: {elevator['escape_velocity']/1000:.2f} km/s") + print(f" 特征能量 C3: {elevator['C3']/1e6:.1f} km²/s²") + print(f" 机动类型: {elevator['maneuver_type']}") + + # ========== 2. 发射场纬度影响分析 (月球任务) ========== + print("\n") + print("=" * 110) + print("LUNAR MISSION: Launch Site Comparison") + print("(No orbital plane change required - only rotation velocity loss)") + print("=" * 110) + + # 月球任务配置（允许倾角与发射场纬度匹配） + mission_lunar = MissionParams( + alpha=mission.alpha, + lunar_orbit_altitude=mission.lunar_orbit_altitude, + delta_v_ground_to_leo_base=mission.delta_v_ground_to_leo_base, + delta_v_tli=mission.delta_v_tli, + delta_v_loi=mission.delta_v_loi, + target_inclination=90.0, # 允许任意倾角 + elevator_length=mission.elevator_length + ) + + latitude_results = compare_launch_sites(engine, mission_lunar, constants) + print_latitude_comparison(latitude_results) + + # 绘制纬度影响图表 (月球任务模式) + print("\n正在生成纬度影响图表 (Lunar Mission)...") + plot_latitude_effects(engine, mission_lunar, constants, lunar_mission=True) + + # ========== 3. 公式总结 ========== + print("\n" + "=" * 70) + print("数学关系总结:") + print("=" * 70) + print(""" + 地球自转速度: + v(φ) = ω × R_E × cos(φ) + + 自转速度损失: + Δv_loss = v(0°) - v(φ) = 465 × (1 - cos(φ)) m/s + + 纬度对燃料的影响: + 由于 ΔV 增加 → R = exp(ΔV/ve) 指数增长 → 燃料比指数增长 + + 关键结论: + - 赤道发射: 自转贡献 465 m/s, 燃料最省 + - 高纬度发射: 损失自转贡献, 且可能需要额外倾角机动 + - 每损失 100 m/s 自转速度 ≈ 增加 ~3% 燃料需求 +""")