主题048：磁纳米机器人控制仿真

场景描述

用一句话简述要解决的工程问题：

如何利用外部磁场精确控制磁纳米机器人在生物体内的运动，实现靶向给药、微创手术等生物医学应用？

---

一、引言

1.1 磁纳米机器人的概念与意义

磁纳米机器人（Magnetic Nanorobots）是指尺寸在纳米至微米量级、能够通过外部磁场进行无线操控的功能化微型装置。这一概念最早由诺贝尔物理学奖得主理查德·费曼在1959年的著名演讲《底部还有很大空间》（There’s Plenty of Room at the Bottom）中提出，他预言人类将能够操控单个原子和分子，制造出微型机器。

进入21世纪，随着纳米技术、微纳加工技术和生物医学工程的飞速发展，磁纳米机器人已经从理论设想逐步走向实际应用。特别是2000年代以来，研究人员成功开发出多种类型的磁纳米机器人，包括磁性纳米颗粒、磁性纳米螺旋、磁性纳米链等，并在靶向给药、细胞操作、微创手术、血栓清除等领域展现出巨大的应用潜力。

磁纳米机器人相比传统医疗手段具有独特优势：

1. 无创或微创操作：通过外部磁场控制，无需开刀即可到达体内深部组织
2. 精准定位：磁场可实现亚微米级的定位精度
3. 实时可控：通过调节磁场参数，可实时改变机器人的运动状态
4. 生物相容性好：磁性材料（如Fe₃O₄）具有良好的生物相容性
5. 多功能集成：可同时携带药物、造影剂、传感器等功能模块
   
   
   
   
   
   
   

1.2 磁纳米机器人的分类

根据驱动机制和运动方式，磁纳米机器人主要分为以下几类：

（1）被动式磁纳米颗粒

这类机器人通常是球形或近似球形的磁性纳米颗粒，通过磁场梯度力进行牵引。它们本身不产生主动运动，完全依赖外部磁场的梯度力进行位置控制。典型的应用是磁性靶向给药，通过外部永磁体或电磁线圈将载药纳米颗粒引导至病灶区域。

（2）螺旋形磁纳米机器人

受自然界细菌（如大肠杆菌）鞭毛运动的启发，研究人员开发了螺旋形磁纳米机器人。这类机器人具有螺旋状结构，在旋转磁场的作用下能够像螺丝一样旋转前进，产生主动推进力。螺旋结构通常采用纳米加工技术制备，材料可以是镍、钴等磁性金属或其合金。

（3）磁性纳米链

由多个磁性纳米颗粒在磁场作用下自组装形成的链状结构。这种结构具有柔性，能够适应复杂的生物环境，如血管网络。通过调节磁场方向和强度，可以控制纳米链的刚度和形状。

（4）柔性磁纳米机器人

采用柔性材料（如水凝胶）与磁性纳米颗粒复合而成，具有更好的生物相容性和变形能力。这类机器人可以设计成各种形状，如蠕虫状、鱼状等，模拟自然界生物的运动方式。

1.3 磁纳米机器人的应用领域

（1）靶向给药

这是磁纳米机器人最具前景的应用之一。传统的化疗药物在杀死癌细胞的同时也会损伤正常细胞，导致严重的副作用。磁纳米机器人可以携带抗癌药物，通过外部磁场精确引导至肿瘤部位，实现药物的定点释放，大大提高治疗效果并降低副作用。

（2）微创手术

磁纳米机器人可以作为微型手术工具，在体内执行精细操作。例如，清除血管中的血栓、修复受损组织、进行细胞级手术等。由于无需开刀，这种方法可以显著降低手术风险和患者恢复时间。

（3）细胞操作

在细胞生物学研究中，磁纳米机器人可以用于单细胞操作，如细胞分离、细胞注射、细胞内物质提取等。相比传统的光学镊子，磁操控具有更深的穿透深度和更大的作用力。

（4）医学成像增强

磁性纳米颗粒本身就是优良的磁共振成像（MRI）造影剂。通过磁纳米机器人技术，可以实现造影剂的精准递送和局部富集，提高成像对比度和分辨率。

（5）生物传感

功能化的磁纳米机器人可以作为移动传感器，在体内特定位置检测生物标志物、pH值、温度等参数，为疾病诊断提供实时信息。

1.4 本教程的学习目标

通过本教程的学习，读者将能够：

1. 理解磁纳米机器人的基本物理原理，包括磁场梯度力、磁扭矩、低雷诺数流体动力学等
2. 掌握磁纳米机器人的运动建模方法，能够建立球形和螺旋形机器人的运动方程
3. 学会使用Python进行磁纳米机器人仿真，包括磁场计算、力计算、轨迹模拟等
4. 了解磁纳米机器人在生物医学中的应用，特别是靶向输运和血管导航
5. 具备分析和优化磁纳米机器人控制策略的能力

---

二、数学/物理模型

2.1 磁场梯度力

磁纳米机器人在非均匀磁场中受到的力称为磁场梯度力，其表达式为：

$$\mathbf{F}=(\mathbf{m}\cdot \nabla )\mathbf{B}$$

其中：

• $\mathbf{F}$ 是磁场梯度力（N）
• $\mathbf{m}$ 是磁纳米机器人的磁矩（A·m²）
• $\mathbf{B}$ 是磁感应强度（T）
• $\nabla$ 是梯度算子

在直角坐标系中，该公式可以展开为：

$$F_i=\sum _{j=x,y,z}{m}_j\frac{\partial B_i}{\partial x_j},i=x,y,z$$

对于球形磁性颗粒，其磁矩可以表示为：

$$\mathbf{m}=V\cdot \mathbf{M}=\frac{4}{3}\pi R^3\cdot \mathbf{M}$$

其中 $V$ 是颗粒体积，$R$ 是半径，$\mathbf{M}$ 是磁化强度。

在低磁场下，磁化强度与磁场成正比（线性磁化）：

$$\mathbf{M}=\chi \mathbf{H}=\frac{\chi }{{\mu }_0(1+\chi )}\mathbf{B}$$

其中 $\chi$ 是磁化率，${\mu }_0$ 是真空磁导率。

2.2 磁扭矩

当磁矩方向与磁场方向不一致时，磁纳米机器人会受到磁扭矩的作用：

$$\mathbf{\tau }=\mathbf{m}\times \mathbf{B}$$

扭矩的大小为：

$$\tau =mB\sin \theta$$

其中 $\theta$ 是磁矩与磁场之间的夹角。

磁扭矩会使磁纳米机器人的磁矩趋向于与磁场方向对齐。在低雷诺数条件下，这种对齐过程非常迅速，可以认为磁矩始终与磁场方向保持一致。

2.3 低雷诺数流体动力学

磁纳米机器人在液体（如血液、水）中运动时，由于尺寸微小，其运动处于低雷诺数（Stokes流） regime。雷诺数的定义为：

$$Re=\frac{\rho vL}{\eta }$$

其中：

• $\rho$ 是流体密度（kg/m³）
• $v$ 是特征速度（m/s）
• $L$ 是特征长度（m）
• $\eta$ 是流体动力粘度（Pa·s）

对于半径为500 nm的纳米颗粒在水中以10 μm/s的速度运动：

$$Re=\frac{1000\times 10\times 10^{-6}\times 500\times 10^{-9}}{10^{-3}}\approx 5\times 10^{-6}\ll 1$$

在如此低的雷诺数下，惯性力可以完全忽略，运动由粘性力主导。根据Stokes定律，球形颗粒在流体中受到的阻力为：

$${\mathbf{F}}_d=-6\pi \eta R\mathbf{v}$$

阻力系数（drag coefficient）为：

$$\zeta =6\pi \eta R$$

在低雷诺数条件下，纳米机器人的运动方程简化为力的平衡方程：

$${\mathbf{F}}_{magnetic}+{\mathbf{F}}_d=0$$

因此，终端速度为：

$$\mathbf{v}=\frac{{\mathbf{F}}_{magnetic}}{6\pi \eta R}$$

这意味着纳米机器人的速度瞬间响应外力变化，没有惯性延迟。

2.4 布朗运动

由于尺寸微小，磁纳米机器人会受到流体分子热运动的影响，产生随机布朗运动。布朗力的大小由涨落-耗散定理决定：

$$⟨F_B^2⟩=\frac{2\zeta k_{B}T}{\Delta t}$$

其中：

• $k_B$ 是玻尔兹曼常数（1.38×10⁻²³ J/K）
• $T$ 是绝对温度（K）
• $\Delta t$ 是时间步长（s）

布朗运动是磁纳米机器人控制的挑战之一，因为它会干扰精确的位置控制。为了克服布朗运动的影响，需要施加足够强的磁场梯度力，使得磁力远大于布朗力的涨落。

2.5 螺旋推进模型

对于螺旋形磁纳米机器人，其推进机制类似于细菌的鞭毛运动。在旋转磁场的作用下，螺旋结构旋转并产生推进力。

螺旋推进的速度可以用以下简化模型描述：

$$v_{swim}=\eta \cdot f\cdot p\cdot \cos \theta$$

其中：

• $\eta$ 是推进效率系数（0 < η < 1）
• $f$ 是旋转频率（Hz）
• $p$ 是螺旋节距（m）
• $\theta$ 是螺旋角

推进效率取决于螺旋的几何形状和流体的粘度。最优的螺旋几何参数需要通过实验或详细的流体力学计算确定。

旋转磁场同步条件：为了使螺旋机器人跟随旋转磁场同步旋转，磁扭矩必须大于粘性阻力矩：

$${\tau }_{mag}>{\tau }_{drag}$$

磁扭矩：${\tau }_{mag}=m{B}_{rot}\sin \varphi$，其中 $\varphi$ 是相位差
粘性阻力矩：${\tau }_{drag}=8\pi \eta r_{helix}^3\omega$，其中 $r_{helix}$ 是螺旋半径，$\omega$ 是角速度

2.6 永磁偶极子磁场

在仿真中，外部磁场通常由永磁体或电磁线圈产生。对于永磁体，可以用磁偶极子模型近似：

$$\mathbf{B}(\mathbf{r})=\frac{{\mu }_0}{4\pi r^3}[3(\mathbf{m}\cdot \hat{\mathbf{r}})\hat{\mathbf{r}}-\mathbf{m}]$$

其中：

• $\mathbf{r}$ 是从磁体到场点的位置矢量
• $r=|\mathbf{r}|$ 是距离
• $\hat{\mathbf{r}}=\mathbf{r}/r$ 是单位矢量
• $\mathbf{m}$ 是永磁体的磁矩

磁场梯度可以通过对磁场表达式求导得到，或者使用数值方法计算。

2.7 圆形线圈磁场

对于圆形电流线圈，轴线上的磁场为：

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

其中：

• $I$ 是电流（A）
• $R$ 是线圈半径（m）
• $z$ 是轴线上的距离（m）

对于偏离轴线的点，需要使用椭圆积分计算，或者使用近似公式。

---

三、环境准备

3.1 依赖库安装

本教程需要以下Python库：

pip install numpy matplotlib scipy pillow

各库的作用：

• NumPy：科学计算基础库，用于数组操作和数值计算
• Matplotlib：数据可视化库，用于绘制静态图表和动画
• SciPy：科学计算库，提供优化、积分等高级功能（本教程中可选）
• Pillow：图像处理库，用于生成GIF动画

3.2 验证安装

运行以下Python代码验证安装：

import numpy as np
import matplotlib.pyplot as plt

print(f"NumPy版本: {np.__version__}")
print("所有依赖库安装成功！")

3.3 推荐开发环境

• Python版本：3.8及以上
• IDE推荐：VS Code、PyCharm、Jupyter Notebook
• 操作系统：Windows、Linux、macOS均可

---

四、完整代码实现

以下是完整的磁纳米机器人控制仿真代码，包含5个案例：

"""
磁纳米机器人控制仿真
===================
本代码模拟磁纳米机器人在梯度磁场中的运动行为
包括：磁场计算、磁力计算、运动轨迹仿真、可视化
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch, Ellipse
from matplotlib.collections import LineCollection
import matplotlib.animation as animation
from matplotlib.colors import LinearSegmentedColormap
import os

# 强制使用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)
ETA_WATER = 1e-3  # 水的动力粘度 (Pa·s)
KB_T = 4.11e-21  # 室温下的热能量 (J, 300K)


# ==================== 磁场计算模块 ====================
def magnetic_field_coil(position, coil_center, coil_radius, current, n_turns=1):
    """
    计算圆形线圈产生的磁场（轴对称近似）
    使用Biot-Savart定律的简化形式

    参数:
        position: 观察点位置 [x, y, z] (m)
        coil_center: 线圈中心 [x, y, z] (m)
        coil_radius: 线圈半径 (m)
        current: 电流 (A)
        n_turns: 线圈匝数

    返回:
        B: 磁感应强度 [Bx, By, Bz] (T)
    """
    r = np.array(position) - np.array(coil_center)
    rho = np.sqrt(r[0]**2 + r[1]**2)  # 到轴的距离
    z = r[2]  # 轴向距离

    # 使用椭圆积分近似计算（简化公式）
    # 基于圆形线圈磁场的标准结果
    a = coil_radius
    I = current * n_turns

    # 避免奇点
    if rho < 1e-10 and abs(z) < 1e-10:
        return np.array([0.0, 0.0, MU_0 * I / (2 * a)])

    # 简化计算（轴附近近似）
    denom = (a**2 + rho**2 + z**2)**1.5
    Bz = MU_0 * I * a**2 / (2 * denom)
    Brho = MU_0 * I * a**2 * rho * z / (2 * denom * (a**2 + z**2))

    # 转换回笛卡尔坐标
    if rho > 1e-10:
        Bx = Brho * r[0] / rho
        By = Brho * r[1] / rho
    else:
        Bx, By = 0.0, 0.0

    return np.array([Bx, By, Bz])


def gradient_field_helmholtz(position, coil_radius, current, separation):
    """
    计算Helmholtz线圈对的磁场梯度
    用于产生均匀的磁场梯度

    参数:
        position: 观察点位置 [x, y, z] (m)
        coil_radius: 线圈半径 (m)
        current: 电流 (A)
        separation: 线圈间距 (m)

    返回:
        B: 磁感应强度 [Bx, By, Bz] (T)
        grad_B: 磁场梯度张量 (T/m)
    """
    z_offset = separation / 2

    # 上线圈（电流正向）
    B1 = magnetic_field_coil(position, [0, 0, z_offset], coil_radius, current)
    # 下线圈（电流反向）
    B2 = magnetic_field_coil(position, [0, 0, -z_offset], coil_radius, -current)

    B = B1 + B2

    # 数值计算梯度
    eps = 1e-6
    grad_B = np.zeros((3, 3))
    for i in range(3):
        pos_plus = position.copy()
        pos_minus = position.copy()
        pos_plus[i] += eps
        pos_minus[i] -= eps

        B_plus = (magnetic_field_coil(pos_plus, [0, 0, z_offset], coil_radius, current) +
                  magnetic_field_coil(pos_plus, [0, 0, -z_offset], coil_radius, -current))
        B_minus = (magnetic_field_coil(pos_minus, [0, 0, z_offset], coil_radius, current) +
                   magnetic_field_coil(pos_minus, [0, 0, -z_offset], coil_radius, -current))

        grad_B[:, i] = (B_plus - B_minus) / (2 * eps)

    return B, grad_B


def field_from_permanent_magnet(position, magnet_pos, magnet_moment):
    """
    计算永磁偶极子产生的磁场

    参数:
        position: 观察点位置 [x, y, z] (m)
        magnet_pos: 磁体位置 [x, y, z] (m)
        magnet_moment: 磁矩 [mx, my, mz] (A·m²)

    返回:
        B: 磁感应强度 [Bx, By, Bz] (T)
        grad_B: 磁场梯度张量 (T/m)
    """
    r_vec = np.array(position) - np.array(magnet_pos)
    r = np.linalg.norm(r_vec)

    if r < 1e-10:
        return np.zeros(3), np.zeros((3, 3))

    m = np.array(magnet_moment)
    r_hat = r_vec / r

    # 偶极子磁场公式: B = (μ₀/4πr³)[3(m·r̂)r̂ - m]
    m_dot_r = np.dot(m, r_hat)
    B = (MU_0 / (4 * np.pi * r**3)) * (3 * m_dot_r * r_hat - m)

    # 计算梯度（数值方法）
    eps = 1e-8
    grad_B = np.zeros((3, 3))
    for i in range(3):
        pos_plus = position.copy()
        pos_minus = position.copy()
        pos_plus[i] += eps
        pos_minus[i] -= eps

        # 计算正负偏移位置的磁场
        r_vec_p = pos_plus - np.array(magnet_pos)
        r_p = np.linalg.norm(r_vec_p)
        r_hat_p = r_vec_p / r_p
        m_dot_r_p = np.dot(m, r_hat_p)
        B_plus = (MU_0 / (4 * np.pi * r_p**3)) * (3 * m_dot_r_p * r_hat_p - m)

        r_vec_m = pos_minus - np.array(magnet_pos)
        r_m = np.linalg.norm(r_vec_m)
        r_hat_m = r_vec_m / r_m
        m_dot_r_m = np.dot(m, r_hat_m)
        B_minus = (MU_0 / (4 * np.pi * r_m**3)) * (3 * m_dot_r_m * r_hat_m - m)

        # 中心差分计算梯度
        grad_B[:, i] = (B_plus - B_minus) / (2 * eps)

    return B, grad_B

代码解析：

1. 位置矢量计算：r_vec = position - magnet_pos 计算从磁体到场点的矢量
2. 距离计算：r = np.linalg.norm(r_vec) 计算欧几里得距离
3. 单位矢量：r_hat = r_vec / r 得到方向单位矢量
4. 点积计算：m_dot_r = np.dot(m, r_hat) 计算磁矩与位置方向的点积
5. 磁场计算：应用磁偶极子公式计算磁场
6. 梯度计算：使用数值差分方法计算磁场梯度张量

5.1.2 Helmholtz线圈磁场梯度

def gradient_field_helmholtz(position, coil_radius, current, separation):
    """
    计算Helmholtz线圈对的磁场梯度
    两个线圈电流方向相反，产生线性梯度场
    """
    z_offset = separation / 2

    # 上线圈（电流正向）
    B1 = magnetic_field_coil(position, [0, 0, z_offset], coil_radius, current)
    # 下线圈（电流反向）
    B2 = magnetic_field_coil(position, [0, 0, -z_offset], coil_radius, -current)

    B = B1 + B2

    # 数值计算梯度
    eps = 1e-6
    grad_B = np.zeros((3, 3))
    for i in range(3):
        pos_plus = position.copy()
        pos_minus = position.copy()
        pos_plus[i] += eps
        pos_minus[i] -= eps

        B_plus = (magnetic_field_coil(pos_plus, [0, 0, z_offset], coil_radius, current) +
                  magnetic_field_coil(pos_plus, [0, 0, -z_offset], coil_radius, -current))
        B_minus = (magnetic_field_coil(pos_minus, [0, 0, z_offset], coil_radius, current) +
                   magnetic_field_coil(pos_minus, [0, 0, -z_offset], coil_radius, -current))

        grad_B[:, i] = (B_plus - B_minus) / (2 * eps)

    return B, grad_B

关键概念：

Helmholtz线圈对（电流反向）可以在中心区域产生近似线性的磁场梯度，这对于精确控制纳米机器人的位置非常有用。

5.2 磁纳米机器人类

5.2.1 基类设计

class MagneticNanorobot:
    """
    磁纳米机器人基类
    模拟在磁场中的受力与运动
    """

    def __init__(self, position, radius, magnetic_moment, mass_density=5000):
        """
        初始化纳米机器人
        """
        self.position = np.array(position, dtype=float)
        self.radius = radius
        self.magnetic_moment_sat = magnetic_moment
        self.mass_density = mass_density

        # 计算体积和质量
        self.volume = (4/3) * np.pi * radius**3
        self.mass = mass_density * self.volume

        # 速度（初始静止）
        self.velocity = np.zeros(3)

        # 磁化方向（初始沿z轴）
        self.magnetization_dir = np.array([0.0, 0.0, 1.0])

        # 磁矩（取决于当前磁化方向）
        self.magnetic_moment = self.magnetic_moment_sat * self.magnetization_dir

        # 阻力系数（Stokes定律）
        self.drag_coefficient = 6 * np.pi * ETA_WATER * radius

设计要点：

1. 物理参数初始化：位置、半径、磁矩、密度
2. 派生参数计算：体积、质量、阻力系数
3. 状态变量：速度、磁化方向
4. 使用NumPy数组：便于矢量运算

5.2.2 磁力计算

def calculate_magnetic_force(self, B, grad_B):
    """
    计算磁场梯度力
    F = (m·∇)B
    """
    # F_i = m_j * ∂B_i/∂x_j = m · (grad_B的第i行)
    F = np.dot(grad_B, self.magnetic_moment)
    return F

这是磁场梯度力的核心计算公式。梯度张量 grad_B 是一个3×3矩阵，其中 grad_B[i, j] = ∂B_i/∂x_j。力矢量通过矩阵-矢量乘法得到。

5.2.3 磁化方向更新

def update_magnetization(self, B, dt):
    """
    更新磁化方向（磁矩与磁场对齐）
    使用指数松弛模型
    """
    B_mag = np.linalg.norm(B)
    if B_mag < 1e-10:
        return

    # 目标方向（沿磁场方向）
    target_dir = B / B_mag

    # 松弛时间（布朗弛豫，典型值~1-100 μs）
    tau_relax = 10e-6  # 10微秒

    # 指数松弛
    decay = np.exp(-dt / tau_relax)
    self.magnetization_dir = decay * self.magnetization_dir + (1 - decay) * target_dir

    # 归一化
    self.magnetization_dir /= np.linalg.norm(self.magnetization_dir)

    # 更新磁矩
    self.magnetic_moment = self.magnetic_moment_sat * self.magnetization_dir

物理意义：

磁纳米颗粒的磁矩会趋向于与外部磁场对齐。这个过程不是瞬时的，而是遵循指数松弛规律。松弛时间取决于颗粒的大小、磁性和周围介质的粘度。

5.2.4 位置更新（低雷诺数近似）

def update_position(self, F, dt):
    """
    更新位置（低雷诺数近似，惯性可忽略）
    v = F / ζ (Stokes定律)
    """
    # 低雷诺数下的终端速度
    self.velocity = F / self.drag_coefficient

    # 更新位置
    self.position += self.velocity * dt

在低雷诺数条件下，惯性力可以忽略，速度瞬间达到终端速度。这是Stokes流的典型特征。

5.2.5 单步更新

def step(self, B, grad_B, dt):
    """
    单步更新：计算力、更新磁化、更新位置
    """
    # 计算磁力
    F_magnetic = self.calculate_magnetic_force(B, grad_B)

    # 添加布朗运动随机力
    F_brownian = np.random.normal(0, np.sqrt(2 * self.drag_coefficient * KB_T / dt), 3)

    # 总力
    F_total = F_magnetic + F_brownian

    # 更新磁化方向
    self.update_magnetization(B, dt)

    # 更新位置
    self.update_position(F_total, dt)

布朗力的幅度由涨落-耗散定理决定，与时间步长的平方根成反比。

5.3 螺旋形机器人子类

class HelicalNanorobot(MagneticNanorobot):
    """
    螺旋形磁纳米机器人
    模拟旋转磁场下的游泳运动
    """

    def __init__(self, position, radius, magnetic_moment, helix_radius, helix_pitch, mass_density=5000):
        super().__init__(position, radius, magnetic_moment, mass_density)

        self.helix_radius = helix_radius
        self.helix_pitch = helix_pitch

        # 螺旋几何参数
        self.helix_length = np.sqrt((2*np.pi*helix_radius)**2 + helix_pitch**2)
        self.tan_theta = helix_pitch / (2 * np.pi * helix_radius)

        # 推进效率系数
        self.propulsion_efficiency = 0.1

螺旋形机器人继承自基类，添加了螺旋几何参数和游泳相关的方法。

---

"""
磁纳米机器人控制仿真
===================
本代码模拟磁纳米机器人在梯度磁场中的运动行为
包括：磁场计算、磁力计算、运动轨迹仿真、可视化
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch, Ellipse
from matplotlib.collections import LineCollection
import matplotlib.animation as animation
from matplotlib.colors import LinearSegmentedColormap
import os

# 强制使用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)
ETA_WATER = 1e-3  # 水的动力粘度 (Pa·s)
KB_T = 4.11e-21  # 室温下的热能量 (J, 300K)


# ==================== 磁场计算模块 ====================
def magnetic_field_coil(position, coil_center, coil_radius, current, n_turns=1):
    """
    计算圆形线圈产生的磁场（轴对称近似）
    使用Biot-Savart定律的简化形式

    参数:
        position: 观察点位置 [x, y, z] (m)
        coil_center: 线圈中心 [x, y, z] (m)
        coil_radius: 线圈半径 (m)
        current: 电流 (A)
        n_turns: 线圈匝数

    返回:
        B: 磁感应强度 [Bx, By, Bz] (T)
    """
    r = np.array(position) - np.array(coil_center)
    rho = np.sqrt(r[0]**2 + r[1]**2)  # 到轴的距离
    z = r[2]  # 轴向距离

    # 使用椭圆积分近似计算（简化公式）
    # 基于圆形线圈磁场的标准结果
    a = coil_radius
    I = current * n_turns

    # 避免奇点
    if rho < 1e-10 and abs(z) < 1e-10:
        return np.array([0.0, 0.0, MU_0 * I / (2 * a)])

    # 简化计算（轴附近近似）
    denom = (a**2 + rho**2 + z**2)**1.5
    Bz = MU_0 * I * a**2 / (2 * denom)
    Brho = MU_0 * I * a**2 * rho * z / (2 * denom * (a**2 + z**2))

    # 转换回笛卡尔坐标
    if rho > 1e-10:
        Bx = Brho * r[0] / rho
        By = Brho * r[1] / rho
    else:
        Bx, By = 0.0, 0.0

    return np.array([Bx, By, Bz])


def gradient_field_helmholtz(position, coil_radius, current, separation):
    """
    计算Helmholtz线圈对的磁场梯度
    用于产生均匀的磁场梯度

    参数:
        position: 观察点位置 [x, y, z] (m)
        coil_radius: 线圈半径 (m)
        current: 电流 (A)
        separation: 线圈间距 (m)

    返回:
        B: 磁感应强度 [Bx, By, Bz] (T)
        grad_B: 磁场梯度张量 (T/m)
    """
    z_offset = separation / 2

    # 上线圈（电流正向）
    B1 = magnetic_field_coil(position, [0, 0, z_offset], coil_radius, current)
    # 下线圈（电流反向）
    B2 = magnetic_field_coil(position, [0, 0, -z_offset], coil_radius, -current)

    B = B1 + B2

    # 数值计算梯度
    eps = 1e-6
    grad_B = np.zeros((3, 3))
    for i in range(3):
        pos_plus = position.copy()
        pos_minus = position.copy()
        pos_plus[i] += eps
        pos_minus[i] -= eps

        B_plus = (magnetic_field_coil(pos_plus, [0, 0, z_offset], coil_radius, current) +
                  magnetic_field_coil(pos_plus, [0, 0, -z_offset], coil_radius, -current))
        B_minus = (magnetic_field_coil(pos_minus, [0, 0, z_offset], coil_radius, current) +
                   magnetic_field_coil(pos_minus, [0, 0, -z_offset], coil_radius, -current))

        grad_B[:, i] = (B_plus - B_minus) / (2 * eps)

    return B, grad_B


def field_from_permanent_magnet(position, magnet_pos, magnet_moment):
    """
    计算永磁偶极子产生的磁场

    参数:
        position: 观察点位置 [x, y, z] (m)
        magnet_pos: 磁体位置 [x, y, z] (m)
        magnet_moment: 磁矩 [mx, my, mz] (A·m²)

    返回:
        B: 磁感应强度 [Bx, By, Bz] (T)
        grad_B: 磁场梯度张量 (T/m)
    """
    r_vec = np.array(position) - np.array(magnet_pos)
    r = np.linalg.norm(r_vec)

    if r < 1e-10:
        return np.zeros(3), np.zeros((3, 3))

    m = np.array(magnet_moment)
    r_hat = r_vec / r

    # 偶极子磁场公式: B = (μ₀/4πr³)[3(m·r̂)r̂ - m]
    m_dot_r = np.dot(m, r_hat)
    B = (MU_0 / (4 * np.pi * r**3)) * (3 * m_dot_r * r_hat - m)

    # 计算梯度（数值方法）
    eps = 1e-8
    grad_B = np.zeros((3, 3))
    for i in range(3):
        pos_plus = position.copy()
        pos_minus = position.copy()
        pos_plus[i] += eps
        pos_minus[i] -= eps

        r_vec_p = pos_plus - np.array(magnet_pos)
        r_p = np.linalg.norm(r_vec_p)
        r_hat_p = r_vec_p / r_p
        m_dot_r_p = np.dot(m, r_hat_p)
        B_plus = (MU_0 / (4 * np.pi * r_p**3)) * (3 * m_dot_r_p * r_hat_p - m)

        r_vec_m = pos_minus - np.array(magnet_pos)
        r_m = np.linalg.norm(r_vec_m)
        r_hat_m = r_vec_m / r_m
        m_dot_r_m = np.dot(m, r_hat_m)
        B_minus = (MU_0 / (4 * np.pi * r_m**3)) * (3 * m_dot_r_m * r_hat_m - m)

        grad_B[:, i] = (B_plus - B_minus) / (2 * eps)

    return B, grad_B


# ==================== 磁纳米机器人类 ====================
class MagneticNanorobot:
    """
    磁纳米机器人基类
    模拟在磁场中的受力与运动
    """

    def __init__(self, position, radius, magnetic_moment, mass_density=5000):
        """
        初始化纳米机器人

        参数:
            position: 初始位置 [x, y, z] (m)
            radius: 机器人半径 (m)
            magnetic_moment: 饱和磁矩 (A·m²)
            mass_density: 材料密度 (kg/m³)
        """
        self.position = np.array(position, dtype=float)
        self.radius = radius
        self.magnetic_moment_sat = magnetic_moment
        self.mass_density = mass_density

        # 计算体积和质量
        self.volume = (4/3) * np.pi * radius**3
        self.mass = mass_density * self.volume

        # 速度（初始静止）
        self.velocity = np.zeros(3)

        # 磁化方向（初始沿z轴）
        self.magnetization_dir = np.array([0.0, 0.0, 1.0])

        # 磁矩（取决于当前磁化方向）
        self.magnetic_moment = self.magnetic_moment_sat * self.magnetization_dir

        # 阻力系数（Stokes定律）
        self.drag_coefficient = 6 * np.pi * ETA_WATER * radius

    def calculate_magnetic_force(self, B, grad_B):
        """
        计算磁场梯度力
        F = (m·∇)B

        参数:
            B: 磁感应强度 (T)
            grad_B: 磁场梯度张量 (T/m)

        返回:
            F: 磁力 [Fx, Fy, Fz] (N)
        """
        # F_i = m_j * ∂B_i/∂x_j = m · (grad_B的第i行)
        F = np.dot(grad_B, self.magnetic_moment)
        return F

    def calculate_magnetic_torque(self, B):
        """
        计算磁扭矩
        τ = m × B

        参数:
            B: 磁感应强度 (T)

        返回:
            tau: 扭矩 [τx, τy, τz] (N·m)
        """
        tau = np.cross(self.magnetic_moment, B)
        return tau

    def update_magnetization(self, B, dt):
        """
        更新磁化方向（磁矩与磁场对齐）
        使用松弛模型

        参数:
            B: 磁感应强度 (T)
            dt: 时间步长 (s)
        """
        B_mag = np.linalg.norm(B)
        if B_mag < 1e-10:
            return

        # 目标方向（沿磁场方向）
        target_dir = B / B_mag

        # 松弛时间（布朗弛豫，典型值~1-100 μs）
        tau_relax = 10e-6  # 10微秒

        # 指数松弛
        decay = np.exp(-dt / tau_relax)
        self.magnetization_dir = decay * self.magnetization_dir + (1 - decay) * target_dir

        # 归一化
        self.magnetization_dir /= np.linalg.norm(self.magnetization_dir)

        # 更新磁矩
        self.magnetic_moment = self.magnetic_moment_sat * self.magnetization_dir

    def update_position(self, F, dt):
        """
        更新位置（低雷诺数近似，惯性可忽略）
        v = F / ζ (Stokes定律)

        参数:
            F: 总力 (N)
            dt: 时间步长 (s)
        """
        # 低雷诺数下的终端速度
        self.velocity = F / self.drag_coefficient

        # 更新位置
        self.position += self.velocity * dt

    def step(self, B, grad_B, dt):
        """
        单步更新

        参数:
            B: 磁感应强度 (T)
            grad_B: 磁场梯度 (T/m)
            dt: 时间步长 (s)
        """
        # 计算磁力
        F_magnetic = self.calculate_magnetic_force(B, grad_B)

        # 添加布朗运动随机力
        F_brownian = np.random.normal(0, np.sqrt(2 * self.drag_coefficient * KB_T / dt), 3)

        # 总力
        F_total = F_magnetic + F_brownian

        # 更新磁化方向
        self.update_magnetization(B, dt)

        # 更新位置
        self.update_position(F_total, dt)


class HelicalNanorobot(MagneticNanorobot):
    """
    螺旋形磁纳米机器人
    模拟旋转磁场下的游泳运动
    """

    def __init__(self, position, radius, magnetic_moment, helix_radius, helix_pitch,
                 mass_density=5000):
        """
        初始化螺旋形纳米机器人

        参数:
            position: 初始位置 [x, y, z] (m)
            radius: 机器人半径 (m)
            magnetic_moment: 饱和磁矩 (A·m²)
            helix_radius: 螺旋半径 (m)
            helix_pitch: 螺旋节距 (m)
            mass_density: 材料密度 (kg/m³)
        """
        super().__init__(position, radius, magnetic_moment, mass_density)

        self.helix_radius = helix_radius
        self.helix_pitch = helix_pitch

        # 螺旋几何参数
        self.helix_length = np.sqrt((2*np.pi*helix_radius)**2 + helix_pitch**2)
        self.tan_theta = helix_pitch / (2 * np.pi * helix_radius)  # 螺旋角正切

        # 旋转角度
        self.rotation_angle = 0.0
        self.angular_velocity = 0.0

        # 推进效率系数（简化模型）
        self.propulsion_efficiency = 0.1

    def calculate_swimming_velocity(self, B_rotating, freq):
        """
        计算旋转磁场下的游泳速度
        基于旋转磁场同步旋转原理

        参数:
            B_rotating: 旋转磁场幅值 (T)
            freq: 旋转频率 (Hz)

        返回:
            v_swim: 游泳速度 [vx, vy, vz] (m/s)
        """
        # 磁扭矩
        tau_mag = np.linalg.norm(self.calculate_magnetic_torque(B_rotating * self.magnetization_dir))

        # 旋转粘性阻力矩（简化）
        tau_drag = 8 * np.pi * ETA_WATER * self.helix_radius**3 * 2 * np.pi * freq

        # 同步条件：磁扭矩 > 阻力矩
        if tau_mag > tau_drag:
            # 游泳速度（沿螺旋轴方向）
            # 简化模型：v = η * f * p * (1 - ω_slip/ω_field)
            v_magnitude = self.propulsion_efficiency * freq * self.helix_pitch
        else:
            # 异步旋转，速度降低
            v_magnitude = self.propulsion_efficiency * freq * self.helix_pitch * (tau_mag / tau_drag)

        # 游泳方向沿磁化方向（螺旋轴方向）
        v_swim = v_magnitude * self.magnetization_dir

        return v_swim

    def step_swimming(self, B, B_rotating, freq, dt):
        """
        游泳步进更新

        参数:
            B: 直流磁场（用于导向）(T)
            B_rotating: 旋转磁场幅值 (T)
            freq: 旋转频率 (Hz)
            dt: 时间步长 (s)
        """
        # 更新磁化方向（与B对齐）
        self.update_magnetization(B, dt)

        # 计算游泳速度
        v_swim = self.calculate_swimming_velocity(B_rotating, freq)

        # 添加布朗运动
        F_brownian = np.random.normal(0, np.sqrt(2 * self.drag_coefficient * KB_T / dt), 3)
        v_brownian = F_brownian / self.drag_coefficient

        # 总速度
        self.velocity = v_swim + v_brownian

        # 更新位置
        self.position += self.velocity * dt

        # 更新旋转角度
        self.rotation_angle += 2 * np.pi * freq * dt


# ==================== 可视化函数 ====================
def plot_magnetic_field_gradient():
    """
    案例1: 磁场梯度可视化
    展示Helmholtz线圈对产生的磁场梯度
    """
    print("正在生成案例1: 磁场梯度可视化...")

    # 参数设置
    coil_radius = 0.05  # 5 cm
    current = 10.0  # 10 A
    separation = 0.06  # 6 cm

    # 创建网格
    x = np.linspace(-0.08, 0.08, 50)
    z = np.linspace(-0.08, 0.08, 50)
    X, Z = np.meshgrid(x, z)

    # 计算磁场
    B_magnitude = np.zeros_like(X)
    Bz = np.zeros_like(X)
    grad_Bz_z = np.zeros_like(X)

    for i in range(len(x)):
        for j in range(len(z)):
            pos = [X[j, i], 0, Z[j, i]]
            B, grad_B = gradient_field_helmholtz(pos, coil_radius, current, separation)
            B_magnitude[j, i] = np.linalg.norm(B)
            Bz[j, i] = B[2]
            grad_Bz_z[j, i] = grad_B[2, 2]  # dBz/dz

    # 创建图形
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))

    # 子图1: 磁场强度分布
    ax1 = axes[0, 0]
    im1 = ax1.contourf(X*100, Z*100, B_magnitude*1000, levels=30, cmap='viridis')
    ax1.set_xlabel('x (cm)', fontsize=12)
    ax1.set_ylabel('z (cm)', fontsize=12)
    ax1.set_title('磁场强度分布 (mT)', fontsize=14)
    plt.colorbar(im1, ax=ax1)

    # 绘制线圈位置
    ax1.axhline(y=separation*100/2, color='red', linestyle='--', linewidth=2, label='线圈位置')
    ax1.axhline(y=-separation*100/2, color='red', linestyle='--', linewidth=2)
    ax1.legend()

    # 子图2: Bz分量分布
    ax2 = axes[0, 1]
    im2 = ax2.contourf(X*100, Z*100, Bz*1000, levels=30, cmap='RdBu_r')
    ax2.set_xlabel('x (cm)', fontsize=12)
    ax2.set_ylabel('z (cm)', fontsize=12)
    ax2.set_title('Bz分量分布 (mT)', fontsize=14)
    plt.colorbar(im2, ax=ax2)
    ax2.axhline(y=separation*100/2, color='black', linestyle='--', linewidth=1)
    ax2.axhline(y=-separation*100/2, color='black', linestyle='--', linewidth=1)

    # 子图3: 磁场梯度分布
    ax3 = axes[1, 0]
    im3 = ax3.contourf(X*100, Z*100, grad_Bz_z, levels=30, cmap='coolwarm')
    ax3.set_xlabel('x (cm)', fontsize=12)
    ax3.set_ylabel('z (cm)', fontsize=12)
    ax3.set_title('磁场梯度 dBz/dz (T/m)', fontsize=14)
    plt.colorbar(im3, ax=ax3)
    ax3.axhline(y=separation*100/2, color='black', linestyle='--', linewidth=1)
    ax3.axhline(y=-separation*100/2, color='black', linestyle='--', linewidth=1)

    # 子图4: 中心轴线上的场分布
    ax4 = axes[1, 1]
    center_idx = len(x) // 2
    ax4.plot(Z[:, center_idx]*100, Bz[:, center_idx]*1000, 'b-', linewidth=2, label='Bz (mT)')
    ax4_twin = ax4.twinx()
    ax4_twin.plot(Z[:, center_idx]*100, grad_Bz_z[:, center_idx], 'r--', linewidth=2, label='dBz/dz (T/m)')

    ax4.set_xlabel('z (cm)', fontsize=12)
    ax4.set_ylabel('Bz (mT)', fontsize=12, color='blue')
    ax4_twin.set_ylabel('dBz/dz (T/m)', fontsize=12, color='red')
    ax4.set_title('中心轴线场分布', fontsize=14)
    ax4.axhline(y=0, color='gray', linestyle='-', alpha=0.3)
    ax4.axvline(x=separation*100/2, color='green', linestyle=':', alpha=0.5, label='线圈位置')
    ax4.axvline(x=-separation*100/2, color='green', linestyle=':', alpha=0.5)
    ax4.legend(loc='upper left')
    ax4_twin.legend(loc='upper right')

    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig1_magnetic_field_gradient.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("案例1完成: fig1_magnetic_field_gradient.png")


def simulate_nanorobot_trajectory():
    """
    案例2: 单纳米机器人运动轨迹仿真
    模拟纳米机器人在梯度磁场中的运动
    """
    print("正在生成案例2: 纳米机器人运动轨迹仿真...")

    # 创建纳米机器人
    # 参数：半径500nm，磁矩10^-15 A·m²
    robot = MagneticNanorobot(
        position=[0.03, 0, 0],  # 初始位置 (30mm, 0, 0)
        radius=500e-9,  # 500 nm
        magnetic_moment=1e-15  # 1 fA·m²
    )

    # 永磁体参数（产生梯度场）
    magnet_pos = np.array([0, 0, 0])
    magnet_moment = np.array([0, 0, 1e-3])  # 1 mA·m²

    # 仿真参数
    dt = 1e-4  # 100 μs
    n_steps = 5000
    save_interval = 50

    # 存储轨迹
    trajectory = []
    forces = []
    times = []

    for step in range(n_steps):
        # 计算磁场
        B, grad_B = field_from_permanent_magnet(robot.position, magnet_pos, magnet_moment)

        # 记录数据
        if step % save_interval == 0:
            trajectory.append(robot.position.copy())
            F = robot.calculate_magnetic_force(B, grad_B)
            forces.append(np.linalg.norm(F))
            times.append(step * dt)

        # 更新机器人状态
        robot.step(B, grad_B, dt)

    trajectory = np.array(trajectory)

    # 创建可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))

    # 子图1: 3D轨迹
    ax1 = axes[0, 0]
    scatter = ax1.scatter(trajectory[:, 0]*1000, trajectory[:, 2]*1000,
                          c=np.arange(len(trajectory)), cmap='viridis', s=10)
    ax1.plot(trajectory[:, 0]*1000, trajectory[:, 2]*1000, 'k--', alpha=0.3, linewidth=0.5)
    ax1.plot(trajectory[0, 0]*1000, trajectory[0, 2]*1000, 'go', markersize=10, label='起点')
    ax1.plot(trajectory[-1, 0]*1000, trajectory[-1, 2]*1000, 'r*', markersize=15, label='终点')
    ax1.plot(magnet_pos[0]*1000, magnet_pos[2]*1000, 'bs', markersize=12, label='磁体')
    ax1.set_xlabel('x (mm)', fontsize=12)
    ax1.set_ylabel('z (mm)', fontsize=12)
    ax1.set_title('纳米机器人运动轨迹 (xz平面)', fontsize=14)
    ax1.legend()
    ax1.set_aspect('equal')
    plt.colorbar(scatter, ax=ax1, label='时间步')

    # 子图2: 距离随时间变化
    ax2 = axes[0, 1]
    distances = np.linalg.norm(trajectory - magnet_pos, axis=1)
    ax2.plot(np.array(times)*1000, distances*1000, 'b-', linewidth=2)
    ax2.set_xlabel('时间 (ms)', fontsize=12)
    ax2.set_ylabel('到磁体距离 (mm)', fontsize=12)
    ax2.set_title('纳米机器人与磁体距离变化', fontsize=14)
    ax2.grid(True, alpha=0.3)

    # 子图3: 磁力随时间变化
    ax3 = axes[1, 0]
    ax3.plot(np.array(times)*1000, np.array(forces)*1e12, 'r-', linewidth=2)
    ax3.set_xlabel('时间 (ms)', fontsize=12)
    ax3.set_ylabel('磁力 (pN)', fontsize=12)
    ax3.set_title('磁场梯度力变化', fontsize=14)
    ax3.grid(True, alpha=0.3)

    # 子图4: 速度分析
    ax4 = axes[1, 1]
    velocities = np.linalg.norm(np.diff(trajectory, axis=0), axis=1) / (save_interval * dt)
    ax4.plot(np.array(times[1:])*1000, velocities*1e3, 'g-', linewidth=2)
    ax4.set_xlabel('时间 (ms)', fontsize=12)
    ax4.set_ylabel('速度 (mm/s)', fontsize=12)
    ax4.set_title('纳米机器人速度变化', fontsize=14)
    ax4.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig2_nanorobot_trajectory.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("案例2完成: fig2_nanorobot_trajectory.png")

    return trajectory, times


def simulate_helical_swimming():
    """
    案例3: 螺旋形磁纳米机器人游泳仿真
    模拟旋转磁场驱动的螺旋机器人游泳
    """
    print("正在生成案例3: 螺旋形磁纳米机器人游泳仿真...")

    # 创建螺旋形纳米机器人
    robot = HelicalNanorobot(
        position=[0, 0, -0.02],  # 初始位置
        radius=200e-9,  # 200 nm
        magnetic_moment=5e-16,  # 0.5 fA·m²
        helix_radius=1e-6,  # 1 μm 螺旋半径
        helix_pitch=5e-6  # 5 μm 节距
    )

    # 仿真参数
    dt = 1e-5  # 10 μs
    n_steps = 3000
    save_interval = 30

    # 旋转磁场参数
    B_dc = np.array([0, 0, 0.01])  # 10 mT 直流导向场
    B_rot_amp = 0.005  # 5 mT 旋转场幅值
    freq = 50  # 50 Hz 旋转频率

    # 存储数据
    trajectory = []
    times = []
    swimming_speeds = []

    for step in range(n_steps):
        # 生成旋转磁场（绕z轴旋转）
        angle = 2 * np.pi * freq * step * dt
        B_rot = B_rot_amp * np.array([np.cos(angle), np.sin(angle), 0])
        B_total = B_dc + B_rot

        # 记录数据
        if step % save_interval == 0:
            trajectory.append(robot.position.copy())
            times.append(step * dt)
            v_swim = robot.calculate_swimming_velocity(B_rot_amp, freq)
            swimming_speeds.append(np.linalg.norm(v_swim))

        # 更新机器人状态
        robot.step_swimming(B_total, B_rot_amp, freq, dt)

    trajectory = np.array(trajectory)

    # 创建可视化
    fig, axes = plt.subplots(2, 2, figsize=(14, 12))

    # 子图1: 3D轨迹
    ax1 = axes[0, 0]
    ax1.plot(trajectory[:, 0]*1000, trajectory[:, 2]*1000, 'b-', linewidth=2)
    ax1.plot(trajectory[0, 0]*1000, trajectory[0, 2]*1000, 'go', markersize=10, label='起点')
    ax1.plot(trajectory[-1, 0]*1000, trajectory[-1, 2]*1000, 'r*', markersize=15, label='终点')
    ax1.set_xlabel('x (mm)', fontsize=12)
    ax1.set_ylabel('z (mm)', fontsize=12)
    ax1.set_title('螺旋机器人游泳轨迹 (xz平面)', fontsize=14)
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # 子图2: 深度随时间变化
    ax2 = axes[0, 1]
    ax2.plot(np.array(times)*1000, trajectory[:, 2]*1000, 'purple', linewidth=2)
    ax2.set_xlabel('时间 (ms)', fontsize=12)
    ax2.set_ylabel('z位置 (mm)', fontsize=12)
    ax2.set_title('游泳深度变化', fontsize=14)
    ax2.grid(True, alpha=0.3)

    # 子图3: 游泳速度
    ax3 = axes[1, 0]
    ax3.plot(np.array(times)*1000, np.array(swimming_speeds)*1e6, 'orange', linewidth=2)
    ax3.set_xlabel('时间 (ms)', fontsize=12)
    ax3.set_ylabel('游泳速度 (μm/s)', fontsize=12)
    ax3.set_title('游泳速度变化', fontsize=14)
    ax3.grid(True, alpha=0.3)

    # 子图4: 轨迹3D示意
    ax4 = axes[1, 1]
    # 绘制螺旋示意图
    t_helix = np.linspace(0, 4*np.pi, 200)
    x_helix = robot.helix_radius * 1e6 * np.cos(t_helix)
    z_helix = robot.helix_pitch * 1e6 * t_helix / (2*np.pi)
    ax4.plot(x_helix, z_helix, 'b-', linewidth=2, label='螺旋结构')
    ax4.fill_between(x_helix, z_helix - robot.radius*1e6, z_helix + robot.radius*1e6,
                     alpha=0.3, color='blue')
    ax4.set_xlabel('x (μm)', fontsize=12)
    ax4.set_ylabel('z (μm)', fontsize=12)
    ax4.set_title('螺旋机器人结构示意图', fontsize=14)
    ax4.legend()
    ax4.set_aspect('equal')
    ax4.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig3_helical_swimming.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("案例3完成: fig3_helical_swimming.png")

    return trajectory, times


def simulate_multi_robot_control():
    """
    案例4: 多机器人协同控制仿真
    模拟多个纳米机器人在同一磁场下的运动
    """
    print("正在生成案例4: 多机器人协同控制仿真...")

    # 创建多个纳米机器人
    n_robots = 5
    robots = []
    colors = plt.cm.tab10(np.linspace(0, 1, n_robots))

    for i in range(n_robots):
        angle = 2 * np.pi * i / n_robots
        r_init = 0.025  # 25 mm
        pos = [r_init * np.cos(angle), r_init * np.sin(angle), 0]
        robot = MagneticNanorobot(
            position=pos,
            radius=400e-9 + i * 100e-9,  # 不同大小
            magnetic_moment=(0.8 + i * 0.1) * 1e-15
        )
        robots.append(robot)

    # 永磁体参数
    magnet_pos = np.array([0, 0, 0])
    magnet_moment = np.array([0, 0, 2e-3])

    # 仿真参数
    dt = 5e-5
    n_steps = 3000
    save_interval = 30

    # 存储所有机器人的轨迹
    all_trajectories = [[] for _ in range(n_robots)]
    times = []

    for step in range(n_steps):
        for i, robot in enumerate(robots):
            B, grad_B = field_from_permanent_magnet(robot.position, magnet_pos, magnet_moment)

            if step % save_interval == 0:
                all_trajectories[i].append(robot.position.copy())

            robot.step(B, grad_B, dt)

        if step % save_interval == 0:
            times.append(step * dt)

    # 创建可视化
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))

    # 子图1: 所有机器人的轨迹
    ax1 = axes[0]
    for i in range(n_robots):
        traj = np.array(all_trajectories[i])
        ax1.plot(traj[:, 0]*1000, traj[:, 1]*1000, color=colors[i],
                linewidth=2, label=f'机器人 {i+1}')
        ax1.plot(traj[0, 0]*1000, traj[0, 1]*1000, 'o', color=colors[i], markersize=8)
        ax1.plot(traj[-1, 0]*1000, traj[-1, 1]*1000, '*', color=colors[i], markersize=12)

    ax1.plot(magnet_pos[0]*1000, magnet_pos[1]*1000, 'ks', markersize=15, label='磁体')
    ax1.set_xlabel('x (mm)', fontsize=12)
    ax1.set_ylabel('y (mm)', fontsize=12)
    ax1.set_title('多机器人协同运动轨迹 (xy平面)', fontsize=14)
    ax1.legend(loc='upper right', fontsize=9)
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)

    # 子图2: 距离磁体的距离随时间变化
    ax2 = axes[1]
    for i in range(n_robots):
        traj = np.array(all_trajectories[i])
        distances = np.linalg.norm(traj - magnet_pos, axis=1)
        ax2.plot(np.array(times)*1000, distances*1000, color=colors[i],
                linewidth=2, label=f'机器人 {i+1}')

    ax2.set_xlabel('时间 (ms)', fontsize=12)
    ax2.set_ylabel('到磁体距离 (mm)', fontsize=12)
    ax2.set_title('各机器人与磁体距离变化', fontsize=14)
    ax2.legend(loc='upper right', fontsize=9)
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig4_multi_robot_control.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("案例4完成: fig4_multi_robot_control.png")


def simulate_vessel_navigation():
    """
    案例5: 血管环境中的靶向输运仿真
    模拟纳米机器人在血管网络中的导航
    """
    print("正在生成案例5: 血管环境中的靶向输运仿真...")

    # 定义血管网络（简化Y型分叉）
    def is_in_vessel(pos, vessel_radius=0.005):
        """检查点是否在血管内"""
        x, y, z = pos

        # 主血管（沿z轴）
        if z > 0 and np.sqrt(x**2 + y**2) < vessel_radius:
            return True

        # 分支1（沿x正方向，45度）
        if z <= 0:
            # 转换到分支坐标系
            x_branch = (x - z) / np.sqrt(2)
            y_branch = y
            if x > 0 and np.sqrt(x_branch**2 + y_branch**2) < vessel_radius:
                return True

        # 分支2（沿x负方向，45度）
        if z <= 0:
            x_branch = (-x - z) / np.sqrt(2)
            y_branch = y
            if x < 0 and np.sqrt(x_branch**2 + y_branch**2) < vessel_radius:
                return True

        return False

    # 创建纳米机器人
    robot = MagneticNanorobot(
        position=[0, 0, 0.03],  # 从主血管上端开始
        radius=300e-9,
        magnetic_moment=1e-15
    )

    # 移动磁体（引导机器人）
    magnet_positions = []
    magnet_moment = np.array([0, 0, 1.5e-3])

    # 仿真参数
    dt = 5e-5
    n_steps = 4000
    save_interval = 40

    trajectory = []
    magnet_traj = []
    times = []

    for step in range(n_steps):
        # 磁体沿预定路径移动（从主血管到分支1）
        if step < 1500:
            # 主血管内下降
            t = step / 1500
            magnet_pos = np.array([0, 0, 0.03 - 0.04 * t])
        elif step < 3000:
            # 进入分支1
            t = (step - 1500) / 1500
            magnet_pos = np.array([0.02 * t, 0, -0.01 - 0.02 * t])
        else:
            # 继续沿分支前进
            t = (step - 3000) / 1000
            magnet_pos = np.array([0.02 + 0.015 * t, 0, -0.03 - 0.015 * t])

        # 计算磁场
        B, grad_B = field_from_permanent_magnet(robot.position, magnet_pos, magnet_moment)

        # 血管壁约束（简化：如果超出血管，添加恢复力）
        if not is_in_vessel(robot.position):
            # 计算到血管中心的距离
            x, y, z = robot.position
            if z > 0:
                center = np.array([0, 0, z])
            elif x > 0:
                center = np.array([(z + 0.01) + 0.01, 0, z])
            else:
                center = np.array([-(z + 0.01) - 0.01, 0, z])

            to_center = center - robot.position
            dist = np.linalg.norm(to_center[:2])  # xy平面距离
            if dist > 0:
                F_wall = 1e-10 * to_center / dist  # 恢复力
                robot.velocity += F_wall / robot.drag_coefficient * dt

        # 记录数据
        if step % save_interval == 0:
            trajectory.append(robot.position.copy())
            magnet_traj.append(magnet_pos.copy())
            times.append(step * dt)

        # 更新机器人
        robot.step(B, grad_B, dt)

    trajectory = np.array(trajectory)
    magnet_traj = np.array(magnet_traj)

    # 创建可视化
    fig, axes = plt.subplots(1, 2, figsize=(16, 7))

    # 子图1: 3D轨迹投影（xz平面）
    ax1 = axes[0]

    # 绘制血管轮廓
    z_main = np.linspace(0, 0.03, 50)
    x_main_upper = 0.005 * np.ones_like(z_main)
    x_main_lower = -0.005 * np.ones_like(z_main)
    ax1.fill_betweenx(z_main*1000, x_main_lower*1000, x_main_upper*1000, alpha=0.2, color='red', label='血管')

    z_branch = np.linspace(-0.03, 0, 50)
    x_branch1 = 0.005 + 0.03 * (0 - z_branch) / 0.03
    x_branch2 = -0.005 - 0.03 * (0 - z_branch) / 0.03
    ax1.fill_betweenx(z_branch*1000, (0.005 + z_branch)*1000, (0.03 - z_branch)*1000, alpha=0.2, color='red')
    ax1.fill_betweenx(z_branch*1000, (-0.005 + z_branch)*1000, (-0.03 + z_branch)*1000, alpha=0.2, color='red')

    # 绘制轨迹
    ax1.plot(trajectory[:, 0]*1000, trajectory[:, 2]*1000, 'b-', linewidth=2, label='机器人轨迹')
    ax1.plot(trajectory[0, 0]*1000, trajectory[0, 2]*1000, 'go', markersize=10, label='起点')
    ax1.plot(trajectory[-1, 0]*1000, trajectory[-1, 2]*1000, 'r*', markersize=15, label='终点')

    # 绘制磁体轨迹
    ax1.plot(magnet_traj[:, 0]*1000, magnet_traj[:, 2]*1000, 'm--', linewidth=1.5, alpha=0.7, label='磁体轨迹')
    ax1.plot(magnet_traj[-1, 0]*1000, magnet_traj[-1, 2]*1000, 'ms', markersize=10, label='磁体')

    ax1.set_xlabel('x (mm)', fontsize=12)
    ax1.set_ylabel('z (mm)', fontsize=12)
    ax1.set_title('血管网络中的靶向输运 (xz平面)', fontsize=14)
    ax1.legend(loc='upper left', fontsize=9)
    ax1.set_aspect('equal')
    ax1.grid(True, alpha=0.3)

    # 子图2: 距离分析
    ax2 = axes[1]
    dist_to_magnet = np.linalg.norm(trajectory - magnet_traj, axis=1)
    ax2.plot(np.array(times)*1000, dist_to_magnet*1000, 'g-', linewidth=2)
    ax2.set_xlabel('时间 (ms)', fontsize=12)
    ax2.set_ylabel('机器人-磁体距离 (mm)', fontsize=12)
    ax2.set_title('机器人与引导磁体的距离', fontsize=14)
    ax2.grid(True, alpha=0.3)

    # 标记阶段
    ax2.axvline(x=75, color='gray', linestyle='--', alpha=0.5, label='进入分叉')
    ax2.axvline(x=150, color='orange', linestyle='--', alpha=0.5, label='进入分支')
    ax2.legend()

    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig5_vessel_navigation.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("案例5完成: fig5_vessel_navigation.png")


def create_animation():
    """
    创建动画展示纳米机器人运动
    """
    print("正在生成动画...")

    # 创建简单的螺旋游泳动画帧
    robot = HelicalNanorobot(
        position=[0, 0, 0],
        radius=200e-9,
        magnetic_moment=5e-16,
        helix_radius=1e-6,
        helix_pitch=5e-6
    )

    # 生成旋转磁场下的运动
    dt = 1e-5
    n_frames = 100
    positions = []

    B_dc = np.array([0, 0, 0.01])
    B_rot_amp = 0.005
    freq = 50

    for frame in range(n_frames):
        for _ in range(30):  # 每帧30步
            angle = 2 * np.pi * freq * frame * 30 * dt
            B_rot = B_rot_amp * np.array([np.cos(angle), np.sin(angle), 0])
            B_total = B_dc + B_rot
            robot.step_swimming(B_total, B_rot_amp, freq, dt)
        positions.append(robot.position.copy())

    positions = np.array(positions)

    # 创建静态图代替动画（为了兼容性）
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))

    # 轨迹图
    ax1 = axes[0]
    scatter = ax1.scatter(positions[:, 0]*1000, positions[:, 2]*1000,
                          c=np.arange(len(positions)), cmap='plasma', s=30)
    ax1.plot(positions[:, 0]*1000, positions[:, 2]*1000, 'k--', alpha=0.3, linewidth=1)
    ax1.set_xlabel('x (mm)', fontsize=12)
    ax1.set_ylabel('z (mm)', fontsize=12)
    ax1.set_title('螺旋机器人游泳轨迹 (100帧)', fontsize=14)
    plt.colorbar(scatter, ax=ax1, label='帧编号')
    ax1.grid(True, alpha=0.3)

    # 深度变化
    ax2 = axes[1]
    ax2.plot(np.arange(len(positions)), positions[:, 2]*1000, 'purple', linewidth=2)
    ax2.set_xlabel('帧编号', fontsize=12)
    ax2.set_ylabel('z位置 (mm)', fontsize=12)
    ax2.set_title('游泳深度随时间变化', fontsize=14)
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig(f'{output_dir}/fig6_animation_frames.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("动画帧完成: fig6_animation_frames.png")


def create_summary_figure():
    """
    创建综合总结图
    """
    print("正在生成综合总结图...")

    fig = plt.figure(figsize=(16, 10))

    # 创建网格布局
    gs = fig.add_gridspec(3, 3, hspace=0.3, wspace=0.3)

    # 子图1: 磁力公式
    ax1 = fig.add_subplot(gs[0, 0])
    ax1.axis('off')
    ax1.text(0.5, 0.5, r'$\mathbf{F} = (\mathbf{m} \cdot \nabla)\mathbf{B}$' + '\n\n' +
             r'$\boldsymbol{\tau} = \mathbf{m} \times \mathbf{B}$',
             fontsize=16, ha='center', va='center',
             bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.5))
    ax1.set_title('磁力与扭矩公式', fontsize=12)

    # 子图2: 阻力公式
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.axis('off')
    ax2.text(0.5, 0.5, r'$\mathbf{F}_d = 6\pi\eta R \mathbf{v}$' + '\n\n' +
             r'$v_{terminal} = \frac{F}{6\pi\eta R}$',
             fontsize=16, ha='center', va='center',
             bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))
    ax2.set_title('Stokes阻力', fontsize=12)

    # 子图3: 游泳速度
    ax3 = fig.add_subplot(gs[0, 2])
    ax3.axis('off')
    ax3.text(0.5, 0.5, r'$v_{swim} = \eta \cdot f \cdot p$' + '\n\n' +
             r'$\eta$: 效率系数' + '\n' +
             r'$f$: 旋转频率' + '\n' +
             r'$p$: 螺旋节距',
             fontsize=14, ha='center', va='center',
             bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))
    ax3.set_title('螺旋游泳速度', fontsize=12)

    # 子图4: 参数对比表
    ax4 = fig.add_subplot(gs[1, :])
    ax4.axis('off')

    table_data = [
        ['参数', '球形机器人', '螺旋机器人', '典型值'],
        ['半径', '200-500 nm', '100-300 nm', '300 nm'],
        ['磁矩', '0.5-2 fA·m²', '0.2-1 fA·m²', '1 fA·m²'],
        ['游泳速度', '0 (被动)', '1-100 μm/s', '10 μm/s'],
        ['驱动力', '磁场梯度', '旋转磁场', '梯度+旋转'],
        ['应用场景', '靶向给药', '微创手术', '细胞操作']
    ]

    table = ax4.table(cellText=table_data, cellLoc='center', loc='center',
                      colWidths=[0.2, 0.25, 0.25, 0.2])
    table.auto_set_font_size(False)
    table.set_fontsize(11)
    table.scale(1, 2)

    # 设置表头样式
    for i in range(4):
        table[(0, i)].set_facecolor('#4472C4')
        table[(0, i)].set_text_props(weight='bold', color='white')

    ax4.set_title('磁纳米机器人参数对比', fontsize=14, pad=20)

    # 子图5: 应用示意图
    ax5 = fig.add_subplot(gs[2, :2])
    ax5.set_xlim(-1, 1)
    ax5.set_ylim(-0.5, 0.5)
    ax5.axis('off')

    # 绘制血管
    from matplotlib.patches import FancyBboxPatch
    vessel = FancyBboxPatch((-0.9, -0.1), 1.8, 0.2, boxstyle="round,pad=0.02",
                            facecolor='pink', edgecolor='red', alpha=0.5)
    ax5.add_patch(vessel)

    # 绘制机器人
    for i in range(5):
        x = -0.7 + i * 0.35
        circle = Circle((x, 0), 0.03, facecolor='blue', edgecolor='navy')
        ax5.add_patch(circle)
        ax5.arrow(x, 0.05, 0.1, 0, head_width=0.02, head_length=0.02, fc='green', ec='green')

    ax5.text(0, 0.3, '靶向给药应用示意图', fontsize=14, ha='center')
    ax5.text(0, -0.3, '磁纳米机器人在血管中向病灶区域输运药物', fontsize=11, ha='center', style='italic')

    # 子图6: 关键结论
    ax6 = fig.add_subplot(gs[2, 2])
    ax6.axis('off')
    conclusions = [
        '1. 磁场梯度力可实现',
        '   精确位置控制',
        '',
        '2. 旋转磁场驱动螺旋',
        '   机器人主动游泳',
        '',
        '3. 多机器人可同时',
        '   在统一场中控制',
        '',
        '4. 血管导航需要',
        '   实时反馈控制'
    ]
    ax6.text(0.1, 0.95, '\n'.join(conclusions), fontsize=11, va='top',
             family='monospace',
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    ax6.set_title('关键结论', fontsize=12)

    plt.savefig(f'{output_dir}/fig7_summary.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("综合总结图完成: fig7_summary.png")


# ==================== 主程序 ====================
def main():
    """主函数：运行所有仿真案例"""
    print("=" * 60)
    print("磁纳米机器人控制仿真")
    print("=" * 60)
    print(f"输出目录: {output_dir}")
    print()

    # 运行所有案例
    plot_magnetic_field_gradient()
    print()

    simulate_nanorobot_trajectory()
    print()

    simulate_helical_swimming()
    print()

    simulate_multi_robot_control()
    print()

    simulate_vessel_navigation()
    print()

    create_animation()
    print()

    create_summary_figure()
    print()

    print("=" * 60)
    print("所有仿真完成！")
    print("=" * 60)


if __name__ == "__main__":
    main()