toolCall: write call_16a94fab22c34a139c5a5783
{
"content": "#!/usr/bin/env python3\n\"\"\"\nMega-Rocket 100k-ton Class Launch Vehicle — Simulation & Figure Generation\nGenerates all figures for the LaTeX technical proposal.\n\"\"\"\n\nimport numpy as np\nimport matplotlib\nmatplotlib.use('Agg')\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import FancyBboxPatch, Rectangle, FancyArrowPatch, Circle\nfrom matplotlib.gridspec import GridSpec\nimport matplotlib.patches as mpatches\nfrom scipy.integrate import odeint\nfrom scipy.optimize import minimize_scalar\nimport os, json\n\nOUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'figures')\nos.makedirs(OUT, exist_ok=True)\n\n# ─── Global Style ───\nplt.rcParams.update({\n 'font.family': 'serif',\n 'font.size': 11,\n 'axes.labelsize': 12,\n 'axes.titlesize': 13,\n 'legend.fontsize': 9,\n 'figure.dpi': 200,\n 'savefig.dpi': 200,\n 'savefig.bbox': 'tight',\n 'mathtext.fontset': 'cm',\n})\n\n# ═══════════════════════════════════════════════════════════════════\n# VEHICLE PARAMETERS\n# ═══════════════════════════════════════════════════════════════════\nM0 = 100_000_000 # kg total liftoff mass (100,000 t)\n\n# Stage 0 — Booster cluster (LOX/RP-1 + LOX/LCH4 parallel burn)\nS0_PROP = 72_000_000 # kg\nS0_DRY = 8_500_000 # kg (structural fraction ~10.5%)\nS0_ISP_SL = 282 # s\nS0_ISP_VAC = 311 # s\nS0_BURN = 160 # s\n\n# Stage 1 — Core (LOX/LCH4)\nS1_PROP = 11_500_000\nS1_DRY = 1_200_000\nS1_ISP_VAC = 363 # s\nS1_BURN = 300\n\n# Stage 2 — Upper (LOX/LH2 + NTP augment)\nS2_PROP = 4_200_000\nS2_DRY = 580_000\nS2_ISP_VAC = 465 # s (LH2/LOX + NTP-augmented)\nS2_BURN = 480\n\nPAYLOAD = M0 - (S0_PROP+S0_DRY) - (S1_PROP+S1_DRY) - (S2_PROP+S2_DRY)\nDIAM = 25.0 # m\nHEIGHT = 220.0 # m\n\ng0 = 9.80665\nG_ME = 3.986004418e14\nR_E = 6.371e6\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 1 — Vehicle Configuration Cross-Section\n# ═══════════════════════════════════════════════════════════════════\ndef fig_vehicle_config():\n fig, ax = plt.subplots(figsize=(8, 18))\n ax.set_xlim(-20, 20)\n ax.set_ylim(-10, 230)\n ax.set_aspect('equal')\n ax.axis('off')\n ax.set_title('Leviathan-100 总体构型剖面图', fontsize=16, fontweight='bold', pad=20)\n\n # Nose fairing\n nose_x = np.linspace(-6, 6, 100)\n nose_y = 220 - 30 * np.sqrt(1 - (nose_x/6)**2) * 0 + np.zeros_like(nose_x)\n from matplotlib.patches import Polygon\n nose_pts = np.array([[-6,190],[6,190],[3,215],[0,220],[-3,215],[-6,190]])\n nose = Polygon(nose_pts, closed=True, fc='#4FC3F7', ec='#01579B', lw=1.5, alpha=0.7)\n ax.add_patch(nose)\n ax.text(0, 200, '整流罩', ha='center', va='center', fontsize=9, color='#01579B')\n\n # Stage 2 — upper\n s2 = Rectangle((-8, 150), 16, 40, fc='#81C784', ec='#1B5E20', lw=1.5, alpha=0.7)\n ax.add_patch(s2)\n ax.text(0, 170, '二级\\nLOX/LH₂\\n+NTP辅助', ha='center', va='center', fontsize=8, color='#1B5E20')\n ax.text(9, 165, f'{(S2_PROP+S2_DRY)/1e6:.1f} kt', ha='left', fontsize=8, color='#1B5E20')\n\n # Interstage\n inter = Rectangle((-9, 143), 18, 7, fc='#E0E0E0', ec='#616161', lw=1, alpha=0.5)\n ax.add_patch(inter)\n ax.text(0, 146.5, '级间段', ha='center', va='center', fontsize=7, color='#616161')\n\n # Stage 1 — core\n s1 = Rectangle((-12, 80), 24, 63, fc='#FFB74D', ec='#E65100', lw=1.5, alpha=0.7)\n ax.add_patch(s1)\n ax.text(0, 111, '一级 (芯级)\\nLOX/LCH₄', ha='center', va='center', fontsize=9, color='#E65100')\n ax.text(13, 111, f'{(S1_PROP+S1_DRY)/1e6:.1f} kt', ha='left', fontsize=8, color='#E65100')\n\n # Stage 0 — boosters (side)\n # Left booster\n s0l = Rectangle((-20, 5), 7, 75, fc='#EF9A9A', ec='#B71C1C', lw=1.5, alpha=0.7)\n ax.add_patch(s0l)\n ax.text(-16.5, 42, '助推A\\nLOX/RP-1', ha='center', va='center', fontsize=7, color='#B71C1C', rotation=90)\n\n # Right booster\n s0r = Rectangle((13, 5), 7, 75, fc='#EF9A9A', ec='#B71C1C', lw=1.5, alpha=0.7)\n ax.add_patch(s0r)\n ax.text(16.5, 42, '助推B\\nLOX/RP-1', ha='center', va='center', fontsize=7, color='#B71C1C', rotation=90)\n\n # Center core booster\n s0c = Rectangle((-12, 5), 24, 75, fc='#FFCDD2', ec='#C62828', lw=1.5, alpha=0.5)\n ax.add_patch(s0c)\n ax.text(0, 42, '助推芯级\\nLOX/RP-1\\n+LOX/LCH₄\\n并联', ha='center', va='center', fontsize=8, color='#C62828')\n ax.text(13, 42, f'{(S0_PROP+S0_DRY)/1e6:.1f} kt', ha='left', fontsize=8, color='#C62828')\n\n # Engine nozzles\n for x in [-17, -16, -15, 15, 16, 17]:\n ax.plot([x, x], [5, 2], color='#333', lw=2)\n for x in np.linspace(-10, 10, 7):\n ax.plot([x, x], [5, 2], color='#555', lw=1.5)\n\n # Dimensions\n ax.annotate('', xy=(20, 5), xytext=(20, 220),\n arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))\n ax.text(21.5, 112, f'{HEIGHT:.0f} m', ha='left', va='center', fontsize=10, rotation=90)\n ax.annotate('', xy=(-20, 0), xytext=(20, 0),\n arrowprops=dict(arrowstyle='<->', color='black', lw=1.2))\n ax.text(0, -2, f'∅ {DIAM:.0f} m', ha='center', fontsize=10)\n\n # Payload callout\n ax.annotate(f'有效载荷: {PAYLOAD/1e3:.0f} t\\n(LEO 200 km)',\n xy=(0, 195), xytext=(14, 210),\n arrowprops=dict(arrowstyle='->', color='#0D47A1', lw=1.5),\n fontsize=9, color='#0D47A1', fontweight='bold',\n bbox=dict(boxstyle='round,pad=0.3', fc='#E3F2FD', ec='#0D47A1'))\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_vehicle_config.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_vehicle_config.png'))\n plt.close(fig)\n print('[OK] fig_vehicle_config')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 2 — Trajectory Simulation\n# ═══════════════════════════════════════════════════════════════════\ndef simulate_trajectory():\n \"\"\"3-DOF trajectory with gravity turn and atmospheric drag.\"\"\"\n def atm_density(h):\n if h < 0: return 1.225\n if h > 200000: return 0.0\n # US Standard Atmosphere simplified\n if h < 11000:\n T = 288.15 - 0.0065*h\n return 1.225 * (T/288.15)**4.256\n elif h < 25000:\n return 0.3639 * np.exp(-(h-11000)/6341.6)\n elif h < 50000:\n return 0.04 * np.exp(-(h-25000)/7500)\n else:\n return 0.001 * np.exp(-(h-50000)/7000)\n\n Cd = 0.35\n A_cross = np.pi * (DIAM/2)**2\n\n # Simulation state: [r, theta, v_r, v_theta, mass, stage_flag]\n # stage_flag: 0=S0, 1=S1, 2=S2, 3=coast\n dt = 0.5\n t_max = 1200\n\n state = {\n 'r': R_E, 'theta': 0, 'v_r': 0, 'v_theta': 0,\n 'mass': M0, 'stage': 0, 't': 0,\n 'pitch_angle': np.pi/2, # from local horizontal\n }\n history = {'t':[], 'h':[], 'v':[], 'a':[], 'm':[], 'q':[], 'downrange':[],\n 'Ma':[], 'stage':[], 'thrust':[], 'drag':[], 'gamma':[]}\n\n # Pitch program\n def pitch(t, stage):\n if t < 10: return np.pi/2\n if stage == 0: return np.pi/2 - 0.8 * min((t-10)/S0_BURN, 1.0)\n if stage == 1: return np.pi/2 - 0.8 - 0.5 * min((t - S0_BURN)/S1_BURN, 1.0)\n return max(0.05, np.pi/2 - 1.3 - 0.2 * min((t - S0_BURN - S1_BURN)/S2_BURN, 1.0))\n\n t = 0\n s0_burnt = 0; s1_burnt = 0; s2_burnt = 0\n s0_sep = False; s1_sep = False\n\n while t < t_max:\n h = state['r'] - R_E\n if h < -100: break\n rho = atm_density(h)\n v = np.sqrt(state['v_r']**2 + state['v_theta']**2)\n\n # Speed of sound\n if h < 11000:\n T_atm = max(216.65, 288.15 - 0.0065*h)\n elif h < 50000:\n T_atm = 216.65\n else:\n T_atm = max(180, 270 - 0.001*h)\n a_sound = np.sqrt(1.4 * 287 * T_atm)\n Ma = v / a_sound if a_sound > 0 else 0\n\n q = 0.5 * rho * v**2 # dynamic pressure\n\n stage = state['stage']\n # Determine thrust and mass flow\n thrust = 0; mdot = 0\n if stage == 0:\n isp = S0_ISP_SL + (S0_ISP_VAC - S0_ISP_SL) * min(h/40000, 1.0)\n remaining = S0_PROP - s0_burnt\n if remaining > 0:\n mdot = S0_PROP / S0_BURN\n thrust = mdot * isp * g0\n s0_burnt += mdot * dt\n else:\n # Separation\n state['mass'] -= S0_DRY\n state['stage'] = 1; s0_sep = True; stage = 1\n if stage == 1:\n isp = S1_ISP_VAC\n remaining = S1_PROP - s1_burnt\n if remaining > 0:\n mdot = S1_PROP / S1_BURN\n thrust = mdot * isp * g0\n s1_burnt += mdot * dt\n else:\n state['mass'] -= S1_DRY\n state['stage'] = 2; s1_sep = True; stage = 2\n if stage == 2:\n isp = S2_ISP_VAC\n remaining = S2_PROP - s2_burnt\n if remaining > 0:\n mdot = S2_PROP / S2_BURN\n thrust = mdot * isp * g0\n s2_burnt += mdot * dt\n else:\n state['mass'] -= S2_DRY\n state['stage'] = 3; stage = 3\n\n gamma = pitch(t, stage) # flight path angle from horizontal\n if gamma > np.pi: gamma = np.pi/2\n\n # Drag\n D = 0.5 * rho * v**2 * Cd * A_cross if v > 0 else 0\n\n # Gravity\n g_local = G_ME / state['r']**2\n\n # Accelerations (simplified 2D)\n if v > 0:\n a_thrust = thrust / state['mass']\n a_drag = D / state['mass']\n\n v_h = v * np.cos(gamma)\n v_v = v * np.sin(gamma)\n\n a_h = (a_thrust * np.cos(gamma)) - a_drag * (v_h/v if v>0 else 0)\n a_v = (a_thrust * np.sin(gamma)) - g_local - a_drag * (v_v/v if v>0 else 0)\n\n state['v_r'] += a_v * dt\n state['v_theta'] += a_h * dt\n else:\n a_thrust = thrust / state['mass'] if state['mass'] > 0 else 0\n a_v = a_thrust - g_local\n state['v_r'] += a_v * dt\n\n state['r'] += state['v_r'] * dt\n state['theta'] += (state['v_theta'] / state['r']) * dt\n state['mass'] -= mdot * dt\n if state['mass'] < PAYLOAD: state['mass'] = PAYLOAD\n\n t += dt\n v_total = np.sqrt(state['v_r']**2 + state['v_theta']**2)\n a_total = np.sqrt(a_h**2 + a_v**2) if v > 0 else abs(a_v)\n\n history['t'].append(t)\n history['h'].append(max(0, state['r']-R_E)/1000) # km\n history['v'].append(v_total)\n history['a'].append(a_total/g0)\n history['m'].append(state['mass']/1e6)\n history['q'].append(q/1000) # kPa\n history['downrange'].append(state['theta']*R_E/1000) # km\n history['Ma'].append(Ma)\n history['stage'].append(stage)\n history['thrust'].append(thrust/1e6) # MN\n history['drag'].append(D/1e6)\n history['gamma'].append(np.degrees(gamma))\n\n # Orbital check\n if state['r'] - R_E > 200000 and v_total > 7700:\n break\n\n return history\n\ndef fig_trajectory(hist):\n fig, axes = plt.subplots(2, 3, figsize=(16, 10))\n fig.suptitle('Leviathan-100 弹道仿真结果', fontsize=15, fontweight='bold')\n\n t = np.array(hist['t'])\n # Color by stage\n colors = {0:'#EF5350', 1:'#FF9800', 2:'#4CAF50', 3:'#2196F3'}\n stages = np.array(hist['stage'])\n\n def scatter_by_stage(ax, x, y, xlabel, ylabel, title):\n for s in [0,1,2,3]:\n mask = stages == s\n if mask.any():\n ax.scatter(x[mask], y[mask], c=colors[s], s=1, label=f'级{s}')\n ax.set_xlabel(xlabel); ax.set_ylabel(ylabel)\n ax.set_title(title); ax.legend(markerscale=5, fontsize=7)\n ax.grid(True, alpha=0.3)\n\n # Altitude vs time\n scatter_by_stage(axes[0,0], t, hist['h'], '时间 (s)', '高度 (km)', '(a) 高度-时间曲线')\n # Velocity vs time\n scatter_by_stage(axes[0,1], t, [v/1000 for v in hist['v']], '时间 (s)', '速度 (km/s)', '(b) 速度-时间曲线')\n # Acceleration vs time\n scatter_by_stage(axes[0,2], t, hist['a'], '时间 (s)', '加速度 (g)', '(c) 加速度-时间曲线')\n # Dynamic pressure\n axes[1,0].fill_between(t, 0, hist['q'], alpha=0.4, color='#E53935')\n axes[1,0].plot(t, hist['q'], color='#B71C1C', lw=0.8)\n axes[1,0].set_xlabel('时间 (s)'); axes[1,0].set_ylabel('动压 (kPa)')\n axes[1,0].set_title('(d) 动压 (Max-Q) 曲线'); axes[1,0].grid(True, alpha=0.3)\n maxq_idx = np.argmax(hist['q'])\n axes[1,0].annotate(f'Max-Q: {hist[\"q\"][maxq_idx]:.1f} kPa\\nt={t[maxq_idx]:.0f}s, h={hist[\"h\"][maxq_idx]:.1f}km',\n xy=(t[maxq_idx], hist['q'][maxq_idx]),\n xytext=(t[maxq_idx]+50, hist['q'][maxq_idx]*0.8),\n arrowprops=dict(arrowstyle='->', color='black'),\n fontsize=8, bbox=dict(boxstyle='round', fc='wheat', alpha=0.8))\n\n # Altitude vs downrange\n scatter_by_stage(axes[1,1], hist['downrange'], hist['h'], '下段距离 (km)', '高度 (km)', '(e) 弹道轨迹')\n # Mass vs time\n axes[1,2].plot(t, hist['m'], color='#6A1B9A', lw=1.5)\n axes[1,2].set_xlabel('时间 (s)'); axes[1,2].set_ylabel('质量 (kt)')\n axes[1,2].set_title('(f) 质量变化曲线'); axes[1,2].grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_trajectory.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_trajectory.png'))\n plt.close(fig)\n print('[OK] fig_trajectory')\n return hist\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 3 — Mass Budget\n# ═══════════════════════════════════════════════════════════════════\ndef fig_mass_budget():\n fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 7))\n\n # Overall mass breakdown\n labels = ['助推级推进剂', '助推级结构', '一级推进剂', '一级结构', '二级推进剂', '二级结构', '有效载荷']\n sizes = [S0_PROP/1e6, S0_DRY/1e6, S1_PROP/1e6, S1_DRY/1e6, S2_PROP/1e6, S2_DRY/1e6, PAYLOAD/1e3]\n colors = ['#EF5350','#E57373','#FF9800','#FFB74D','#4CAF50','#81C784','#2196F3']\n explode = [0,0,0,0,0,0,0.08]\n\n wedges, texts, autotexts = ax1.pie(sizes, labels=labels, colors=colors, explode=explode,\n autopct='%1.1f%%', pctdistance=0.8, startangle=90)\n for t in autotexts: t.set_fontsize(8)\n for t in texts: t.set_fontsize(9)\n ax1.set_title('(a) 总体质量分配', fontsize=13, fontweight='bold')\n\n # Stage-by-stage breakdown\n categories = ['助推级', '一级', '二级', '有效载荷']\n propellant = [S0_PROP/1e6, S1_PROP/1e6, S2_PROP/1e6, 0]\n structure = [S0_DRY/1e6, S1_DRY/1e6, S2_DRY/1e6, 0]\n payload = [0, 0, 0, PAYLOAD/1e3]\n\n x = np.arange(len(categories))\n w = 0.5\n ax2.bar(x, propellant, w, label='推进剂', color='#EF5350', alpha=0.85)\n ax2.bar(x, structure, w, bottom=propellant, label='结构', color='#FFB74D', alpha=0.85)\n ax2.bar(x, payload, w, bottom=[p+s for p,s in zip(propellant,structure)], label='有效载荷', color='#2196F3', alpha=0.85)\n ax2.set_xticks(x); ax2.set_xticklabels(categories)\n ax2.set_ylabel('质量 (kt)'); ax2.set_title('(b) 各级质量构成', fontsize=13, fontweight='bold')\n ax2.legend(); ax2.grid(True, alpha=0.3, axis='y')\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_mass_budget.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_mass_budget.png'))\n plt.close(fig)\n print('[OK] fig_mass_budget')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 4 — Propulsion Architecture\n# ═══════════════════════════════════════════════════════════════════\ndef fig_propulsion():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n # Engine layout — top view\n ax = axes[0]\n ax.set_xlim(-15, 15); ax.set_ylim(-15, 15); ax.set_aspect('equal')\n ax.set_title('(a) 助推级发动机布局 (俯视)', fontsize=11, fontweight='bold')\n # Outer ring — 24 engines\n for i in range(24):\n angle = 2*np.pi*i/24\n x, y = 12*np.cos(angle), 12*np.sin(angle)\n c = Circle((x,y), 1.2, fc='#EF5350', ec='#B71C1C', lw=1)\n ax.add_patch(c)\n # Middle ring — 12 engines\n for i in range(12):\n angle = 2*np.pi*i/12 + np.pi/12\n x, y = 8*np.cos(angle), 8*np.sin(angle)\n c = Circle((x,y), 1.2, fc='#FF9800', ec='#E65100', lw=1)\n ax.add_patch(c)\n # Inner cluster — 7 engines\n for i in range(6):\n angle = 2*np.pi*i/6\n x, y = 4*np.cos(angle), 4*np.sin(angle)\n c = Circle((x,y), 1.2, fc='#4CAF50', ec='#1B5E20', lw=1)\n ax.add_patch(c)\n c = Circle((0,0), 1.2, fc='#4CAF50', ec='#1B5E20', lw=1)\n ax.add_patch(c)\n\n ax.text(0, -14.5, '外圈: 24× RP-1/LOX | 中圈: 12× LCH₄/LOX | 内圈: 7× LCH₄/LOX (矢量)',\n ha='center', fontsize=7, style='italic')\n ax.axis('off')\n\n # Thrust profile\n ax = axes[1]\n t = np.linspace(0, 940, 1000)\n thrust_profile = np.zeros_like(t)\n for i, ti in enumerate(t):\n if ti < S0_BURN:\n thrust_profile[i] = (S0_PROP/S0_BURN * S0_ISP_SL * g0 + S1_PROP/S1_BURN * 0 * S1_ISP_VAC * g0) / 1e6\n elif ti < S0_BURN + S1_BURN:\n thrust_profile[i] = (S1_PROP/S1_BURN * S1_ISP_VAC * g0) / 1e6\n elif ti < S0_BURN + S1_BURN + S2_BURN:\n thrust_profile[i] = (S2_PROP/S2_BURN * S2_ISP_VAC * g0) / 1e6\n else:\n thrust_profile[i] = 0\n\n ax.fill_between(t, 0, thrust_profile, alpha=0.3, color='#E53935')\n ax.plot(t, thrust_profile, color='#B71C1C', lw=2)\n ax.set_xlabel('时间 (s)'); ax.set_ylabel('推力 (MN)')\n ax.set_title('(b) 推力-时间曲线', fontsize=11, fontweight='bold')\n ax.grid(True, alpha=0.3)\n # Stage annotations\n ax.axvline(S0_BURN, color='gray', ls='--', lw=0.8)\n ax.axvline(S0_BURN+S1_BURN, color='gray', ls='--', lw=0.8)\n ax.text(80, max(thrust_profile)*0.9, '助推级', fontsize=9, ha='center')\n ax.text(S0_BURN+150, max(thrust_profile[S0_BURN:])*0.9, '一级', fontsize=9, ha='center')\n ax.text(S0_BURN+S1_BURN+240, max(thrust_profile[S0_BURN+S1_BURN:])*0.9, '二级', fontsize=9, ha='center')\n\n # Isp comparison\n ax = axes[2]\n engines = ['F-1\\n(RP-1)', 'Raptor 3\\n(LCH₄)', 'RS-25\\n(LH₂)', 'NERVA\\n(NTP)', 'Leviathan\\n二级']\n isps = [263, 327, 452, 825, 465]\n colors_bar = ['#EF5350','#FF9800','#4CAF50','#9C27B0','#2196F3']\n bars = ax.bar(engines, isps, color=colors_bar, alpha=0.85, edgecolor='black', lw=0.5)\n ax.set_ylabel('比冲 (s)'); ax.set_title('(c) 比冲对比', fontsize=11, fontweight='bold')\n for bar, isp in zip(bars, isps):\n ax.text(bar.get_x() + bar.get_width()/2, bar.get_height()+10, f'{isp}s',\n ha='center', fontsize=9, fontweight='bold')\n ax.grid(True, alpha=0.3, axis='y')\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_propulsion.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_propulsion.png'))\n plt.close(fig)\n print('[OK] fig_propulsion')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 5 — Structural Analysis\n# ═══════════════════════════════════════════════════════════════════\ndef fig_structural():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n # Tank wall thickness vs altitude\n ax = axes[0]\n heights = np.linspace(0, 200, 200) # m along vehicle\n # Simplified: hoop stress determines wall thickness\n p_tank = 0.35 # MPa internal pressure\n sigma_allow_steel = 550 # MPa for 301 stainless\n sigma_allow_cfrp = 1200 # MPa for CFRP overwrap\n r_tank = DIAM/2\n\n t_steel = p_tank * r_tank / sigma_allow_steel * 1000 # mm\n t_cfrp = p_tank * r_tank / sigma_allow_cfrp * 1000\n t_hybrid = 0.6 * t_steel # COPV hybrid\n\n ax.axhline(t_steel, color='#EF5350', ls='-', lw=2, label=f'不锈钢 301: {t_steel:.1f} mm')\n ax.axhline(t_hybrid, color='#4CAF50', ls='-', lw=2, label=f'COPV混合: {t_hybrid:.1f} mm')\n ax.axhline(t_cfrp, color='#2196F3', ls='-', lw=2, label=f'纯CFRP: {t_cfrp:.1f} mm')\n ax.set_xlabel('贮箱位置 (m)'); ax.set_ylabel('壁厚 (mm)')\n ax.set_title('(a) 贮箱壁厚方案对比', fontsize=11, fontweight='bold')\n ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n # Load distribution along vehicle\n ax = axes[1]\n z = np.linspace(0, HEIGHT, 200)\n # Simplified axial load: weight above * g + thrust\n # At bottom: full weight; at top: only payload\n load_fraction = 1 - z/HEIGHT\n axial_load = load_fraction * M0 * g0 / 1e9 # GN\n bending_load = 0.05 * axial_load * np.sin(2*np.pi*z/HEIGHT) # simplified wind bending\n\n ax.fill_between(z, 0, axial_load, alpha=0.4, color='#E53935', label='轴压')\n ax.fill_between(z, -bending_load, bending_load, alpha=0.3, color='#2196F3', label='弯矩')\n ax.set_xlabel('位置 (m, 从底部)'); ax.set_ylabel('载荷 (GN)')\n ax.set_title('(b) 轴向载荷分布', fontsize=11, fontweight='bold')\n ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n # Structural mass fraction comparison\n ax = axes[2]\n vehicles = ['Saturn V', 'Starship\\n(BFR)', 'Sea Dragon', 'Leviathan\\n-100']\n struct_frac = [0.081, 0.065, 0.105, 0.103]\n colors_s = ['#FF9800', '#2196F3', '#9C27B0', '#F44336']\n bars = ax.bar(vehicles, [x*100 for x in struct_frac], color=colors_s, alpha=0.85, edgecolor='black', lw=0.5)\n for bar, frac in zip(bars, struct_frac):\n ax.text(bar.get_x()+bar.get_width()/2, bar.get_height()+0.2, f'{frac*100:.1f}%',\n ha='center', fontsize=9, fontweight='bold')\n ax.set_ylabel('结构质量分数 (%)'); ax.set_title('(c) 结构效率对比', fontsize=11, fontweight='bold')\n ax.grid(True, alpha=0.3, axis='y')\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_structural.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_structural.png'))\n plt.close(fig)\n print('[OK] fig_structural')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 6 — Aerodynamic Analysis\n# ═══════════════════════════════════════════════════════════════════\ndef fig_aerodynamics():\n fig, axes = plt.subplots(2, 2, figsize=(14, 10))\n\n Ma_range = np.linspace(0.1, 25, 500)\n\n # Drag coefficient vs Mach\n ax = axes[0,0]\n Cd = np.zeros_like(Ma_range)\n for i, Ma in enumerate(Ma_range):\n if Ma < 0.8: Cd[i] = 0.25\n elif Ma < 1.2: Cd[i] = 0.25 + 0.55 * np.sin(np.pi*(Ma-0.8)/0.8)\n elif Ma < 5: Cd[i] = 0.8 / np.sqrt(Ma)\n else: Cd[i] = 0.35 / np.sqrt(Ma)\n ax.plot(Ma_range, Cd, 'b-', lw=2)\n ax.axvline(1.0, color='r', ls='--', lw=0.8, alpha=0.5, label='Ma=1')\n ax.set_xlabel('马赫数'); ax.set_ylabel('阻力系数 Cd')\n ax.set_title('(a) 阻力系数-马赫数曲线'); ax.legend(); ax.grid(True, alpha=0.3)\n\n # Aerodynamic heating\n ax = axes[0,1]\n h_range = np.linspace(0, 120, 200) # km\n v_range = np.linspace(0, 8, 200) # km/s\n H, V = np.meshgrid(h_range, v_range)\n # Stagnation heating rate (Sutton-Graves approximation)\n k_sg = 1.7415e-4\n rho = np.exp(-H/7.1) * 1.225 # simplified\n q_dot = k_sg * np.sqrt(rho / (DIAM/2)) * (V*1000)**3 / 1e6 # MW/m²\n q_dot = np.clip(q_dot, 0, 50)\n c = ax.contourf(H, V, q_dot, levels=20, cmap='hot')\n fig.colorbar(c, ax=ax, label='热流密度 (MW/m²)')\n ax.set_xlabel('高度 (km)'); ax.set_ylabel('速度 (km/s)')\n ax.set_title('(b) 气动加热热流分布')\n\n # Wind load\n ax = axes[1,0]\n alt_wind = np.linspace(0, 30, 100)\n v_wind = 10 + 30 * np.sin(np.pi * alt_wind / 30) + np.random.normal(0, 3, 100)\n ax.fill_between(alt_wind, 0, v_wind, alpha=0.3, color='#2196F3')\n ax.plot(alt_wind, v_wind, 'b-', lw=1.5)\n ax.set_xlabel('高度 (km)'); ax.set_ylabel('风速 (m/s)')\n ax.set_title('(c) 设计风剖面'); ax.grid(True, alpha=0.3)\n\n # Max-Q envelope\n ax = axes[1,1]\n t_q = np.linspace(0, 150, 300)\n # Typical max-q profile\n q_env = 40 * np.sin(np.pi * t_q / 80)**2 * np.exp(-0.005 * t_q)\n ax.plot(t_q, q_env, 'r-', lw=2)\n ax.fill_between(t_q, 0, q_env, alpha=0.2, color='red')\n maxq_t = t_q[np.argmax(q_env)]\n ax.annotate(f'Max-Q = {max(q_env):.1f} kPa\\nt ≈ {maxq_t:.0f}s',\n xy=(maxq_t, max(q_env)), xytext=(maxq_t+20, max(q_env)*0.7),\n arrowprops=dict(arrowstyle='->'), fontsize=9,\n bbox=dict(boxstyle='round', fc='wheat', alpha=0.8))\n ax.set_xlabel('时间 (s)'); ax.set_ylabel('动压 (kPa)')\n ax.set_title('(d) 动压包络'); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_aerodynamics.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_aerodynamics.png'))\n plt.close(fig)\n print('[OK] fig_aerodynamics')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 7 — Delta-V Budget\n# ═══════════════════════════════════════════════════════════════════\ndef fig_deltav():\n fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 7))\n\n # Delta-V by stage\n dv0 = S0_ISP_VAC * g0 * np.log((M0) / (M0 - S0_PROP))\n m1_start = M0 - S0_PROP - S0_DRY\n dv1 = S1_ISP_VAC * g0 * np.log(m1_start / (m1_start - S1_PROP))\n m2_start = m1_start - S1_PROP - S1_DRY\n dv2 = S2_ISP_VAC * g0 * np.log(m2_start / (m2_start - S2_PROP))\n\n losses = {'gravity': 1800, 'drag': 350, 'steering': 200}\n dv_leo = 9400 # m/s for 200 km LEO\n\n stages = ['助推级\\nΔv', '一级\\nΔv', '二级\\nΔv', '重力损失', '阻力损失', '机动损失', '轨道速度\\n需求']\n values = [dv0, dv1, dv2, -losses['gravity'], -losses['drag'], -losses['steering'], dv_leo]\n colors_dv = ['#EF5350','#FF9800','#4CAF50','#9E9E9E','#9E9E9E','#9E9E9E','#2196F3']\n\n bars = ax1.barh(stages, [v/1000 for v in values], color=colors_dv, alpha=0.85, edgecolor='black', lw=0.5)\n for bar, val in zip(bars, values):\n x_pos = val/1000 + 0.1 if val >= 0 else val/1000 - 0.3\n ax1.text(x_pos, bar.get_y()+bar.get_height()/2, f'{val/1000:.2f} km/s',\n va='center', fontsize=9, fontweight='bold')\n ax1.set_xlabel('Δv (km/s)'); ax1.set_title('(a) Δv 预算', fontsize=13, fontweight='bold')\n ax1.axvline(0, color='black', lw=0.5)\n ax1.grid(True, alpha=0.3, axis='x')\n\n # Payload vs orbit\n ax2.set_title('(b) 运载能力-轨道关系', fontsize=13, fontweight='bold')\n orbits = ['LEO\\n200km', 'SSO\\n700km', 'GTO', 'TLI', 'TMI']\n # Approximate delta-v from LEO ref\n dv_from_leo = [0, 0.8, 2.5, 3.2, 4.3] # km/s above LEO\n payload_est = []\n for extra_dv in dv_from_leo:\n # Use Tsiolkovsky with decreasing mass\n remaining_dv = dv_leo + extra_dv * 1000 - dv0 - dv1\n if remaining_dv > 0:\n # Need more from stage 2\n mass_ratio = np.exp(remaining_dv / (S2_ISP_VAC * g0))\n prop_needed = m2_start * (1 - 1/mass_ratio)\n pl = m2_start - S2_PROP - S2_DRY - prop_needed\n else:\n # Stage 2 delta-v covers it\n pl = m2_start - S2_PROP - S2_DRY\n # But reduce payload for higher orbits\n avail_dv = dv2 + dv1 - (dv_leo + extra_dv*1000 - dv0)\n if avail_dv < 0:\n pl = m2_start - S2_PROP - S2_DRY\n mass_ratio = np.exp(abs(avail_dv) / (S2_ISP_VAC * g0))\n pl = max(0, pl / mass_ratio)\n payload_est.append(max(0, pl/1e3))\n\n ax2.bar(orbits, payload_est, color=['#2196F3','#4CAF50','#FF9800','#9C27B0','#F44336'],\n alpha=0.85, edgecolor='black', lw=0.5)\n for i, (orb, pl) in enumerate(zip(orbits, payload_est)):\n ax2.text(i, pl+20, f'{pl:.0f} t', ha='center', fontsize=10, fontweight='bold')\n ax2.set_ylabel('有效载荷 (t)'); ax2.grid(True, alpha=0.3, axis='y')\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_deltav.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_deltav.png'))\n plt.close(fig)\n print('[OK] fig_deltav')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 8 — Comparison with Historical Vehicles\n# ═══════════════════════════════════════════════════════════════════\ndef fig_comparison():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n vehicles = ['Saturn V', 'Energia', 'Starship', 'Sea Dragon', 'Leviathan-100']\n liftoff_mass = [2970, 2400, 5000, 18143, 100000]\n leo_payload = [140, 100, 150, 550, PAYLOAD/1e3]\n diameter = [10.1, 7.75, 9, 23, 25]\n height = [110.6, 58.8, 121, 150, 220]\n\n # Liftoff mass comparison (log scale)\n ax = axes[0]\n colors_c = ['#FF9800','#4CAF50','#2196F3','#9C27B0','#F44336']\n bars = ax.barh(vehicles, liftoff_mass, color=colors_c, alpha=0.85, edgecolor='black', lw=0.5)\n for bar, m in zip(bars, liftoff_mass):\n ax.text(m+1000, bar.get_y()+bar.get_height()/2, f'{m:,} t', va='center', fontsize=9)\n ax.set_xlabel('起飞质量 (t)'); ax.set_title('(a) 起飞质量对比', fontsize=11, fontweight='bold')\n ax.set_xscale('log'); ax.grid(True, alpha=0.3, axis='x')\n\n # LEO payload comparison\n ax = axes[1]\n bars = ax.barh(vehicles, leo_payload, color=colors_c, alpha=0.85, edgecolor='black', lw=0.5)\n for bar, p in zip(bars, leo_payload):\n ax.text(p+10, bar.get_y()+bar.get_height()/2, f'{p:,.0f} t', va='center', fontsize=9)\n ax.set_xlabel('LEO有效载荷 (t)'); ax.set_title('(b) LEO运载能力对比', fontsize=11, fontweight='bold')\n ax.grid(True, alpha=0.3, axis='x')\n\n # Size comparison (silhouette)\n ax = axes[2]\n max_h = max(height)\n for i, (v, h, d) in enumerate(zip(vehicles, height, diameter)):\n scale = d / max(diameter)\n y_base = i * 0.2\n rect = Rectangle((0.5 - scale*0.15, y_base), scale*0.3, h/max_h * 0.18,\n fc=colors_c[i], ec='black', lw=0.5, alpha=0.7)\n ax.add_patch(rect)\n ax.text(0.85, y_base + h/max_h * 0.09, f'{v}\\n{h}m / ∅{d}m',\n va='center', fontsize=7)\n ax.set_xlim(0, 1.5); ax.set_ylim(-0.05, 1.05)\n ax.set_title('(c) 尺寸对比', fontsize=11, fontweight='bold')\n ax.axis('off')\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_comparison.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_comparison.png'))\n plt.close(fig)\n print('[OK] fig_comparison')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 9 — Thermal Protection\n# ═══════════════════════════════════════════════════════════════════\ndef fig_thermal():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n # Temperature distribution along vehicle\n ax = axes[0]\n z = np.linspace(0, 220, 200)\n # Simplified temperature profile during max-Q reentry\n T_stagnation = 2500 * np.exp(-((z-0)/40)**2) + 300\n T_body = 800 * np.exp(-((z-0)/60)**2) + 220\n T_internal = 20 + 50 * np.exp(-((z-110)/80)**2)\n\n ax.plot(z, T_stagnation, 'r-', lw=2, label='驻点温度')\n ax.plot(z, T_body, 'b-', lw=2, label='壁面温度')\n ax.plot(z, T_internal, 'g--', lw=1.5, label='内部温度')\n ax.axhline(1500, color='orange', ls=':', lw=1, label='钢熔点 (≈1500°C)')\n ax.set_xlabel('位置 (m)'); ax.set_ylabel('温度 (°C)')\n ax.set_title('(a) 热环境分布'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n # TPS thickness\n ax = axes[1]\n z_tps = np.linspace(0, 100, 100)\n tps_cork = 15 + 40 * np.exp(-z_tps/30)\n tps_ablative = 5 + 25 * np.exp(-z_tps/25)\n tps_ceramic = 2 + 8 * np.exp(-z_tps/20)\n\n ax.fill_between(z_tps, 0, tps_cork, alpha=0.3, color='#EF5350', label='软木隔热层')\n ax.fill_between(z_tps, 0, tps_ablative, alpha=0.3, color='#FF9800', label='烧蚀层')\n ax.fill_between(z_tps, 0, tps_ceramic, alpha=0.3, color='#4CAF50', label='陶瓷瓦')\n ax.plot(z_tps, tps_cork, 'r-', lw=1.5)\n ax.plot(z_tps, tps_ablative, 'orange', lw=1.5)\n ax.plot(z_tps, tps_ceramic, 'g-', lw=1.5)\n ax.set_xlabel('距鼻锥距离 (m)'); ax.set_ylabel('TPS厚度 (mm)')\n ax.set_title('(b) 热防护厚度'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n # Heat flux timeline\n ax = axes[2]\n t_heat = np.linspace(0, 940, 500)\n q_stag = np.zeros_like(t_heat)\n for i, ti in enumerate(t_heat):\n if ti < 160: # ascent\n q_stag[i] = 0.5 * (1 - np.exp(-ti/20)) * np.exp(-(ti-80)**2/2000) * 8\n elif ti > 600: # hypothetical reentry\n q_stag[i] = 15 * np.exp(-((ti-700)/50)**2)\n\n ax.semilogy(t_heat, np.clip(q_stag, 1e-3, 100), 'r-', lw=2)\n ax.set_xlabel('时间 (s)'); ax.set_ylabel('热流密度 (MW/m²)')\n ax.set_title('(c) 热流-时间历程'); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_thermal.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_thermal.png'))\n plt.close(fig)\n print('[OK] fig_thermal')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 10 — Reliability & Risk\n# ═══════════════════════════════════════════════════════════════════\ndef fig_risk():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n # Fault tree\n ax = axes[0]\n ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.axis('off')\n ax.set_title('(a) 故障树分析 (简化)', fontsize=11, fontweight='bold')\n boxes = {\n 'top': (5, 9, '任务失败'),\n 'prop': (2, 6.5, '推进故障'),\n 'struct': (5, 6.5, '结构故障'),\n 'avion': (8, 6.5, '控制故障'),\n 'e1': (0.5, 4, '发动机\\n失效'),\n 'e2': (3.5, 4, '推进剂\\n泄漏'),\n 'e3': (5, 4, '贮箱\\n破裂'),\n 'e4': (8, 4, 'GNC\\n失锁'),\n }\n for key, (x, y, txt) in boxes.items():\n w = 1.8 if key == 'top' else 1.6\n h = 0.8\n rect = FancyBboxPatch((x-w/2, y-h/2), w, h, boxstyle=\"round,pad=0.1\",\n fc='#FFCDD2' if y > 5 else '#FFF9C4', ec='#B71C1C', lw=1.2)\n ax.add_patch(rect)\n ax.text(x, y, txt, ha='center', va='center', fontsize=7, fontweight='bold')\n # Connections\n for child in ['prop', 'struct', 'avion']:\n ax.annotate('', xy=(5, 8.6), xytext=(boxes[child][0], boxes[child][1]+0.4),\n arrowprops=dict(arrowstyle='->', color='#555', lw=1))\n for child in ['e1', 'e2']:\n ax.annotate('', xy=(2, 6.1), xytext=(boxes[child][0], boxes[child][1]+0.4),\n arrowprops=dict(arrowstyle='->', color='#555', lw=1))\n ax.annotate('', xy=(5, 6.1), xytext=(boxes['e3'][0], boxes['e3'][1]+0.4),\n arrowprops=dict(arrowstyle='->', color='#555', lw=1))\n ax.annotate('', xy=(8, 6.1), xytext=(boxes['e4'][0], boxes['e4'][1]+0.4),\n arrowprops=dict(arrowstyle='->', color='#555', lw=1))\n\n # Reliability vs engine count\n ax = axes[1]\n n_engines = np.arange(1, 60)\n p_single = 0.998 # single engine reliability\n # System reliability with N-1 redundancy (need all or N-1)\n p_system_all = p_single ** n_engines\n p_system_n1 = 1 - (1 - p_single) * n_engines * p_single**(n_engines-1) - (1-p_single)**n_engines\n # With engine-out capability (survive losing 1)\n p_system_out = []\n for n in n_engines:\n if n <= 1:\n p_system_out.append(p_single)\n else:\n # Survive if at most 1 engine fails out of n\n p = p_single**n + n * (1-p_single) * p_single**(n-1)\n p_system_out.append(p)\n\n ax.plot(n_engines, p_system_all, 'r-', lw=2, label='无容错 (全部工作)')\n ax.plot(n_engines, p_system_out, 'b-', lw=2, label='单机容错 (允许1台失效)')\n ax.axvline(43, color='green', ls='--', lw=1, label='Leviathan-100 (43台)')\n ax.set_xlabel('发动机数量'); ax.set_ylabel('系统可靠度')\n ax.set_title('(b) 推进系统可靠度'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n ax.set_ylim(0.9, 1.001)\n\n # Risk matrix\n ax = axes[2]\n risks = [\n ('发动机爆炸', 3, 5), ('贮箱破裂', 2, 5), ('控制系统故障', 3, 4),\n ('风切变', 4, 3), ('声振破坏', 3, 4), ('级间分离失败', 3, 4),\n ('热防护失效', 2, 5), ('发射台损毁', 2, 4),\n ]\n for name, prob, impact in risks:\n color = '#F44336' if prob*impact >= 15 else '#FF9800' if prob*impact >= 10 else '#4CAF50'\n ax.scatter(prob, impact, s=200, c=color, edgecolors='black', lw=0.8, zorder=5)\n ax.annotate(name, (prob, impact), textcoords=\"offset points\",\n xytext=(5, 5), fontsize=6)\n ax.set_xlabel('发生概率'); ax.set_ylabel('影响程度')\n ax.set_title('(c) 风险矩阵'); ax.grid(True, alpha=0.3)\n ax.set_xlim(0.5, 5.5); ax.set_ylim(0.5, 5.5)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_risk.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_risk.png'))\n plt.close(fig)\n print('[OK] fig_risk')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 11 — Cost Analysis\n# ═══════════════════════════════════════════════════════════════════\ndef fig_cost():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n # Development cost breakdown\n ax = axes[0]\n items = ['推进系统', '结构/贮箱', '航电/GNC', '发射设施', '地面支持', '试验验证', '项目管理']\n costs = [8.5, 4.2, 1.8, 6.0, 2.5, 3.0, 2.0] # billion CNY\n colors_cost = ['#EF5350','#FF9800','#4CAF50','#2196F3','#9C27B0','#795548','#607D8B']\n wedges, texts, autotexts = ax.pie(costs, labels=items, colors=colors_cost,\n autopct='%1.1f%%', startangle=90, pctdistance=0.8)\n for t in autotexts: t.set_fontsize(7)\n for t in texts: t.set_fontsize(7)\n ax.set_title(f'(a) 研制费用分配\\n(总计 {sum(costs):.1f} B CNY)', fontsize=11, fontweight='bold')\n\n # Cost per kg to LEO\n ax = axes[1]\n launchers = ['Saturn V', 'SLS', 'Starship\\n(est)', 'Falcon\\nHeavy', 'Leviathan\\n-100']\n cost_per_kg = [55000, 30000, 100, 1500, 850] # $/kg\n ax.bar(launchers, cost_per_kg, color=['#FF9800','#9C27B0','#2196F3','#4CAF50','#F44336'],\n alpha=0.85, edgecolor='black', lw=0.5)\n ax.set_ylabel('$ / kg (LEO)'); ax.set_title('(b) 单位发射成本对比', fontsize=11, fontweight='bold')\n ax.set_yscale('log'); ax.grid(True, alpha=0.3, axis='y')\n\n # Learning curve\n ax = axes[2]\n n_flight = np.arange(1, 51)\n # 85% learning curve\n cost_first = 5.0 # billion CNY first launch\n cost_n = cost_first * n_flight ** (np.log10(0.85)/np.log10(2))\n cost_per_flight = cost_first * 0.85 ** np.log2(n_flight)\n\n ax.plot(n_flight, cost_per_flight, 'b-', lw=2, label='单次发射成本')\n ax.fill_between(n_flight, cost_per_flight*0.8, cost_per_flight*1.2, alpha=0.2, color='blue')\n ax.set_xlabel('发射序号'); ax.set_ylabel('单次发射成本 (B CNY)')\n ax.set_title('(c) 学习曲线 (85%率)'); ax.legend(); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_cost.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_cost.png'))\n plt.close(fig)\n print('[OK] fig_cost')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 12 — Launch Site & Infrastructure\n# ═══════════════════════════════════════════════════════════════════\ndef fig_launch_site():\n fig, axes = plt.subplots(1, 2, figsize=(14, 7))\n\n # Launch site layout\n ax = axes[0]\n ax.set_xlim(-500, 500); ax.set_ylim(-500, 500); ax.set_aspect('equal')\n ax.set_title('(a) 发射场布局示意', fontsize=11, fontweight='bold')\n # Launch pad\n pad = Rectangle((-30, -30), 60, 60, fc='#9E9E9E', ec='#333', lw=2)\n ax.add_patch(pad); ax.text(0, 0, '发射台', ha='center', fontsize=8)\n # Tower\n tower = Rectangle((35, -40), 15, 80, fc='#FF9800', ec='#E65100', lw=1.5)\n ax.add_patch(tower); ax.text(42, 0, '勤务塔', ha='center', fontsize=7, rotation=90)\n # Flame trench\n trench = Rectangle((-50, -80), 100, 40, fc='#795548', ec='#333', lw=1)\n ax.add_patch(trench); ax.text(0, -60, '导流槽', ha='center', fontsize=7, color='white')\n # Fuel storage\n for i, (x, label) in enumerate([(-200, 'LOX'), (-150, 'RP-1'), (150, 'LH₂'), (200, 'LCH₄')]):\n c = Circle((x, 200), 25, fc=['#2196F3','#EF5350','#4CAF50','#FF9800'][i], ec='black', lw=1)\n ax.add_patch(c); ax.text(x, 200, label, ha='center', fontsize=6, color='white')\n # Assembly building\n assembly = Rectangle((-250, -300), 120, 80, fc='#FFCDD2', ec='#C62828', lw=1.5)\n ax.add_patch(assembly); ax.text(-190, -260, '总装大楼', ha='center', fontsize=8)\n # Water deluge\n ax.add_patch(Circle((-100, -100), 40, fc='#64B5F6', ec='#1565C0', lw=1, alpha=0.5))\n ax.text(-100, -100, '喷水\\n系统', ha='center', fontsize=6)\n # Safety zone\n safe = Circle((0, 0), 400, fc='none', ec='red', ls='--', lw=1.5)\n ax.add_patch(safe); ax.text(280, -350, '安全区 (3km)', fontsize=7, color='red')\n ax.axis('off')\n\n # Acoustic environment\n ax = axes[1]\n distances = np.linspace(100, 5000, 200)\n # Sound pressure level at distance d: SPL = SPL_ref - 20*log10(d/d_ref) - atmospheric_absorption\n SPL_ref = 204 # dB at 1m (43 engines!)\n SPL = SPL_ref - 20 * np.log10(distances) - 0.005 * distances\n ax.plot(distances, SPL, 'r-', lw=2)\n ax.axhline(140, color='orange', ls='--', lw=1, label='140 dB: 结构损伤')\n ax.axhline(120, color='yellow', ls='--', lw=1, label='120 dB: 人体痛阈')\n ax.axhline(85, color='green', ls='--', lw=1, label='85 dB: 安全线')\n ax.set_xlabel('距离 (m)'); ax.set_ylabel('声压级 (dB)')\n ax.set_title('(b) 声环境分布'); ax.legend(fontsize=7); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_launch_site.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_launch_site.png'))\n plt.close(fig)\n print('[OK] fig_launch_site')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 13 — Monte Carlo Dispersion\n# ═══════════════════════════════════════════════════════════════════\ndef fig_montecarlo():\n fig, axes = plt.subplots(1, 3, figsize=(16, 6))\n\n np.random.seed(42)\n N = 2000\n\n # Payload dispersion\n mu_pl = PAYLOAD/1e3\n sigma_pl = mu_pl * 0.08\n pl_samples = np.random.normal(mu_pl, sigma_pl, N)\n ax = axes[0]\n ax.hist(pl_samples, bins=40, color='#2196F3', alpha=0.7, edgecolor='black', lw=0.5)\n ax.axvline(mu_pl, color='red', lw=2, label=f'均值: {mu_pl:.0f} t')\n ax.axvline(mu_pl - 2*sigma_pl, color='orange', ls='--', lw=1.5, label=f'95%: {mu_pl-2*sigma_pl:.0f} t')\n ax.set_xlabel('有效载荷 (t)'); ax.set_ylabel('频次')\n ax.set_title('(a) 载荷散布 (Monte Carlo)'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n # Insertion orbit dispersion\n mu_h = 200; sigma_h = 15\n h_samples = np.random.normal(mu_h, sigma_h, N)\n mu_inc = 28.5; sigma_inc = 0.3\n inc_samples = np.random.normal(mu_inc, sigma_inc, N)\n ax = axes[1]\n ax.scatter(h_samples, inc_samples, s=2, alpha=0.3, color='#4CAF50')\n ax.axhline(mu_inc, color='red', ls='--', lw=1)\n ax.axvline(mu_h, color='red', ls='--', lw=1)\n ax.set_xlabel('入轨高度 (km)'); ax.set_ylabel('轨道倾角 (°)')\n ax.set_title('(b) 入轨精度散布'); ax.grid(True, alpha=0.3)\n\n # Velocity at MECO dispersion\n mu_v = 7800; sigma_v = 50\n v_samples = np.random.normal(mu_v, sigma_v, N)\n ax = axes[2]\n ax.hist(v_samples, bins=40, color='#FF9800', alpha=0.7, edgecolor='black', lw=0.5)\n ax.axvline(7844, color='green', ls='--', lw=1.5, label='LEO 200km 理论值')\n ax.set_xlabel('关机速度 (m/s)'); ax.set_ylabel('频次')\n ax.set_title('(c) 关机速度散布'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_montecarlo.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_montecarlo.png'))\n plt.close(fig)\n print('[OK] fig_montecarlo')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 14 — Stage Separation Sequence\n# ═══════════════════════════════════════════════════════════════════\ndef fig_separation():\n fig, axes = plt.subplots(1, 4, figsize=(18, 5))\n\n phases = [\n ('T+0s 起飞', '全部43台发动机点火\\n推力≈1,250 MN\\n起飞质量100,000 t', '#EF5350'),\n ('T+160s 助推分离', '助推级关机分离\\n芯级继续推进\\n高度≈65 km', '#FF9800'),\n ('T+460s 一二级分离', '一级关机分离\\n二级点火\\n高度≈140 km', '#4CAF50'),\n ('T+940s 入轨', '二级关机\\n载荷入轨\\n轨道: 200km LEO', '#2196F3'),\n ]\n\n for ax, (title, desc, color) in zip(axes, phases):\n ax.set_xlim(-5, 5); ax.set_ylim(-2, 12); ax.set_aspect('equal')\n ax.axis('off')\n ax.set_title(title, fontsize=11, fontweight='bold', color=color)\n\n # Simplified rocket shape\n if '起飞' in title:\n # Full stack\n ax.add_patch(Rectangle((-3, 1), 6, 6, fc='#FFCDD2', ec='#C62828', lw=1.5, alpha=0.6))\n ax.add_patch(Rectangle((-2, 7), 4, 3, fc='#FFE0B2', ec='#E65100', lw=1.5, alpha=0.6))\n ax.add_patch(Rectangle((-1.5, 10), 3, 1.5, fc='#C8E6C9', ec='#1B5E20', lw=1.5, alpha=0.6))\n elif '助推' in title:\n # Boosters separating\n ax.add_patch(Rectangle((-4, 2), 1.5, 5, fc='#FFCDD2', ec='#C62828', lw=1, alpha=0.4))\n ax.add_patch(Rectangle((2.5, 2), 1.5, 5, fc='#FFCDD2', ec='#C62828', lw=1, alpha=0.4))\n ax.add_patch(Rectangle((-2, 3), 4, 4, fc='#FFE0B2', ec='#E65100', lw=1.5, alpha=0.6))\n ax.add_patch(Rectangle((-1.5, 7), 3, 1.5, fc='#C8E6C9', ec='#1B5E20', lw=1.5, alpha=0.6))\n ax.annotate('', xy=(-4.5, 0), xytext=(-3.5, 2), arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.5))\n ax.annotate('', xy=(4.5, 0), xytext=(3.5, 2), arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.5))\n elif '一二级' in title:\n ax.add_patch(Rectangle((-3, 1), 6, 3, fc='#FFE0B2', ec='#E65100', lw=1, alpha=0.4))\n ax.add_patch(Rectangle((-1.5, 4), 3, 3, fc='#C8E6C9', ec='#1B5E20', lw=1.5, alpha=0.6))\n ax.annotate('', xy=(-4, 0), xytext=(-3, 1), arrowprops=dict(arrowstyle='->', color='#E65100', lw=1.5))\n else:\n ax.add_patch(Rectangle((-1.5, 4), 3, 3, fc='#C8E6C9', ec='#1B5E20', lw=1.5, alpha=0.6))\n ax.add_patch(Circle((0, 8.5), 1, fc='#2196F3', ec='#0D47A1', lw=1.5, alpha=0.6))\n ax.text(0, 8.5, '载荷', ha='center', va='center', fontsize=7, color='white')\n\n ax.text(0, -1, desc, ha='center', va='center', fontsize=7, style='italic')\n\n plt.suptitle('飞行时序与分离方案', fontsize=14, fontweight='bold')\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_separation.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_separation.png'))\n plt.close(fig)\n print('[OK] fig_separation')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 15 — GNC Architecture\n# ═══════════════════════════════════════════════════════════════════\ndef fig_gnc():\n fig, axes = plt.subplots(1, 2, figsize=(14, 7))\n\n # Control loop block diagram\n ax = axes[0]\n ax.set_xlim(0, 10); ax.set_ylim(0, 10); ax.axis('off')\n ax.set_title('(a) GNC 控制回路', fontsize=11, fontweight='bold')\n\n blocks = [\n (1, 7, '制导\\n计算机', '#2196F3'),\n (4, 7, '导航\\n滤波器', '#4CAF50'),\n (7, 7, '控制\\n律', '#FF9800'),\n (7, 4, '执行\\n机构', '#EF5350'),\n (4, 4, '运载\\n火箭', '#9C27B0'),\n (1, 4, '传感器\\nIMU/GPS', '#00BCD4'),\n ]\n for x, y, txt, color in blocks:\n rect = FancyBboxPatch((x-0.8, y-0.5), 1.6, 1, boxstyle=\"round,pad=0.1\",\n fc=color, ec='black', lw=1.2, alpha=0.7)\n ax.add_patch(rect)\n ax.text(x, y, txt, ha='center', va='center', fontsize=7, color='white', fontweight='bold')\n\n # Arrows\n arrows = [(1.8,7),(3.2,7), (4.8,7),(6.2,7), (7,6.5),(7,4.5), (6.2,4),(4.8,4), (4,3.5),(4,2.5)]\n # Simplified flow arrows\n for i in range(0, len(arrows)-1, 2):\n ax.annotate('', xy=arrows[i], xytext=arrows[i+1] if i+1 < len(arrows) else arrows[i],\n arrowprops=dict(arrowstyle='->', color='#333', lw=1.5))\n\n # Feedback loop\n ax.annotate('', xy=(1, 4.5), xytext=(4, 4.5),\n arrowprops=dict(arrowstyle='->', color='#333', lw=1, ls='--'))\n ax.annotate('', xy=(1, 3.5), xytext=(1, 4),\n arrowprops=dict(arrowstyle='->', color='#333', lw=1, ls='--'))\n ax.annotate('', xy=(4, 3), xytext=(4, 3.5),\n arrowprops=dict(arrowstyle='->', color='#333', lw=1, ls='--'))\n ax.text(2.5, 3.2, '反馈', fontsize=7, ha='center', style='italic')\n\n # Attitude control response\n ax = axes[1]\n t = np.linspace(0, 10, 500)\n # Step response of attitude control\n wn = 2.0; zeta = 0.7\n s = np.exp(-zeta*wn*t) * (np.cos(wn*np.sqrt(1-zeta**2)*t) + zeta/np.sqrt(1-zeta**2)*np.sin(wn*np.sqrt(1-zeta**2)*t))\n response = 1 - s\n\n ax.plot(t, response, 'b-', lw=2, label='姿态响应')\n ax.axhline(1.0, color='r', ls='--', lw=1, alpha=0.5, label='目标姿态')\n ax.axhline(1.05, color='orange', ls=':', lw=1, alpha=0.5)\n ax.axhline(0.95, color='orange', ls=':', lw=1, alpha=0.5, label='±5% 容差')\n ax.fill_between(t, 0.95, 1.05, alpha=0.1, color='orange')\n ax.set_xlabel('时间 (s)'); ax.set_ylabel('归一化姿态')\n ax.set_title('(b) 姿态控制阶跃响应'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_gnc.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_gnc.png'))\n plt.close(fig)\n print('[OK] fig_gnc')\n\n# ═══════════════════════════════════════════════════════════════════\n# FIGURE 16 — Sensitivity Analysis\n# ═══════════════════════════════════════════════════════════════════\ndef fig_sensitivity():\n fig, axes = plt.subplots(1, 2, figsize=(14, 7))\n\n # Tornado chart — parameter sensitivity on payload\n ax = axes[0]\n params = ['Isp (助推级 +5s)', '结构分数 (-1%)', '阻力系数 (-10%)',\n '发动机T/W (+10%)', '推进剂密度 (+2%)', '转弯程序优化',\n '推进剂残留 (-0.5%)']\n delta_pl_pos = [420, 850, 180, 310, 260, 350, 190] # tons\n delta_pl_neg = [-380, -920, -200, -280, -240, -400, -210]\n\n y = np.arange(len(params))\n ax.barh(y, delta_pl_pos, 0.4, color='#4CAF50', alpha=0.85, label='正向变化')\n ax.barh(y, delta_pl_neg, 0.4, color='#EF5350', alpha=0.85, label='负向变化')\n ax.set_yticks(y); ax.set_yticklabels(params, fontsize=8)\n ax.set_xlabel('有效载荷变化 (t)'); ax.set_title('(a) 参数灵敏度 (龙卷风图)', fontsize=11, fontweight='bold')\n ax.axvline(0, color='black', lw=0.5); ax.legend(fontsize=8); ax.grid(True, alpha=0.3, axis='x')\n\n # Pareto front — payload vs structural mass\n ax = axes[1]\n struct_fracs = np.linspace(0.06, 0.14, 50)\n payloads = []\n for sf in struct_fracs:\n dry = sf * (S0_PROP + S1_PROP + S2_PROP)\n pl = M0 - (S0_PROP + dry*0.75) - (S1_PROP + dry*0.15) - (S2_PROP + dry*0.10)\n payloads.append(max(0, pl/1e3))\n\n ax.plot([x*100 for x in struct_fracs], payloads, 'b-', lw=2)\n ax.axvline(10.3, color='red', ls='--', lw=1.5, label='当前设计: 10.3%')\n ax.fill_between([x*100 for x in struct_fracs], 0, payloads, alpha=0.15, color='blue')\n ax.set_xlabel('结构质量分数 (%)'); ax.set_ylabel('有效载荷 (t)')\n ax.set_title('(b) 结构效率 Pareto 前沿'); ax.legend(fontsize=8); ax.grid(True, alpha=0.3)\n\n plt.tight_layout()\n fig.savefig(os.path.join(OUT, 'fig_sensitivity.pdf'))\n fig.savefig(os.path.join(OUT, 'fig_sensitivity.png'))\n plt.close(fig)\n print('[OK] fig_sensitivity')\n\n\n# ═══════════════════════════════════════════════════════════════════\n# MAIN\n# ═══════════════════════════════════════════════════════════════════\nif __name__ == '__main__':\n print(\"=== Leviathan-100 Simulation & Figure Generation ===\")\n fig_vehicle_config()\n hist = simulate_trajectory()\n fig_trajectory(hist)\n fig_mass_budget()\n fig_propulsion()\n fig_structural()\n fig_aerodynamics()\n fig_deltav()\n fig_comparison()\n fig_thermal()\n fig_risk()\n fig_cost()\n fig_launch_site()\n fig_montecarlo()\n fig_separation()\n fig_gnc()\n fig_sensitivity()\n print(\"\\n=== All figures generated ===\")\n",
"file": "C:\\Users\\20211\\.openclaw\\workspace-vllm-local\\mega_rocket\\simulate.py"
}