p1: fig

p1/combined_scenario.py

@@ -515,7 +515,7 @@ def print_summary(df: pd.DataFrame):
 
 def plot_three_scenarios_comparison(
     year_min: float = 100,
-    year_max: float = 300,
+    year_max: float = 250,
     save_path: str = '/Volumes/Files/code/mm/20260130_b/three_scenarios_comparison.png'
 ):
     """
@@ -579,74 +579,45 @@ def plot_three_scenarios_comparison(
         else:
             combined_energy.append(np.nan)
     
-    # 创建图表
+    # 创建图表 - 使用0.618黄金比例，缩小尺寸以放大字体
-    fig, ax = plt.subplots(figsize=(14, 8))
+    fig_width = 8
+    fig_height = fig_width * 0.618
+    fig, ax = plt.subplots(figsize=(fig_width, fig_height))
+    
+    # 中偏低饱和度配色
+    color_rocket = '#B85450'    # 暗砖红色
+    color_elevator = '#5B8266'  # 灰绿色
+    color_combined = '#4A6FA5'  # 灰蓝色
     
     # 绘制三条折线
-    ax.plot(years, rocket_energy, 'r-', linewidth=2.5, label='Rocket Only', marker='', markersize=4)
+    ax.plot(years, rocket_energy, color=color_rocket, 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, elevator_energy, color=color_elevator, 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.plot(years, combined_energy, color=color_combined, 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=combined_min_years, color=color_combined, 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=elevator_only_min_years, color=color_elevator, linestyle=':', alpha=0.7, linewidth=1.5)
-    ax.axvline(x=rocket_only_min_years, color='red', linestyle=':', alpha=0.7, linewidth=1.5)
+    ax.axvline(x=rocket_only_min_years, color=color_rocket, linestyle=':', alpha=0.7, linewidth=1.5)
     
-    # 添加标注
+    # 添加标注 - 使用低饱和度背景色，避免遮挡曲线
-    y_pos = ax.get_ylim()[1] * 0.95
+    y_max = ax.get_ylim()[1]
     ax.annotate(f'Combined\nMin: {combined_min_years:.0f}y', 
-                xy=(combined_min_years, y_pos), fontsize=9, ha='center',
+                xy=(combined_min_years, y_max * 0.55), fontsize=9, ha='center',
-                bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
+                bbox=dict(boxstyle='round', facecolor='#C5D5E8', edgecolor=color_combined, alpha=0.9))
     ax.annotate(f'Elevator\nMin: {elevator_only_min_years:.0f}y', 
-                xy=(elevator_only_min_years, y_pos * 0.85), fontsize=9, ha='center',
+                xy=(elevator_only_min_years, y_max * 0.72), fontsize=9, ha='center',
-                bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
+                bbox=dict(boxstyle='round', facecolor='#C8DBC8', edgecolor=color_elevator, alpha=0.9))
     ax.annotate(f'Rocket\nMin: {rocket_only_min_years:.0f}y', 
-                xy=(rocket_only_min_years, y_pos * 0.75), fontsize=9, ha='center',
+                xy=(rocket_only_min_years, y_max * 0.55), fontsize=9, ha='center',
-                bbox=dict(boxstyle='round', facecolor='lightsalmon', alpha=0.8))
+                bbox=dict(boxstyle='round', facecolor='#E8C8C5', edgecolor=color_rocket, alpha=0.9))
     
-    # 填充区域表示节能效果
+    ax.set_xlabel('Completion Time (years)', fontsize=11)
-    # Combined比Rocket节省的能量
+    ax.set_ylabel('Total Energy Consumption (PJ)', fontsize=11)
-    rocket_arr = np.array(rocket_energy)
+    ax.legend(loc='upper left', fontsize=9)
-    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)
     
-    # 添加数据表格
+    plt.tight_layout()
-    # 选择几个关键年份的数据
-    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}")
@@ -700,6 +671,6 @@ if __name__ == "__main__":
         print(f"{name:<25} {data['years']:>12.0f} {data['energy_PJ']:>15.0f} {savings:>17.0f}%")
     print("=" * 90)
     
-    # 生成三方案对比折线图 (100-300年)
+    # 生成三方案对比折线图 (100-250年)
-    print("\n\n正在生成三方案对比折线图 (100-300年)...")
+    print("\n\n正在生成三方案对比折线图 (100-250年)...")
-    plot_three_scenarios_comparison(year_min=100, year_max=300)
+    plot_three_scenarios_comparison(year_min=100, year_max=250)

p1/pareto_optimization.py

@@ -796,6 +796,60 @@ def plot_combined_range_analysis(
     return fig
 
 
+def plot_marginal_benefit_only(
+    df: pd.DataFrame,
+    knee_analysis: Dict,
+    save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/marginal_benefit.png'
+):
+    """
+    仅绘制边际收益图（右下角图的简化版）
+    - 长宽比 0.618
+    - 中低饱和度色系
+    - 无标题
+    - 无 Elevator share 线
+    - Y轴范围 450-600
+    """
+    # 使用0.618黄金比例，缩小尺寸以放大字体
+    fig_width = 8
+    fig_height = fig_width * 0.618
+    fig, ax = plt.subplots(figsize=(fig_width, fig_height))
+    
+    T = df['years'].values
+    E = df['energy_PJ'].values
+    
+    # 边际节省
+    dE_dT = np.gradient(E, T)
+    marginal_savings = -dE_dT  # 每多1年节省的能量
+    
+    # 中低饱和度色系
+    color_main = '#4A6FA5'  # 灰蓝色
+    color_fill = '#A8C0D8'  # 浅灰蓝
+    color_vline = '#B85450'  # 暗砖红
+    
+    # 边际收益曲线
+    ax.plot(T, marginal_savings, color=color_main, linewidth=2.5, label='Marginal Saving (PJ/year)')
+    ax.fill_between(T, marginal_savings, alpha=0.25, color=color_fill)
+    
+    # 标记推荐点
+    rec = knee_analysis['recommended']
+    ax.axvline(x=rec['years'], color=color_vline, linestyle=':', linewidth=2,
+               label=f'Recommended: {rec["years"]:.0f}y')
+    
+    ax.set_xlabel('Completion Time (years)', fontsize=11)
+    ax.set_ylabel('Marginal Energy Saving (PJ/year)', fontsize=11)
+    ax.legend(loc='upper right', fontsize=9)
+    ax.grid(True, alpha=0.3)
+    
+    # 设置Y轴范围 450-600
+    ax.set_ylim(450, 600)
+    
+    plt.tight_layout()
+    plt.savefig(save_path, dpi=150, bbox_inches='tight')
+    print(f"边际收益图已保存至: {save_path}")
+    
+    return fig
+
+
 def plot_combined_decision(
     df: pd.DataFrame,
     knee_analysis: Dict,

p1/richards_curve_2010.py

@@ -103,43 +103,49 @@ def main():
     print(f"K / Physical limit = {K_fixed/physical_max:.2f}x")
     
     # ========== Create Visualization ==========
-    fig, ax = plt.subplots(figsize=(12, 7))
+    # 缩小图片尺寸（比例不变），使字体相对更大
+    fig, ax = plt.subplots(figsize=(8, 4.67))
+    
+    # 中低饱和度配色
+    color_data = '#5D6D7E'       # 灰蓝色 - 数据点
+    color_curve = '#52796F'      # 暗绿色 - S曲线
+    color_target = '#7B9EA8'     # 灰蓝绿色 - 2050标记
     
     # Historical data points
-    ax.scatter(years, launches, color='#2C3E50', s=100, alpha=0.9, 
+    ax.scatter(years, launches, color=color_data, s=80, alpha=0.85, 
                label='Historical Data (2010-2025)', zorder=4, edgecolor='white', linewidth=1)
     
     # Generate smooth S-curve
-    years_smooth = np.linspace(start_year, 2080, 500)
+    years_smooth = np.linspace(start_year, 2055, 500)
     t_smooth = years_smooth - start_year
     pred_smooth = richards(t_smooth, K_fixed, r, t0, v)
     
     # S-curve prediction
-    ax.plot(years_smooth, pred_smooth, color='#27AE60', lw=3, 
+    ax.plot(years_smooth, pred_smooth, color=color_curve, lw=2.5, 
             label=f'Richards Model (K={K_fixed}, R²={r_squared:.3f})', zorder=2)
     
     # K=4298 saturation line
-    ax.axhline(K_fixed, color='#27AE60', ls=':', lw=2, alpha=0.7,
+    ax.axhline(K_fixed, color=color_curve, ls=':', lw=1.5, alpha=0.6,
                label=f'K = {K_fixed}')
     
     # Mark 2050 line only
-    ax.axvline(2050, color='#3498DB', ls=':', lw=2, alpha=0.8)
+    ax.axvline(2050, color=color_target, ls=':', lw=1.5, alpha=0.7)
-    ax.text(2051, K_fixed*0.85, '2050\n(Target)', fontsize=10, color='#3498DB')
+    ax.text(2050.5, K_fixed*0.83, '2050\n(Target)', fontsize=9, color=color_target)
     
     # Only show 2050 prediction point
     t_2050 = 2050 - start_year
     pred_2050 = richards(t_2050, K_fixed, r, t0, v)
-    ax.scatter([2050], [pred_2050], color='#3498DB', s=80, marker='D', zorder=4)
+    ax.scatter([2050], [pred_2050], color=color_target, s=60, marker='D', zorder=4)
     ax.annotate(f'{pred_2050:.0f}', xy=(2050, pred_2050), 
-                xytext=(2051, pred_2050+150),
+                xytext=(2050.5, pred_2050+180),
-                fontsize=9, color='#3498DB', fontweight='bold')
+                fontsize=9, color=color_target, fontweight='bold')
     
     # Formatting
-    ax.set_xlabel('Year', fontsize=12)
+    ax.set_xlabel('Year', fontsize=11)
-    ax.set_ylabel('Annual Launches', fontsize=12)
+    ax.set_ylabel('Annual Launches', fontsize=11)
-    ax.legend(loc='upper left', fontsize=10)
+    ax.legend(loc='upper left', fontsize=9)
-    ax.grid(True, alpha=0.3)
+    ax.grid(True, alpha=0.25)
-    ax.set_xlim(2008, 2080)
+    ax.set_xlim(2010, 2055)
     ax.set_ylim(0, K_fixed * 1.15)
     
     plt.tight_layout()

p1/sensitivity_analysis.py

@@ -585,6 +585,96 @@ def plot_sensitivity_analysis(
     return fig
 
 
+def plot_tornado_chart_standalone(
+    payload_df: pd.DataFrame,
+    elevator_df: pd.DataFrame, 
+    frequency_df: pd.DataFrame,
+    alpha_df: pd.DataFrame,
+    save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/tornado_chart.png'
+):
+    """
+    单独绘制 Tornado Chart（敏感性汇总图）
+    - XY轴反转（垂直柱状图）
+    - 长宽比 0.618
+    - 无标题
+    - 有图例
+    - 低饱和度色系
+    """
+    # 使用0.618黄金比例
+    fig_width = 8
+    fig_height = fig_width * 0.618
+    fig, ax = plt.subplots(figsize=(fig_width, fig_height))
+    
+    # 计算敏感性指标
+    baseline_T = 139
+    
+    sensitivities = [
+        ('Payload\n(100→150t)', 
+         (payload_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100,
+         (payload_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100),
+        ('Elevator Cap.\n(80%→120%)', 
+         (elevator_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100,
+         (elevator_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100),
+        ('Launch Freq.\n(0.5→2/day)', 
+         (frequency_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100,
+         (frequency_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100),
+        ('Struct. Coef.\n(0.06→0.14)', 
+         (alpha_df['optimal_T_lambda504'].iloc[0] - baseline_T) / baseline_T * 100,
+         (alpha_df['optimal_T_lambda504'].iloc[-1] - baseline_T) / baseline_T * 100),
+    ]
+    
+    # 低饱和度配色
+    color_positive = '#7A9CBF'  # 灰蓝色（参数增加导致T*增加）
+    color_negative = '#C4816B'  # 灰橙红色（参数增加导致T*减少）
+    
+    x_pos = np.arange(len(sensitivities))
+    bar_width = 0.6
+    
+    # 绘制柱状图（XY反转，使用垂直柱状图）
+    for i, (name, low, high) in enumerate(sensitivities):
+        bar_color = color_positive if high > low else color_negative
+        bar_height = high - low
+        bar_bottom = min(low, high)
+        
+        ax.bar(i, bar_height, bottom=bar_bottom, width=bar_width, 
+               color=bar_color, alpha=0.8, edgecolor='#555555', linewidth=0.8)
+        
+        # 标记端点
+        ax.plot(i, low, 'o', color='#444444', markersize=7, zorder=5)
+        ax.plot(i, high, '^', color='#444444', markersize=7, zorder=5)
+        ax.plot([i, i], [low, high], '-', color='#444444', linewidth=1.5, zorder=4)
+    
+    # 基准线
+    ax.axhline(y=0, color='#333333', linestyle='-', linewidth=1.2)
+    
+    # 设置轴
+    ax.set_xticks(x_pos)
+    ax.set_xticklabels([s[0] for s in sensitivities], fontsize=10)
+    ax.set_ylabel('Change in Optimal Timeline T* (%)', fontsize=11)
+    ax.grid(True, alpha=0.3, axis='y')
+    
+    # 添加图例
+    from matplotlib.patches import Patch
+    from matplotlib.lines import Line2D
+    legend_elements = [
+        Patch(facecolor=color_negative, edgecolor='#555555', alpha=0.8, 
+              label='Parameter ↑ → T* ↓'),
+        Patch(facecolor=color_positive, edgecolor='#555555', alpha=0.8, 
+              label='Parameter ↑ → T* ↑'),
+        Line2D([0], [0], marker='o', color='#444444', linestyle='None', 
+               markersize=6, label='Low value'),
+        Line2D([0], [0], marker='^', color='#444444', linestyle='None', 
+               markersize=6, label='High value'),
+    ]
+    ax.legend(handles=legend_elements, loc='upper right', fontsize=9, framealpha=0.9)
+    
+    plt.tight_layout()
+    plt.savefig(save_path, dpi=150, bbox_inches='tight')
+    print(f"Tornado chart saved to: {save_path}")
+    
+    return fig
+
+
 def plot_tradeoff_comparison(
     save_path: str = '/Volumes/Files/code/mm/20260130_b/p1/sensitivity_tradeoff.png'
 ):

p1/specific_energy_comparison.py

@@ -79,29 +79,29 @@ class LaunchSite:
 # 预定义的发射场
 LAUNCH_SITES = {
     # 赤道参考
-    '赤道': LaunchSite("Equator (参考)", 0.0),
+    '赤道': LaunchSite("Equator", 0.0),
     
     # 法属圭亚那 (接近赤道)
-    'Kourou': LaunchSite("Kourou (French Guiana)", 5.2),
+    'Kourou': LaunchSite("French Guiana", 5.2),
     
     # 印度
-    'SDSC': LaunchSite("Satish Dhawan (India)", 13.7),
+    'SDSC': LaunchSite("Satish Dhawan", 13.7),
     
     # 美国
-    'Texas': LaunchSite("Boca Chica (Texas)", 26.0),
+    'Texas': LaunchSite("Texas", 26.0),
-    'Florida': LaunchSite("Cape Canaveral (Florida)", 28.5),
+    'Florida': LaunchSite("Florida", 28.5),
-    'California': LaunchSite("Vandenberg (California)", 34.7),
+    'California': LaunchSite("California", 34.7),
-    'Virginia': LaunchSite("Wallops (Virginia)", 37.8),
+    'Virginia': LaunchSite("Virginia", 37.8),
-    'Alaska': LaunchSite("Kodiak (Alaska)", 57.4),
+    'Alaska': LaunchSite("Alaska", 57.4),
     
     # 中国
-    'Taiyuan': LaunchSite("Taiyuan (China)", 38.8),
+    'Taiyuan': LaunchSite("Taiyuan", 38.8),
     
     # 新西兰
-    'Mahia': LaunchSite("Mahia (New Zealand)", -39.3),  # 南半球
+    'Mahia': LaunchSite("Mahia", -39.3),  # 南半球
     
     # 哈萨克斯坦
-    'Baikonur': LaunchSite("Baikonur (Kazakhstan)", 45.6),
+    'Baikonur': LaunchSite("Kazakhstan", 45.6),
 }
 
 
@@ -535,8 +535,12 @@ def plot_latitude_effects(
         mission_plot = mission
         title_suffix = " (Target: Equatorial Orbit)"
     
-    # Golden ratio aspect: width=10, height=10*0.618
+    # Golden ratio aspect: 缩小尺寸使字体相对更大
-    fig, ax = plt.subplots(figsize=(10, 10 * 0.618))
+    fig, ax = plt.subplots(figsize=(8, 8 * 0.618))
+    
+    # 中低饱和度配色
+    color_curve = '#52796F'      # 暗绿色 - 曲线
+    color_point = '#8E7B9A'      # 灰紫色 - 数据点
     
     # Continuous latitude range
     latitudes = np.linspace(0, 65, 100)
@@ -549,29 +553,32 @@ def plot_latitude_effects(
         fuel_ratios.append(result['fuel_ratio'])
     
     # Plot continuous curve
-    ax.plot(latitudes, fuel_ratios, 'g-', linewidth=2, label='Fuel Ratio vs Latitude')
+    ax.plot(latitudes, fuel_ratios, color=color_curve, linewidth=2, label='Fuel Ratio vs Latitude')
     
     # Mark launch sites with custom offsets to avoid overlap with curve
     # Offsets: (x, y) in points - positive y moves text up, positive x moves right
+    # 纬度参考: Equator(0), French Guiana(5.2), Satish Dhawan(13.7), Texas(26), 
+    #          Florida(28.5), California(34.7), Virginia(37.8), Taiyuan(38.8), 
+    #          Mahia(39.3), Kazakhstan(45.6), Alaska(57.4)
     label_offsets = {
         '赤道': (5, -15),          # Below curve
-        'Kourou': (5, 8),          # Above curve
+        'Kourou': (-75, 15),       # Left above (French Guiana)
-        'SDSC': (5, -15),          # Below curve
+        'SDSC': (8, -15),          # Below right (Satish Dhawan)
-        'Texas': (-50, 8),         # Left and above
+        'Texas': (-38, -15),       # Left below
-        'Florida': (5, 8),         # Above curve
+        'Florida': (8, -15),       # Right below
-        'California': (-60, -8),   # Left
+        'California': (8, -15),    # Right below
-        'Virginia': (5, 8),        # Above curve
+        'Virginia': (8, 12),       # Right above
-        'Taiyuan': (-45, -12),     # Left and below
+        'Taiyuan': (-52, 12),      # Left above (避开California)
-        'Mahia': (5, 8),           # Above curve
+        'Mahia': (8, -18),         # Right below
-        'Alaska': (5, -15),        # Below curve
+        'Alaska': (8, 8),          # Right above
-        'Baikonur': (5, 8),        # Above curve
+        'Baikonur': (8, 8),        # Right above (Kazakhstan)
     }
     
     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)
-        ax.plot(abs_lat, result['fuel_ratio'], 'mo', markersize=8)
+        ax.plot(abs_lat, result['fuel_ratio'], 'o', color=color_point, markersize=8)
         # Simplified label
         label = site.name.split('(')[0].strip()
         # Get custom offset for this site, default to (5, 8)
@@ -583,7 +590,7 @@ def plot_latitude_effects(
     ax.set_ylabel('Fuel / Payload Mass Ratio', fontsize=12)
     # ax.set_title(f'Fuel Requirement vs Launch Latitude{title_suffix}', fontsize=13)
     ax.grid(True, alpha=0.3)
-    ax.set_xlim(0, 65)
+    ax.set_xlim(-2, 65)
     
     plt.tight_layout()
     plt.savefig(save_path, dpi=150, bbox_inches='tight')