主题044：无线充电系统磁场

引言

无线充电技术（Wireless Power Transfer, WPT）作为现代电力电子领域的重要突破，正在深刻改变我们的生活方式。从智能手机的便捷充电到电动汽车的无线充电，从医疗植入设备的供电到工业机器人的能量传输，无线充电技术展现出广阔的应用前景。

无线充电的核心原理是利用电磁感应或磁共振实现电能的无线传输。当两个线圈靠近时，发射线圈产生的交变磁场会在接收线圈中感应出电动势，从而实现能量传输。这种技术的效率取决于线圈的设计、耦合系数、工作频率以及补偿网络的选择。

本教程将系统讲解无线充电系统的磁场建模与仿真方法，包括：

• 磁耦合基本原理与互感计算
• 线圈设计与磁场分布分析
• 补偿网络拓扑与谐振条件
• 传输效率优化与参数分析
• 偏移容忍性与多线圈系统

通过Python数值仿真，我们将深入理解无线充电系统的工作特性，掌握从理论计算到工程设计的完整流程。

---

第一部分：无线充电基础理论

1.1 电磁感应原理

无线充电技术基于法拉第电磁感应定律。当穿过闭合回路的磁通量发生变化时，回路中会产生感应电动势：

$$\mathcal{E}=-\frac{d\Phi }{dt}$$

其中：

• $\mathcal{E}$：感应电动势 [V]
• $\Phi$：磁通量 [Wb]
• $t$：时间 [s]

1.1.1 互感现象

当两个线圈靠近时，一个线圈中的电流产生的磁场会穿过另一个线圈，形成互感耦合。设线圈1中的电流为$I_1$，在线圈2中产生的磁链为${\Psi }_{21}$，则互感$M$定义为：

$$M=\frac{{\Psi }_{21}}{I_1}=\frac{N_2{\Phi }_{21}}{I_1}$$

互感具有以下重要性质：

• 互易性：$M_{12}=M_{21}=M$
• 对称性：互感与线圈的几何形状、相对位置和介质有关
• 可正可负：取决于线圈的绕向和相对位置

1.1.2 耦合系数

耦合系数$k$描述了两个线圈之间磁耦合的紧密程度：

$$k=\frac{M}{\sqrt{L_1{L}_2}}$$

其中：

• $L_1,L_2$：两个线圈的自感 [H]
• $M$：互感 [H]
• $k$：耦合系数 [0-1]

耦合系数的物理意义：

• $k=0$：无耦合（线圈相距很远或正交）
• $k=1$：全耦合（理想变压器）
• $0<k<1$：部分耦合（实际无线充电系统）

影响耦合系数的因素：

1. 线圈间距：距离越大，耦合越弱
2. 线圈对齐：横向偏移会显著降低耦合
3. 线圈尺寸：尺寸匹配度影响耦合效率
4. 介质材料：磁芯可增强耦合

1.2 等效电路模型

1.2.1 基本T型等效电路

无线电能传输系统可以用T型等效电路表示：

        L1 - M          M          L2 - M
    ━━━┯━━━━/\/\/\━━━━┯━━━━/\/\/\━━━━┯━━━
       │               │               │
      ═╪═ C1          ═╪═            ═╪═ C2
       │               │               │
    ━━━┷━━━━━━━━━━━━━━━┷━━━━━━━━━━━━━━━┷━━━
       │                               │
      ═╪═ R1                          ═╪═ R2 + RL
       │                               │
    ━━━┷━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┷━━━

其中：

• $L_1,L_2$：发射和接收线圈自感
• $M$：互感
• $R_1,R_2$：线圈电阻
• $C_1,C_2$：补偿电容
• $R_L$：负载电阻

1.2.2 反射阻抗

从发射端看，接收回路的影响可以等效为反射阻抗：

$$Z_r=\frac{(\omega M{)}^2}{Z_2}$$

其中$Z_2$是接收回路的总阻抗。反射阻抗的实部代表有功功率传输，虚部代表无功功率交换。

1.3 谐振与补偿

1.3.1 谐振条件

为了提高传输效率，通常使系统工作在谐振状态。谐振条件为：

$$\omega =\frac{1}{\sqrt{L_1{C}_1}}=\frac{1}{\sqrt{L_2{C}_2}}$$

谐振时：

• 电抗相互抵消，输入阻抗最小
• 无功功率为零，功率因数为1
• 传输效率达到最大

1.3.2 补偿网络类型

根据补偿电容的连接方式，可分为四种基本拓扑：

1. 串联-串联（SS）补偿

• 发射端和接收端都采用串联补偿
• 特点：输入阻抗随负载变化小，适合恒压输出
• 应用：电池充电、LED驱动

2. 串联-并联（SP）补偿

• 发射端串联，接收端并联
• 特点：输出电压稳定，适合恒压负载
• 应用：手机无线充电

3. 并联-串联（PS）补偿

• 发射端并联，接收端串联
• 特点：输入电流稳定，适合恒流源
• 应用：恒流充电

4. 并联-并联（PP）补偿

• 发射端和接收端都采用并联补偿
• 特点：输入阻抗高，适合高电压应用
• 应用：高功率传输

---

第二部分：线圈设计与磁场分析

2.1 线圈类型选择

2.1.1 圆形线圈

圆形线圈是最常用的无线充电线圈形式，具有对称性好、磁场分布均匀的特点。

自感计算公式：

$$L={\mu }_0{N}^{2}R\left[\ln \left(\frac{8R}{a}\right)-2\right]$$

其中：

• $N$：匝数
• $R$：线圈半径 [m]
• $a$：导线半径 [m]

优点：

• 磁场分布轴对称
• 易于设计和制造
• 耦合系数计算简单

缺点：

• 横向偏移容忍性一般
• 空间利用率不高

2.1.2 方形线圈

方形线圈在平面电路中更容易实现，常用于PCB集成。

自感计算公式：

$$L=\frac{{\mu }_0{N}^{2}s}{\pi }\left[\ln \left(\frac{s}{a}\right)+0.25\right]$$

其中$s$为边长。

优点：

• 便于PCB布线
• 空间利用率高
• 易于阵列扩展

缺点：

• 角落处磁场集中
• 耦合系数略低于圆形

2.1.3 螺旋线圈

螺旋线圈通过平面螺旋结构实现多匝绕组，适合薄型应用。

特点：

• 单层结构，厚度小
• 匝间电容较大
• 自谐振频率受限

2.2 磁场分布计算

2.2.1 单匝圆形线圈的磁场

根据毕奥-萨伐尔定律，单匝圆形线圈在轴线上产生的磁场为：

$$B_z=\frac{{\mu }_{0}I{R}^2}{2(R^2+z^2{)}^{3/2}}$$

其中：

• $I$：电流 [A]
• $R$：线圈半径 [m]
• $z$：轴向距离 [m]

在任意位置的磁场需要通过数值积分计算。

2.2.2 多匝线圈的磁场

多匝线圈的磁场是各匝磁场的叠加。对于紧密绕制的线圈，可以近似为：

$$B\approx N\cdot B_{single}$$

其中$N$为匝数，$B_{single}$为单匝磁场。

2.2.3 磁场可视化

通过数值计算可以得到磁场在空间中的分布：

• 磁力线：表示磁场方向的曲线
• 等磁密线：磁场强度相等的连线
• 矢量场：各点的磁场方向和大小

2.3 互感计算方法

2.3.1 纽曼公式

互感的严格计算使用纽曼（Neumann）公式：

$$M=\frac{{\mu }_0}{4\pi }{\oint }_{C_1}{\oint }_{C_2}\frac{d{\vec{l}}_1\cdot d{\vec{l}}_2}{|{\vec{r}}_1-{\vec{r}}_2|}$$

这是一个双重线积分，需要通过数值方法求解。

2.3.2 简化计算公式

对于同轴圆形线圈，互感可以用简化公式近似：

$$M\approx \frac{{\mu }_0\pi N_1{N}_2{R}_1^2{R}_2^2}{2(R_1^2+d^2{)}^{3/2}}$$

其中$d$为线圈间距。这个公式在$d\gg R$时不准确，但对于近距离耦合可用。

2.3.3 数值计算方法

实际工程中通常采用数值方法计算互感：

1. 将线圈离散为多个电流元
2. 计算每对电流元之间的互感
3. 对所有电流元对求和

---

第三部分：传输效率分析

3.1 效率计算公式

3.1.1 串联-串联（SS）补偿

对于SS补偿拓扑，传输效率为：

$$\eta =\frac{R_L}{R_2+R_L}\cdot \frac{(\omega M{)}^2/(R_2+R_L)}{R_1+(\omega M{)}^2/(R_2+R_L)}$$

当$R_1=R_2=R$且$R_L=R$时，效率简化为：

$$\eta =\frac{(\omega M{)}^2}{(2R{)}^2+(\omega M{)}^2}$$

3.1.2 品质因数与效率

定义线圈的品质因数：

$$Q_1=\frac{\omega L_1}{R_1},Q_2=\frac{\omega L_2}{R_2}$$

传输效率可以表示为：

$$\eta =\frac{k^2{Q}_1{Q}_2}{1+k^2{Q}_1{Q}_2}$$

这表明：

• 耦合系数$k$越大，效率越高
• 品质因数$Q$越高，效率越高
• 当$k^2{Q}_1{Q}_2\gg 1$时，效率接近100%

3.2 影响效率的因素

3.2.1 线圈间距

线圈间距$d$对效率的影响：

• 距离增加 → 耦合系数减小 → 效率降低
• 近距离（$d<R$）：效率高，但对准要求严格
• 远距离（$d>2R$）：效率低，但容忍度大

经验公式：

$$k\propto \frac{1}{d^3}$$

3.2.2 横向偏移

横向偏移$\Delta x$对效率的影响：

• 偏移增加 → 耦合系数减小 → 效率降低
• 偏移超过线圈半径时，效率急剧下降

3.2.3 工作频率

工作频率$f$对效率的影响：

• 频率增加 → 感抗增加 → $Q$值增加 → 效率提高
• 但频率过高会增加开关损耗和辐射损耗
• 常用频率范围：100 kHz - 10 MHz

3.3 效率优化策略

3.3.1 提高耦合系数

方法：

1. 减小线圈间距
2. 增加线圈尺寸
3. 使用磁芯增强磁场
4. 优化线圈形状（如DD线圈）

3.3.2 提高品质因数

方法：

1. 使用低电阻导线（利兹线）
2. 增加线圈匝数
3. 提高工作频率
4. 使用高$Q$值电容

3.3.3 阻抗匹配

通过调整负载阻抗或补偿网络，实现最大功率传输：

$$R_{L,opt}=R_2\sqrt{1+k^2{Q}_1{Q}_2}$$

---

第四部分：Python仿真实现

4.1 代码架构设计

本仿真代码采用面向对象设计，主要包含以下类：

1. CoilParameters：线圈参数数据类
2. CoilModel：线圈磁场模型
3. WPTSystemParameters：WPT系统参数
4. WPTSystemModel：WPT系统模型

4.2 核心算法实现

4.2.1 线圈几何建模

def get_coil_points(self, num_points: int = 100) -> np.ndarray:
    """
    获取线圈几何点坐标
    
    根据线圈类型（圆形、方形、螺旋）生成线圈上的点坐标，
    用于后续的磁场计算和互感计算。
    """
    params = self.params
    all_points = []
    
    if params.coil_type == CoilType.CIRCULAR:
        # 圆形线圈：生成同心圆
        for n in range(params.N):
            # 计算当前匝的半径
            r = params.R_inner + (params.R_outer - params.R_inner) * n / (params.N - 1)
            
            theta = np.linspace(0, 2*np.pi, num_points, endpoint=False)
            x = r * np.cos(theta)
            y = r * np.sin(theta)
            z = np.zeros_like(theta)
            
            points = np.column_stack([x, y, z])
            all_points.append(points)

这段代码实现了圆形线圈的几何建模。对于多匝线圈，每匝的半径从内半径线性变化到外半径，形成同心圆结构。

4.2.2 磁场计算（Biot-Savart定律）

def magnetic_field_single_turn(self, r_turn: float, I: float, 
                                x: np.ndarray, y: np.ndarray, z: np.ndarray):
    """
    计算单匝圆形线圈产生的磁场
    
    使用Biot-Savart定律数值积分计算磁场：
    dB = (μ₀I / 4π) (dl × r) / r³
    """
    num_segments = 100
    theta = np.linspace(0, 2*np.pi, num_segments, endpoint=False)
    d_theta = 2 * np.pi / num_segments
    
    # 线圈上的点
    x_c = r_turn * np.cos(theta)
    y_c = r_turn * np.sin(theta)
    z_c = np.zeros_like(theta)
    
    # 电流元 dl
    dx_c = -r_turn * np.sin(theta) * d_theta
    dy_c = r_turn * np.cos(theta) * d_theta
    
    # 初始化磁场
    Bx = np.zeros_like(x)
    By = np.zeros_like(y)
    Bz = np.zeros_like(z)
    
    # 计算每个电流元对场点的贡献
    for i in range(num_segments):
        # 从电流元指向场点的向量
        rx = x - x_c[i]
        ry = y - y_c[i]
        rz = z - z_c[i]
        
        # 距离
        r_mag = np.sqrt(rx**2 + ry**2 + rz**2)
        r_mag = np.maximum(r_mag, 1e-10)  # 避免除零
        
        # 叉积 dl × r
        cross_x = dy_c[i] * rz - dz_c[i] * ry
        cross_y = dz_c[i] * rx - dx_c[i] * rz
        cross_z = dx_c[i] * ry - dy_c[i] * rx
        
        # Biot-Savart定律
        dB = self.mu0 * I / (4 * np.pi * r_mag**3)
        
        Bx += dB * cross_x
        By += dB * cross_y
        Bz += dB * cross_z

这段代码实现了基于Biot-Savart定律的磁场数值计算。将线圈离散为多个电流元，计算每个电流元对场点的贡献并叠加。

4.2.3 互感计算

def mutual_inductance_numerical(self) -> float:
    """
    数值计算互感
    
    使用Neumann公式: M = (μ₀/4π) ∮∮ (dl₁·dl₂)/r
    """
    # 获取线圈点
    tx_points = self.tx_model.get_coil_points(num_points=30)
    rx_points = self.rx_model.get_coil_points(num_points=30)
    
    M = 0.0
    
    # 数值积分
    for i in range(len(tx_points) - 1):
        dl1 = tx_points[i+1] - tx_points[i]
        mid1 = (tx_points[i+1] + tx_points[i]) / 2
        
        for j in range(len(rx_points) - 1):
            dl2 = rx_points[j+1] - rx_points[j]
            mid2 = (rx_points[j+1] + rx_points[j]) / 2
            
            # 距离
            r_vec = mid2 - mid1
            r = np.linalg.norm(r_vec)
            
            if r > 1e-10:
                # 点积 dl₁·dl₂
                dot = np.dot(dl1, dl2)
                M += self.mu0 / (4 * np.pi) * dot / r
    
    return M

这段代码使用纽曼公式数值计算两个线圈之间的互感。通过双重线积分，计算发射线圈和接收线圈之间的磁耦合。

4.2.4 耦合系数计算

def coupling_coefficient(self) -> Tuple[float, float, float, float]:
    """
    计算耦合系数
    
    k = M / √(L₁ × L₂)
    """
    M = self.mutual_inductance_analytical()
    
    L1 = self.tx_model.self_inductance()
    L2 = self.rx_model.self_inductance()
    
    k = M / np.sqrt(L1 * L2)
    
    # 限制在合理范围
    k = np.clip(k, 0, 1)
    
    return k, M, L1, L2

耦合系数是衡量两个线圈磁耦合紧密程度的关键参数，取值范围为0到1。

4.2.5 传输效率计算

def transmission_efficiency(self, k: float, L1: float, L2: float) -> Dict:
    """
    计算传输效率
    
    对于SS补偿拓扑：
    η = P_out / P_in = (ωM)²RL / [(R₂+RL)²(R₁ + (ωM)²/(R₂+RL))]
    """
    params = self.params
    
    f = params.f_operating
    omega = 2 * np.pi * f
    R_L = params.R_load
    
    # 互感
    M = k * np.sqrt(L1 * L2)
    
    if params.compensation == CompensationType.SS:
        # 串联-串联补偿
        # 反射阻抗
        Z_reflected = (omega * M)**2 / (R2 + R_L)
        
        # 输入阻抗
        Z_in = R1 + Z_reflected
        
        # 输入功率
        I1 = params.V_in / Z_in
        P_in = params.V_in * I1
        
        # 次级电流
        I2 = omega * M * I1 / (R2 + R_L)
        
        # 输出功率
        P_out = I2**2 * R_L
        
        # 效率
        eta = P_out / P_in if P_in > 0 else 0
    
    return {
        'efficiency': eta,
        'mutual_inductance': M,
        'Q1': omega * L1 / R1,
        'Q2': omega * L2 / (R2 + R_L)
    }

这段代码计算了SS补偿拓扑下的传输效率，包括反射阻抗、输入阻抗、输入输出功率等中间参数。

4.3 可视化分析

代码实现了多种可视化功能：

1. 磁场分布图：展示线圈周围的磁场强度和磁力线分布
2. 效率特性图：效率随距离、偏移、负载、频率的变化
3. 补偿网络对比：四种补偿拓扑的效率对比
4. 线圈结构对比：不同形状线圈的磁场特性

4.4 GIF动画生成

代码生成3个GIF动画：

1. 磁场随距离变化：展示线圈间距变化时磁场分布的变化
2. 横向偏移影响：展示接收线圈横向偏移对耦合的影响
3. 能量传输动画：展示电流波形和功率传输的动态过程

---

第五部分：仿真结果分析

5.1 系统参数分析

从仿真结果可以看到典型无线充电系统的参数：

参数	数值	说明
发射线圈自感 L1	19.04 μH	10匝圆形线圈
接收线圈自感 L2	19.04 μH	与发射端对称
互感 M	13.68 μH	距离20mm时
耦合系数 k	0.7184	强耦合
补偿电容 C1/C2	133.07 nF	100kHz谐振
传输效率 η	97.67%	理想条件下
品质因数 Q1	119.60	发射端

5.2 磁场分布特征

5.2.1 轴向磁场分布

沿轴线（Z方向）的磁场分布特点：

• 线圈中心处磁场最强
• 随距离增加，磁场快速衰减
• 在接收线圈位置（Z=20mm），磁场约为发射端的30-40%

5.2.2 径向磁场分布

在垂直于轴线的平面内：

• 线圈中心区域磁场均匀
• 边缘处磁场快速下降
• 横向偏移会显著降低耦合

5.3 效率特性分析

5.3.1 距离对效率的影响

效率随距离的变化规律：

• 近距离（<10mm）：效率>95%
• 中等距离（20mm）：效率约90%
• 远距离（>50mm）：效率<50%

这种衰减主要由耦合系数的立方反比关系决定。

5.3.2 偏移对效率的影响

横向偏移对效率的影响：

• 无偏移：效率最高
• 偏移<10mm：效率下降约10%
• 偏移>30mm：效率急剧下降

这表明无线充电对对准精度有一定要求。

5.4 补偿网络对比

不同补偿拓扑的效率特性：

SS补偿：

• 效率随距离变化相对平缓
• 适合恒压输出应用
• 对负载变化不敏感

SP补偿：

• 输出电压稳定
• 适合电池充电
• 轻载时效率较低

PS补偿：

• 输入电流稳定
• 适合恒流源
• 控制较复杂

PP补偿：

• 输入阻抗高
• 适合高功率应用
• 对参数敏感

---

第六部分：工程设计实例

6.1 实例1：手机无线充电设计

设计规格：

• 输出功率：$P_{out}=10$ W
• 输出电压：$V_{out}=5$ V
• 工作频率：$f=100$ kHz
• 线圈间距：$d=3-5$ mm
• 目标效率：$\eta >75\%$

设计步骤：

1. 选择线圈尺寸：

   • 发射线圈直径：50 mm
   • 接收线圈直径：40 mm
   • 匝数：10匝

2. 计算自感：

   发射线圈（R=25mm, N=10）：
   $$L_1={\mu }_0{N}^{2}R\left[\ln \left(\frac{8R}{a}\right)-2\right]\approx 25\text{ μH}$$

   接收线圈（R=20mm, N=10）：
   $$L_2\approx 20\text{ μH}$$

3. 计算互感（d=5mm）：
   $$M\approx \frac{{\mu }_0\pi N_1{N}_2{R}_1^2{R}_2^2}{2(R_1^2+d^2{)}^{3/2}}\approx 8\text{ μH}$$

4. 计算耦合系数：
   $$k=\frac{M}{\sqrt{L_1{L}_2}}=\frac{8}{\sqrt{25\times 20}}\approx 0.36$$

5. 设计补偿电容：
   $$C_1=\frac{1}{(2\pi f{)}^2{L}_1}=\frac{1}{(2\pi \times 100\text{k}{)}^2\times 25\text{μ}}\approx 101\text{ nF}$$

   $$C_2=\frac{1}{(2\pi f{)}^2{L}_2}\approx 127\text{ nF}$$

6. 计算负载电阻：
   $$R_L=\frac{V_{out}^2}{P_{out}}=\frac{5^2}{10}=2.5\text{ Ω}$$

7. 估算效率：

   假设线圈电阻$R_1=R_2=0.2$ Ω：
   $$Q_1=\frac{\omega L_1}{R_1}=\frac{2\pi \times 100\text{k}\times 25\text{μ}}{0.2}\approx 78.5$$

   $$\eta =\frac{k^2{Q}_1{Q}_2}{1+k^2{Q}_1{Q}_2}\approx 80\%$$

8. 热设计：

   损耗功率：
   $$P_{loss}=P_{out}\times \frac{1-\eta }{\eta }=10\times \frac{0.2}{0.8}=2.5\text{ W}$$

   需要确保散热良好，温升不超过40°C。

6.2 实例2：电动汽车无线充电

设计规格：

• 输出功率：$P_{out}=3.3$ kW
• 输出电压：$V_{out}=400$ V
• 工作频率：$f=85$ kHz (SAE J2954标准)
• 气隙：$d=100-200$ mm
• 目标效率：$\eta >90\%$

设计要点：

1. 线圈设计：

   • 采用大直径线圈（400-600mm）
   • 使用利兹线减小趋肤效应
   • 铁氧体磁芯增强耦合

2. 补偿网络：

   • 采用双边LCC补偿
   • 实现恒流/恒压输出
   • 提高系统鲁棒性

3. 安全设计：

   • 活体检测
   • 金属异物检测
   • 过温保护

6.3 实例3：多设备无线充电

设计挑战：

• 同时给多个设备充电
• 各设备位置不固定
• 功率分配优化

解决方案：

1. 发射端阵列：

   • 多个小线圈阵列
   • 动态选择激活线圈
   • 实现区域覆盖

2. 控制策略：

   • 检测接收端位置
   • 动态调整工作频率
   • 优化功率分配

---

第七部分：常见问题与解决方案

7.1 效率低下

现象：

• 传输效率低于预期
• 发热严重

原因：

• 线圈对准不良
• 补偿网络失谐
• 线圈电阻过大
• 工作频率偏离谐振点

解决方案：

1. 增加对准引导机制
2. 使用自适应调谐
3. 采用利兹线降低电阻
4. 实现频率跟踪

7.2 EMI/EMC问题

现象：

• 干扰其他电子设备
• 辐射超标

原因：

• 工作频率选择不当
• 屏蔽不足
• 谐波含量高

解决方案：

1. 选择ISM频段（如6.78MHz, 13.56MHz）
2. 增加电磁屏蔽
3. 优化驱动波形
4. 使用滤波器

7.3 偏移容忍性不足

现象：

• 稍微移动设备就断充
• 充电位置要求严格

原因：

• 耦合系数对偏移敏感
• 线圈尺寸小
• 无位置自适应

解决方案：

1. 增大线圈尺寸
2. 采用多线圈阵列
3. 实现位置检测和自适应
4. 使用磁导引

7.4 热管理问题

现象：

• 充电器和设备发热
• 效率下降

原因：

• 损耗功率大
• 散热不良
• 环境温度高

解决方案：

1. 优化线圈设计降低损耗
2. 增加散热片或风扇
3. 限制最大输出功率
4. 实现温度监控

7.5 金属异物发热

现象：

• 金属物体在充电区域发热
• 安全隐患

原因：

• 涡流效应
• 磁场穿透金属

解决方案：

1. 金属异物检测（FOD）
2. 降低工作频率
3. 增加屏蔽层
4. 功率限制

---

第八部分：进阶话题

8.1 磁共振无线电能传输

8.1.1 磁共振原理

磁共振无线电能传输（MCR-WPT）利用谐振原理实现中距离高效传输：

• 发射和接收线圈都调谐到相同频率
• 通过强耦合磁共振实现能量传输
• 可以实现数米距离的高效传输

8.1.2 强耦合区域

当耦合系数$k$大于临界耦合系数$k_c$时，系统进入强耦合区域：

$$k_c=\frac{1}{\sqrt{Q_1{Q}_2}}$$

在强耦合区域：

• 传输效率对距离变化不敏感
• 存在频率分裂现象
• 可以实现高效率中距离传输

8.2 多线圈系统

8.2.1 中继线圈

在发射和接收线圈之间增加中继线圈：

• 扩展传输距离
• 绕过障碍物
• 实现多跳传输

8.2.2 多接收端

一个发射端同时给多个接收端供电：

• 需要功率分配控制
• 各接收端相互影响
• 适合物联网应用

8.3 新型线圈结构

8.3.1 DD线圈

双D形（DD）线圈：

• 提高横向偏移容忍性
• 磁场分布更均匀
• 适合电动汽车充电

8.3.2 极化线圈

利用极化特性：

• 减少相邻系统干扰
• 提高空间复用率
• 适合密集部署

8.4 智能控制技术

8.4.1 自适应调谐

实时调整补偿网络：

• 跟踪谐振频率
• 适应环境变化
• 保持高效率

8.4.2 通信与控制

发射与接收端通信：

• 功率需求协商
• 状态监测
• 安全保护

---

"""
主题044：无线充电系统磁场仿真
=====================================
本代码实现无线充电系统的磁场建模与仿真，包括：
1. 发射/接收线圈磁场分布计算
2. 互感与耦合系数计算
3. 传输效率分析
4. 磁场可视化与动画生成
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')  # 使用无头后端
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch
from matplotlib.collections import LineCollection
import json
from dataclasses import dataclass, asdict
from typing import Tuple, Optional, List, Dict
from enum import Enum
import warnings
warnings.filterwarnings('ignore')

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


class CoilType(Enum):
    """线圈类型"""
    CIRCULAR = "圆形"
    SQUARE = "方形"
    RECTANGULAR = "矩形"
    SPIRAL = "螺旋"


class CompensationType(Enum):
    """补偿网络类型"""
    SS = "串联-串联"  # 发射串联，接收串联
    SP = "串联-并联"  # 发射串联，接收并联
    PS = "并联-串联"  # 发射并联，接收串联
    PP = "并联-并联"  # 发射并联，接收并联


@dataclass
class CoilParameters:
    """线圈参数"""
    # 几何参数
    coil_type: CoilType = CoilType.CIRCULAR
    N: int = 10                    # 匝数
    R_outer: float = 0.05          # 外半径 [m]
    R_inner: float = 0.01          # 内半径 [m] (用于螺旋线圈)
    width: float = 0.05            # 宽度 [m] (用于矩形线圈)
    length: float = 0.05           # 长度 [m] (用于矩形线圈)
    wire_diameter: float = 0.001   # 线径 [m]
    
    # 电气参数
    L: float = 50e-6               # 自感 [H]
    R: float = 0.1                 # 电阻 [Ω]
    f_resonant: float = 100e3      # 谐振频率 [Hz]
    
    # 位置参数
    position: Tuple[float, float, float] = (0, 0, 0)  # (x, y, z) [m]
    orientation: Tuple[float, float] = (0, 0)         # (theta, phi) [rad]


@dataclass
class WPTSystemParameters:
    """无线电能传输系统参数"""
    # 发射线圈
    tx_coil: CoilParameters = None
    
    # 接收线圈
    rx_coil: CoilParameters = None
    
    # 系统参数
    f_operating: float = 100e3     # 工作频率 [Hz]
    V_in: float = 12.0             # 输入电压 [V]
    R_load: float = 10.0           # 负载电阻 [Ω]
    
    # 补偿网络
    compensation: CompensationType = CompensationType.SS
    
    # 耦合参数
    distance: float = 0.02         # 线圈间距 [m]
    misalignment_x: float = 0.0    # X方向偏移 [m]
    misalignment_y: float = 0.0    # Y方向偏移 [m]
    
    # 材料参数
    mu0: float = 4 * np.pi * 1e-7  # 真空磁导率
    
    def __post_init__(self):
        if self.tx_coil is None:
            self.tx_coil = CoilParameters()
        if self.rx_coil is None:
            self.rx_coil = CoilParameters(
                position=(0, 0, self.distance)
            )


class CoilModel:
    """线圈磁场模型"""
    
    def __init__(self, params: CoilParameters):
        self.params = params
        self.mu0 = 4 * np.pi * 1e-7
        
    def get_coil_points(self, num_points: int = 100) -> np.ndarray:
        """
        获取线圈几何点坐标
        
        参数:
            num_points: 每匝采样点数
        返回:
            points: 形状为 (N*num_points, 3) 的坐标数组
        """
        params = self.params
        all_points = []
        
        if params.coil_type == CoilType.CIRCULAR:
            # 圆形线圈
            for n in range(params.N):
                # 计算当前匝的半径
                if params.N == 1:
                    r = params.R_outer
                else:
                    r = params.R_inner + (params.R_outer - params.R_inner) * n / (params.N - 1)
                
                theta = np.linspace(0, 2*np.pi, num_points, endpoint=False)
                x = r * np.cos(theta)
                y = r * np.sin(theta)
                z = np.zeros_like(theta)
                
                points = np.column_stack([x, y, z])
                all_points.append(points)
                
        elif params.coil_type == CoilType.SQUARE:
            # 方形线圈
            for n in range(params.N):
                if params.N == 1:
                    s = 2 * params.R_outer
                else:
                    s = 2 * (params.R_inner + (params.R_outer - params.R_inner) * n / (params.N - 1))
                
                # 四条边
                side_points = num_points // 4
                
                # 底边
                x1 = np.linspace(-s/2, s/2, side_points)
                y1 = np.full_like(x1, -s/2)
                z1 = np.zeros_like(x1)
                
                # 右边
                y2 = np.linspace(-s/2, s/2, side_points)
                x2 = np.full_like(y2, s/2)
                z2 = np.zeros_like(y2)
                
                # 顶边
                x3 = np.linspace(s/2, -s/2, side_points)
                y3 = np.full_like(x3, s/2)
                z3 = np.zeros_like(x3)
                
                # 左边
                y4 = np.linspace(s/2, -s/2, side_points)
                x4 = np.full_like(y4, -s/2)
                z4 = np.zeros_like(y4)
                
                x = np.concatenate([x1, x2, x3, x4])
                y = np.concatenate([y1, y2, y3, y4])
                z = np.concatenate([z1, z2, z3, z4])
                
                points = np.column_stack([x, y, z])
                all_points.append(points)
                
        elif params.coil_type == CoilType.SPIRAL:
            # 螺旋线圈
            theta = np.linspace(0, 2*np.pi*params.N, num_points*params.N)
            r = params.R_inner + (params.R_outer - params.R_inner) * theta / (2*np.pi*params.N)
            
            x = r * np.cos(theta)
            y = r * np.sin(theta)
            z = np.zeros_like(theta)
            
            points = np.column_stack([x, y, z])
            all_points.append(points)
        
        # 合并所有点并应用位置偏移
        if all_points:
            points = np.vstack(all_points)
            # 应用位置偏移
            points[:, 0] += params.position[0]
            points[:, 1] += params.position[1]
            points[:, 2] += params.position[2]
            return points
        else:
            return np.array([])
    
    def magnetic_field_single_turn(self, r_turn: float, I: float, 
                                    x: np.ndarray, y: np.ndarray, z: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """
        计算单匝圆形线圈产生的磁场 (Biot-Savart定律)
        
        参数:
            r_turn: 线圈半径 [m]
            I: 电流 [A]
            x, y, z: 场点坐标 [m]
        返回:
            Bx, By, Bz: 磁场分量 [T]
        """
        # 数值积分计算磁场
        num_segments = 100
        theta = np.linspace(0, 2*np.pi, num_segments, endpoint=False)
        d_theta = 2 * np.pi / num_segments
        
        # 线圈上的点
        x_c = r_turn * np.cos(theta)
        y_c = r_turn * np.sin(theta)
        z_c = np.zeros_like(theta)
        
        # 电流元 dl
        dx_c = -r_turn * np.sin(theta) * d_theta
        dy_c = r_turn * np.cos(theta) * d_theta
        dz_c = np.zeros_like(theta)
        
        # 初始化磁场
        Bx = np.zeros_like(x)
        By = np.zeros_like(y)
        Bz = np.zeros_like(z)
        
        # 计算每个电流元对场点的贡献
        for i in range(num_segments):
            # 从电流元指向场点的向量
            rx = x - x_c[i]
            ry = y - y_c[i]
            rz = z - z_c[i]
            
            # 距离
            r_mag = np.sqrt(rx**2 + ry**2 + rz**2)
            r_mag = np.maximum(r_mag, 1e-10)  # 避免除零
            
            # 叉积 dl × r
            cross_x = dy_c[i] * rz - dz_c[i] * ry
            cross_y = dz_c[i] * rx - dx_c[i] * rz
            cross_z = dx_c[i] * ry - dy_c[i] * rx
            
            # Biot-Savart定律: dB = (mu0 * I / 4pi) * (dl × r) / r³
            dB = self.mu0 * I / (4 * np.pi * r_mag**3)
            
            Bx += dB * cross_x
            By += dB * cross_y
            Bz += dB * cross_z
        
        return Bx, By, Bz
    
    def magnetic_field(self, I: float, x: np.ndarray, y: np.ndarray, 
                       z: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """
        计算整个线圈的磁场
        
        参数:
            I: 电流 [A]
            x, y, z: 场点坐标 [m]
        返回:
            Bx, By, Bz: 磁场分量 [T]
        """
        params = self.params
        Bx = np.zeros_like(x)
        By = np.zeros_like(y)
        Bz = np.zeros_like(z)
        
        if params.coil_type == CoilType.CIRCULAR:
            # 多匝圆形线圈
            for n in range(params.N):
                if params.N == 1:
                    r = params.R_outer
                else:
                    r = params.R_inner + (params.R_outer - params.R_inner) * n / (params.N - 1)
                
                Bx_n, By_n, Bz_n = self.magnetic_field_single_turn(r, I, x, y, z)
                Bx += Bx_n
                By += By_n
                Bz += Bz_n
        else:
            # 对于其他类型，使用数值积分
            points = self.get_coil_points(num_points=50)
            if len(points) > 1:
                for i in range(len(points) - 1):
                    # 电流元
                    dl = points[i+1] - points[i]
                    mid_point = (points[i+1] + points[i]) / 2
                    
                    # 从电流元指向场点的向量
                    rx = x - mid_point[0]
                    ry = y - mid_point[1]
                    rz = z - mid_point[2]
                    
                    # 距离
                    r_mag = np.sqrt(rx**2 + ry**2 + rz**2)
                    r_mag = np.maximum(r_mag, 1e-10)
                    
                    # 叉积 dl × r
                    cross_x = dl[1] * rz - dl[2] * ry
                    cross_y = dl[2] * rx - dl[0] * rz
                    cross_z = dl[0] * ry - dl[1] * rx
                    
                    # Biot-Savart
                    dB = self.mu0 * I / (4 * np.pi * r_mag**3)
                    
                    Bx += dB * cross_x
                    By += dB * cross_y
                    Bz += dB * cross_z
        
        return Bx, By, Bz
    
    def self_inductance(self) -> float:
        """
        计算线圈自感 (简化公式)
        
        返回:
            L: 自感 [H]
        """
        params = self.params
        
        if params.coil_type == CoilType.CIRCULAR:
            # 圆形线圈自感近似公式
            # L ≈ μ0 * N² * R * [ln(8R/a) - 2]
            R_avg = (params.R_outer + params.R_inner) / 2
            a = params.wire_diameter / 2
            
            L_single = self.mu0 * R_avg * (np.log(8 * R_avg / a) - 2)
            L = params.N**2 * L_single
            
        elif params.coil_type == CoilType.SQUARE:
            # 方形线圈
            s = 2 * params.R_outer
            a = params.wire_diameter / 2
            
            L = self.mu0 * params.N**2 * s / np.pi * (np.log(s / a) + 0.25)
            
        else:
            # 使用参数中的值
            L = params.L
            
        return L


class WPTSystemModel:
    """无线电能传输系统模型"""
    
    def __init__(self, params: WPTSystemParameters):
        self.params = params
        self.tx_model = CoilModel(params.tx_coil)
        self.rx_model = CoilModel(params.rx_coil)
        self.mu0 = 4 * np.pi * 1e-7
        
    def mutual_inductance_numerical(self) -> float:
        """
        数值计算互感
        
        使用Neumann公式: M = (μ0/4π) ∮∮ (dl1·dl2)/r
        
        返回:
            M: 互感 [H]
        """
        # 获取线圈点
        tx_points = self.tx_model.get_coil_points(num_points=30)
        rx_points = self.rx_model.get_coil_points(num_points=30)
        
        if len(tx_points) == 0 or len(rx_points) == 0:
            return 0.0
        
        M = 0.0
        
        # 数值积分
        for i in range(len(tx_points) - 1):
            dl1 = tx_points[i+1] - tx_points[i]
            mid1 = (tx_points[i+1] + tx_points[i]) / 2
            
            for j in range(len(rx_points) - 1):
                dl2 = rx_points[j+1] - rx_points[j]
                mid2 = (rx_points[j+1] + rx_points[j]) / 2
                
                # 距离
                r_vec = mid2 - mid1
                r = np.linalg.norm(r_vec)
                
                if r > 1e-10:
                    # 点积 dl1·dl2
                    dot = np.dot(dl1, dl2)
                    M += self.mu0 / (4 * np.pi) * dot / r
        
        return M
    
    def mutual_inductance_analytical(self) -> float:
        """
        解析计算圆形线圈互感 (简化公式)
        
        适用于同轴圆形线圈
        
        返回:
            M: 互感 [H]
        """
        params = self.params
        
        # 提取线圈参数
        R1 = params.tx_coil.R_outer
        R2 = params.rx_coil.R_outer
        d = params.distance
        N1 = params.tx_coil.N
        N2 = params.rx_coil.N
        
        # 使用椭圆积分近似
        # 对于同轴线圈，可以使用Maxwell公式
        
        # 简化公式 (当 d >> R 时不准确，但对于近距离耦合可用)
        # M ≈ μ0 * π * N1 * N2 * R1² * R2² / (2 * (R1² + d²)^(3/2))
        
        # 更准确的近似
        k_squared = 4 * R1 * R2 / ((R1 + R2)**2 + d**2)
        
        # 简化计算
        M = self.mu0 * N1 * N2 * np.sqrt(R1 * R2) * (
            (2 / k_squared - 1) * np.log(4 / np.sqrt(1 - k_squared + 1e-10)) - 2 / k_squared + 1
        )
        
        # 如果结果不合理，使用简化公式
        if M < 0 or not np.isfinite(M):
            M = self.mu0 * np.pi * N1 * N2 * R1**2 * R2**2 / (2 * (R1**2 + d**2)**1.5)
        
        return M
    
    def coupling_coefficient(self, use_numerical: bool = False) -> float:
        """
        计算耦合系数
        
        k = M / sqrt(L1 * L2)
        
        参数:
            use_numerical: 是否使用数值方法计算互感
        返回:
            k: 耦合系数 [0-1]
        """
        if use_numerical:
            M = self.mutual_inductance_numerical()
        else:
            M = self.mutual_inductance_analytical()
        
        L1 = self.tx_model.self_inductance()
        L2 = self.rx_model.self_inductance()
        
        k = M / np.sqrt(L1 * L2)
        
        # 限制在合理范围
        k = np.clip(k, 0, 1)
        
        return k, M, L1, L2
    
    def compensation_capacitors(self, L1: float, L2: float) -> Tuple[float, float]:
        """
        计算补偿电容
        
        参数:
            L1: 发射线圈自感 [H]
            L2: 接收线圈自感 [H]
        返回:
            C1: 发射端补偿电容 [F]
            C2: 接收端补偿电容 [F]
        """
        f = self.params.f_operating
        omega = 2 * np.pi * f
        
        # 串联补偿: 1/(omega²L)
        C1 = 1 / (omega**2 * L1)
        C2 = 1 / (omega**2 * L2)
        
        return C1, C2
    
    def transmission_efficiency(self, k: float, L1: float, L2: float, 
                                 R1: Optional[float] = None, 
                                 R2: Optional[float] = None) -> Dict:
        """
        计算传输效率
        
        参数:
            k: 耦合系数
            L1, L2: 自感 [H]
            R1, R2: 线圈电阻 [Ω]
        返回:
            效率相关参数字典
        """
        params = self.params
        
        if R1 is None:
            R1 = params.tx_coil.R
        if R2 is None:
            R2 = params.rx_coil.R
        
        f = params.f_operating
        omega = 2 * np.pi * f
        R_L = params.R_load
        
        # 互感
        M = k * np.sqrt(L1 * L2)
        
        # 补偿电容
        C1, C2 = self.compensation_capacitors(L1, L2)
        
        # 根据补偿类型计算效率
        if params.compensation == CompensationType.SS:
            # 串联-串联补偿
            # 反射阻抗
            Z_reflected = (omega * M)**2 / (R2 + R_L)
            
            # 输入阻抗
            Z_in = R1 + Z_reflected
            
            # 输入功率
            I1 = params.V_in / Z_in
            P_in = params.V_in * I1
            
            # 次级电流
            I2 = omega * M * I1 / (R2 + R_L)
            
            # 输出功率
            P_out = I2**2 * R_L
            
            # 效率
            eta = P_out / P_in if P_in > 0 else 0
            
        elif params.compensation == CompensationType.SP:
            # 串联-并联补偿
            # 简化的效率计算
            Q1 = omega * L1 / R1
            Q2 = omega * L2 / (R2 + R_L)
            eta = (k**2 * Q1 * Q2) / (1 + k**2 * Q1 * Q2)
            
        else:
            # 其他类型使用近似公式
            Q1 = omega * L1 / R1
            Q2 = omega * L2 / (R2 + R_L)
            eta = (k**2 * Q1 * Q2) / (1 + k**2 * Q1 * Q2 + (Q2/Q1))
        
        return {
            'efficiency': eta,
            'mutual_inductance': M,
            'C1': C1,
            'C2': C2,
            'Q1': omega * L1 / R1,
            'Q2': omega * L2 / (R2 + R_L)
        }
    
    def efficiency_vs_distance(self, distances: np.ndarray) -> np.ndarray:
        """
        计算不同距离下的效率
        
        参数:
            distances: 距离数组 [m]
        返回:
            efficiencies: 效率数组
        """
        efficiencies = []
        
        for d in distances:
            # 临时修改距离
            original_distance = self.params.distance
            self.params.distance = d
            self.params.rx_coil.position = (0, 0, d)
            
            # 计算耦合系数和效率
            k, M, L1, L2 = self.coupling_coefficient()
            result = self.transmission_efficiency(k, L1, L2)
            efficiencies.append(result['efficiency'])
            
            # 恢复原始距离
            self.params.distance = original_distance
            self.params.rx_coil.position = (0, 0, original_distance)
        
        return np.array(efficiencies)
    
    def efficiency_vs_misalignment(self, offsets: np.ndarray) -> np.ndarray:
        """
        计算不同偏移下的效率
        
        参数:
            offsets: 偏移数组 [m]
        返回:
            efficiencies: 效率数组
        """
        efficiencies = []
        
        for offset in offsets:
            # 临时修改偏移
            self.params.rx_coil.position = (offset, 0, self.params.distance)
            
            # 计算耦合系数和效率
            k, M, L1, L2 = self.coupling_coefficient()
            result = self.transmission_efficiency(k, L1, L2)
            efficiencies.append(result['efficiency'])
            
            # 恢复
            self.params.rx_coil.position = (0, 0, self.params.distance)
        
        return np.array(efficiencies)


def create_visualizations():
    """创建可视化图表"""
    
    print("=" * 70)
    print("主题044：无线充电系统磁场仿真")
    print("=" * 70)
    
    # 创建系统参数
    tx_coil = CoilParameters(
        coil_type=CoilType.CIRCULAR,
        N=10,
        R_outer=0.05,
        R_inner=0.02,
        L=50e-6,
        R=0.1,
        position=(0, 0, 0)
    )
    
    rx_coil = CoilParameters(
        coil_type=CoilType.CIRCULAR,
        N=10,
        R_outer=0.05,
        R_inner=0.02,
        L=50e-6,
        R=0.1,
        position=(0, 0, 0.02)
    )
    
    params = WPTSystemParameters(
        tx_coil=tx_coil,
        rx_coil=rx_coil,
        f_operating=100e3,
        V_in=12.0,
        R_load=10.0,
        distance=0.02,
        compensation=CompensationType.SS
    )
    
    model = WPTSystemModel(params)
    
    # 计算基本参数
    print("\n1. 计算系统基本参数...")
    k, M, L1, L2 = model.coupling_coefficient()
    C1, C2 = model.compensation_capacitors(L1, L2)
    result = model.transmission_efficiency(k, L1, L2)
    
    print(f"   发射线圈自感 L1 = {L1*1e6:.2f} μH")
    print(f"   接收线圈自感 L2 = {L2*1e6:.2f} μH")
    print(f"   互感 M = {M*1e6:.2f} μH")
    print(f"   耦合系数 k = {k:.4f}")
    print(f"   发射端补偿电容 C1 = {C1*1e9:.2f} nF")
    print(f"   接收端补偿电容 C2 = {C2*1e9:.2f} nF")
    print(f"   传输效率 η = {result['efficiency']*100:.2f}%")
    print(f"   品质因数 Q1 = {result['Q1']:.2f}, Q2 = {result['Q2']:.2f}")
    
    # 图1: 线圈磁场分布
    print("\n2. 生成线圈磁场分布图...")
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 创建网格
    x = np.linspace(-0.1, 0.1, 50)
    z = np.linspace(-0.02, 0.06, 50)
    X, Z = np.meshgrid(x, z)
    Y = np.zeros_like(X)
    
    # 计算磁场
    I_tx = 1.0  # 发射线圈电流
    Bx, By, Bz = model.tx_model.magnetic_field(I_tx, X, Y, Z)
    B_magnitude = np.sqrt(Bx**2 + By**2 + Bz**2)
    
    # 子图1: 磁场强度分布
    ax = axes[0, 0]
    im = ax.contourf(X*1000, Z*1000, B_magnitude*1e6, levels=50, cmap='viridis')
    plt.colorbar(im, ax=ax, label='|B| [μT]')
    
    # 绘制线圈
    theta = np.linspace(0, 2*np.pi, 100)
    r_tx = tx_coil.R_outer * 1000
    r_rx = rx_coil.R_outer * 1000
    ax.plot(r_tx * np.cos(theta), np.zeros_like(theta), 'r-', linewidth=3, label='发射线圈')
    ax.plot(r_rx * np.cos(theta), np.ones_like(theta)*params.distance*1000, 'b-', linewidth=3, label='接收线圈')
    
    ax.set_xlabel('X [mm]', fontsize=11)
    ax.set_ylabel('Z [mm]', fontsize=11)
    ax.set_title('磁场强度分布 (XZ平面)', fontsize=12)
    ax.legend()
    ax.set_aspect('equal')
    
    # 子图2: 磁力线
    ax = axes[0, 1]
    
    # 绘制磁力线
    ax.streamplot(X*1000, Z*1000, Bx*1e6, Bz*1e6, 
                  color=B_magnitude*1e6, cmap='plasma', density=1.5, linewidth=1)
    
    # 绘制线圈
    ax.plot(r_tx * np.cos(theta), np.zeros_like(theta), 'w-', linewidth=3)
    ax.plot(r_rx * np.cos(theta), np.ones_like(theta)*params.distance*1000, 'w-', linewidth=3)
    
    ax.set_xlabel('X [mm]', fontsize=11)
    ax.set_ylabel('Z [mm]', fontsize=11)
    ax.set_title('磁力线分布', fontsize=12)
    ax.set_aspect('equal')
    
    # 子图3: 轴向磁场分布
    ax = axes[1, 0]
    
    z_axis = np.linspace(-0.02, 0.06, 200)
    x_center = np.zeros_like(z_axis)
    y_center = np.zeros_like(z_axis)
    
    Bx_axis, By_axis, Bz_axis = model.tx_model.magnetic_field(I_tx, x_center, y_center, z_axis)
    B_axis = np.sqrt(Bx_axis**2 + By_axis**2 + Bz_axis**2)
    
    ax.plot(z_axis*1000, B_axis*1e6, 'b-', linewidth=2, label='|B|')
    ax.axvline(x=0, color='r', linestyle='--', alpha=0.5, label='发射线圈')
    ax.axvline(x=params.distance*1000, color='g', linestyle='--', alpha=0.5, label='接收线圈')
    ax.set_xlabel('Z [mm]', fontsize=11)
    ax.set_ylabel('|B| [μT]', fontsize=11)
    ax.set_title('轴向磁场分布', fontsize=12)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图4: 径向磁场分布
    ax = axes[1, 1]
    
    x_radial = np.linspace(-0.1, 0.1, 200)
    z_tx = np.zeros_like(x_radial)
    z_rx = np.full_like(x_radial, params.distance)
    y_radial = np.zeros_like(x_radial)
    
    Bx_tx, By_tx, Bz_tx = model.tx_model.magnetic_field(I_tx, x_radial, y_radial, z_tx)
    B_tx = np.sqrt(Bx_tx**2 + By_tx**2 + Bz_tx**2)
    
    Bx_rx, By_rx, Bz_rx = model.tx_model.magnetic_field(I_tx, x_radial, y_radial, z_rx)
    B_rx = np.sqrt(Bx_rx**2 + By_rx**2 + Bz_rx**2)
    
    ax.plot(x_radial*1000, B_tx*1e6, 'r-', linewidth=2, label='发射线圈平面')
    ax.plot(x_radial*1000, B_rx*1e6, 'b-', linewidth=2, label='接收线圈平面')
    ax.set_xlabel('X [mm]', fontsize=11)
    ax.set_ylabel('|B| [μT]', fontsize=11)
    ax.set_title('径向磁场分布', fontsize=12)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig1_magnetic_field_distribution.png', dpi=150, bbox_inches='tight')
    print("   fig1_magnetic_field_distribution.png 已保存")
    plt.close()
    
    # 图2: 效率特性分析
    print("\n3. 生成效率特性分析图...")
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    # 子图1: 效率vs距离
    ax = axes[0, 0]
    distances = np.linspace(0.005, 0.1, 100)
    efficiencies = model.efficiency_vs_distance(distances)
    
    ax.plot(distances*1000, efficiencies*100, 'b-', linewidth=2)
    ax.axvline(x=params.distance*1000, color='r', linestyle='--', alpha=0.5, 
               label=f'设计距离 = {params.distance*1000:.1f} mm')
    ax.set_xlabel('线圈间距 [mm]', fontsize=11)
    ax.set_ylabel('传输效率 [%]', fontsize=11)
    ax.set_title('传输效率 vs 线圈间距', fontsize=12)
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 100)
    
    # 子图2: 效率vs横向偏移
    ax = axes[0, 1]
    offsets = np.linspace(0, 0.08, 100)
    efficiencies_offset = model.efficiency_vs_misalignment(offsets)
    
    ax.plot(offsets*1000, efficiencies_offset*100, 'g-', linewidth=2)
    ax.set_xlabel('横向偏移 [mm]', fontsize=11)
    ax.set_ylabel('传输效率 [%]', fontsize=11)
    ax.set_title('传输效率 vs 横向偏移', fontsize=12)
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 100)
    
    # 子图3: 耦合系数vs距离
    ax = axes[1, 0]
    k_values = []
    for d in distances:
        original_distance = model.params.distance
        model.params.distance = d
        model.params.rx_coil.position = (0, 0, d)
        k, _, _, _ = model.coupling_coefficient()
        k_values.append(k)
        model.params.distance = original_distance
        model.params.rx_coil.position = (0, 0, original_distance)
    
    ax.plot(distances*1000, k_values, 'r-', linewidth=2)
    ax.set_xlabel('线圈间距 [mm]', fontsize=11)
    ax.set_ylabel('耦合系数 k', fontsize=11)
    ax.set_title('耦合系数 vs 线圈间距', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    # 子图4: 互感vs距离
    ax = axes[1, 1]
    M_values = []
    for d in distances:
        original_distance = model.params.distance
        model.params.distance = d
        model.params.rx_coil.position = (0, 0, d)
        _, M, _, _ = model.coupling_coefficient()
        M_values.append(M*1e6)
        model.params.distance = original_distance
        model.params.rx_coil.position = (0, 0, original_distance)
    
    ax.plot(distances*1000, M_values, 'm-', linewidth=2)
    ax.set_xlabel('线圈间距 [mm]', fontsize=11)
    ax.set_ylabel('互感 M [μH]', fontsize=11)
    ax.set_title('互感 vs 线圈间距', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig2_efficiency_analysis.png', dpi=150, bbox_inches='tight')
    print("   fig2_efficiency_analysis.png 已保存")
    plt.close()
    
    # 图3: 补偿网络对比
    print("\n4. 生成补偿网络对比图...")
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))
    
    compensations = [
        CompensationType.SS,
        CompensationType.SP,
        CompensationType.PS,
        CompensationType.PP
    ]
    
    colors = ['blue', 'red', 'green', 'purple']
    
    # 子图1: 不同补偿网络的效率对比
    ax = axes[0, 0]
    
    for comp, color in zip(compensations, colors):
        model.params.compensation = comp
        efficiencies_comp = model.efficiency_vs_distance(distances)
        ax.plot(distances*1000, efficiencies_comp*100, 
                color=color, linewidth=2, label=comp.value)
    
    ax.set_xlabel('线圈间距 [mm]', fontsize=11)
    ax.set_ylabel('传输效率 [%]', fontsize=11)
    ax.set_title('不同补偿网络的效率对比', fontsize=12)
    ax.legend()
    ax.grid(True, alpha=0.3)
    ax.set_ylim(0, 100)
    
    # 恢复默认补偿
    model.params.compensation = params.compensation
    
    # 子图2: 负载特性
    ax = axes[0, 1]
    
    R_loads = np.linspace(1, 50, 100)
    efficiencies_load = []
    
    for R_L in R_loads:
        model.params.R_load = R_L
        k, M, L1, L2 = model.coupling_coefficient()
        result = model.transmission_efficiency(k, L1, L2)
        efficiencies_load.append(result['efficiency']*100)
    
    model.params.R_load = params.R_load
    
    ax.plot(R_loads, efficiencies_load, 'b-', linewidth=2)
    ax.axvline(x=params.R_load, color='r', linestyle='--', alpha=0.5,
               label=f'设计负载 = {params.R_load} Ω')
    ax.set_xlabel('负载电阻 [Ω]', fontsize=11)
    ax.set_ylabel('传输效率 [%]', fontsize=11)
    ax.set_title('效率 vs 负载电阻', fontsize=12)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图3: 频率特性
    ax = axes[1, 0]
    
    frequencies = np.linspace(50e3, 200e3, 200)
    efficiencies_freq = []
    
    for f in frequencies:
        model.params.f_operating = f
        k, M, L1, L2 = model.coupling_coefficient()
        result = model.transmission_efficiency(k, L1, L2)
        efficiencies_freq.append(result['efficiency']*100)
    
    model.params.f_operating = params.f_operating
    
    ax.plot(frequencies/1e3, efficiencies_freq, 'g-', linewidth=2)
    ax.axvline(x=params.f_operating/1e3, color='r', linestyle='--', alpha=0.5,
               label=f'谐振频率 = {params.f_operating/1e3:.0f} kHz')
    ax.set_xlabel('频率 [kHz]', fontsize=11)
    ax.set_ylabel('传输效率 [%]', fontsize=11)
    ax.set_title('效率 vs 工作频率', fontsize=12)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 子图4: 品质因数影响
    ax = axes[1, 1]
    
    Q_factors = np.linspace(10, 200, 100)
    efficiencies_Q = []
    
    for Q in Q_factors:
        # 假设Q1 = Q2 = Q
        R1 = 2 * np.pi * params.f_operating * L1 / Q
        R2 = 2 * np.pi * params.f_operating * L2 / Q
        
        k, M, L1, L2 = model.coupling_coefficient()
        result = model.transmission_efficiency(k, L1, L2, R1, R2)
        efficiencies_Q.append(result['efficiency']*100)
    
    ax.plot(Q_factors, efficiencies_Q, 'm-', linewidth=2)
    ax.set_xlabel('品质因数 Q', fontsize=11)
    ax.set_ylabel('传输效率 [%]', fontsize=11)
    ax.set_title('效率 vs 品质因数', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig3_compensation_comparison.png', dpi=150, bbox_inches='tight')
    print("   fig3_compensation_comparison.png 已保存")
    plt.close()
    
    # 图4: 线圈结构对比
    print("\n5. 生成线圈结构对比图...")
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    coil_types = [
        (CoilType.CIRCULAR, "圆形线圈"),
        (CoilType.SQUARE, "方形线圈"),
        (CoilType.SPIRAL, "螺旋线圈")
    ]
    
    for idx, (coil_type, name) in enumerate(coil_types):
        ax = axes[idx // 2, idx % 2]
        
        # 创建线圈模型
        test_coil = CoilParameters(
            coil_type=coil_type,
            N=10,
            R_outer=0.05,
            R_inner=0.02,
            position=(0, 0, 0)
        )
        coil_model = CoilModel(test_coil)
        points = coil_model.get_coil_points(num_points=50)
        
        # 绘制线圈
        if len(points) > 0:
            ax.plot(points[:, 0]*1000, points[:, 1]*1000, 'b-', linewidth=2)
            ax.scatter(points[::10, 0]*1000, points[::10, 1]*1000, 
                      c='red', s=20, zorder=5)
        
        ax.set_xlabel('X [mm]', fontsize=11)
        ax.set_ylabel('Y [mm]', fontsize=11)
        ax.set_title(name, fontsize=12)
        ax.set_aspect('equal')
        ax.grid(True, alpha=0.3)
        
        # 计算并显示自感
        L = coil_model.self_inductance()
        ax.text(0.05, 0.95, f'自感 L = {L*1e6:.1f} μH', 
               transform=ax.transAxes, fontsize=10,
               verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    # 第4个子图：不同匝数的对比
    ax = axes[1, 1]
    
    N_values = [5, 10, 15, 20]
    L_values = []
    
    for N in N_values:
        test_coil = CoilParameters(
            coil_type=CoilType.CIRCULAR,
            N=N,
            R_outer=0.05,
            R_inner=0.02
        )
        coil_model = CoilModel(test_coil)
        L = coil_model.self_inductance()
        L_values.append(L*1e6)
    
    ax.plot(N_values, L_values, 'bo-', linewidth=2, markersize=8)
    ax.set_xlabel('匝数 N', fontsize=11)
    ax.set_ylabel('自感 L [μH]', fontsize=11)
    ax.set_title('自感 vs 匝数', fontsize=12)
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('fig4_coil_comparison.png', dpi=150, bbox_inches='tight')
    print("   fig4_coil_comparison.png 已保存")
    plt.close()
    
    return model, params


def create_gif_animations():
    """创建GIF动画"""
    
    print("\n6. 生成GIF动画...")
    
    from matplotlib.animation import FuncAnimation
    from PIL import Image
    import io
    
    # 创建系统
    tx_coil = CoilParameters(
        coil_type=CoilType.CIRCULAR,
        N=10,
        R_outer=0.05,
        R_inner=0.02,
        L=50e-6,
        R=0.1,
        position=(0, 0, 0)
    )
    
    rx_coil = CoilParameters(
        coil_type=CoilType.CIRCULAR,
        N=10,
        R_outer=0.05,
        R_inner=0.02,
        L=50e-6,
        R=0.1,
        position=(0, 0, 0.02)
    )
    
    params = WPTSystemParameters(
        tx_coil=tx_coil,
        rx_coil=rx_coil,
        f_operating=100e3,
        V_in=12.0,
        R_load=10.0,
        distance=0.02,
        compensation=CompensationType.SS
    )
    
    model = WPTSystemModel(params)
    
    # 动画1: 磁场随距离变化
    print("  生成磁场随距离变化动画...")
    
    distances_anim = np.linspace(0.01, 0.08, 30)
    frames1 = []
    
    for d in distances_anim:
        fig, axes = plt.subplots(1, 2, figsize=(14, 6))
        
        # 更新接收线圈位置
        model.params.distance = d
        model.params.rx_coil.position = (0, 0, d)
        
        # 计算耦合系数
        k, M, L1, L2 = model.coupling_coefficient()
        
        # 创建网格
        x = np.linspace(-0.1, 0.1, 40)
        z = np.linspace(-0.02, max(0.1, d + 0.05), 40)
        X, Z = np.meshgrid(x, z)
        Y = np.zeros_like(X)
        
        # 计算磁场
        Bx, By, Bz = model.tx_model.magnetic_field(1.0, X, Y, Z)
        B_magnitude = np.sqrt(Bx**2 + By**2 + Bz**2)
        
        # 左图：磁场强度
        ax = axes[0]
        im = ax.contourf(X*1000, Z*1000, B_magnitude*1e6, levels=30, cmap='viridis', vmin=0, vmax=50)
        plt.colorbar(im, ax=ax, label='|B| [μT]')
        
        theta = np.linspace(0, 2*np.pi, 100)
        ax.plot(50 * np.cos(theta), np.zeros_like(theta), 'r-', linewidth=3, label='发射线圈')
        ax.plot(50 * np.cos(theta), np.ones_like(theta)*d*1000, 'b-', linewidth=3, label='接收线圈')
        
        ax.set_xlabel('X [mm]', fontsize=11)
        ax.set_ylabel('Z [mm]', fontsize=11)
        ax.set_title(f'磁场分布 (距离 = {d*1000:.1f} mm)', fontsize=12)
        ax.legend(loc='upper right')
        ax.set_aspect('equal')
        
        # 右图：参数显示
        ax = axes[1]
        ax.axis('off')
        
        info_text = f"""
        无线充电系统参数
        ====================
        
        线圈间距: {d*1000:.1f} mm
        耦合系数: {k:.4f}
        互感: {M*1e6:.2f} μH
        
        发射线圈自感: {L1*1e6:.1f} μH
        接收线圈自感: {L2*1e6:.1f} μH
        
        工作频率: {params.f_operating/1e3:.0f} kHz
        """
        
        ax.text(0.1, 0.5, info_text, transform=ax.transAxes, fontsize=12,
               verticalalignment='center', fontfamily='monospace',
               bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
        
        plt.tight_layout()
        
        # 保存到内存
        buf = io.BytesIO()
        plt.savefig(buf, format='png', dpi=100, bbox_inches='tight')
        buf.seek(0)
        frames1.append(Image.open(buf))
        plt.close()
    
    # 保存GIF
    frames1[0].save('anim1_field_vs_distance.gif',
                   save_all=True, append_images=frames1[1:],
                   duration=200, loop=0)
    print("    anim1_field_vs_distance.gif 已保存")
    
    # 动画2: 横向偏移对效率的影响
    print("  生成横向偏移动画...")
    
    offsets_anim = np.linspace(0, 0.08, 30)
    frames2 = []
    
    for offset in offsets_anim:
        fig, axes = plt.subplots(1, 2, figsize=(14, 6))
        
        # 更新接收线圈位置
        model.params.rx_coil.position = (offset, 0, params.distance)
        
        # 计算耦合系数和效率
        k, M, L1, L2 = model.coupling_coefficient()
        result = model.transmission_efficiency(k, L1, L2)
        
        # 左图：线圈位置示意
        ax = axes[0]
        
        # 绘制发射线圈
        theta = np.linspace(0, 2*np.pi, 100)
        ax.plot(50 * np.cos(theta), 50 * np.sin(theta), 'r-', linewidth=3, label='发射线圈')
        
        # 绘制接收线圈（偏移）
        ax.plot(50 * np.cos(theta) + offset*1000, 50 * np.sin(theta), 'b-', linewidth=3, label='接收线圈')
        
        # 绘制磁场耦合示意
        ax.annotate('', xy=(offset*1000, 0), xytext=(0, 0),
                   arrowprops=dict(arrowstyle='<->', color='green', lw=2))
        ax.text(offset*1000/2, 5, f'偏移 = {offset*1000:.1f} mm', 
               ha='center', fontsize=10, color='green')
        
        ax.set_xlabel('X [mm]', fontsize=11)
        ax.set_ylabel('Y [mm]', fontsize=11)
        ax.set_title('线圈位置示意 (俯视图)', fontsize=12)
        ax.legend()
        ax.set_aspect('equal')
        ax.set_xlim(-80, 80)
        ax.set_ylim(-80, 80)
        ax.grid(True, alpha=0.3)
        
        # 右图：参数显示
        ax = axes[1]
        ax.axis('off')
        
        info_text = f"""
        横向偏移影响分析
        ====================
        
        横向偏移: {offset*1000:.1f} mm
        耦合系数: {k:.4f}
        传输效率: {result['efficiency']*100:.1f}%
        
        互感: {M*1e6:.2f} μH
        
        发射端Q值: {result['Q1']:.1f}
        接收端Q值: {result['Q2']:.1f}
        """
        
        ax.text(0.1, 0.5, info_text, transform=ax.transAxes, fontsize=12,
               verticalalignment='center', fontfamily='monospace',
               bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
        
        plt.tight_layout()
        
        buf = io.BytesIO()
        plt.savefig(buf, format='png', dpi=100, bbox_inches='tight')
        buf.seek(0)
        frames2.append(Image.open(buf))
        plt.close()
    
    frames2[0].save('anim2_misalignment_effect.gif',
                   save_all=True, append_images=frames2[1:],
                   duration=200, loop=0)
    print("    anim2_misalignment_effect.gif 已保存")
    
    # 动画3: 电流波形与能量传输
    print("  生成能量传输动画...")
    
    t = np.linspace(0, 3e-5, 100)  # 3个周期
    f = params.f_operating
    omega = 2 * np.pi * f
    
    frames3 = []
    
    for i in range(0, len(t), 3):  # 33帧
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        
        idx = int(i)
        t_current = t[:idx+1]
        
        # 发射端电流
        I1 = np.sin(omega * t_current)
        # 接收端电流（相位差）
        k, M, L1, L2 = model.coupling_coefficient()
        I2 = k * np.sin(omega * t_current - np.pi/2)
        
        # 子图1: 电流波形
        ax = axes[0, 0]
        ax.plot(t_current*1e6, I1, 'r-', linewidth=2, label='发射端电流')
        ax.plot(t_current*1e6, I2, 'b-', linewidth=2, label='接收端电流')
        ax.set_xlabel('时间 [μs]', fontsize=11)
        ax.set_ylabel('电流 [A]', fontsize=11)
        ax.set_title('线圈电流波形', fontsize=12)
        ax.legend()
        ax.grid(True, alpha=0.3)
        ax.set_xlim(0, 30)
        
        # 子图2: 瞬时功率
        ax = axes[0, 1]
        P_in = params.V_in * I1
        P_out = I2**2 * params.R_load
        ax.plot(t_current*1e6, P_in, 'r-', linewidth=2, label='输入功率')
        ax.plot(t_current*1e6, P_out, 'b-', linewidth=2, label='输出功率')
        ax.set_xlabel('时间 [μs]', fontsize=11)
        ax.set_ylabel('功率 [W]', fontsize=11)
        ax.set_title('瞬时功率', fontsize=12)
        ax.legend()
        ax.grid(True, alpha=0.3)
        ax.set_xlim(0, 30)
        
        # 子图3: 磁场矢量图
        ax = axes[1, 0]
        
        # 简化显示
        theta = np.linspace(0, 2*np.pi, 8, endpoint=False)
        for j, th in enumerate(theta):
            r = 0.03
            x_pos = r * np.cos(th)
            y_pos = r * np.sin(th)
            
            # 磁场方向（简化）
            phase = omega * t_current[-1] if len(t_current) > 0 else 0
            B_mag = 20 * np.sin(phase + th)
            
            color = 'red' if B_mag > 0 else 'blue'
            ax.arrow(x_pos*1000, y_pos*1000, 0, B_mag*0.5, 
                    head_width=2, head_length=1, fc=color, ec=color, alpha=0.7)
        
        # 绘制线圈
        theta_coil = np.linspace(0, 2*np.pi, 100)
        ax.plot(50 * np.cos(theta_coil), 50 * np.sin(theta_coil), 'k-', linewidth=2)
        
        ax.set_xlabel('X [mm]', fontsize=11)
        ax.set_ylabel('Y [mm]', fontsize=11)
        ax.set_title('磁场矢量分布', fontsize=12)
        ax.set_aspect('equal')
        ax.set_xlim(-70, 70)
        ax.set_ylim(-70, 70)
        ax.grid(True, alpha=0.3)
        
        # 子图4: 状态信息
        ax = axes[1, 1]
        ax.axis('off')
        
        if len(P_out) > 0 and len(P_in) > 0 and np.mean(P_in) > 0:
            eta = np.mean(P_out) / np.mean(P_in) * 100
        else:
            eta = 0
        
        info_text = f"""
        能量传输状态
        ====================
        
        时间: {t_current[-1]*1e6:.1f} μs
        
        发射电流: {I1[-1]:.2f} A
        接收电流: {I2[-1]:.2f} A
        
        输入功率: {P_in[-1]:.2f} W
        输出功率: {P_out[-1]:.2f} W
        
        瞬时效率: {eta:.1f}%
        
        耦合系数: {k:.3f}
        """
        
        ax.text(0.1, 0.5, info_text, transform=ax.transAxes, fontsize=12,
               verticalalignment='center', fontfamily='monospace',
               bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
        
        plt.tight_layout()
        
        buf = io.BytesIO()
        plt.savefig(buf, format='png', dpi=100, bbox_inches='tight')
        buf.seek(0)
        frames3.append(Image.open(buf))
        plt.close()
    
    frames3[0].save('anim3_power_transfer.gif',
                   save_all=True, append_images=frames3[1:],
                   duration=150, loop=0)
    print("    anim3_power_transfer.gif 已保存")
    
    print("\n  所有GIF动画生成完成！")


def save_results(model: WPTSystemModel, params: WPTSystemParameters):
    """保存计算结果到JSON文件"""
    
    k, M, L1, L2 = model.coupling_coefficient()
    C1, C2 = model.compensation_capacitors(L1, L2)
    result = model.transmission_efficiency(k, L1, L2)
    
    results = {
        'system_parameters': {
            'operating_frequency_Hz': params.f_operating,
            'input_voltage_V': params.V_in,
            'load_resistance_Ohm': params.R_load,
            'coil_distance_mm': params.distance * 1000,
            'compensation_type': params.compensation.value
        },
        'coil_parameters': {
            'tx_turns': params.tx_coil.N,
            'tx_radius_mm': params.tx_coil.R_outer * 1000,
            'tx_inductance_uH': L1 * 1e6,
            'tx_resistance_Ohm': params.tx_coil.R,
            'rx_turns': params.rx_coil.N,
            'rx_radius_mm': params.rx_coil.R_outer * 1000,
            'rx_inductance_uH': L2 * 1e6,
            'rx_resistance_Ohm': params.rx_coil.R
        },
        'coupling_parameters': {
            'mutual_inductance_uH': M * 1e6,
            'coupling_coefficient': k
        },
        'compensation_parameters': {
            'C1_nF': C1 * 1e9,
            'C2_nF': C2 * 1e9
        },
        'performance': {
            'transmission_efficiency_percent': result['efficiency'] * 100,
            'Q_factor_tx': result['Q1'],
            'Q_factor_rx': result['Q2']
        }
    }
    
    with open('wpt_results.json', 'w', encoding='utf-8') as f:
        json.dump(results, f, indent=2, ensure_ascii=False)
    
    print("\n  结果已保存到: wpt_results.json")


if __name__ == '__main__':
    # 创建可视化
    model, params = create_visualizations()
    
    # 创建GIF动画
    create_gif_animations()
    
    # 保存结果
    save_results(model, params)
    
    print("\n" + "=" * 70)
    print("主题044仿真完成！")
    print("=" * 70)
    print("\n生成的文件:")
    print("  - fig1_magnetic_field_distribution.png: 磁场分布")
    print("  - fig2_efficiency_analysis.png: 效率特性分析")
    print("  - fig3_compensation_comparison.png: 补偿网络对比")
    print("  - fig4_coil_comparison.png: 线圈结构对比")
    print("  - anim1_field_vs_distance.gif: 磁场随距离变化")
    print("  - anim2_misalignment_effect.gif: 横向偏移影响")
    print("  - anim3_power_transfer.gif: 能量传输动画")
    print("  - wpt_results.json: 结果数据")
    print("=" * 70)