静磁场仿真-主题058_磁疗设备仿真-磁疗设备仿真

kkchenjj

24人浏览 · 2026-04-08 20:27:29

kkchenjj · 2026-04-08 20:27:29 发布

主题058：磁疗设备仿真

场景描述

磁疗是一种利用磁场作用于人体以达到治疗目的的物理治疗方法，广泛应用于疼痛管理、炎症消退、骨折愈合、神经系统疾病辅助治疗等领域。本教程将系统讲解磁疗设备的工作原理、磁场分布特性、生物效应分析和安全性评估，通过Python仿真实现静磁贴片、脉冲磁场治疗仪的设计与优化。

在这里插入图片描述

数学/物理模型

1. 永磁体磁场模型

永磁体产生的磁场可以用磁偶极子模型近似：

$B⃗(r⃗)=μ04πr3[3(m⃗⋅r^)r^−m⃗]\vec{B}(\vec{r}) = \frac{\mu_0}{4\pi r^3} [3(\vec{m} \cdot \hat{r})\hat{r} - \vec{m}]$

其中：

$B⃗\vec{B}$ ：磁感应强度 (T)
$μ0=4π×10−7\mu_0 = 4\pi \times 10^{-7}$ H/m：真空磁导率
$m⃗\vec{m}$ ：磁矩 (A·m²)
$r⃗\vec{r}$ ：距离矢量
$r$ ：距离

磁矩与永磁体参数的关系：

$\frac{B_r V}{\mu_0}$

其中：

$B_r$ ：剩磁 (T)
$V$ ：永磁体体积 (m³)

2. 圆形线圈轴向磁场

圆形线圈轴线上的磁场：

$Bz(z)=μ0NIR22(R2+z2)3/2B_z(z) = \frac{\mu_0 N I R^2}{2(R^2 + z^2)^{3/2}}$

其中：

$N$ ：线圈匝数
$I$ ：电流 (A)
$R$ ：线圈半径 (m)
$z$ ：轴向距离 (m)

3. 磁场在组织中的衰减

磁场在人体组织中的衰减遵循指数规律：

$B_0 e^{-d/\delta}$

其中：

$B_0$ ：表面磁场 (T)
$d$ ：组织深度 (m)
$δ\delta$ ：穿透深度 (m)

4. 趋肤深度

时变磁场的趋肤深度：

$δ=2ωμσ=1πfμσ\delta = \sqrt{\frac{2}{\omega \mu \sigma}} = \sqrt{\frac{1}{\pi f \mu \sigma}}$

其中：

$f$ ：频率 (Hz)
$μ=μ0μr\mu = \mu_0 \mu_r$ ：磁导率
$σ\sigma$ ：电导率 (S/m)

5. 感应电场

时变磁场产生的感应电场（法拉第定律）：

$∇×E⃗=−∂B⃗∂t\nabla \times \vec{E} = -\frac{\partial \vec{B}}{\partial t}$

对于均匀变化的磁场：

$\frac{r}{2} \frac{dB}{dt}$

环境准备

依赖库安装

pip install numpy matplotlib scipy pillow

导入库

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, Ellipse
import os

# 强制使用Agg后端
plt.switch_backend('Agg')

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

完整代码实现

案例1：静磁治疗磁场分布

def static_magnetic_therapy():
    """
    静磁治疗磁场分布
    模拟永磁体贴片在人体表面的磁场分布
    """
    # 永磁体贴片参数
    magnet_length = 0.02      # 磁贴长度 (m) = 20mm
    magnet_width = 0.01       # 磁贴宽度 (m) = 10mm
    magnet_thickness = 0.003  # 磁贴厚度 (m) = 3mm
    B_remanent = 0.3          # 剩磁 (T)
    
    # 创建计算网格
    x = np.linspace(-0.05, 0.05, 200)
    y = np.linspace(-0.03, 0.03, 120)
    X, Y = np.meshgrid(x, y)
    
    # 磁偶极子磁场计算
    def magnet_field(x_pos, y_pos, magnet_x, magnet_y, magnet_m):
        r_x = x_pos - magnet_x
        r_y = y_pos - magnet_y
        r = np.sqrt(r_x**2 + r_y**2)
        r = np.maximum(r, 1e-6)
        
        Bx = MU_0 * magnet_m / (4 * np.pi) * (3 * r_x * r_y) / r**5
        By = MU_0 * magnet_m / (4 * np.pi) * (3 * r_y**2 - r**2) / r**5
        
        return Bx, By
    
    # 计算磁场
    magnet_volume = magnet_length * magnet_width * magnet_thickness
    magnet_m = B_remanent * magnet_volume / MU_0
    
    Bx, By = magnet_field(X, Y, 0, 0, magnet_m)
    B_magnitude = np.sqrt(Bx**2 + By**2)
    
    # 深度衰减计算
    depths = np.linspace(0, 0.05, 100)
    B_at_depths = [B_remanent * np.exp(-d / 0.01) * 1000 for d in depths]
    
    # 可视化...

案例2：脉冲磁场产生与传播

def pulsed_magnetic_field():
    """
    脉冲磁场产生与传播
    模拟脉冲磁场治疗仪的工作原理
    """
    # 线圈参数
    N_turns = 100      # 匝数
    I_peak = 10        # 峰值电流 (A)
    R_coil = 0.1       # 线圈半径 (m)
    L_coil = 0.05      # 线圈长度 (m)
    
    # 脉冲参数
    pulse_frequency = 10   # 脉冲频率 (Hz)
    pulse_duration = 0.01  # 脉冲宽度 (s)
    
    # 时间
    t = np.linspace(0, 0.5, 1000)
    
    # 脉冲电流波形
    def pulse_current(t, f, duration, I_peak):
        period = 1 / f
        current = np.zeros_like(t)
        for i, time in enumerate(t):
            phase = time % period
            if phase < duration:
                current[i] = I_peak * np.sin(np.pi * phase / duration)
        return current
    
    current_waveform = pulse_current(t, pulse_frequency, pulse_duration, I_peak)
    
    # 轴向磁场计算
    def coil_field_on_axis(z, N, I, R, L):
        n = N / L
        if abs(z) < L/2:
            return MU_0 * n * I * 0.8
        else:
            return MU_0 * n * I * 0.8 * np.exp(-(abs(z) - L/2) / R)
    
    z_positions = np.linspace(-0.2, 0.2, 100)
    B_static = [coil_field_on_axis(z, N_turns, I_peak, R_coil, L_coil) * 1000 
                for z in z_positions]
    
    # 可视化...

案例3：磁场在人体组织中的穿透

def tissue_penetration():
    """
    磁场在人体组织中的穿透
    分析不同频率和强度的磁场在人体组织中的衰减
    """
    # 人体组织电磁参数
    tissues = {
        '皮肤': {'sigma': 0.1, 'epsilon_r': 1000, 'mu_r': 1.0},
        '脂肪': {'sigma': 0.04, 'epsilon_r': 20, 'mu_r': 1.0},
        '肌肉': {'sigma': 0.35, 'epsilon_r': 100, 'mu_r': 1.0},
        '骨骼': {'sigma': 0.02, 'epsilon_r': 30, 'mu_r': 1.0},
    }
    
    # 频率范围
    frequencies = np.logspace(0, 6, 100)
    
    # 趋肤深度计算
    def skin_depth(sigma, mu, f):
        if f == 0:
            return np.inf
        return np.sqrt(2 / (2 * np.pi * f * mu * sigma))
    
    # 计算各组织的趋肤深度
    skin_depths = {name: [] for name in tissues}
    
    for f in frequencies:
        for name, params in tissues.items():
            mu = MU_0 * params['mu_r']
            delta = skin_depth(params['sigma'], mu, f)
            skin_depths[name].append(delta * 1000)
    
    # 深度-场强分布
    depths = np.linspace(0, 100, 200)
    B_static = 100 * np.exp(-depths / 200)
    B_low_freq = 100 * np.exp(-depths / 150)
    B_mid_freq = 100 * np.exp(-depths / 50)
    B_high_freq = 100 * np.exp(-depths / 10)
    
    # 可视化...

案例4：治疗区域优化设计

def treatment_zone_optimization():
    """
    治疗区域优化设计
    优化磁疗设备的磁场分布以获得均匀的治疗区域
    """
    x = np.linspace(-0.1, 0.1, 200)
    y = np.linspace(-0.1, 0.1, 200)
    X, Y = np.meshgrid(x, y)
    
    # 方案1: 单线圈
    def single_coil_field(X, Y, I, N, R):
        r = np.sqrt(X**2 + Y**2)
        r = np.maximum(r, 1e-6)
        Bz = MU_0 * N * I * R**2 / (2 * (R**2 + r**2)**(3/2))
        return Bz * 1000
    
    # 方案2: 双线圈
    def double_coil_field(X, Y, I, N, R, separation):
        r1 = np.sqrt(X**2 + (Y - separation/2)**2)
        r1 = np.maximum(r1, 1e-6)
        B1 = MU_0 * N * I * R**2 / (2 * (R**2 + r1**2)**(3/2))
        
        r2 = np.sqrt(X**2 + (Y + separation/2)**2)
        r2 = np.maximum(r2, 1e-6)
        B2 = MU_0 * N * I * R**2 / (2 * (R**2 + r2**2)**(3/2))
        
        return (B1 + B2) * 1000
    
    # 方案3: 四线圈阵列
    def quad_coil_field(X, Y, I, N, R):
        positions = [(0.03, 0.03), (-0.03, 0.03), (0.03, -0.03), (-0.03, -0.03)]
        B_total = np.zeros_like(X)
        
        for pos in positions:
            r = np.sqrt((X - pos[0])**2 + (Y - pos[1])**2)
            r = np.maximum(r, 1e-6)
            B = MU_0 * N * I * R**2 / (2 * (R**2 + r**2)**(3/2))
            B_total += B
        
        return B_total * 1000
    
    # 参数
    I, N, R = 5, 50, 0.03
    separation = 0.06
    
    # 计算各方案磁场
    B_single = single_coil_field(X, Y, I, N, R)
    B_double = double_coil_field(X, Y, I, N, R, separation)
    B_quad = quad_coil_field(X, Y, I, N, R)
    
    # 计算均匀性
    center_mask = (np.abs(X) < 0.02) & (np.abs(Y) < 0.02)
    uniformity_single = np.std(B_single[center_mask]) / np.mean(B_single[center_mask]) * 100
    uniformity_double = np.std(B_double[center_mask]) / np.mean(B_double[center_mask]) * 100
    uniformity_quad = np.std(B_quad[center_mask]) / np.mean(B_quad[center_mask]) * 100
    
    # 可视化...

案例5：安全性评估

def safety_assessment():
    """
    安全性评估
    评估磁疗设备的磁场暴露安全性
    """
    # ICNIRP安全标准
    B_public_limit = 40       # mT (公众暴露)
    B_occupational_limit = 2000  # mT (职业暴露)
    
    # 频率范围
    frequencies = np.logspace(0, 6, 100)
    
    # ICNIRP参考水平
    B_ref_head = np.where(frequencies < 300, 
                          30 / (frequencies / 1000),
                          0.1)
    B_ref_limbs = B_ref_head * 3
    
    # 设备磁场衰减
    distances = np.linspace(0.01, 0.5, 100)
    
    I_device, N_device, R_device = 10, 100, 0.05
    
    B_at_distances = []
    for d in distances:
        B = MU_0 * N_device * I_device * R_device**2 / (2 * (R_device**2 + d**2)**(3/2))
        B_at_distances.append(B * 1000)
    
    # 安全距离计算
    safe_distance_public = None
    safe_distance_occupational = None
    
    for i, B in enumerate(B_at_distances):
        if safe_distance_public is None and B < B_public_limit:
            safe_distance_public = distances[i] * 1000
        if safe_distance_occupational is None and B < B_occupational_limit:
            safe_distance_occupational = distances[i] * 1000
    
    # 可视化...

代码深度解析

磁偶极子磁场计算

def magnet_field(x_pos, y_pos, magnet_x, magnet_y, magnet_m):
    """
    计算磁偶极子产生的磁场
    
    参数:
        x_pos, y_pos: 计算点坐标
        magnet_x, magnet_y: 磁偶极子位置
        magnet_m: 磁矩 (A·m²)
    
    返回:
        Bx, By: 磁场x和y分量
    
    原理:
        磁偶极子磁场公式：
        B = (μ₀/4πr³)[3(m·r̂)r̂ - m]
        
        对于z方向磁矩，简化为：
        Bx = (μ₀m/4π)(3xy)/r⁵
        By = (μ₀m/4π)(3y²-r²)/r⁵
    """
    r_x = x_pos - magnet_x
    r_y = y_pos - magnet_y
    r = np.sqrt(r_x**2 + r_y**2)
    r = np.maximum(r, 1e-6)  # 避免除零
    
    # 磁场分量
    Bx = MU_0 * magnet_m / (4 * np.pi) * (3 * r_x * r_y) / r**5
    By = MU_0 * magnet_m / (4 * np.pi) * (3 * r_y**2 - r**2) / r**5
    
    return Bx, By

关键说明：

使用磁偶极子近似适用于远场计算
近场需要更复杂的模型（如有限元方法）
添加小量避免除零错误

趋肤深度计算

def skin_depth(sigma, mu, f):
    """
    计算趋肤深度
    
    参数:
        sigma: 电导率 (S/m)
        mu: 磁导率 (H/m)
        f: 频率 (Hz)
    
    返回:
        趋肤深度 (m)
    
    原理:
        电磁波在导体中的穿透深度：
        δ = √(2/ωμσ) = √(1/πfμσ)
        
        物理意义：
        - 频率越高，趋肤深度越小
        - 电导率越高，趋肤深度越小
        - 静磁场(f=0)无趋肤效应
    """
    if f == 0:
        return np.inf
    return np.sqrt(2 / (2 * np.pi * f * mu * sigma))

物理意义：

趋肤深度表示电磁波衰减到表面值的1/e时的深度
静磁场无趋肤效应，可深入人体内部
高频磁场主要集中在组织表面

脉冲电流生成

def pulse_current(t, f, duration, I_peak):
    """
    生成正弦脉冲电流波形
    
    参数:
        t: 时间数组
        f: 脉冲频率 (Hz)
        duration: 脉冲宽度 (s)
        I_peak: 峰值电流 (A)
    
    返回:
        电流波形数组
    
    波形特征:
        - 周期性脉冲
        - 每个脉冲为正弦半波
        - 占空比 = duration × f
    """
    period = 1 / f
    current = np.zeros_like(t)
    
    for i, time in enumerate(t):
        phase = time % period  # 当前周期内的相位
        if phase < duration:
            # 正弦脉冲: sin(π × phase/duration)
            current[i] = I_peak * np.sin(np.pi * phase / duration)
    
    return current

波形特点：

使用正弦半波减少高频谐波
占空比影响平均功率和治疗效果
频率通常在1-100Hz范围内

均匀性评估

# 计算磁场均匀性
center_mask = (np.abs(X) < 0.02) & (np.abs(Y) < 0.02)
uniformity = np.std(B_field[center_mask]) / np.mean(B_field[center_mask]) * 100

评估方法：

选择中心区域进行评估
使用变异系数(CV)作为均匀性指标
CV越小，均匀性越好
医疗级设备通常要求CV < 10%

"""
磁疗设备仿真
============
本代码模拟磁疗设备的磁场分布和生物效应分析，包括：
1. 静磁治疗磁场分布
2. 脉冲磁场产生与传播
3. 磁场在人体组织中的穿透
4. 治疗区域优化设计
5. 安全性评估
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch, Arc, FancyArrowPatch, Ellipse
import matplotlib.animation as animation
from PIL import Image
import os
from scipy.special import ellipk, ellipe
from scipy.integrate import odeint

# 强制使用Agg后端，避免GUI窗口
plt.switch_backend('Agg')

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

# 创建输出目录
output_dir = os.path.dirname(os.path.abspath(__file__))

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

# ==================== 案例1: 静磁治疗磁场分布 ====================
def static_magnetic_therapy():
    """
    静磁治疗磁场分布
    模拟永磁体贴片在人体表面的磁场分布
    """
    print("=" * 60)
    print("案例1: 静磁治疗磁场分布")
    print("=" * 60)
    
    # 永磁体贴片参数
    magnet_length = 0.02  # 磁贴长度 (m) = 20mm
    magnet_width = 0.01   # 磁贴宽度 (m) = 10mm
    magnet_thickness = 0.003  # 磁贴厚度 (m) = 3mm
    B_remanent = 0.3  # 剩磁 (T) - 医用磁贴典型值
    
    # 创建计算网格
    x = np.linspace(-0.05, 0.05, 200)  # -5cm to 5cm
    y = np.linspace(-0.03, 0.03, 120)  # -3cm to 3cm
    X, Y = np.meshgrid(x, y)
    
    # 计算永磁体产生的磁场 (简化模型)
    # 使用磁偶极子近似
    def magnet_field(x_pos, y_pos, magnet_x, magnet_y, magnet_m):
        """计算单个磁偶极子的磁场"""
        r_x = x_pos - magnet_x
        r_y = y_pos - magnet_y
        r = np.sqrt(r_x**2 + r_y**2)
        r = np.maximum(r, 1e-6)
        
        # 磁场分量 (2D近似)
        Bx = MU_0 * magnet_m / (4 * np.pi) * (3 * r_x * r_y) / r**5
        By = MU_0 * magnet_m / (4 * np.pi) * (3 * r_y**2 - r**2) / r**5
        
        return Bx, By
    
    # 磁贴中心位置
    magnet_x = 0
    magnet_y = 0
    
    # 磁矩 (与体积和剩磁成正比)
    magnet_volume = magnet_length * magnet_width * magnet_thickness
    magnet_m = B_remanent * magnet_volume / MU_0
    
    # 计算磁场
    Bx, By = magnet_field(X, Y, magnet_x, magnet_y, magnet_m)
    B_magnitude = np.sqrt(Bx**2 + By**2)
    
    # 计算人体表面不同深度的磁场
    depths = np.linspace(0, 0.05, 100)  # 0 to 5cm深度
    B_at_depths = []
    
    for d in depths:
        # 简化模型：磁场随深度指数衰减
        B_surface = B_remanent
        B_depth = B_surface * np.exp(-d / 0.01)  # 衰减长度1cm
        B_at_depths.append(B_depth * 1000)  # 转换为mT
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 磁贴结构示意图
    ax1 = axes[0, 0]
    ax1.set_xlim(-0.04, 0.04)
    ax1.set_ylim(-0.02, 0.04)
    ax1.set_aspect('equal')
    ax1.axis('off')
    ax1.set_title('医用磁贴结构示意图', fontsize=12, fontweight='bold')
    
    # 绘制人体皮肤
    skin = Rectangle((-0.04, -0.005), 0.08, 0.005, 
                     facecolor='#FFDBAC', edgecolor='black', linewidth=2)
    ax1.add_patch(skin)
    ax1.text(0, -0.0025, '皮肤', ha='center', va='center', fontsize=10)
    
    # 绘制磁贴
    magnet = Rectangle((-magnet_length/2, 0), magnet_length, magnet_thickness,
                       facecolor='red', edgecolor='black', linewidth=2, alpha=0.7)
    ax1.add_patch(magnet)
    ax1.text(0, magnet_thickness/2, 'N', ha='center', va='center', 
            fontsize=12, fontweight='bold', color='white')
    
    # 绘制磁力线
    for y_offset in np.linspace(-0.008, 0.008, 5):
        ax1.annotate('', xy=(0.025, y_offset), xytext=(magnet_length/2, y_offset/2),
                    arrowprops=dict(arrowstyle='->', color='blue', lw=1, alpha=0.5))
        ax1.annotate('', xy=(-0.025, y_offset), xytext=(-magnet_length/2, y_offset/2),
                    arrowprops=dict(arrowstyle='->', color='blue', lw=1, alpha=0.5))
    
    # 标注
    ax1.annotate('磁场方向', xy=(0.03, 0.015), fontsize=9, color='blue')
    
    # 图2: 磁场分布热力图
    ax2 = axes[0, 1]
    
    # 限制显示范围
    mask = (X >= -0.04) & (X <= 0.04) & (Y >= -0.02) & (Y <= 0.03)
    B_display = np.where(mask, B_magnitude * 1000, np.nan)  # 转换为mT
    
    im = ax2.contourf(X * 1000, Y * 1000, B_display, levels=20, cmap='jet')
    plt.colorbar(im, ax=ax2, label='磁场强度 (mT)')
    
    # 绘制磁贴位置
    magnet_rect = Rectangle((-magnet_length/2 * 1000, 0), 
                           magnet_length * 1000, magnet_thickness * 1000,
                           facecolor='red', edgecolor='black', linewidth=2, alpha=0.5)
    ax2.add_patch(magnet_rect)
    
    ax2.set_xlabel('x (mm)', fontsize=11)
    ax2.set_ylabel('y (mm)', fontsize=11)
    ax2.set_title('磁场分布热力图', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    
    # 图3: 磁场随深度衰减
    ax3 = axes[1, 0]
    ax3.plot(depths * 1000, B_at_depths, 'b-', linewidth=2)
    ax3.fill_between(depths * 1000, 0, B_at_depths, alpha=0.3, color='blue')
    
    # 标记治疗区域
    therapeutic_range = (np.array(B_at_depths) >= 5) & (np.array(B_at_depths) <= 100)
    if np.any(therapeutic_range):
        therapeutic_depths = depths[therapeutic_range] * 1000
        ax3.axvspan(therapeutic_depths[0], therapeutic_depths[-1], 
                   alpha=0.2, color='green', label='有效治疗区')
    
    ax3.axhline(y=5, color='green', linestyle='--', alpha=0.5, label='最小治疗场强')
    ax3.axhline(y=100, color='red', linestyle='--', alpha=0.5, label='安全上限')
    
    ax3.set_xlabel('深度 (mm)', fontsize=11)
    ax3.set_ylabel('磁场强度 (mT)', fontsize=11)
    ax3.set_title('磁场随组织深度衰减', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4: 不同磁贴尺寸的对比
    ax4 = axes[1, 1]
    
    sizes = [10, 15, 20, 25, 30]  # mm
    B_surface_values = []
    penetration_depths = []
    
    for size in sizes:
        # 简化模型：磁场与尺寸成正比
        B_surf = B_remanent * (size / 20) * 1000  # mT
        B_surface_values.append(B_surf)
        
        # 穿透深度 (场强降至5mT)
        if B_surf > 5:
            pen_depth = -0.01 * np.log(5 / B_surf) * 1000  # mm
            penetration_depths.append(pen_depth)
        else:
            penetration_depths.append(0)
    
    ax4_twin = ax4.twinx()
    
    line1 = ax4.plot(sizes, B_surface_values, 'bo-', linewidth=2, markersize=8, label='表面场强')
    line2 = ax4_twin.plot(sizes, penetration_depths, 'rs-', linewidth=2, markersize=8, label='穿透深度')
    
    ax4.set_xlabel('磁贴尺寸 (mm)', fontsize=11)
    ax4.set_ylabel('表面磁场 (mT)', fontsize=11, color='blue')
    ax4_twin.set_ylabel('穿透深度 (mm)', fontsize=11, color='red')
    ax4.set_title('磁贴尺寸对治疗效果的影响', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    # 合并图例
    lines = line1 + line2
    labels = [l.get_label() for l in lines]
    ax4.legend(lines, labels, loc='center right')
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig1_static_magnetic.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 静磁治疗分析完成")
    print(f"  表面磁场: {B_remanent*1000:.1f} mT")
    print(f"  1cm深度场强: {B_at_depths[20]:.1f} mT")
    
    return B_magnitude, B_at_depths

# ==================== 案例2: 脉冲磁场产生与传播 ====================
def pulsed_magnetic_field():
    """
    脉冲磁场产生与传播
    模拟脉冲磁场治疗仪的工作原理
    """
    print("\n" + "=" * 60)
    print("案例2: 脉冲磁场产生与传播")
    print("=" * 60)
    
    # 线圈参数
    N_turns = 100  # 匝数
    I_peak = 10  # 峰值电流 (A)
    R_coil = 0.1  # 线圈半径 (m)
    L_coil = 0.05  # 线圈长度 (m)
    
    # 脉冲参数
    pulse_frequency = 10  # 脉冲频率 (Hz)
    pulse_duration = 0.01  # 脉冲宽度 (s)
    
    # 时间
    t = np.linspace(0, 0.5, 1000)  # 0.5秒
    
    # 脉冲电流波形
    def pulse_current(t, f, duration, I_peak):
        """生成脉冲电流波形"""
        period = 1 / f
        current = np.zeros_like(t)
        for i, time in enumerate(t):
            phase = time % period
            if phase < duration:
                # 正弦脉冲
                current[i] = I_peak * np.sin(np.pi * phase / duration)
        return current
    
    current_waveform = pulse_current(t, pulse_frequency, pulse_duration, I_peak)
    
    # 线圈中心轴向磁场
    def coil_field_on_axis(z, N, I, R, L):
        """计算螺线管轴线磁场"""
        n = N / L
        # 简化模型
        if abs(z) < L/2:
            return MU_0 * n * I * 0.8  # 内部近似均匀
        else:
            return MU_0 * n * I * 0.8 * np.exp(-(abs(z) - L/2) / R)
    
    # 不同位置的磁场
    z_positions = np.linspace(-0.2, 0.2, 100)
    B_static = [coil_field_on_axis(z, N_turns, I_peak, R_coil, L_coil) * 1000 
                for z in z_positions]  # mT
    
    # 脉冲期间的磁场变化
    t_pulse = np.linspace(0, pulse_duration, 100)
    current_pulse = I_peak * np.sin(np.pi * t_pulse / pulse_duration)
    B_pulse = [MU_0 * N_turns / L_coil * I * 1000 for I in current_pulse]
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 脉冲磁场治疗仪结构
    ax1 = axes[0, 0]
    ax1.set_xlim(-0.15, 0.15)
    ax1.set_ylim(-0.1, 0.1)
    ax1.set_aspect('equal')
    ax1.axis('off')
    ax1.set_title('脉冲磁场治疗仪结构', fontsize=12, fontweight='bold')
    
    # 绘制线圈
    coil_outer = Circle((0, 0), R_coil + 0.02, fill=False, edgecolor='blue', linewidth=3)
    coil_inner = Circle((0, 0), R_coil, fill=False, edgecolor='blue', linewidth=3)
    ax1.add_patch(coil_outer)
    ax1.add_patch(coil_inner)
    
    # 绘制线圈匝
    for angle in np.linspace(0, 2*np.pi, 12, endpoint=False):
        x_pos = R_coil * np.cos(angle)
        y_pos = R_coil * np.sin(angle)
        ax1.plot([x_pos, x_pos], [y_pos-0.01, y_pos+0.01], 'b-', linewidth=2)
    
    # 绘制治疗区域
    treatment_zone = Circle((0, 0), 0.08, fill=True, facecolor='lightblue', 
                           edgecolor='blue', linewidth=1, alpha=0.3, linestyle='--')
    ax1.add_patch(treatment_zone)
    ax1.text(0, 0, '治疗区', ha='center', va='center', fontsize=10, color='blue')
    
    # 绘制患者示意
    patient = Ellipse((0, -0.06), 0.2, 0.04, facecolor='#FFDBAC', edgecolor='black')
    ax1.add_patch(patient)
    ax1.text(0, -0.06, '患者', ha='center', va='center', fontsize=9)
    
    # 图2: 脉冲电流波形
    ax2 = axes[0, 1]
    ax2.plot(t * 1000, current_waveform, 'b-', linewidth=2)
    ax2.fill_between(t * 1000, 0, current_waveform, alpha=0.3, color='blue')
    
    # 标记脉冲
    for i in range(3):
        pulse_start = i * 1000 / pulse_frequency
        ax2.axvline(x=pulse_start, color='red', linestyle='--', alpha=0.3)
    
    ax2.set_xlabel('时间 (ms)', fontsize=11)
    ax2.set_ylabel('电流 (A)', fontsize=11)
    ax2.set_title(f'脉冲电流波形 (频率: {pulse_frequency}Hz)', fontsize=12, fontweight='bold')
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, 300)
    
    # 图3: 轴向磁场分布
    ax3 = axes[1, 0]
    ax3.plot(z_positions * 1000, B_static, 'r-', linewidth=2)
    ax3.fill_between(z_positions * 1000, 0, B_static, alpha=0.3, color='red')
    
    # 标记线圈位置
    ax3.axvspan(-L_coil/2 * 1000, L_coil/2 * 1000, alpha=0.2, color='blue', label='线圈区域')
    
    ax3.set_xlabel('轴向位置 (mm)', fontsize=11)
    ax3.set_ylabel('磁场强度 (mT)', fontsize=11)
    ax3.set_title('线圈轴向磁场分布', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4: 脉冲期间的磁场变化
    ax4 = axes[1, 1]
    ax4.plot(t_pulse * 1000, B_pulse, 'g-', linewidth=2)
    ax4.fill_between(t_pulse * 1000, 0, B_pulse, alpha=0.3, color='green')
    
    # 标记峰值
    max_B_idx = np.argmax(B_pulse)
    ax4.plot(t_pulse[max_B_idx] * 1000, B_pulse[max_B_idx], 'r*', markersize=15)
    ax4.annotate(f'峰值: {B_pulse[max_B_idx]:.1f} mT', 
                xy=(t_pulse[max_B_idx] * 1000, B_pulse[max_B_idx]),
                xytext=(t_pulse[max_B_idx] * 1000 + 1, B_pulse[max_B_idx] - 2),
                arrowprops=dict(arrowstyle='->', color='red'),
                fontsize=10, color='red')
    
    ax4.set_xlabel('时间 (ms)', fontsize=11)
    ax4.set_ylabel('磁场强度 (mT)', fontsize=11)
    ax4.set_title('单个脉冲期间的磁场变化', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig2_pulsed_field.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 脉冲磁场分析完成")
    print(f"  峰值磁场: {max(B_pulse):.1f} mT")
    print(f"  中心场强: {B_static[len(B_static)//2]:.1f} mT")
    
    return B_static, B_pulse

# ==================== 案例3: 磁场在人体组织中的穿透 ====================
def tissue_penetration():
    """
    磁场在人体组织中的穿透
    分析不同频率和强度的磁场在人体组织中的衰减
    """
    print("\n" + "=" * 60)
    print("案例3: 磁场在人体组织中的穿透")
    print("=" * 60)
    
    # 人体组织电磁参数
    tissues = {
        '皮肤': {'sigma': 0.1, 'epsilon_r': 1000, 'mu_r': 1.0, 'color': '#FFDBAC'},
        '脂肪': {'sigma': 0.04, 'epsilon_r': 20, 'mu_r': 1.0, 'color': '#FFE4B5'},
        '肌肉': {'sigma': 0.35, 'epsilon_r': 100, 'mu_r': 1.0, 'color': '#CD5C5C'},
        '骨骼': {'sigma': 0.02, 'epsilon_r': 30, 'mu_r': 1.0, 'color': '#F5F5DC'},
    }
    
    # 频率范围 (Hz)
    frequencies = np.logspace(0, 6, 100)  # 1Hz to 1MHz
    
    # 趋肤深度计算
    def skin_depth(sigma, mu, f):
        """计算趋肤深度"""
        if f == 0:
            return np.inf
        return np.sqrt(2 / (2 * np.pi * f * mu * sigma))
    
    # 计算各组织的趋肤深度
    skin_depths = {name: [] for name in tissues}
    
    for f in frequencies:
        for name, params in tissues.items():
            mu = MU_0 * params['mu_r']
            delta = skin_depth(params['sigma'], mu, f)
            skin_depths[name].append(delta * 1000)  # 转换为mm
    
    # 深度-场强分布
    depths = np.linspace(0, 100, 200)  # 0 to 100mm
    
    # 静磁场 (0Hz) - 无趋肤效应
    B_static = 100 * np.exp(-depths / 200)  # 缓慢衰减
    
    # 低频脉冲 (10Hz)
    B_low_freq = 100 * np.exp(-depths / 150)
    
    # 中频 (1kHz)
    B_mid_freq = 100 * np.exp(-depths / 50)
    
    # 高频 (100kHz)
    B_high_freq = 100 * np.exp(-depths / 10)
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: 人体组织分层模型
    ax1 = axes[0, 0]
    ax1.set_xlim(0, 0.1)
    ax1.set_ylim(0, 0.12)
    ax1.axis('off')
    ax1.set_title('人体组织分层模型', fontsize=12, fontweight='bold')
    
    # 绘制各层组织
    layer_thickness = [0.002, 0.01, 0.03, 0.05]  # m
    layer_names = ['皮肤', '脂肪', '肌肉', '骨骼']
    y_pos = 0.1
    
    for thickness, name in zip(layer_thickness, layer_names):
        rect = Rectangle((0.02, y_pos - thickness), 0.06, thickness,
                        facecolor=tissues[name]['color'], edgecolor='black', linewidth=1)
        ax1.add_patch(rect)
        ax1.text(0.05, y_pos - thickness/2, name, ha='center', va='center', fontsize=10)
        y_pos -= thickness
    
    # 绘制磁场方向
    for y in np.linspace(0.02, 0.09, 5):
        ax1.annotate('', xy=(0.09, y), xytext=(0.01, y),
                    arrowprops=dict(arrowstyle='->', color='blue', lw=1.5))
    ax1.text(0.05, 0.115, '磁场方向', ha='center', fontsize=10, color='blue')
    
    # 图2: 趋肤深度 vs 频率
    ax2 = axes[0, 1]
    
    for name, depths_list in skin_depths.items():
        ax2.loglog(frequencies, depths_list, linewidth=2, label=name, 
                  color=tissues[name]['color'], marker='o', markersize=3)
    
    ax2.set_xlabel('频率 (Hz)', fontsize=11)
    ax2.set_ylabel('趋肤深度 (mm)', fontsize=11)
    ax2.set_title('不同组织的趋肤深度', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3, which='both')
    
    # 图3: 不同频率的穿透深度
    ax3 = axes[1, 0]
    
    ax3.plot(depths, B_static, 'b-', linewidth=2, label='静磁 (0Hz)')
    ax3.plot(depths, B_low_freq, 'g-', linewidth=2, label='低频 (10Hz)')
    ax3.plot(depths, B_mid_freq, 'orange', linewidth=2, label='中频 (1kHz)')
    ax3.plot(depths, B_high_freq, 'r-', linewidth=2, label='高频 (100kHz)')
    
    # 标记治疗有效区
    ax3.axhline(y=5, color='green', linestyle='--', alpha=0.5, label='最小有效场强')
    
    ax3.set_xlabel('深度 (mm)', fontsize=11)
    ax3.set_ylabel('相对磁场强度 (%)', fontsize=11)
    ax3.set_title('不同频率磁场的穿透能力', fontsize=12, fontweight='bold')
    ax3.legend()
    ax3.grid(True, alpha=0.3)
    
    # 图4: 组织参数汇总
    ax4 = axes[1, 1]
    ax4.axis('off')
    ax4.set_title('人体组织电磁参数', fontsize=12, fontweight='bold')
    
    # 参数表格
    tissue_data = [
        ['组织', '电导率 (S/m)', '相对介电常数', '典型厚度'],
        ['皮肤', '0.1', '1000', '2mm'],
        ['脂肪', '0.04', '20', '10mm'],
        ['肌肉', '0.35', '100', '30mm'],
        ['骨骼', '0.02', '30', '50mm'],
    ]
    
    table = ax4.table(cellText=tissue_data[1:], colLabels=tissue_data[0],
                     cellLoc='center', loc='center', bbox=[0, 0.2, 1, 0.7])
    table.auto_set_font_size(False)
    table.set_fontsize(10)
    table.scale(1, 2)
    
    # 设置表头样式
    for i in range(4):
        table[(0, i)].set_facecolor('#4CAF50')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    # 设置行颜色
    for i in range(1, 5):
        for j in range(4):
            table[(i, j)].set_facecolor(tissues[tissue_data[i][0]]['color'])
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig3_tissue_penetration.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 组织穿透分析完成")
    print(f"  肌肉组织1kHz趋肤深度: {skin_depths['肌肉'][50]:.1f} mm")
    print(f"  静磁场100mm衰减至: {B_static[-1]:.1f}%")
    
    return skin_depths

# ==================== 案例4: 治疗区域优化设计 ====================
def treatment_zone_optimization():
    """
    治疗区域优化设计
    优化磁疗设备的磁场分布以获得均匀的治疗区域
    """
    print("\n" + "=" * 60)
    print("案例4: 治疗区域优化设计")
    print("=" * 60)
    
    # 创建计算网格
    x = np.linspace(-0.1, 0.1, 200)
    y = np.linspace(-0.1, 0.1, 200)
    X, Y = np.meshgrid(x, y)
    
    # 方案1: 单线圈
    def single_coil_field(X, Y, I, N, R):
        """单圆形线圈磁场 (简化模型)"""
        r = np.sqrt(X**2 + Y**2)
        r = np.maximum(r, 1e-6)
        # 轴向磁场近似
        Bz = MU_0 * N * I * R**2 / (2 * (R**2 + r**2)**(3/2))
        return Bz * 1000  # mT
    
    # 方案2: 双线圈 (Helmholtz-like)
    def double_coil_field(X, Y, I, N, R, separation):
        """双线圈磁场"""
        # 线圈1
        r1 = np.sqrt(X**2 + (Y - separation/2)**2)
        r1 = np.maximum(r1, 1e-6)
        B1 = MU_0 * N * I * R**2 / (2 * (R**2 + r1**2)**(3/2))
        
        # 线圈2
        r2 = np.sqrt(X**2 + (Y + separation/2)**2)
        r2 = np.maximum(r2, 1e-6)
        B2 = MU_0 * N * I * R**2 / (2 * (R**2 + r2**2)**(3/2))
        
        return (B1 + B2) * 1000  # mT
    
    # 方案3: 四线圈阵列
    def quad_coil_field(X, Y, I, N, R):
        """四线圈阵列磁场"""
        positions = [(0.03, 0.03), (-0.03, 0.03), (0.03, -0.03), (-0.03, -0.03)]
        B_total = np.zeros_like(X)
        
        for pos in positions:
            r = np.sqrt((X - pos[0])**2 + (Y - pos[1])**2)
            r = np.maximum(r, 1e-6)
            B = MU_0 * N * I * R**2 / (2 * (R**2 + r**2)**(3/2))
            B_total += B
        
        return B_total * 1000  # mT
    
    # 参数
    I = 5  # A
    N = 50  # 匝
    R = 0.03  # m
    separation = 0.06  # m
    
    # 计算各方案磁场
    B_single = single_coil_field(X, Y, I, N, R)
    B_double = double_coil_field(X, Y, I, N, R, separation)
    B_quad = quad_coil_field(X, Y, I, N, R)
    
    # 计算均匀性 (中心区域)
    center_mask = (np.abs(X) < 0.02) & (np.abs(Y) < 0.02)
    
    uniformity_single = np.std(B_single[center_mask]) / np.mean(B_single[center_mask]) * 100
    uniformity_double = np.std(B_double[center_mask]) / np.mean(B_double[center_mask]) * 100
    uniformity_quad = np.std(B_quad[center_mask]) / np.mean(B_quad[center_mask]) * 100
    
    fig, axes = plt.subplots(2, 3, figsize=(16, 10))
    
    # 方案1: 单线圈
    ax1 = axes[0, 0]
    im1 = ax1.contourf(X * 1000, Y * 1000, B_single, levels=20, cmap='jet')
    plt.colorbar(im1, ax=ax1, label='B (mT)')
    circle1 = Circle((0, 0), R * 1000, fill=False, edgecolor='white', linewidth=2, linestyle='--')
    ax1.add_patch(circle1)
    ax1.set_title(f'单线圈\n均匀性: {uniformity_single:.1f}%', fontsize=11, fontweight='bold')
    ax1.set_xlabel('x (mm)')
    ax1.set_ylabel('y (mm)')
    ax1.set_aspect('equal')
    
    # 方案2: 双线圈
    ax2 = axes[0, 1]
    im2 = ax2.contourf(X * 1000, Y * 1000, B_double, levels=20, cmap='jet')
    plt.colorbar(im2, ax=ax2, label='B (mT)')
    circle2a = Circle((0, separation/2 * 1000), R * 1000, fill=False, edgecolor='white', linewidth=2)
    circle2b = Circle((0, -separation/2 * 1000), R * 1000, fill=False, edgecolor='white', linewidth=2)
    ax2.add_patch(circle2a)
    ax2.add_patch(circle2b)
    ax2.set_title(f'双线圈\n均匀性: {uniformity_double:.1f}%', fontsize=11, fontweight='bold')
    ax2.set_xlabel('x (mm)')
    ax2.set_ylabel('y (mm)')
    ax2.set_aspect('equal')
    
    # 方案3: 四线圈
    ax3 = axes[0, 2]
    im3 = ax3.contourf(X * 1000, Y * 1000, B_quad, levels=20, cmap='jet')
    plt.colorbar(im3, ax=ax3, label='B (mT)')
    positions = [(30, 30), (-30, 30), (30, -30), (-30, -30)]
    for pos in positions:
        circle = Circle(pos, R * 1000, fill=False, edgecolor='white', linewidth=2)
        ax3.add_patch(circle)
    ax3.set_title(f'四线圈阵列\n均匀性: {uniformity_quad:.1f}%', fontsize=11, fontweight='bold')
    ax3.set_xlabel('x (mm)')
    ax3.set_ylabel('y (mm)')
    ax3.set_aspect('equal')
    
    # 中心线磁场分布对比
    ax4 = axes[1, 0]
    center_line = Y[:, len(x)//2] * 1000  # mm
    
    ax4.plot(center_line, B_single[:, len(x)//2], 'b-', linewidth=2, label='单线圈')
    ax4.plot(center_line, B_double[:, len(x)//2], 'r-', linewidth=2, label='双线圈')
    ax4.plot(center_line, B_quad[:, len(x)//2], 'g-', linewidth=2, label='四线圈')
    
    ax4.axhline(y=10, color='green', linestyle='--', alpha=0.5, label='治疗下限')
    ax4.axhline(y=50, color='red', linestyle='--', alpha=0.5, label='治疗上限')
    
    ax4.set_xlabel('y位置 (mm)', fontsize=11)
    ax4.set_ylabel('磁场强度 (mT)', fontsize=11)
    ax4.set_title('中心线磁场分布对比', fontsize=12, fontweight='bold')
    ax4.legend()
    ax4.grid(True, alpha=0.3)
    
    # 均匀性对比
    ax5 = axes[1, 1]
    
    schemes = ['单线圈', '双线圈', '四线圈']
    uniformities = [uniformity_single, uniformity_double, uniformity_quad]
    colors = ['blue', 'red', 'green']
    
    bars = ax5.bar(schemes, uniformities, color=colors, alpha=0.7)
    
    for bar, uni in zip(bars, uniformities):
        height = bar.get_height()
        ax5.text(bar.get_x() + bar.get_width()/2., height + 0.5,
                f'{uni:.1f}%', ha='center', va='bottom', fontsize=10)
    
    ax5.set_ylabel('不均匀度 (%)', fontsize=11)
    ax5.set_title('治疗区域均匀性对比', fontsize=12, fontweight='bold')
    ax5.grid(True, alpha=0.3, axis='y')
    
    # 优化建议
    ax6 = axes[1, 2]
    ax6.axis('off')
    ax6.set_title('线圈配置优化建议', fontsize=12, fontweight='bold')
    
    recommendations = [
        '1. 小面积精准治疗：单线圈',
        '   - 优点: 结构简单，成本低',
        '   - 缺点: 均匀性差',
        '',
        '2. 中等面积治疗：双线圈',
        '   - 优点: 均匀性较好',
        '   - 适用: 关节、腰部治疗',
        '',
        '3. 大面积均匀治疗：四线圈',
        '   - 优点: 均匀性最佳',
        '   - 适用: 背部、全身治疗',
        '',
        '优化要点:',
        '- 线圈间距 ≈ 线圈半径',
        '- 采用反向串联抵消杂散场',
        '- 增加铁芯提高场强',
    ]
    
    y_pos = 0.9
    for line in recommendations:
        if line.startswith(('1.', '2.', '3.')):
            ax6.text(0.05, y_pos, line, fontsize=10, fontweight='bold', 
                    transform=ax6.transAxes, color='blue')
        elif line.startswith('优化要点:'):
            ax6.text(0.05, y_pos, line, fontsize=10, fontweight='bold', 
                    transform=ax6.transAxes, color='red')
        else:
            ax6.text(0.05, y_pos, line, fontsize=9, transform=ax6.transAxes)
        y_pos -= 0.06
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig4_treatment_optimization.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 治疗区域优化完成")
    print(f"  单线圈均匀性: {uniformity_single:.1f}%")
    print(f"  双线圈均匀性: {uniformity_double:.1f}%")
    print(f"  四线圈均匀性: {uniformity_quad:.1f}%")
    
    return uniformity_single, uniformity_double, uniformity_quad

# ==================== 案例5: 安全性评估 ====================
def safety_assessment():
    """
    安全性评估
    评估磁疗设备的磁场暴露安全性
    """
    print("\n" + "=" * 60)
    print("案例5: 磁疗设备安全性评估")
    print("=" * 60)
    
    # 国际安全标准 (ICNIRP guidelines)
    # 静磁场暴露限值
    B_public_limit = 40  # mT (公众暴露)
    B_occupational_limit = 2000  # mT (职业暴露)
    
    # 时变磁场限值 (取决于频率)
    frequencies = np.logspace(0, 6, 100)  # 1Hz to 1MHz
    
    # ICNIRP参考水平 (简化模型)
    # 对于头部/躯干
    B_ref_head = np.where(frequencies < 300, 
                          30 / (frequencies / 1000),  # 低频: 与频率成反比
                          0.1)  # 高频: 固定值
    
    # 对于四肢
    B_ref_limbs = B_ref_head * 3  # 四肢限值是头部的3倍
    
    # 计算不同距离处的磁场衰减
    distances = np.linspace(0.01, 0.5, 100)  # 1cm to 50cm
    
    # 设备参数
    I_device = 10  # A
    N_device = 100  # 匝
    R_device = 0.05  # m
    
    # 轴向磁场
    B_at_distances = []
    for d in distances:
        B = MU_0 * N_device * I_device * R_device**2 / (2 * (R_device**2 + d**2)**(3/2))
        B_at_distances.append(B * 1000)  # mT
    
    # 安全距离计算
    safe_distance_public = None
    safe_distance_occupational = None
    
    for i, B in enumerate(B_at_distances):
        if safe_distance_public is None and B < B_public_limit:
            safe_distance_public = distances[i] * 1000  # mm
        if safe_distance_occupational is None and B < B_occupational_limit:
            safe_distance_occupational = distances[i] * 1000  # mm
    
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))
    
    # 图1: ICNIRP安全限值
    ax1 = axes[0, 0]
    
    ax1.loglog(frequencies, B_ref_head, 'b-', linewidth=2, label='头部/躯干')
    ax1.loglog(frequencies, B_ref_limbs, 'r-', linewidth=2, label='四肢')
    
    # 标记治疗常用频率
    ax1.axvline(x=10, color='green', linestyle='--', alpha=0.5, label='磁疗频率(10Hz)')
    
    ax1.set_xlabel('频率 (Hz)', fontsize=11)
    ax1.set_ylabel('磁场限值 (mT)', fontsize=11)
    ax1.set_title('ICNIRP磁场暴露参考水平', fontsize=12, fontweight='bold')
    ax1.legend()
    ax1.grid(True, alpha=0.3, which='both')
    
    # 图2: 设备磁场衰减与安全距离
    ax2 = axes[0, 1]
    
    ax2.semilogy(distances * 1000, B_at_distances, 'b-', linewidth=2)
    ax2.fill_between(distances * 1000, 0, B_at_distances, alpha=0.3, color='blue')
    
    # 标记安全限值
    ax2.axhline(y=B_public_limit, color='red', linestyle='--', linewidth=2, 
               label=f'公众限值: {B_public_limit}mT')
    ax2.axhline(y=B_occupational_limit, color='orange', linestyle='--', linewidth=2,
               label=f'职业限值: {B_occupational_limit}mT')
    
    # 标记安全距离
    if safe_distance_public:
        ax2.axvline(x=safe_distance_public, color='red', linestyle=':', alpha=0.7)
        ax2.text(safe_distance_public + 10, 100, f'公众安全距离\n{safe_distance_public:.0f}mm',
                fontsize=9, color='red')
    
    ax2.set_xlabel('距离 (mm)', fontsize=11)
    ax2.set_ylabel('磁场强度 (mT)', fontsize=11)
    ax2.set_title('设备磁场衰减与安全距离', fontsize=12, fontweight='bold')
    ax2.legend()
    ax2.grid(True, alpha=0.3)
    
    # 图3: 安全区域示意图
    ax3 = axes[1, 0]
    ax3.set_xlim(-0.3, 0.3)
    ax3.set_ylim(-0.3, 0.3)
    ax3.set_aspect('equal')
    ax3.axis('off')
    ax3.set_title('磁疗设备安全区域示意图', fontsize=12, fontweight='bold')
    
    # 绘制设备
    device = Circle((0, 0), 0.05, facecolor='gray', edgecolor='black', linewidth=2)
    ax3.add_patch(device)
    ax3.text(0, 0, '设备', ha='center', va='center', fontsize=10, color='white')
    
    # 绘制治疗区域
    treatment = Circle((0, 0), 0.08, fill=True, facecolor='lightgreen', 
                       edgecolor='green', linewidth=2, alpha=0.3)
    ax3.add_patch(treatment)
    ax3.text(0, 0.065, '治疗区', ha='center', fontsize=9, color='green')
    
    # 绘制公众安全边界
    if safe_distance_public:
        public_safe = Circle((0, 0), safe_distance_public / 1000, fill=False, 
                            edgecolor='red', linewidth=2, linestyle='--')
        ax3.add_patch(public_safe)
        ax3.text(safe_distance_public/1000 + 0.02, 0, '公众安全边界', 
                fontsize=9, color='red', rotation=90, va='center')
    
    # 绘制警告区
    warning_zone = Circle((0, 0), 0.15, fill=True, facecolor='yellow', 
                         edgecolor='orange', linewidth=1, alpha=0.2)
    ax3.add_patch(warning_zone)
    ax3.text(0.1, 0.1, '警告区\n(职业人员可进入)', fontsize=8, color='orange')
    
    # 图4: 安全评估汇总
    ax4 = axes[1, 1]
    ax4.axis('off')
    ax4.set_title('磁疗设备安全评估汇总', fontsize=12, fontweight='bold')
    
    # 评估结果
    safety_data = [
        ['评估项目', '数值', '状态'],
        ['设备表面场强', f'{B_at_distances[0]:.1f} mT', '超标'],
        ['治疗区场强', f'{B_at_distances[7]:.1f} mT', '正常'],
        ['公众安全距离', f'{safe_distance_public:.0f} mm' if safe_distance_public else 'N/A', '需标识'],
        ['职业安全距离', f'{safe_distance_occupational:.0f} mm' if safe_distance_occupational else 'N/A', '正常'],
        ['10Hz限值', '3000 mT', '符合'],
    ]
    
    table = ax4.table(cellText=safety_data[1:], colLabels=safety_data[0],
                     cellLoc='center', loc='center', bbox=[0, 0.3, 1, 0.6])
    table.auto_set_font_size(False)
    table.set_fontsize(10)
    table.scale(1, 2)
    
    # 设置表头样式
    for i in range(3):
        table[(0, i)].set_facecolor('#4CAF50')
        table[(0, i)].set_text_props(weight='bold', color='white')
    
    # 设置状态列颜色
    for i in range(1, 6):
        status = safety_data[i][2]
        if status == '超标':
            table[(i, 2)].set_facecolor('#FFCCCC')
        elif status == '正常':
            table[(i, 2)].set_facecolor('#CCFFCC')
        elif status == '需标识':
            table[(i, 2)].set_facecolor('#FFFFCC')
    
    # 安全建议
    y_pos = 0.25
    ax4.text(0.05, y_pos, '安全建议:', fontsize=11, fontweight='bold', 
            transform=ax4.transAxes, color='red')
    y_pos -= 0.05
    suggestions = [
        '1. 设备表面需加屏蔽罩',
        '2. 设置明显的安全警示标识',
        '3. 限制非治疗人员靠近',
        '4. 定期检查设备漏磁',
        '5. 操作人员需接受安全培训',
    ]
    for suggestion in suggestions:
        ax4.text(0.05, y_pos, suggestion, fontsize=9, transform=ax4.transAxes)
        y_pos -= 0.04
    
    plt.tight_layout()
    plt.savefig(os.path.join(output_dir, 'fig5_safety_assessment.png'), dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"✓ 安全性评估完成")
    print(f"  公众安全距离: {safe_distance_public:.0f} mm" if safe_distance_public else "  公众安全距离: 需计算")
    print(f"  设备表面场强: {B_at_distances[0]:.1f} mT")
    
    return safe_distance_public, B_at_distances[0]

# ==================== 主程序 ====================
if __name__ == "__main__":
    print("\n" + "=" * 70)
    print("磁疗设备仿真")
    print("=" * 70)
    
    # 运行所有案例
    static_magnetic_therapy()
    pulsed_magnetic_field()
    tissue_penetration()
    treatment_zone_optimization()
    safety_assessment()
    
    print("\n" + "=" * 70)
    print("所有仿真案例已完成！")
    print("=" * 70)
    print("\n生成的文件:")
    print("  - fig1_static_magnetic.png: 静磁治疗磁场分布")
    print("  - fig2_pulsed_field.png: 脉冲磁场产生与传播")
    print("  - fig3_tissue_penetration.png: 磁场在人体组织中的穿透")
    print("  - fig4_treatment_optimization.png: 治疗区域优化设计")
    print("  - fig5_safety_assessment.png: 安全性评估")

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

OpenClaw源码剖析之Signal模块：Signal消息反应系统如何用300行代码征服AI智能体（AI Agent)？OpenClaw流量密码Agent全解析

本文深入拆解了现象级开源AI智能体框架 OpenClaw 中的 Signal 模块，揭示其如何仅凭约300行精炼代码，便构建出一个高效、可靠的消息反应系统。文章详细剖析了该模块的核心机制：如何利用 Signal 协议的安全特性建立私密通信通道，通过事件驱动模型实时捕获并解析用户指令，并无缝对接 OpenClaw 的核心执行引擎与技能（Skills）插件体系。重点阐述了其轻量级架构如何在保障端到端加