主题049：电磁发射系统仿真

引言

电磁发射技术是一种利用电磁力将物体加速到高速的新型发射技术。与传统的化学发射相比，电磁发射具有初速高、可控性好、无化学污染等优点，在军事、航天、科研等领域具有广阔的应用前景。

电磁发射技术的发展历程

电磁发射技术的概念最早可以追溯到19世纪末，但直到20世纪70年代，随着脉冲功率技术和超导技术的发展，电磁发射才真正进入实用化研究阶段。

重要里程碑：

• 1901年：挪威科学家克里斯蒂安·伯克兰首次提出电磁炮概念
• 1944年：德国科学家约阿希姆·汉斯勒设计了最早的电磁炮原型
• 1970年代：澳大利亚国立大学建成第一台实用化电磁轨道炮
• 1980年代：美国启动电磁炮军事应用研究计划
• 2000年代：美国海军电磁轨道炮项目取得重大突破
• 2010年代：电磁发射技术向多领域扩展应用

电磁发射技术的分类

根据加速原理的不同，电磁发射技术主要分为以下几类：

1. 电磁轨道炮（Railgun）：

   • 利用洛伦兹力沿导轨加速电枢
   • 结构简单，加速能力强
   • 主要用于高速动能武器

2. 电磁线圈炮（Coilgun）：

   • 利用线圈产生的磁场吸引或排斥弹丸
   • 无接触发射，磨损小
   • 适合发射精密载荷

3. 重接炮（Reconnection Gun）：

   • 利用磁场重联原理加速
   • 效率较高
   • 研究阶段

4. 感应线圈炮（Induction Coilgun）：

   • 利用感应电流产生推力
   • 无需电接触
   • 适合发射非导电材料
     
     
     
     
     
     
     

电磁发射的应用领域

军事应用：

• 舰载电磁轨道炮：射程超过200公里，速度可达7马赫
• 电磁弹射器：用于航空母舰舰载机发射
• 电磁迫击炮：提高射速和精度

航天应用：

• 电磁助推发射：将载荷加速到第一宇宙速度
• 电磁轨道发射系统：降低航天发射成本
• 月球/火星资源发射：利用电磁发射将资源送入轨道

民用应用：

• 高速碰撞实验：材料冲击性能研究
• 地质勘探：高速穿地弹
• 消防灭火：远距离灭火弹发射

本主题内容概述

本主题将系统介绍电磁发射系统的原理和设计方法，包括：

• 电磁轨道炮的工作原理和数学模型
• 等离子体电枢的物理特性
• 发射效率的影响因素和优化方法
• 电磁线圈炮的设计原理
• 使用Python进行电磁发射系统动态仿真

---

电磁发射技术概述

电磁发射的基本原理

电磁发射利用电磁力（洛伦兹力）对物体进行加速。根据电磁学理论，载流导体在磁场中受到的力为：

$$\mathbf{F}=I\mathbf{L}\times \mathbf{B}$$

其中：

• $I$ 是电流
• $\mathbf{L}$ 是导体长度矢量
• $\mathbf{B}$ 是磁感应强度

对于电磁轨道炮，电枢（armature）在两根平行导轨之间滑动，电流从一根导轨流入，经过电枢，从另一根导轨流出。电流产生的磁场与电流相互作用，产生沿导轨方向的推力。

电磁发射的优势

1. 高初速

化学发射受限于火药燃气膨胀速度，通常不超过2 km/s。而电磁发射理论上可以达到10 km/s以上，甚至接近轨道速度（7.9 km/s）。

2. 可控性好

电磁发射的推力可以通过调节电流精确控制，实现：

• 精确的速度控制
• 可编程的加速度曲线
• 重复性好

3. 无化学污染

电磁发射不需要推进剂，因此：

• 无燃烧产物
• 无烟雾和火光
• 发射装置清洁

4. 安全性高

• 无易燃易爆化学品
• 弹药储存安全
• 发射过程可控

5. 成本低（潜在）

虽然初始投资高，但：

• 弹药成本低（无推进剂）
• 维护成本低
• 可重复使用

电磁发射的技术挑战

1. 脉冲功率技术

电磁发射需要巨大的瞬时功率（吉瓦级），对电源系统提出极高要求：

• 大容量储能系统
• 高速开关技术
• 脉冲成形网络

2. 导轨和电枢的烧蚀

高速滑动接触产生：

• 电弧烧蚀
• 摩擦磨损
• 材料熔化

3. 热管理

能量转换过程中的损耗转化为热量：

• 导轨发热
• 电枢过热
• 需要有效冷却

4. 结构强度

巨大的电磁力和惯性力：

• 导轨需要高强度支撑
• 整体结构抗冲击
• 振动控制

---

电磁轨道炮原理

电磁轨道炮的结构

电磁轨道炮由以下主要部分组成：

1. 导轨（Rails）：

   • 两根平行的导电轨道
   • 通常采用铜或铝材料
   • 长度从几米到几十米

2. 电枢（Armature）：

   • 在导轨之间滑动
   • 传导电流
   • 推动弹丸

3. 弹丸（Projectile）：

   • 被加速的物体
   • 通常包含在电枢内或与电枢一体化

4. 电源系统：

   • 脉冲功率源
   • 储能装置（电容器、电感器、旋转机械等）
   • 开关系统

电磁轨道炮的工作原理

电流分布

当电源接通时，电流从一根导轨流入，经过电枢，从另一根导轨流出。电流在导轨和电枢中产生的磁场可以用毕奥-萨伐尔定律计算。

洛伦兹力产生

电枢中的电流与导轨电流产生的磁场相互作用，产生洛伦兹力：

$$F=\frac{1}{2}{L}^{\prime }{I}^2$$

其中：

• $L^{\prime }$ 是导轨的单位长度电感（H/m）
• $I$ 是电流（A）

这个力推动电枢沿导轨加速。

速度计算

根据能量守恒，电枢获得的动能为：

$$\frac{1}{2}m{v}^2=\int Fdx=\frac{1}{2}{L}^{\prime }{I}^{2}x$$

因此出口速度为：

$$v=\sqrt{\frac{L^{\prime }{I}^{2}x}{m}}$$

其中：

• $m$ 是电枢和弹丸的总质量
• $x$ 是导轨长度
• $v$ 是出口速度

导轨电感计算

对于两根半径为$a$、中心间距为$d$的平行圆导轨，单位长度电感为：

$$L^{\prime }=\frac{{\mu }_0}{\pi }\ln \left(\frac{d-a}{a}\right)$$

其中${\mu }_0=4\pi \times 10^{-7}$ H/m是真空磁导率。

实际应用中，为了提高电感梯度（从而增加推力），通常采用：

• 扁平导轨（增大等效间距）
• 增强型结构（添加铁磁材料）

电流波形的影响

电流波形对发射性能有重要影响：

恒定电流

最简单的电流波形，推力恒定：

$$F=\frac{1}{2}{L}^{\prime }{I}_0^2=\text{常数}$$

加速度恒定，速度线性增加：

$$v(t)=\frac{F}{m}t=\frac{L^{\prime }{I}_0^2}{2m}t$$

衰减电流

实际系统中，由于电源内阻和导轨电阻，电流通常会衰减。假设电流按指数衰减：

$$I(t)=I_0{e}^{-t/\tau }$$

推力相应衰减：

$$F(t)=\frac{1}{2}{L}^{\prime }{I}_0^2{e}^{-2t/\tau }$$

最优电流波形

为了实现最大效率或最小烧蚀，可以设计特定的电流波形。例如，为了保持恒定加速度：

$$F=ma=\text{常数}$$

需要：

$$I(t)=\sqrt{\frac{2ma}{L^{\prime }}}=\text{常数}$$

但这在实际中难以实现，因为随着速度增加，反电动势增大，维持恒定电流需要更高的电压。

---

等离子体电枢理论

电枢的类型

根据电枢的材料和工作原理，可以分为：

1. 固体电枢

• 材料：铜、铝、复合材料
• 优点：结构简单，接触稳定
• 缺点：高速时产生电弧，烧蚀严重
• 适用速度：< 2 km/s

2. 等离子体电枢

• 原理：利用电弧等离子体传导电流
• 优点：无机械接触，适合超高速
• 缺点：电弧不稳定，烧蚀导轨
• 适用速度：2-10 km/s

3. 混合电枢

• 结构：固体电枢+等离子体尾焰
• 优点：结合两者优势
• 缺点：结构复杂

等离子体物理基础

等离子体的定义

等离子体是物质的第四态，由电离的气体组成，包含：

• 自由电子
• 正离子
• 中性原子/分子

等离子体的基本参数

1. 电子密度 $n_e$：单位体积内的电子数
2. 电子温度 $T_e$：电子的热运动能量
3. 电离度 $\alpha$：电离原子占总原子的比例

萨哈方程

描述热平衡状态下电离度的方程：

$$\frac{n_e{n}_i}{n_n}=\frac{2Z_i}{Z_n}{\left(\frac{2\pi m_e{k}_{B}T}{h^2}\right)}^{3/2}{e}^{-E_i/k_{B}T}$$

其中：

• $n_e$、$n_i$、$n_n$ 分别是电子、离子、中性粒子的数密度
• $Z_i$、$Z_n$ 是配分函数
• $E_i$ 是电离能
• $k_B$ 是玻尔兹曼常数

等离子体电枢的导电机制

电弧的形成

当固体电枢与导轨之间的接触失效（由于高速或烧蚀），电流通过电弧传导：

1. 触发电弧：间隙击穿，产生初始电离
2. 电弧维持：焦耳热维持高温，持续电离
3. 电弧运动：洛伦兹力和气流推动电弧前进

电弧电压

电弧电压包括：

$$V_{arc}=V_{cathode}+V_{anode}+E_{column}\cdot l_{arc}$$

其中：

• $V_{cathode}$、$V_{anode}$ 是阴极和阳极压降（通常10-20V）
• $E_{column}$ 是弧柱电场强度（通常10-100 V/cm）
• $l_{arc}$ 是弧柱长度

电弧电阻

电弧的电阻率与温度和压力有关：

$${\rho }_{arc}=\frac{1}{{\sigma }_{arc}}$$

其中电导率${\sigma }_{arc}$可以用Spitzer公式估算：

$${\sigma }_{arc}=\frac{3}{16\sqrt{\pi }}\frac{(k_B{T}_e{)}^{3/2}}{m_e^{1/2}Z{e}^2\ln \Lambda }$$

其中$\ln \Lambda$是库仑对数。

等离子体与导轨的相互作用

烧蚀机制

等离子体电枢对导轨的烧蚀主要包括：

1. 热烧蚀：

   • 高温等离子体加热导轨表面
   • 材料熔化、汽化
   • 热流密度可达$10^9$ W/m²

2. 电弧烧蚀：

   • 电弧根部的高电流密度
   • 局部熔化、蒸发
   • 材料喷溅

3. 化学烧蚀：

   • 等离子体与导轨材料反应
   • 氧化、氮化

烧蚀的影响因素

• 电流密度：烧蚀率 ∝ $J^n$（n≈2-3）
• 电弧速度：速度越高，烧蚀越严重
• 导轨材料：铜、铝、复合材料
• 冷却条件：水冷、气冷

减少烧蚀的方法

1. 材料优化：

   • 使用耐高温材料（钨、钼）
   • 复合材料（铜-钨、铜-钼）
   • 表面涂层

2. 结构优化：

   • 增加导轨截面积（降低电流密度）
   • 优化导轨形状

3. 操作优化：

   • 控制电流波形
   • 预加速（降低初始电弧强度）

---

电磁轨道炮设计

导轨设计

导轨材料选择

导轨材料需要满足以下要求：

• 高导电率（减少电阻损耗）
• 高强度（承受电磁力）
• 耐高温（抵抗烧蚀）
• 良好的耐磨性

常用材料比较：

材料	导电率 (%IACS)	强度 (MPa)	熔点 (°C)	特点
铜	100	200-400	1083	导电性好，强度低
铝	61	100-300	660	轻量，强度低
铜合金	80-90	400-800	1000-1100	平衡性能
铜-钨	50-70	800-1200	>2000	耐高温

导轨几何设计

导轨截面形状影响：

• 电流分布
• 电感梯度
• 机械强度
• 散热能力

常见截面形状：

• 矩形：简单，易加工
• 圆形：电流分布均匀
• 楔形：增加电感梯度
• 增强型：添加铁磁材料提高电感

导轨长度和间距

导轨长度决定加速距离：

$$L_{rail}=\frac{m{v}_{exit}^2}{L^{\prime }{I}^2}$$

导轨间距影响：

• 电感梯度
• 机械稳定性
• 电枢设计

典型参数：

• 长度：3-10 m
• 间距：2-5 cm

电枢设计

固体电枢设计

固体电枢的关键参数：

1. 质量：影响加速度和出口速度
2. 形状：影响接触和稳定性
3. 材料：影响导电性和耐磨性

常见形状：

• C形：简单，易制造
• U形：增加接触面积
• 多指形：分散电流

等离子体电枢设计

等离子体电枢需要：

• 稳定的电弧
• 足够的导电截面
• 对导轨的最小烧蚀

设计考虑：

• 初始触发方式
• 电弧约束结构
• 气体/材料注入

电源系统设计

储能系统类型

1. 电容器储能：

   • 原理：电容器储存电能
   • 优点：技术成熟，响应快
   • 缺点：能量密度低
   • 应用：小型电磁炮

2. 电感储能：

   • 原理：电感储存磁能
   • 优点：能量密度高
   • 缺点：开关技术复杂
   • 应用：大型系统

3. 旋转机械储能：

   • 原理：飞轮或发电机储能
   • 优点：储能密度高，可重复
   • 缺点：机械复杂
   • 应用：连续发射

4. 混合型储能：

   • 组合多种储能方式
   • 优化性能

脉冲成形网络

脉冲成形网络（PFN）用于将储能系统的能量转换为适合电磁炮的电流脉冲。

基本类型：

• LC网络：产生振荡电流
• 传输线：产生方波脉冲
• 主动控制：实时调节电流

开关技术

电磁炮需要高速、大电流开关：

1. 机械开关：

   • 优点：简单可靠
   • 缺点：速度慢，寿命短

2. 爆炸丝开关：

   • 原理：金属丝爆炸断开
   • 优点：速度快
   • 缺点：一次性使用

3. 半导体开关：

   • 类型：IGBT、GTO、IGCT
   • 优点：可控，可重复
   • 缺点：成本高，需要保护

4. 超导开关：

   • 原理：超导态-正常态转变
   • 优点：无损耗，速度快
   • 缺点：需要低温

---

发射效率分析

能量转换过程

电磁发射系统的能量转换过程：

1. 储能阶段：

   • 电源向储能装置充电
   • 效率：${\eta }_{charge}$

2. 放电阶段：

   • 储能装置向导轨放电
   • 效率：${\eta }_{discharge}$

3. 加速阶段：

   • 电磁能转化为弹丸动能
   • 效率：${\eta }_{acceleration}$

总效率：

$${\eta }_{total}={\eta }_{charge}\times {\eta }_{discharge}\times {\eta }_{acceleration}$$

能量损耗机制

1. 电阻损耗

导轨和电枢的电阻产生焦耳热：

$$P_{resistive}=I^{2}R$$

总电阻损耗能量：

$$E_{resistive}=\int I^{2}Rdt$$

2. 电弧损耗

等离子体电枢的电弧电压产生损耗：

$$P_{arc}=V_{arc}I$$

3. 磁滞损耗

如果使用铁磁材料增强磁场，会产生磁滞损耗：

$$E_{hysteresis}=\oint HdB$$

4. 机械损耗

摩擦、振动等机械损耗：

$$E_{mechanical}=\int F_{friction}dx$$

效率优化方法

1. 提高电感梯度

电感梯度$L^{\prime }$直接影响推力：

$$F=\frac{1}{2}{L}^{\prime }{I}^2$$

提高$L^{\prime }$的方法：

• 优化导轨几何形状
• 使用铁磁材料增强
• 增加导轨间距（权衡机械强度）

2. 降低电阻

• 使用高导电率材料
• 增大导轨截面积
• 优化电流分布

3. 优化电流波形

根据速度要求设计最优电流波形：

• 初始阶段：大电流，快速加速
• 后期阶段：维持电流，克服反电动势

4. 减少电弧损耗

• 优化电枢设计
• 控制电弧电压
• 使用固体电枢（低速时）

效率计算实例

假设一个电磁轨道炮的参数：

• 导轨长度：$L=6$ m
• 电感梯度：$L^{\prime }=0.5$ μH/m
• 峰值电流：$I_{max}=3$ MA
• 弹丸质量：$m=10$ kg
• 出口速度：$v=2000$ m/s

动能输出：

$$E_{kinetic}=\frac{1}{2}m{v}^2=\frac{1}{2}\times 10\times (2000{)}^2=20\text{ MJ}$$

输入电能（假设电流恒定）：

$$E_{input}=\frac{1}{2}{L}^{\prime }{I}^{2}L=\frac{1}{2}\times 0.5\times 10^{-6}\times (3\times 10^6{)}^2\times 6=13.5\text{ MJ}$$

理论效率：

$${\eta }_{theoretical}=\frac{E_{kinetic}}{E_{input}}=\frac{20}{13.5}=148\%$$

这个结果表明，仅靠电感储能不足以达到目标速度，需要考虑：

• 电源持续供能
• 更高的电流
• 更长的导轨

实际系统中，考虑各种损耗，典型效率为10-30%。

---

电磁线圈炮原理

电磁线圈炮的结构

电磁线圈炮由一系列同轴线圈组成，弹丸（通常是铁磁材料或感应环）被线圈的磁场吸引或排斥而加速。

主要组成部分：

1. 驱动线圈：产生磁场
2. 弹丸：被加速的物体
3. 位置传感器：检测弹丸位置
4. 控制系统：控制线圈通电时序

电磁线圈炮的工作原理

吸引型线圈炮

当线圈通电时，产生磁场吸引铁磁弹丸：

$$F=\frac{1}{2}\frac{dL}{dx}{I}^2$$

其中$dL/dx$是线圈电感随弹丸位置的变化率。

排斥型线圈炮（感应型）

线圈通电时，在导电弹丸中感应出涡流：

$$\epsilon =-\frac{d\Phi }{dt}$$

感应电流产生的磁场与线圈磁场相互作用，产生排斥力：

$$F=\frac{{\mu }_0{I}_1{I}_{2}A}{2\pi d}$$

多级线圈炮

单级线圈的加速能力有限，通常采用多级加速：

1. 时序控制：

   • 弹丸到达线圈中心前通电（吸引）
   • 弹丸通过中心后断电（避免减速）

2. 同步方法：

   • 位置传感器触发
   • 预测控制
   • 闭环反馈

线圈炮与轨道炮的比较

特性	线圈炮	轨道炮
接触方式	无接触	滑动接触
适用速度	中低速(<2km/s)	高速(>2km/s)
效率	较高(30-50%)	较低(10-30%)
磨损	小	大
结构复杂度	高	低
成本	高	低

---

电磁发射系统仿真

以下Python程序实现了电磁轨道炮的动态仿真，包括：

1. 电磁力计算
2. 运动方程求解
3. 电路方程求解
4. 能量分析
5. 可视化

"""
电磁发射系统仿真
主题049：电磁发射系统仿真
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("电磁发射系统仿真程序")
print("主题049：电磁发射系统仿真")
print("=" * 70)

# ==================== 1. 物理常数和参数 ====================
mu_0 = 4 * np.pi * 1e-7  # 真空磁导率 (H/m)

print("\n【仿真参数设置】")

# ==================== 2. 电磁轨道炮参数 ====================
# 导轨参数
rail_length = 6.0  # 导轨长度 (m)
rail_spacing = 0.05  # 导轨间距 (m)
rail_radius = 0.02  # 导轨半径 (m)
rail_resistivity = 1.7e-8  # 铜电阻率 (Ω·m)

# 计算单位长度电感 (简化模型)
L_prime = mu_0 / np.pi * np.log((rail_spacing - rail_radius) / rail_radius)
print(f"\n【导轨参数】")
print(f"  导轨长度: {rail_length} m")
print(f"  导轨间距: {rail_spacing * 100:.1f} cm")
print(f"  单位长度电感 L': {L_prime * 1e6:.2f} μH/m")

# 电枢参数
armature_mass = 0.5  # 电枢质量 (kg)
projectile_mass = 2.0  # 弹丸质量 (kg)
total_mass = armature_mass + projectile_mass  # 总质量 (kg)
armature_resistance = 1e-4  # 电枢电阻 (Ω)

print(f"\n【电枢和弹丸参数】")
print(f"  电枢质量: {armature_mass} kg")
print(f"  弹丸质量: {projectile_mass} kg")
print(f"  总质量: {total_mass} kg")

# 电源参数
C_bank = 10e-3  # 电容组容量 (F)
V_charge = 10e3  # 充电电压 (V)
L_pulse = 50e-6  # 脉冲电感 (H)
R_circuit = 5e-3  # 电路电阻 (Ω)

print(f"\n【电源参数】")
print(f"  电容组容量: {C_bank * 1e3:.1f} mF")
print(f"  充电电压: {V_charge / 1e3:.1f} kV")
print(f"  储能: {0.5 * C_bank * V_charge**2 / 1e6:.2f} MJ")

# ==================== 3. 电磁轨道炮仿真 ====================
def railgun_simulation(rail_length, L_prime, total_mass, C, V0, L_pulse, R_total, dt=1e-5):
    """
    电磁轨道炮动态仿真
    
    参数:
    rail_length: 导轨长度 (m)
    L_prime: 单位长度电感 (H/m)
    total_mass: 总质量 (kg)
    C: 电容 (F)
    V0: 初始电压 (V)
    L_pulse: 脉冲电感 (H)
    R_total: 总电阻 (Ω)
    dt: 时间步长 (s)
    
    返回:
    t_array: 时间数组
    x_array: 位置数组
    v_array: 速度数组
    I_array: 电流数组
    F_array: 推力数组
    """
    
    # 初始化
    t = 0
    x = 0  # 位置
    v = 0  # 速度
    I = 0  # 电流
    V_cap = V0  # 电容电压
    
    # 存储数组
    t_list = [t]
    x_list = [x]
    v_list = [v]
    I_list = [I]
    F_list = [0]
    V_list = [V_cap]
    
    # 仿真循环
    max_steps = 1000000
    step = 0
    
    while x < rail_length and step < max_steps:
        # 计算电感
        L_rail = L_prime * x if x > 0 else L_prime * 0.001  # 避免除零
        L_total = L_pulse + L_rail
        
        # 计算反电动势
        back_emf = L_prime * I * v if I > 0 else 0
        
        # 电路方程: V_cap = L_total * dI/dt + R_total * I + back_emf
        # 简化的电流更新
        if x < rail_length:
            dI_dt = (V_cap - R_total * I - back_emf) / L_total
            I_new = I + dI_dt * dt
        else:
            I_new = 0
        
        # 电容放电
        dV_dt = -I / C
        V_cap = V_cap + dV_dt * dt
        
        # 计算电磁力
        F = 0.5 * L_prime * I**2
        
        # 运动方程
        a = F / total_mass
        v_new = v + a * dt
        x_new = x + v * dt
        
        # 更新
        t = t + dt
        I = I_new
        v = v_new
        x = x_new
        
        # 存储
        t_list.append(t)
        x_list.append(x)
        v_list.append(v)
        I_list.append(I)
        F_list.append(F)
        V_list.append(V_cap)
        
        step += 1
        
        # 检查电容电压
        if V_cap < 100:
            break
    
    return (np.array(t_list), np.array(x_list), np.array(v_list), 
            np.array(I_list), np.array(F_list), np.array(V_list))

# 运行仿真
print("\n【运行电磁轨道炮仿真...】")
t_rg, x_rg, v_rg, I_rg, F_rg, V_rg = railgun_simulation(
    rail_length, L_prime, total_mass, C_bank, V_charge, 
    L_pulse, R_circuit + armature_resistance
)

print(f"  仿真步数: {len(t_rg)}")
print(f"  仿真时间: {t_rg[-1] * 1e3:.2f} ms")
print(f"  出口位置: {x_rg[-1]:.3f} m")
print(f"  出口速度: {v_rg[-1]:.1f} m/s ({v_rg[-1]/343:.1f} 马赫)")
print(f"  出口电流: {I_rg[-1]/1e3:.1f} kA")

# 计算性能指标
kinetic_energy = 0.5 * total_mass * v_rg[-1]**2
initial_energy = 0.5 * C_bank * V_charge**2
efficiency = kinetic_energy / initial_energy * 100

print(f"\n【性能指标】")
print(f"  初始电能: {initial_energy / 1e6:.2f} MJ")
print(f"  出口动能: {kinetic_energy / 1e6:.2f} MJ")
print(f"  能量效率: {efficiency:.1f}%")
print(f"  平均推力: {np.mean(F_rg)/1e3:.1f} kN")
print(f"  峰值推力: {np.max(F_rg)/1e3:.1f} kN")

# ==================== 4. 等离子体电枢仿真 ====================
def plasma_armature_simulation(t_array, I_array, rail_spacing, dt=1e-5):
    """
    等离子体电枢仿真
    
    参数:
    t_array: 时间数组
    I_array: 电流数组
    rail_spacing: 导轨间距 (m)
    dt: 时间步长 (s)
    
    返回:
    T_plasma: 等离子体温度数组 (K)
    R_plasma: 等离子体电阻数组 (Ω)
    V_arc: 电弧电压数组 (V)
    """
    
    n_steps = len(t_array)
    T_plasma = np.zeros(n_steps)
    R_plasma = np.zeros(n_steps)
    V_arc = np.zeros(n_steps)
    
    # 初始条件
    T_plasma[0] = 5000  # 初始温度 (K)
    
    # 等离子体参数
    m_plasma = 1e-6  # 等离子体质量 (kg)
    c_v = 1000  # 比热容 (J/kg/K)
    
    for i in range(1, n_steps):
        I = I_array[i]
        
        # 等离子体电阻 (简化模型)
        # 电阻率随温度变化
        rho_plasma = 1e-4 * (10000 / T_plasma[i-1])**1.5  # Ω·m
        A_plasma = 1e-4  # 等离子体截面积 (m²)
        R_plasma[i] = rho_plasma * rail_spacing / A_plasma
        
        # 电弧电压
        V_arc[i] = I * R_plasma[i] + 20  # 加20V电极压降
        
        # 焦耳加热
        P_joule = I**2 * R_plasma[i]
        
        # 温度更新 (简化)
        dT = P_joule * dt / (m_plasma * c_v)
        T_plasma[i] = T_plasma[i-1] + dT
        
        # 辐射冷却 (简化)
        T_plasma[i] = T_plasma[i] - 100 * dt  # 冷却项
        
        # 限制温度范围
        T_plasma[i] = np.clip(T_plasma[i], 3000, 50000)
    
    return T_plasma, R_plasma, V_arc

print("\n【运行等离子体电枢仿真...】")
T_plasma, R_plasma, V_arc = plasma_armature_simulation(t_rg, I_rg, rail_spacing)

print(f"  等离子体温度范围: {np.min(T_plasma):.0f} K - {np.max(T_plasma):.0f} K")
print(f"  平均电弧电压: {np.mean(V_arc):.1f} V")
print(f"  电弧能量损耗: {np.trapz(V_arc * I_rg, t_rg) / 1e3:.1f} kJ")

# ==================== 5. 多级线圈炮仿真 ====================
def coilgun_simulation(n_stages, coil_length, coil_spacing, B_max, 
                       projectile_mass, projectile_radius, dt=1e-4):
    """
    多级线圈炮仿真
    
    参数:
    n_stages: 级数
    coil_length: 线圈长度 (m)
    coil_spacing: 线圈间距 (m)
    B_max: 峰值磁场 (T)
    projectile_mass: 弹丸质量 (kg)
    projectile_radius: 弹丸半径 (m)
    dt: 时间步长 (s)
    
    返回:
    t_array: 时间数组
    x_array: 位置数组
    v_array: 速度数组
    stage_active: 各级激活状态
    """
    
    # 线圈位置
    coil_positions = np.arange(n_stages) * (coil_length + coil_spacing) + coil_length / 2
    
    # 初始化
    t = 0
    x = 0
    v = 0
    
    t_list = [t]
    x_list = [x]
    v_list = [v]
    active_stages = [[]]
    
    max_steps = 500000
    step = 0
    
    while x < coil_positions[-1] + coil_length and step < max_steps:
        # 计算当前激活的线圈
        active = []
        F_total = 0
        
        for i, pos in enumerate(coil_positions):
            # 线圈激活条件: 弹丸接近线圈但未通过中心
            if pos - coil_length < x < pos:
                active.append(i)
                
                # 计算磁场力 (简化模型)
                # 力与磁场梯度和弹丸磁化强度成正比
                distance_to_center = pos - x
                B_field = B_max * np.exp(-distance_to_center**2 / (2 * (coil_length/4)**2))
                
                # 磁力 (简化)
                mu_eff = 4 * np.pi * 1e-7 * 1000  # 有效磁导率
                V_projectile = 4/3 * np.pi * projectile_radius**3
                M = B_field / mu_eff  # 磁化强度
                F_magnetic = mu_eff * M * np.gradient(B_field, distance_to_center if distance_to_center != 0 else 0.001) * V_projectile
                
                F_total += max(0, F_magnetic)  # 只保留吸引力
        
        # 运动方程
        a = F_total / projectile_mass
        v_new = v + a * dt
        x_new = x + v * dt
        
        # 更新
        t = t + dt
        v = v_new
        x = x_new
        
        # 存储
        t_list.append(t)
        x_list.append(x)
        v_list.append(v)
        active_stages.append(active.copy())
        
        step += 1
    
    return np.array(t_list), np.array(x_list), np.array(v_list), active_stages

print("\n【运行多级线圈炮仿真...】")
n_stages = 10
coil_length = 0.1  # 10 cm
coil_spacing = 0.05  # 5 cm
B_max = 10.0  # 10 T
projectile_radius = 0.015  # 1.5 cm

t_cg, x_cg, v_cg, stages_cg = coilgun_simulation(
    n_stages, coil_length, coil_spacing, B_max,
    projectile_mass, projectile_radius
)

print(f"  级数: {n_stages}")
print(f"  仿真时间: {t_cg[-1] * 1e3:.2f} ms")
print(f"  出口速度: {v_cg[-1]:.1f} m/s")
print(f"  出口动能: {0.5 * projectile_mass * v_cg[-1]**2 / 1e3:.1f} kJ")

# ==================== 6. 可视化 ====================
print("\n【生成可视化图表】")

# 图1: 电磁轨道炮性能
fig1, axes1 = plt.subplots(2, 2, figsize=(14, 10))

# 位置和速度
ax1 = axes1[0, 0]
ax1.plot(t_rg * 1e3, x_rg, 'b-', linewidth=2, label='位置')
ax1.set_xlabel('时间 (ms)', fontsize=11)
ax1.set_ylabel('位置 (m)', fontsize=11, color='b')
ax1.tick_params(axis='y', labelcolor='b')
ax1.set_title('弹丸位置和速度', fontsize=12)
ax1.grid(True, alpha=0.3)

ax1_twin = ax1.twinx()
ax1_twin.plot(t_rg * 1e3, v_rg, 'r--', linewidth=2, label='速度')
ax1_twin.set_ylabel('速度 (m/s)', fontsize=11, color='r')
ax1_twin.tick_params(axis='y', labelcolor='r')
ax1_twin.axhline(y=343, color='g', linestyle=':', alpha=0.5, label='音速')

# 电流和推力
ax2 = axes1[0, 1]
ax2.plot(t_rg * 1e3, I_rg / 1e3, 'b-', linewidth=2, label='电流')
ax2.set_xlabel('时间 (ms)', fontsize=11)
ax2.set_ylabel('电流 (kA)', fontsize=11, color='b')
ax2.tick_params(axis='y', labelcolor='b')
ax2.set_title('电流和推力', fontsize=12)
ax2.grid(True, alpha=0.3)

ax2_twin = ax2.twinx()
ax2_twin.plot(t_rg * 1e3, F_rg / 1e3, 'r--', linewidth=2, label='推力')
ax2_twin.set_ylabel('推力 (kN)', fontsize=11, color='r')
ax2_twin.tick_params(axis='y', labelcolor='r')

# 电容电压和能量
ax3 = axes1[1, 0]
ax3.plot(t_rg * 1e3, V_rg / 1e3, 'b-', linewidth=2)
ax3.set_xlabel('时间 (ms)', fontsize=11)
ax3.set_ylabel('电容电压 (kV)', fontsize=11)
ax3.set_title('电容放电曲线', fontsize=12)
ax3.grid(True, alpha=0.3)

# 能量转换
ax4 = axes1[1, 1]
E_cap = 0.5 * C_bank * V_rg**2 / 1e6  # MJ
E_kin = 0.5 * total_mass * v_rg**2 / 1e6  # MJ
E_loss = E_cap[0] - E_cap - E_kin  # MJ

ax4.plot(t_rg * 1e3, E_cap, 'b-', linewidth=2, label='电容储能')
ax4.plot(t_rg * 1e3, E_kin, 'r-', linewidth=2, label='弹丸动能')
ax4.plot(t_rg * 1e3, E_loss, 'g-', linewidth=2, label='能量损耗')
ax4.set_xlabel('时间 (ms)', fontsize=11)
ax4.set_ylabel('能量 (MJ)', fontsize=11)
ax4.set_title('能量转换过程', fontsize=12)
ax4.legend(loc='best')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('railgun_performance.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: railgun_performance.png")

# 图2: 等离子体电枢特性
fig2, axes2 = plt.subplots(2, 2, figsize=(14, 10))

# 等离子体温度
ax1 = axes2[0, 0]
ax1.plot(t_rg * 1e3, T_plasma / 1e3, 'r-', linewidth=2)
ax1.set_xlabel('时间 (ms)', fontsize=11)
ax1.set_ylabel('温度 (×1000 K)', fontsize=11)
ax1.set_title('等离子体温度', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.axhline(y=20, color='orange', linestyle='--', alpha=0.5, label='典型工作温度')
ax1.legend()

# 电弧电压
ax2 = axes2[0, 1]
ax2.plot(t_rg * 1e3, V_arc, 'b-', linewidth=2)
ax2.set_xlabel('时间 (ms)', fontsize=11)
ax2.set_ylabel('电弧电压 (V)', fontsize=11)
ax2.set_title('电弧电压', fontsize=12)
ax2.grid(True, alpha=0.3)

# 等离子体电阻
ax3 = axes2[1, 0]
ax3.semilogy(t_rg * 1e3, R_plasma * 1e3, 'g-', linewidth=2)
ax3.set_xlabel('时间 (ms)', fontsize=11)
ax3.set_ylabel('等离子体电阻 (mΩ)', fontsize=11)
ax3.set_title('等离子体电阻', fontsize=12)
ax3.grid(True, alpha=0.3)

# 电弧功率
P_arc = V_arc * I_rg / 1e6  # MW
ax4 = axes2[1, 1]
ax4.plot(t_rg * 1e3, P_arc, 'purple', linewidth=2)
ax4.fill_between(t_rg * 1e3, P_arc, alpha=0.3)
ax4.set_xlabel('时间 (ms)', fontsize=11)
ax4.set_ylabel('电弧功率 (MW)', fontsize=11)
ax4.set_title('电弧功率', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('plasma_armature.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: plasma_armature.png")

# 图3: 多级线圈炮性能
fig3, axes3 = plt.subplots(1, 2, figsize=(14, 5))

# 位置和速度
ax1 = axes3[0]
ax1.plot(t_cg * 1e3, x_cg, 'b-', linewidth=2, label='位置')
ax1.set_xlabel('时间 (ms)', fontsize=11)
ax1.set_ylabel('位置 (m)', fontsize=11, color='b')
ax1.tick_params(axis='y', labelcolor='b')
ax1.set_title('线圈炮弹丸运动', fontsize=12)
ax1.grid(True, alpha=0.3)

# 标记线圈位置
coil_positions = np.arange(n_stages) * (coil_length + coil_spacing) + coil_length / 2
for i, pos in enumerate(coil_positions):
    ax1.axvline(x=t_cg[np.argmin(np.abs(x_cg - pos))] * 1e3, 
                color='gray', linestyle='--', alpha=0.3)
    if i % 2 == 0:
        ax1.text(t_cg[np.argmin(np.abs(x_cg - pos))] * 1e3, 
                max(x_cg) * 0.9, f'{i+1}', fontsize=8, ha='center')

ax1_twin = ax1.twinx()
ax1_twin.plot(t_cg * 1e3, v_cg, 'r--', linewidth=2, label='速度')
ax1_twin.set_ylabel('速度 (m/s)', fontsize=11, color='r')
ax1_twin.tick_params(axis='y', labelcolor='r')

# 加速度
a_cg = np.gradient(v_cg, t_cg)
ax2 = axes3[1]
ax2.plot(t_cg * 1e3, a_cg / 9.81 / 1e3, 'g-', linewidth=2)
ax2.set_xlabel('时间 (ms)', fontsize=11)
ax2.set_ylabel('加速度 (×1000 g)', fontsize=11)
ax2.set_title('线圈炮加速度', fontsize=12)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('coilgun_performance.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: coilgun_performance.png")

# 图4: 电磁发射系统结构示意图
fig4, ax = plt.subplots(1, 1, figsize=(14, 6))

# 绘制导轨
rail_y = 0.5
ax.plot([0, rail_length * 100], [rail_y, rail_y], 'b-', linewidth=8, label='导轨')
ax.plot([0, rail_length * 100], [rail_y + rail_spacing * 100, rail_y + rail_spacing * 100], 
        'b-', linewidth=8)

# 绘制电枢
armature_x = rail_length * 100 * 0.6
armature = FancyBboxPatch((armature_x - 2, rail_y + 2), 4, rail_spacing * 100 - 4,
                          boxstyle="round,pad=0.5", 
                          fill=True, facecolor='orange', edgecolor='red', linewidth=2)
ax.add_patch(armature)

# 绘制弹丸
projectile = FancyBboxPatch((armature_x + 2, rail_y + rail_spacing * 100 / 2 - 1.5), 
                            8, 3,
                            boxstyle="round,pad=0.3", 
                            fill=True, facecolor='gray', edgecolor='black', linewidth=1)
ax.add_patch(projectile)

# 绘制电流方向
ax.annotate('', xy=(5, rail_y - 3), xytext=(5, rail_y - 8),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.text(5, rail_y - 10, '电流', ha='center', fontsize=10, color='red')

ax.annotate('', xy=(armature_x, rail_y + rail_spacing * 100 + 5), 
            xytext=(armature_x, rail_y + rail_spacing * 100 + 10),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))

ax.annotate('', xy=(rail_length * 100 - 5, rail_y + rail_spacing * 100 + 3), 
            xytext=(rail_length * 100 - 5, rail_y + rail_spacing * 100 + 8),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))

# 绘制磁场
for i in range(5):
    x_pos = 20 + i * rail_length * 100 / 5
    circle = Circle((x_pos, rail_y + rail_spacing * 100 / 2), 3, 
                    fill=False, color='green', linewidth=1.5)
    ax.add_patch(circle)
    ax.annotate('', xy=(x_pos + 2, rail_y + rail_spacing * 100 / 2), 
                xytext=(x_pos - 2, rail_y + rail_spacing * 100 / 2),
                arrowprops=dict(arrowstyle='->', color='green', lw=1.5))

ax.text(rail_length * 100 / 2, rail_y + rail_spacing * 100 + 15, 
        '磁场 ⊗ (进入纸面)', ha='center', fontsize=10, color='green')

# 绘制推力
ax.annotate('', xy=(armature_x + 15, rail_y + rail_spacing * 100 / 2), 
            xytext=(armature_x + 5, rail_y + rail_spacing * 100 / 2),
            arrowprops=dict(arrowstyle='->', color='purple', lw=3))
ax.text(armature_x + 20, rail_y + rail_spacing * 100 / 2, 
        '推力 F', fontsize=11, color='purple', va='center')

# 标注
ax.text(armature_x, rail_y + rail_spacing * 100 / 2, '电枢', 
        ha='center', va='center', fontsize=9, color='darkred')
ax.text(armature_x + 6, rail_y + rail_spacing * 100 / 2, '弹丸', 
        ha='center', va='center', fontsize=8, color='white')

ax.set_xlim(-10, rail_length * 100 + 20)
ax.set_ylim(-10, rail_y + rail_spacing * 100 + 25)
ax.set_aspect('equal')
ax.axis('off')
ax.set_title('电磁轨道炮结构示意图', fontsize=14)

plt.tight_layout()
plt.savefig('railgun_schematic.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: railgun_schematic.png")

# 图5: 效率分析
fig5, axes5 = plt.subplots(1, 2, figsize=(14, 5))

# 能量分配饼图
ax1 = axes5[0]
E_final_kin = kinetic_energy / 1e6  # MJ
E_resistive = np.trapz(I_rg**2 * (R_circuit + armature_resistance), t_rg) / 1e6  # MJ
E_remaining = 0.5 * C_bank * V_rg[-1]**2 / 1e6  # MJ
E_arc_total = np.trapz(V_arc * I_rg, t_rg) / 1e6  # MJ
E_other = initial_energy / 1e6 - E_final_kin - E_resistive - E_remaining - E_arc_total

energy_labels = ['弹丸动能', '电阻损耗', '电弧损耗', '剩余电能', '其他损耗']
energy_values = [max(0, E_final_kin), max(0, E_resistive), max(0, E_arc_total), 
                 max(0, E_remaining), max(0, E_other)]
colors = ['green', 'red', 'orange', 'blue', 'gray']

wedges, texts, autotexts = ax1.pie(energy_values, labels=energy_labels, autopct='%1.1f%%',
                                     colors=colors, startangle=90)
ax1.set_title(f'能量分配 (总能量: {initial_energy/1e6:.2f} MJ)', fontsize=12)

# 效率对比
ax2 = axes5[1]
systems = ['电磁轨道炮\n(本仿真)', '电磁轨道炮\n(典型)', '线圈炮\n(典型)', '化学炮\n(典型)']
efficiencies = [efficiency, 20, 40, 30]
colors_bar = ['blue', 'lightblue', 'orange', 'gray']

bars = ax2.bar(systems, efficiencies, color=colors_bar)
ax2.set_ylabel('效率 (%)', fontsize=11)
ax2.set_title('不同发射系统效率对比', fontsize=12)
ax2.grid(True, alpha=0.3, axis='y')
ax2.set_ylim(0, 50)

for bar, val in zip(bars, efficiencies):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
             f'{val:.1f}%', ha='center', fontsize=10)

plt.tight_layout()
plt.savefig('efficiency_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: efficiency_analysis.png")

# 图6: 参数敏感性分析
fig6, axes6 = plt.subplots(2, 2, figsize=(14, 10))

# 电流对速度的影响
ax1 = axes6[0, 0]
I_range = np.linspace(1e6, 5e6, 50)  # 1-5 MA
v_exit_current = []
for I_test in I_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        rail_length, L_prime, total_mass, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_current.append(v_test[-1])

ax1.plot(I_range / 1e6, v_exit_current, 'b-', linewidth=2)
ax1.set_xlabel('峰值电流 (MA)', fontsize=11)
ax1.set_ylabel('出口速度 (m/s)', fontsize=11)
ax1.set_title('电流对出口速度的影响', fontsize=12)
ax1.grid(True, alpha=0.3)

# 质量对速度的影响
ax2 = axes6[0, 1]
m_range = np.linspace(1, 10, 50)  # 1-10 kg
v_exit_mass = []
for m_test in m_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        rail_length, L_prime, m_test, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_mass.append(v_test[-1])

ax2.plot(m_range, v_exit_mass, 'r-', linewidth=2)
ax2.set_xlabel('弹丸质量 (kg)', fontsize=11)
ax2.set_ylabel('出口速度 (m/s)', fontsize=11)
ax2.set_title('弹丸质量对出口速度的影响', fontsize=12)
ax2.grid(True, alpha=0.3)

# 导轨长度对速度的影响
ax3 = axes6[1, 0]
L_range = np.linspace(2, 12, 50)  # 2-12 m
v_exit_length = []
for L_test in L_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        L_test, L_prime, total_mass, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_length.append(v_test[-1])

ax3.plot(L_range, v_exit_length, 'g-', linewidth=2)
ax3.set_xlabel('导轨长度 (m)', fontsize=11)
ax3.set_ylabel('出口速度 (m/s)', fontsize=11)
ax3.set_title('导轨长度对出口速度的影响', fontsize=12)
ax3.grid(True, alpha=0.3)

# 电感梯度对速度的影响
ax4 = axes6[1, 1]
Lprime_range = np.linspace(0.3e-6, 1.0e-6, 50)  # 0.3-1.0 μH/m
v_exit_Lprime = []
for Lp_test in Lprime_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        rail_length, Lp_test, total_mass, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_Lprime.append(v_test[-1])

ax4.plot(Lprime_range * 1e6, v_exit_Lprime, 'purple', linewidth=2)
ax4.set_xlabel('电感梯度 (μH/m)', fontsize=11)
ax4.set_ylabel('出口速度 (m/s)', fontsize=11)
ax4.set_title('电感梯度对出口速度的影响', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('parameter_sensitivity.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: parameter_sensitivity.png")

# ==================== 7. 生成GIF动画 ====================
print("\n【生成GIF动画】")

# 轨道炮发射动画
fig_anim, ax_anim = plt.subplots(1, 1, figsize=(12, 4))

# 准备数据
n_frames = 100
indices = np.linspace(0, len(t_rg)-1, n_frames, dtype=int)

def update_railgun(frame):
    ax_anim.clear()
    idx = indices[frame]
    
    # 绘制导轨
    ax_anim.plot([0, rail_length * 100], [0, 0], 'b-', linewidth=6)
    ax_anim.plot([0, rail_length * 100], [rail_spacing * 100, rail_spacing * 100], 'b-', linewidth=6)
    
    # 绘制电枢和弹丸
    x_pos = x_rg[idx] * 100
    armature = FancyBboxPatch((x_pos - 1, 2), 2, rail_spacing * 100 - 4,
                              boxstyle="round,pad=0.3", 
                              fill=True, facecolor='orange', edgecolor='red', linewidth=1.5)
    ax_anim.add_patch(armature)
    
    projectile = FancyBboxPatch((x_pos + 1, rail_spacing * 100 / 2 - 1), 
                                4, 2,
                                boxstyle="round,pad=0.2", 
                                fill=True, facecolor='gray', edgecolor='black', linewidth=1)
    ax_anim.add_patch(projectile)
    
    # 绘制电流（用箭头密度表示）
    current_normalized = I_rg[idx] / np.max(I_rg)
    n_arrows = int(5 + 10 * current_normalized)
    for i in range(n_arrows):
        x_arrow = 5 + i * (rail_length * 100 - 10) / n_arrows
        ax_anim.annotate('', xy=(x_arrow + 2, -3), xytext=(x_arrow, -3),
                        arrowprops=dict(arrowstyle='->', color='red', lw=1.5))
    
    # 标注
    ax_anim.set_xlim(-5, rail_length * 100 + 10)
    ax_anim.set_ylim(-8, rail_spacing * 100 + 8)
    ax_anim.set_aspect('equal')
    ax_anim.axis('off')
    ax_anim.set_title(f'电磁轨道炮发射过程\n时间: {t_rg[idx]*1e3:.2f} ms, '
                     f'速度: {v_rg[idx]:.0f} m/s, 电流: {I_rg[idx]/1e3:.0f} kA', 
                     fontsize=12)
    
    return []

anim_rg = FuncAnimation(fig_anim, update_railgun, frames=n_frames, interval=100, blit=False)
anim_rg.save('railgun_animation.gif', writer='pillow', fps=10, dpi=100)
plt.close()
print("  已保存: railgun_animation.gif")

# ==================== 8. 输出总结 ====================
print("\n" + "=" * 70)
print("电磁发射系统仿真完成总结")
print("=" * 70)

print("\n【电磁轨道炮性能】")
print(f"  导轨长度: {rail_length} m")
print(f"  电感梯度: {L_prime * 1e6:.2f} μH/m")
print(f"  总质量: {total_mass} kg")
print(f"  出口速度: {v_rg[-1]:.1f} m/s ({v_rg[-1]/343:.1f} 马赫)")
print(f"  出口动能: {kinetic_energy / 1e6:.2f} MJ")
print(f"  能量效率: {efficiency:.1f}%")

print("\n【等离子体电枢特性】")
print(f"  温度范围: {np.min(T_plasma)/1e3:.1f} - {np.max(T_plasma)/1e3:.1f} ×1000 K")
print(f"  平均电弧电压: {np.mean(V_arc):.1f} V")
print(f"  电弧能量损耗: {np.trapz(V_arc * I_rg, t_rg) / 1e3:.1f} kJ")

print("\n【多级线圈炮性能】")
print(f"  级数: {n_stages}")
print(f"  出口速度: {v_cg[-1]:.1f} m/s")
print(f"  出口动能: {0.5 * projectile_mass * v_cg[-1]**2 / 1e3:.1f} kJ")

print("\n【生成的文件】")
print("  1. railgun_performance.png - 轨道炮性能")
print("  2. plasma_armature.png - 等离子体电枢特性")
print("  3. coilgun_performance.png - 线圈炮性能")
print("  4. railgun_schematic.png - 轨道炮结构示意图")
print("  5. efficiency_analysis.png - 效率分析")
print("  6. parameter_sensitivity.png - 参数敏感性分析")
print("  7. railgun_animation.gif - 轨道炮发射动画")

print("\n" + "=" * 70)
print("仿真完成!")
print("=" * 70)

---

案例分析

案例1：美国海军电磁轨道炮项目

项目背景：
美国海军研究实验室（NRL）自2005年起开展电磁轨道炮研究，目标是开发一种射程超过200公里的舰载武器。

系统规格：

• 导轨长度：约10米
• 弹丸质量：约10-20 kg
• 出口速度：约2500 m/s（7马赫）
• 射程：超过200公里
• 储能：约30 MJ

技术挑战：

1. 导轨烧蚀：高速滑动产生严重烧蚀，限制发射次数
2. 脉冲功率：需要吉瓦级瞬时功率
3. 热管理：大量热量需要有效散发
4. 系统集成：将电磁炮集成到舰艇上

解决方案：

• 使用铜-钨复合材料导轨
• 主动冷却系统
• 先进的脉冲功率模块
• 模块化设计

项目进展：
2010年代进行了多次成功测试，但由于技术挑战和成本问题，项目优先级在2020年代被调整。

案例2：电磁弹射器（EMALS）

应用背景：
电磁弹射器用于航空母舰上弹射舰载机，取代传统的蒸汽弹射器。

系统特点：

• 加速距离：约100米
• 弹射质量：2-45吨（不同机型）
• 出口速度：约70-300 km/h
• 能量：约100-122 MJ

技术优势：

1. 可控性好：可以精确调节弹射力，适应不同机型
2. 可靠性高：无蒸汽泄漏，维护简单
3. 能量效率高：约90%（蒸汽弹射约5%）
4. 体积小：比蒸汽弹射节省空间

系统组成：

• 直线感应电机
• 储能系统（飞轮）
• 功率转换系统
• 控制系统

应用：
已部署在美国福特级航空母舰上。

案例3：电磁助推发射系统

概念：
利用电磁发射将航天器或载荷加速到高速，然后火箭发动机点火进入轨道。

系统设想：

• 轨道长度：数公里
• 出口速度：2-3 km/s
• 载荷质量：数吨
• 应用：小型卫星发射

优势：

• 降低火箭燃料消耗
• 提高有效载荷比例
• 可重复使用
• 发射成本低

挑战：

• 超高速下的空气阻力
• 巨大的加速度（数百g）
• 系统建设成本高

研究进展：
NASA和其他机构进行过多次研究，但尚未实现实用化。

---

附录：Python代码使用说明

运行环境要求

• Python 3.7+
• NumPy
• Matplotlib

安装依赖

pip install numpy matplotlib

运行程序

python electromagnetic_launch_simulation.py

输出文件

程序将生成以下文件：

1. railgun_performance.png - 轨道炮性能分析图
2. plasma_armature.png - 等离子体电枢特性图
3. coilgun_performance.png - 线圈炮性能图
4. railgun_schematic.png - 轨道炮结构示意图
5. efficiency_analysis.png - 效率分析图
6. parameter_sensitivity.png - 参数敏感性分析图
7. railgun_animation.gif - 轨道炮发射动画

参数修改

可以通过修改代码中的参数来研究不同配置：

• rail_length - 导轨长度
• L_prime - 电感梯度
• total_mass - 弹丸质量
• C_bank - 电容组容量
• V_charge - 充电电压

---

"""
电磁发射系统仿真
主题049：电磁发射系统仿真
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, FancyBboxPatch, Circle
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Agg')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 70)
print("电磁发射系统仿真程序")
print("主题049：电磁发射系统仿真")
print("=" * 70)

# ==================== 1. 物理常数和参数 ====================
mu_0 = 4 * np.pi * 1e-7  # 真空磁导率 (H/m)

print("\n【仿真参数设置】")

# ==================== 2. 电磁轨道炮参数 ====================
# 导轨参数
rail_length = 6.0  # 导轨长度 (m)
rail_spacing = 0.05  # 导轨间距 (m)
rail_radius = 0.02  # 导轨半径 (m)
rail_resistivity = 1.7e-8  # 铜电阻率 (Ω·m)

# 计算单位长度电感 (简化模型)
L_prime = mu_0 / np.pi * np.log((rail_spacing - rail_radius) / rail_radius)
print(f"\n【导轨参数】")
print(f"  导轨长度: {rail_length} m")
print(f"  导轨间距: {rail_spacing * 100:.1f} cm")
print(f"  单位长度电感 L': {L_prime * 1e6:.2f} μH/m")

# 电枢参数
armature_mass = 0.5  # 电枢质量 (kg)
projectile_mass = 2.0  # 弹丸质量 (kg)
total_mass = armature_mass + projectile_mass  # 总质量 (kg)
armature_resistance = 1e-4  # 电枢电阻 (Ω)

print(f"\n【电枢和弹丸参数】")
print(f"  电枢质量: {armature_mass} kg")
print(f"  弹丸质量: {projectile_mass} kg")
print(f"  总质量: {total_mass} kg")

# 电源参数
C_bank = 10e-3  # 电容组容量 (F)
V_charge = 10e3  # 充电电压 (V)
L_pulse = 50e-6  # 脉冲电感 (H)
R_circuit = 5e-3  # 电路电阻 (Ω)

print(f"\n【电源参数】")
print(f"  电容组容量: {C_bank * 1e3:.1f} mF")
print(f"  充电电压: {V_charge / 1e3:.1f} kV")
print(f"  储能: {0.5 * C_bank * V_charge**2 / 1e6:.2f} MJ")

# ==================== 3. 电磁轨道炮仿真 ====================
def railgun_simulation(rail_length, L_prime, total_mass, C, V0, L_pulse, R_total, dt=1e-5):
    """
    电磁轨道炮动态仿真
    """
    
    # 初始化
    t = 0
    x = 0  # 位置
    v = 0  # 速度
    I = 0  # 电流
    V_cap = V0  # 电容电压
    
    # 存储数组
    t_list = [t]
    x_list = [x]
    v_list = [v]
    I_list = [I]
    F_list = [0]
    V_list = [V_cap]
    
    # 仿真循环
    max_steps = 1000000
    step = 0
    
    while x < rail_length and step < max_steps:
        # 计算电感
        L_rail = L_prime * x if x > 0 else L_prime * 0.001  # 避免除零
        L_total = L_pulse + L_rail
        
        # 计算反电动势
        back_emf = L_prime * I * v if I > 0 else 0
        
        # 电路方程
        if x < rail_length:
            dI_dt = (V_cap - R_total * I - back_emf) / L_total
            I_new = I + dI_dt * dt
        else:
            I_new = 0
        
        # 电容放电
        dV_dt = -I / C
        V_cap = V_cap + dV_dt * dt
        
        # 计算电磁力
        F = 0.5 * L_prime * I**2
        
        # 运动方程
        a = F / total_mass
        v_new = v + a * dt
        x_new = x + v * dt
        
        # 更新
        t = t + dt
        I = I_new
        v = v_new
        x = x_new
        
        # 存储
        t_list.append(t)
        x_list.append(x)
        v_list.append(v)
        I_list.append(I)
        F_list.append(F)
        V_list.append(V_cap)
        
        step += 1
        
        # 检查电容电压
        if V_cap < 100:
            break
    
    return (np.array(t_list), np.array(x_list), np.array(v_list), 
            np.array(I_list), np.array(F_list), np.array(V_list))

# 运行仿真
print("\n【运行电磁轨道炮仿真...】")
t_rg, x_rg, v_rg, I_rg, F_rg, V_rg = railgun_simulation(
    rail_length, L_prime, total_mass, C_bank, V_charge, 
    L_pulse, R_circuit + armature_resistance
)

print(f"  仿真步数: {len(t_rg)}")
print(f"  仿真时间: {t_rg[-1] * 1e3:.2f} ms")
print(f"  出口位置: {x_rg[-1]:.3f} m")
print(f"  出口速度: {v_rg[-1]:.1f} m/s ({v_rg[-1]/343:.1f} 马赫)")
print(f"  出口电流: {I_rg[-1]/1e3:.1f} kA")

# 计算性能指标
kinetic_energy = 0.5 * total_mass * v_rg[-1]**2
initial_energy = 0.5 * C_bank * V_charge**2
efficiency = kinetic_energy / initial_energy * 100

print(f"\n【性能指标】")
print(f"  初始电能: {initial_energy / 1e6:.2f} MJ")
print(f"  出口动能: {kinetic_energy / 1e6:.2f} MJ")
print(f"  能量效率: {efficiency:.1f}%")
print(f"  平均推力: {np.mean(F_rg)/1e3:.1f} kN")
print(f"  峰值推力: {np.max(F_rg)/1e3:.1f} kN")

# ==================== 4. 等离子体电枢仿真 ====================
def plasma_armature_simulation(t_array, I_array, rail_spacing, dt=1e-5):
    """
    等离子体电枢仿真
    """
    
    n_steps = len(t_array)
    T_plasma = np.zeros(n_steps)
    R_plasma = np.zeros(n_steps)
    V_arc = np.zeros(n_steps)
    
    # 初始条件
    T_plasma[0] = 5000  # 初始温度 (K)
    
    # 等离子体参数
    m_plasma = 1e-6  # 等离子体质量 (kg)
    c_v = 1000  # 比热容 (J/kg/K)
    
    for i in range(1, n_steps):
        I = I_array[i]
        
        # 等离子体电阻 (简化模型)
        rho_plasma = 1e-4 * (10000 / T_plasma[i-1])**1.5  # Ω·m
        A_plasma = 1e-4  # 等离子体截面积 (m²)
        R_plasma[i] = rho_plasma * rail_spacing / A_plasma
        
        # 电弧电压
        V_arc[i] = I * R_plasma[i] + 20  # 加20V电极压降
        
        # 焦耳加热
        P_joule = I**2 * R_plasma[i]
        
        # 温度更新 (简化)
        dT = P_joule * dt / (m_plasma * c_v)
        T_plasma[i] = T_plasma[i-1] + dT
        
        # 辐射冷却 (简化)
        T_plasma[i] = T_plasma[i] - 100 * dt  # 冷却项
        
        # 限制温度范围
        T_plasma[i] = np.clip(T_plasma[i], 3000, 50000)
    
    return T_plasma, R_plasma, V_arc

print("\n【运行等离子体电枢仿真...】")
T_plasma, R_plasma, V_arc = plasma_armature_simulation(t_rg, I_rg, rail_spacing)

print(f"  等离子体温度范围: {np.min(T_plasma):.0f} K - {np.max(T_plasma):.0f} K")
print(f"  平均电弧电压: {np.mean(V_arc):.1f} V")
print(f"  电弧能量损耗: {np.trapezoid(V_arc * I_rg, t_rg) / 1e3:.1f} kJ")

# ==================== 5. 多级线圈炮仿真 ====================
def coilgun_simulation(n_stages, coil_length, coil_spacing, B_max, 
                       projectile_mass, projectile_radius, dt=1e-4):
    """
    多级线圈炮仿真
    """
    
    # 线圈位置
    coil_positions = np.arange(n_stages) * (coil_length + coil_spacing) + coil_length / 2
    
    # 初始化
    t = 0
    x = 0
    v = 0
    
    t_list = [t]
    x_list = [x]
    v_list = [v]
    active_stages = [[]]
    
    max_steps = 500000
    step = 0
    
    while x < coil_positions[-1] + coil_length and step < max_steps:
        # 计算当前激活的线圈
        active = []
        F_total = 0
        
        for i, pos in enumerate(coil_positions):
            # 线圈激活条件
            if pos - coil_length < x < pos:
                active.append(i)
                
                # 计算磁场力 (简化模型)
                distance_to_center = pos - x
                B_field = B_max * np.exp(-distance_to_center**2 / (2 * (coil_length/4)**2))
                
                # 磁力 (简化)
                mu_eff = 4 * np.pi * 1e-7 * 1000  # 有效磁导率
                V_projectile = 4/3 * np.pi * projectile_radius**3
                M = B_field / mu_eff  # 磁化强度
                grad_B = -B_field * distance_to_center / (coil_length/4)**2
                F_magnetic = mu_eff * M * grad_B * V_projectile
                
                F_total += max(0, F_magnetic)  # 只保留吸引力
        
        # 运动方程
        a = F_total / projectile_mass
        v_new = v + a * dt
        x_new = x + v * dt
        
        # 更新
        t = t + dt
        v = v_new
        x = x_new
        
        # 存储
        t_list.append(t)
        x_list.append(x)
        v_list.append(v)
        active_stages.append(active.copy())
        
        step += 1
    
    return np.array(t_list), np.array(x_list), np.array(v_list), active_stages

print("\n【运行多级线圈炮仿真...】")
n_stages = 10
coil_length = 0.1  # 10 cm
coil_spacing = 0.05  # 5 cm
B_max = 10.0  # 10 T
projectile_radius = 0.015  # 1.5 cm

t_cg, x_cg, v_cg, stages_cg = coilgun_simulation(
    n_stages, coil_length, coil_spacing, B_max,
    projectile_mass, projectile_radius
)

print(f"  级数: {n_stages}")
print(f"  仿真时间: {t_cg[-1] * 1e3:.2f} ms")
print(f"  出口速度: {v_cg[-1]:.1f} m/s")
print(f"  出口动能: {0.5 * projectile_mass * v_cg[-1]**2 / 1e3:.1f} kJ")

# ==================== 6. 可视化 ====================
print("\n【生成可视化图表】")

# 图1: 电磁轨道炮性能
fig1, axes1 = plt.subplots(2, 2, figsize=(14, 10))

# 位置和速度
ax1 = axes1[0, 0]
ax1.plot(t_rg * 1e3, x_rg, 'b-', linewidth=2, label='位置')
ax1.set_xlabel('时间 (ms)', fontsize=11)
ax1.set_ylabel('位置 (m)', fontsize=11, color='b')
ax1.tick_params(axis='y', labelcolor='b')
ax1.set_title('弹丸位置和速度', fontsize=12)
ax1.grid(True, alpha=0.3)

ax1_twin = ax1.twinx()
ax1_twin.plot(t_rg * 1e3, v_rg, 'r--', linewidth=2, label='速度')
ax1_twin.set_ylabel('速度 (m/s)', fontsize=11, color='r')
ax1_twin.tick_params(axis='y', labelcolor='r')
ax1_twin.axhline(y=343, color='g', linestyle=':', alpha=0.5, label='音速')

# 电流和推力
ax2 = axes1[0, 1]
ax2.plot(t_rg * 1e3, I_rg / 1e3, 'b-', linewidth=2, label='电流')
ax2.set_xlabel('时间 (ms)', fontsize=11)
ax2.set_ylabel('电流 (kA)', fontsize=11, color='b')
ax2.tick_params(axis='y', labelcolor='b')
ax2.set_title('电流和推力', fontsize=12)
ax2.grid(True, alpha=0.3)

ax2_twin = ax2.twinx()
ax2_twin.plot(t_rg * 1e3, F_rg / 1e3, 'r--', linewidth=2, label='推力')
ax2_twin.set_ylabel('推力 (kN)', fontsize=11, color='r')
ax2_twin.tick_params(axis='y', labelcolor='r')

# 电容电压和能量
ax3 = axes1[1, 0]
ax3.plot(t_rg * 1e3, V_rg / 1e3, 'b-', linewidth=2)
ax3.set_xlabel('时间 (ms)', fontsize=11)
ax3.set_ylabel('电容电压 (kV)', fontsize=11)
ax3.set_title('电容放电曲线', fontsize=12)
ax3.grid(True, alpha=0.3)

# 能量转换
ax4 = axes1[1, 1]
E_cap = 0.5 * C_bank * V_rg**2 / 1e6  # MJ
E_kin = 0.5 * total_mass * v_rg**2 / 1e6  # MJ
E_loss = E_cap[0] - E_cap - E_kin  # MJ

ax4.plot(t_rg * 1e3, E_cap, 'b-', linewidth=2, label='电容储能')
ax4.plot(t_rg * 1e3, E_kin, 'r-', linewidth=2, label='弹丸动能')
ax4.plot(t_rg * 1e3, E_loss, 'g-', linewidth=2, label='能量损耗')
ax4.set_xlabel('时间 (ms)', fontsize=11)
ax4.set_ylabel('能量 (MJ)', fontsize=11)
ax4.set_title('能量转换过程', fontsize=12)
ax4.legend(loc='best')
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('railgun_performance.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: railgun_performance.png")

# 图2: 等离子体电枢特性
fig2, axes2 = plt.subplots(2, 2, figsize=(14, 10))

# 等离子体温度
ax1 = axes2[0, 0]
ax1.plot(t_rg * 1e3, T_plasma / 1e3, 'r-', linewidth=2)
ax1.set_xlabel('时间 (ms)', fontsize=11)
ax1.set_ylabel('温度 (×1000 K)', fontsize=11)
ax1.set_title('等离子体温度', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.axhline(y=20, color='orange', linestyle='--', alpha=0.5, label='典型工作温度')
ax1.legend()

# 电弧电压
ax2 = axes2[0, 1]
ax2.plot(t_rg * 1e3, V_arc, 'b-', linewidth=2)
ax2.set_xlabel('时间 (ms)', fontsize=11)
ax2.set_ylabel('电弧电压 (V)', fontsize=11)
ax2.set_title('电弧电压', fontsize=12)
ax2.grid(True, alpha=0.3)

# 等离子体电阻
ax3 = axes2[1, 0]
ax3.semilogy(t_rg * 1e3, R_plasma * 1e3, 'g-', linewidth=2)
ax3.set_xlabel('时间 (ms)', fontsize=11)
ax3.set_ylabel('等离子体电阻 (mΩ)', fontsize=11)
ax3.set_title('等离子体电阻', fontsize=12)
ax3.grid(True, alpha=0.3)

# 电弧功率
P_arc = V_arc * I_rg / 1e6  # MW
ax4 = axes2[1, 1]
ax4.plot(t_rg * 1e3, P_arc, 'purple', linewidth=2)
ax4.fill_between(t_rg * 1e3, P_arc, alpha=0.3)
ax4.set_xlabel('时间 (ms)', fontsize=11)
ax4.set_ylabel('电弧功率 (MW)', fontsize=11)
ax4.set_title('电弧功率', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('plasma_armature.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: plasma_armature.png")

# 图3: 多级线圈炮性能
fig3, axes3 = plt.subplots(1, 2, figsize=(14, 5))

# 位置和速度
ax1 = axes3[0]
ax1.plot(t_cg * 1e3, x_cg, 'b-', linewidth=2, label='位置')
ax1.set_xlabel('时间 (ms)', fontsize=11)
ax1.set_ylabel('位置 (m)', fontsize=11, color='b')
ax1.tick_params(axis='y', labelcolor='b')
ax1.set_title('线圈炮弹丸运动', fontsize=12)
ax1.grid(True, alpha=0.3)

# 标记线圈位置
coil_positions = np.arange(n_stages) * (coil_length + coil_spacing) + coil_length / 2
for i, pos in enumerate(coil_positions):
    ax1.axvline(x=t_cg[np.argmin(np.abs(x_cg - pos))] * 1e3, 
                color='gray', linestyle='--', alpha=0.3)
    if i % 2 == 0:
        ax1.text(t_cg[np.argmin(np.abs(x_cg - pos))] * 1e3, 
                max(x_cg) * 0.9, f'{i+1}', fontsize=8, ha='center')

ax1_twin = ax1.twinx()
ax1_twin.plot(t_cg * 1e3, v_cg, 'r--', linewidth=2, label='速度')
ax1_twin.set_ylabel('速度 (m/s)', fontsize=11, color='r')
ax1_twin.tick_params(axis='y', labelcolor='r')

# 加速度
a_cg = np.gradient(v_cg, t_cg)
ax2 = axes3[1]
ax2.plot(t_cg * 1e3, a_cg / 9.81 / 1e3, 'g-', linewidth=2)
ax2.set_xlabel('时间 (ms)', fontsize=11)
ax2.set_ylabel('加速度 (×1000 g)', fontsize=11)
ax2.set_title('线圈炮加速度', fontsize=12)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('coilgun_performance.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: coilgun_performance.png")

# 图4: 电磁发射系统结构示意图
fig4, ax = plt.subplots(1, 1, figsize=(14, 6))

# 绘制导轨
rail_y = 0.5
ax.plot([0, rail_length * 100], [rail_y, rail_y], 'b-', linewidth=8, label='导轨')
ax.plot([0, rail_length * 100], [rail_y + rail_spacing * 100, rail_y + rail_spacing * 100], 
        'b-', linewidth=8)

# 绘制电枢
armature_x = rail_length * 100 * 0.6
armature = FancyBboxPatch((armature_x - 2, rail_y + 2), 4, rail_spacing * 100 - 4,
                          boxstyle="round,pad=0.5", 
                          fill=True, facecolor='orange', edgecolor='red', linewidth=2)
ax.add_patch(armature)

# 绘制弹丸
projectile = FancyBboxPatch((armature_x + 2, rail_y + rail_spacing * 100 / 2 - 1.5), 
                            8, 3,
                            boxstyle="round,pad=0.3", 
                            fill=True, facecolor='gray', edgecolor='black', linewidth=1)
ax.add_patch(projectile)

# 绘制电流方向
ax.annotate('', xy=(5, rail_y - 3), xytext=(5, rail_y - 8),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.text(5, rail_y - 10, '电流', ha='center', fontsize=10, color='red')

ax.annotate('', xy=(armature_x, rail_y + rail_spacing * 100 + 5), 
            xytext=(armature_x, rail_y + rail_spacing * 100 + 10),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))

ax.annotate('', xy=(rail_length * 100 - 5, rail_y + rail_spacing * 100 + 3), 
            xytext=(rail_length * 100 - 5, rail_y + rail_spacing * 100 + 8),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))

# 绘制磁场
for i in range(5):
    x_pos = 20 + i * rail_length * 100 / 5
    circle = Circle((x_pos, rail_y + rail_spacing * 100 / 2), 3, 
                    fill=False, color='green', linewidth=1.5)
    ax.add_patch(circle)
    ax.annotate('', xy=(x_pos + 2, rail_y + rail_spacing * 100 / 2), 
                xytext=(x_pos - 2, rail_y + rail_spacing * 100 / 2),
                arrowprops=dict(arrowstyle='->', color='green', lw=1.5))

ax.text(rail_length * 100 / 2, rail_y + rail_spacing * 100 + 15, 
        '磁场 ⊗ (进入纸面)', ha='center', fontsize=10, color='green')

# 绘制推力
ax.annotate('', xy=(armature_x + 15, rail_y + rail_spacing * 100 / 2), 
            xytext=(armature_x + 5, rail_y + rail_spacing * 100 / 2),
            arrowprops=dict(arrowstyle='->', color='purple', lw=3))
ax.text(armature_x + 20, rail_y + rail_spacing * 100 / 2, 
        '推力 F', fontsize=11, color='purple', va='center')

# 标注
ax.text(armature_x, rail_y + rail_spacing * 100 / 2, '电枢', 
        ha='center', va='center', fontsize=9, color='darkred')
ax.text(armature_x + 6, rail_y + rail_spacing * 100 / 2, '弹丸', 
        ha='center', va='center', fontsize=8, color='white')

ax.set_xlim(-10, rail_length * 100 + 20)
ax.set_ylim(-10, rail_y + rail_spacing * 100 + 25)
ax.set_aspect('equal')
ax.axis('off')
ax.set_title('电磁轨道炮结构示意图', fontsize=14)

plt.tight_layout()
plt.savefig('railgun_schematic.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: railgun_schematic.png")

# 图5: 效率分析
fig5, axes5 = plt.subplots(1, 2, figsize=(14, 5))

# 能量分配饼图
ax1 = axes5[0]
E_final_kin = kinetic_energy / 1e6  # MJ
E_resistive = np.trapezoid(I_rg**2 * (R_circuit + armature_resistance), t_rg) / 1e6  # MJ
E_remaining = 0.5 * C_bank * V_rg[-1]**2 / 1e6  # MJ
E_arc_total = np.trapezoid(V_arc * I_rg, t_rg) / 1e6  # MJ
E_other = initial_energy / 1e6 - E_final_kin - E_resistive - E_remaining - E_arc_total

energy_labels = ['弹丸动能', '电阻损耗', '电弧损耗', '剩余电能', '其他损耗']
energy_values = [max(0, E_final_kin), max(0, E_resistive), max(0, E_arc_total), 
                 max(0, E_remaining), max(0, E_other)]
colors = ['green', 'red', 'orange', 'blue', 'gray']

wedges, texts, autotexts = ax1.pie(energy_values, labels=energy_labels, autopct='%1.1f%%',
                                     colors=colors, startangle=90)
ax1.set_title(f'能量分配 (总能量: {initial_energy/1e6:.2f} MJ)', fontsize=12)

# 效率对比
ax2 = axes5[1]
systems = ['电磁轨道炮\n(本仿真)', '电磁轨道炮\n(典型)', '线圈炮\n(典型)', '化学炮\n(典型)']
efficiencies = [efficiency, 20, 40, 30]
colors_bar = ['blue', 'lightblue', 'orange', 'gray']

bars = ax2.bar(systems, efficiencies, color=colors_bar)
ax2.set_ylabel('效率 (%)', fontsize=11)
ax2.set_title('不同发射系统效率对比', fontsize=12)
ax2.grid(True, alpha=0.3, axis='y')
ax2.set_ylim(0, 50)

for bar, val in zip(bars, efficiencies):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 1, 
             f'{val:.1f}%', ha='center', fontsize=10)

plt.tight_layout()
plt.savefig('efficiency_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: efficiency_analysis.png")

# 图6: 参数敏感性分析
fig6, axes6 = plt.subplots(2, 2, figsize=(14, 10))

# 电流对速度的影响
ax1 = axes6[0, 0]
I_range = np.linspace(1e6, 5e6, 50)  # 1-5 MA
v_exit_current = []
for I_test in I_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        rail_length, L_prime, total_mass, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_current.append(v_test[-1])

ax1.plot(I_range / 1e6, v_exit_current, 'b-', linewidth=2)
ax1.set_xlabel('峰值电流 (MA)', fontsize=11)
ax1.set_ylabel('出口速度 (m/s)', fontsize=11)
ax1.set_title('电流对出口速度的影响', fontsize=12)
ax1.grid(True, alpha=0.3)

# 质量对速度的影响
ax2 = axes6[0, 1]
m_range = np.linspace(1, 10, 50)  # 1-10 kg
v_exit_mass = []
for m_test in m_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        rail_length, L_prime, m_test, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_mass.append(v_test[-1])

ax2.plot(m_range, v_exit_mass, 'r-', linewidth=2)
ax2.set_xlabel('弹丸质量 (kg)', fontsize=11)
ax2.set_ylabel('出口速度 (m/s)', fontsize=11)
ax2.set_title('弹丸质量对出口速度的影响', fontsize=12)
ax2.grid(True, alpha=0.3)

# 导轨长度对速度的影响
ax3 = axes6[1, 0]
L_range = np.linspace(2, 12, 50)  # 2-12 m
v_exit_length = []
for L_test in L_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        L_test, L_prime, total_mass, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_length.append(v_test[-1])

ax3.plot(L_range, v_exit_length, 'g-', linewidth=2)
ax3.set_xlabel('导轨长度 (m)', fontsize=11)
ax3.set_ylabel('出口速度 (m/s)', fontsize=11)
ax3.set_title('导轨长度对出口速度的影响', fontsize=12)
ax3.grid(True, alpha=0.3)

# 电感梯度对速度的影响
ax4 = axes6[1, 1]
Lprime_range = np.linspace(0.3e-6, 1.0e-6, 50)  # 0.3-1.0 μH/m
v_exit_Lprime = []
for Lp_test in Lprime_range:
    t_test, x_test, v_test, I_test_arr, F_test, V_test = railgun_simulation(
        rail_length, Lp_test, total_mass, C_bank, V_charge, 
        L_pulse, R_circuit + armature_resistance
    )
    v_exit_Lprime.append(v_test[-1])

ax4.plot(Lprime_range * 1e6, v_exit_Lprime, 'purple', linewidth=2)
ax4.set_xlabel('电感梯度 (μH/m)', fontsize=11)
ax4.set_ylabel('出口速度 (m/s)', fontsize=11)
ax4.set_title('电感梯度对出口速度的影响', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('parameter_sensitivity.png', dpi=150, bbox_inches='tight')
plt.close()
print("  已保存: parameter_sensitivity.png")

# ==================== 7. 生成GIF动画 ====================
print("\n【生成GIF动画】")

# 轨道炮发射动画
fig_anim, ax_anim = plt.subplots(1, 1, figsize=(12, 4))

# 准备数据
n_frames = 100
indices = np.linspace(0, len(t_rg)-1, n_frames, dtype=int)

def update_railgun(frame):
    ax_anim.clear()
    idx = indices[frame]
    
    # 绘制导轨
    ax_anim.plot([0, rail_length * 100], [0, 0], 'b-', linewidth=6)
    ax_anim.plot([0, rail_length * 100], [rail_spacing * 100, rail_spacing * 100], 'b-', linewidth=6)
    
    # 绘制电枢和弹丸
    x_pos = x_rg[idx] * 100
    armature = FancyBboxPatch((x_pos - 1, 2), 2, rail_spacing * 100 - 4,
                              boxstyle="round,pad=0.3", 
                              fill=True, facecolor='orange', edgecolor='red', linewidth=1.5)
    ax_anim.add_patch(armature)
    
    projectile = FancyBboxPatch((x_pos + 1, rail_spacing * 100 / 2 - 1), 
                                4, 2,
                                boxstyle="round,pad=0.2", 
                                fill=True, facecolor='gray', edgecolor='black', linewidth=1)
    ax_anim.add_patch(projectile)
    
    # 绘制电流（用箭头密度表示）
    current_normalized = I_rg[idx] / np.max(I_rg)
    n_arrows = int(5 + 10 * current_normalized)
    for i in range(n_arrows):
        x_arrow = 5 + i * (rail_length * 100 - 10) / n_arrows
        ax_anim.annotate('', xy=(x_arrow + 2, -3), xytext=(x_arrow, -3),
                        arrowprops=dict(arrowstyle='->', color='red', lw=1.5))
    
    # 标注
    ax_anim.set_xlim(-5, rail_length * 100 + 10)
    ax_anim.set_ylim(-8, rail_spacing * 100 + 8)
    ax_anim.set_aspect('equal')
    ax_anim.axis('off')
    ax_anim.set_title(f'电磁轨道炮发射过程\n时间: {t_rg[idx]*1e3:.2f} ms, '
                     f'速度: {v_rg[idx]:.0f} m/s, 电流: {I_rg[idx]/1e3:.0f} kA', 
                     fontsize=12)
    
    return []

anim_rg = FuncAnimation(fig_anim, update_railgun, frames=n_frames, interval=100, blit=False)
anim_rg.save('railgun_animation.gif', writer='pillow', fps=10, dpi=100)
plt.close()
print("  已保存: railgun_animation.gif")

# ==================== 8. 输出总结 ====================
print("\n" + "=" * 70)
print("电磁发射系统仿真完成总结")
print("=" * 70)

print("\n【电磁轨道炮性能】")
print(f"  导轨长度: {rail_length} m")
print(f"  电感梯度: {L_prime * 1e6:.2f} μH/m")
print(f"  总质量: {total_mass} kg")
print(f"  出口速度: {v_rg[-1]:.1f} m/s ({v_rg[-1]/343:.1f} 马赫)")
print(f"  出口动能: {kinetic_energy / 1e6:.2f} MJ")
print(f"  能量效率: {efficiency:.1f}%")

print("\n【等离子体电枢特性】")
print(f"  温度范围: {np.min(T_plasma)/1e3:.1f} - {np.max(T_plasma)/1e3:.1f} ×1000 K")
print(f"  平均电弧电压: {np.mean(V_arc):.1f} V")
print(f"  电弧能量损耗: {np.trapezoid(V_arc * I_rg, t_rg) / 1e3:.1f} kJ")

print("\n【多级线圈炮性能】")
print(f"  级数: {n_stages}")
print(f"  出口速度: {v_cg[-1]:.1f} m/s")
print(f"  出口动能: {0.5 * projectile_mass * v_cg[-1]**2 / 1e3:.1f} kJ")

print("\n【生成的文件】")
print("  1. railgun_performance.png - 轨道炮性能")
print("  2. plasma_armature.png - 等离子体电枢特性")
print("  3. coilgun_performance.png - 线圈炮性能")
print("  4. railgun_schematic.png - 轨道炮结构示意图")
print("  5. efficiency_analysis.png - 效率分析")
print("  6. parameter_sensitivity.png - 参数敏感性分析")
print("  7. railgun_animation.gif - 轨道炮发射动画")

print("\n" + "=" * 70)
print("仿真完成!")
print("=" * 70)