Basilisk 居中偶极子磁场模型与 IGRF 系数的数学推导

说明：本文推导 Basilisk 中居中偶极子磁场模块 MagneticFieldCenteredDipole 如何利用 IGRF 的三个一阶高斯系数 $g_1^0,g_1^1,h_1^1$ 以及行星半径 $R$ 计算航天器位置处的磁场矢量。

0. 磁偶极子物理基础

在深入推导数学公式之前，我们需要理解磁偶极子的物理本质和在地球物理学中的应用。这一章将从基本概念出发，建立对磁偶极子场的直观理解。

0.1 磁偶极子的基本定义

0.1.1 什么是磁偶极子

磁偶极子（magnetic dipole）是电磁学中的基本概念，表示一个具有两个磁极（南极和北极）的磁源。在数学上，它可以看作一个极小的闭合电流环或一对无限接近的磁单极子对。

形象地说，一个小磁铁（如指南针）、一个小型线圈通电后都可以近似为磁偶极子。它的特征是：

• 磁场线从北极出发，环绕后进入南极，形成闭合回路；
• 在远处观察（距离远大于偶极子尺寸），磁场分布由一个矢量——磁偶极矩 $\mathbf{m}$——完全描述；
• 磁偶极子没有孤立的磁荷（磁单极子至今未被观测到），磁场线总是闭合的。

0.1.2 与电偶极子的比较

特性	电偶极子	磁偶极子
构成	正负电荷对 $+q$ 和 $-q$，间距 $\mathbf{d}$	一对磁极（N极和S极）或小电流环
偶极矩定义	$\mathbf{p}=q\mathbf{d}$	$\mathbf{m}=I\mathbf{A}$（电流 $I$ × 面积矢量 $\mathbf{A}$）
场的势	电势 $\varphi =\frac{1}{4\pi {ϵ}_0}\frac{\mathbf{p}\cdot \mathbf{r}}{r^3}$	磁标势 ${\varphi }_m=\frac{{\mu }_0}{4\pi }\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}$
场的形式	$\mathbf{E}=-\nabla \varphi$	$\mathbf{B}=-\nabla {\varphi }_m$（静磁近似）
衰减规律	$E\propto 1/r^3$	$B\propto 1/r^3$
场线	从正电荷发散到负电荷（开放）	闭合环绕（无磁荷）

相似之处：

• 数学形式完全类比：都是 $1/r^3$ 衰减的矢量场；
• 都有"$3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}-\mathbf{m}$"的标准形式；
• 都可以通过标量势 $\varphi$ 来描述（静场情况）。

区别：

• 物理源不同：电偶极子来自电荷对，磁偶极子来自电流或内禀磁矩（如电子自旋）；
• 场线拓扑不同：电场线开放，磁场线闭合（反映了"没有磁单极子"这一基本事实）；
• 在介质中的行为：电偶极子受电场力矩，磁偶极子受磁场力矩。

0.2 地球磁场的磁偶极子来源

0.2.1 地磁场的观测特征

地球磁场在地表及近地空间具有以下特征：

• 主磁场强度：赤道约 $30,000\text{nT}$，极区约 $60,000\text{nT}$；
• 倾角（inclination）：在赤道为 $0°$（水平），在极区为 $\pm 90°$（垂直）；
• 偏角（declination）：磁北极与地理北极不重合，且随位置变化；
• 长期变化（secular variation）：磁极位置每年移动数十公里，场强缓慢变化。

这些特征与一个倾斜的磁偶极子高度吻合：约 90% 的地表磁场可以用一个居中偶极子解释。

0.2.2 地核发电机理论（简要）

地球磁场的物理来源是地核发电机（geodynamo）过程：

1. 地球内部结构：

   • 外核（约 2900–5100 km 深）：液态铁镍合金，温度约 4000–6000 K；
   • 内核：固态，但外核流体可以自由对流。

2. 发电机机制：

   • 热对流 + 地球自转驱动外核中的导电流体运动；
   • 流体运动切割已存在的磁场线，通过电磁感应产生电流（法拉第定律）；
   • 这些电流反过来增强磁场——形成自激发电机（self-sustaining dynamo）。

3. 为什么是偶极子主导：

   • 地核的球对称几何和科里奥利力（地球自转效应）倾向于产生轴对称的磁场模式；
   • 在球谐展开中，一阶项（偶极子）能量最低、最稳定，因此占主导；
   • 高阶项（四极子、八极子等）也存在，但能量较小（~10%）。

4. 偶极子的倾斜：

   • 理想情况下偶极轴应与自转轴重合，但实际地核流体运动的复杂性导致偶极轴相对地理轴倾斜约 11°；
   • 这个倾角对应于 IGRF 系数中 $g_1^1$ 和 $h_1^1$ 的非零值。

关键结论：地球磁场本质上是一个由地核流体运动产生的、以倾斜偶极子为主的复杂磁场，IGRF/WMM 用球谐展开来精确描述，而一阶项恰好对应这个主导的偶极子分量。

0.3 磁偶极矩矢量的数学描述

0.3.1 磁偶极矩的定义

对于一个小型闭合电流环（面积 $A$，电流 $I$），其磁偶极矩定义为：

$$\mathbf{m}=I\mathbf{A}$$

其中 $\mathbf{A}$ 是面积矢量，方向由右手定则确定（四指沿电流方向，拇指指向 $\mathbf{A}$）。

• 单位：SI 制中为 $\text{A}\cdot {\text{m}}^2$（安培·米²）或等效的 $\text{J/T}$（焦耳/特斯拉）；
• 在地磁学中，常用 $\text{nT}$ 表示磁场，因此偶极矩的"有效单位"是 $\text{T}$（结合参考半径归一化）。

0.3.2 物理意义

• $|\mathbf{m}|$：偶极子的强度，越大磁场越强；
• $\text{方向}(\mathbf{m})$：从南极指向北极（磁场线在外部从北极出发）；
• $\mathbf{m}\cdot \hat{\mathbf{r}}$：偶极矩在径向方向的投影，决定该方向上的磁场贡献。

0.3.3 在 IGRF/Basilisk 中的表示

在地球物理学中，偶极矩通过球谐系数表示。对于 IGRF 一阶项：

$${\mathbf{m}}_P=\left[\begin{array}{c}{g}_1^1\\ h_1^1\\ g_1^0\end{array}\right]\text{（在行星固定坐标系 }\mathcal{P}\text{ 中）}$$

• $g_1^0$：沿自转轴（$z$ 轴）的分量；
• $g_1^1,h_1^1$：赤道面内（$x,y$ 方向）的分量，体现偶极轴的倾斜。

对于地球（IGRF 2020）：

$${\mathbf{m}}_{\text{Earth}}=\left[\begin{array}{c}-2318\\ 5817\\ -30926\end{array}\right]\text{nT}=\left[\begin{array}{c}-2.318\\ 5.817\\ -30.926\end{array}\right]\times 10^{-6}\text{T}$$

• 负号：$g_1^0<0$ 表示磁偶极子的北极在地理南半球（实际上磁南极在地理北极附近，这是历史命名习惯）；
• 归一化：这些系数已经包含了 ${\mu }_0/(4\pi )$ 和参考半径 $R$ 的归一化，可以直接用于计算磁场。

0.4 磁偶极子场的空间分布特征

0.4.1 磁场线的形状

磁偶极子的磁场线具有标志性的"苹果形"或"葫芦形"分布：

        北极 (N)
          ↑
       ╱  |  ╲
      ╱   |   ╲     ← 磁场线环绕
     |    |    |
  ---+----+----+---  ← 赤道面
     |    |    |
      ╲   |   ╱
       ╲  |  ╱
          ↓
        南极 (S)

• 极区：磁场线近乎垂直，从北极发散，向南极汇聚；
• 赤道：磁场线近乎水平，环绕赤道；
• 闭合性：所有磁场线都是闭合的（从北极出发，绕到南极，再通过偶极子内部回到北极）。

0.4.2 强度的 $1/r^3$ 衰减规律

磁偶极子场的强度随距离 $r$ 的变化为：

$$B(r)\propto \frac{1}{r^3}$$

这比单极子场（如点电荷产生的电场 $E\propto 1/r^2$）衰减更快。

物理原因：

• 单极子是"点源"，场强随球面积（$\propto r^2$）稀释，故 $\propto 1/r^2$；
• 偶极子是"一对反向的点源"，在远处它们的贡献相互抵消一部分，净场强多衰减一个 $1/r$ 因子，故 $\propto 1/r^3$。

实际意义：

• 在地球表面（$r\approx R$）磁场约 $30,000\text{nT}$；
• 在 LEO 轨道（$r\approx 1.1R$）已衰减到约 $20,000\text{nT}$；
• 在 GEO 轨道（$r\approx 6.6R$）只剩约 $100\text{nT}$。

0.4.3 方向性特征

磁偶极子场的方向取决于观测点相对偶极轴的位置：

$$\mathbf{B}={\left(\frac{R}{r}\right)}^3\left[3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}-\mathbf{m}\right]$$

1. 在偶极轴上（$\hat{\mathbf{r}}\parallel \mathbf{m}$）：

   • $\hat{\mathbf{r}}\cdot \mathbf{m}=|\mathbf{m}|$
   • $\mathbf{B}=\frac{R^3}{r^3}(3|\mathbf{m}|\hat{\mathbf{r}}-\mathbf{m})=\frac{R^3}{r^3}\cdot 2|\mathbf{m}|\hat{\mathbf{r}}$
   • 磁场沿径向，指向外（北极）或内（南极），强度为偶极矩的 2 倍。

2. 在赤道面上（$\hat{\mathbf{r}}\perp \mathbf{m}$）：

   • $\hat{\mathbf{r}}\cdot \mathbf{m}=0$
   • $\mathbf{B}=-\frac{R^3}{r^3}\mathbf{m}$
   • 磁场反向于偶极矩，沿纬向环绕，强度为偶极矩的 1 倍。

3. 一般位置：磁场既有径向分量，也有切向分量，总体指向"从北极到南极"的磁场线方向。

磁倾角（inclination $I$）：

$$\tan I=2\tan \lambda$$

其中 $\lambda$ 是地理纬度。这是偶极子场的经典关系式（对于倾斜偶极需修正）。

0.5 IGRF 一阶球谐系数与偶极子的对应

0.5.1 为什么一阶项对应偶极子

在球谐展开中：

$$V(r,\theta ,\varphi )=a\sum _{n=1}^{\infty }{\left(\frac{a}{r}\right)}^{n+1}\sum _{m=0}^n(g_n^m\cos m\varphi +h_n^m\sin m\varphi )P_n^m(\cos \theta )$$

• $n=1$ 项对应 偶极子（两极性）；
• $n=2$ 项对应 四极子（四极性）；
• $n=3$ 项对应 八极子，以此类推。

数学上：

• $n=1$ 时，勒让德函数 $P_1^0(\cos \theta )=\cos \theta$，$P_1^1(\cos \theta )=\sin \theta$；
• 这恰好对应于偶极子势 ${\varphi }_m\propto (\mathbf{m}\cdot \mathbf{r})/r^3$ 在球坐标下的展开。

物理上：

• 偶极子是最简单的"有方向性"的磁源（比单极子高一级）；
• 地核发电机的球对称性自然产生以偶极为主的场结构。

0.5.2 为什么地球主磁场可以用一阶近似

能量占比：

• IGRF 完整模型展开到 $n=13$，但各阶能量贡献为：

  • $n=1$（偶极）：~90%
  • $n=2$（四极）：~5%
  • $n=3$（八极）：~2%
  • $n\ge 4$：~3%

• 因此，只用 $g_1^0,g_1^1,h_1^1$ 三个系数，就能在全球尺度上重现约 90% 的地磁场特征。

适用范围：

• 大尺度、远场：航天器轨道（LEO/MEO/GEO）、全球磁场分布；
• 粗略导航：罗经、磁强计定姿（精度要求不高时）。

不适用：

• 近地表局部磁异常（地壳磁化、矿藏）；
• 精密测量（地磁勘探、科学卫星）；
• 极区细节（磁极附近高阶项影响增大）。

0.6 磁偶极子近似的局限性

0.6.1 无法描述的地磁场特征

1. 磁异常（magnetic anomalies）：

   • 定义：局部磁场偏离偶极子预测的区域；
   • 成因：地壳中磁性矿物（如磁铁矿）、火山岩浆、断层构造；
   • 特征：空间尺度小（几十到几百公里），幅度可达数百 nT；
   • 例子：库尔斯克磁异常（俄罗斯）、南大西洋异常区（South Atlantic Anomaly, SAA）。
   • IGRF 处理：需要 $n\ge 2$ 的高阶项，或用局部模型（如 EMM）。

2. 高阶多极场（higher-order multipoles）：

   • 四极子（$n=2$）：在极区和赤道之间的过渡带影响显著；
   • 八极子（$n=3$）及以上：刻画更精细的空间结构；
   • 偶极近似误差：在地表可达 $\pm 3000\text{nT}$（约 10%）。

3. 时间变化（secular variation）：

   • 偶极子模型假设磁矩 $\mathbf{m}$ 恒定，但实际上地核流动导致磁场每年变化几十 nT；
   • IGRF 提供 ${\dot{g}}_n^m,{\dot{h}}_n^m$（年变化率），但偶极近似无法捕捉这种复杂的时变。

4. 非轴对称结构：

   • 真实地磁场不完全轴对称（偶极轴倾斜只是一阶修正）；
   • 经度方向的变化（如磁偏角的复杂分布）需要 $m\ne 0$ 的高阶项。

0.6.2 在 Basilisk 中的实际影响

适用场景（MagneticFieldCenteredDipole 足够）：

• LEO 卫星姿态控制：磁力矩器（magnetorquer）定姿，只需知道磁场方向和量级，10% 误差可接受；
• 初步任务设计：轨道仿真、功率预算、磁环境评估；
• 非地球行星：水星、木星等，只有偶极子数据时。

需要完整模型场景（应使用 MagneticFieldWMM）：

• 精密磁测：磁强计校准、科学任务（如 Swarm 卫星）;
• 近地表应用：无人机导航、地磁勘探；
• SAA 区域：辐射防护、电子设备可靠性分析。

Basilisk 的设计哲学：

• 提供两套工具：
  • MagneticFieldCenteredDipole：简单、快速、适合多行星；
  • MagneticFieldWMM：高精度、地球专用、读取完整 WMM.COF 系数文件。
• 用户根据任务需求选择。

0.7 本章小结

本章从物理基础出发，建立了对磁偶极子的全面理解：

1. 概念：磁偶极子是两极磁源，与电偶极子类比，但磁场线闭合；
2. 来源：地球磁场主要由地核发电机产生，以偶极子为主（~90% 能量）；
3. 数学：用磁矩矢量 $\mathbf{m}$ 描述，在 IGRF 中对应 $(g_1^1,h_1^1,g_1^0)$；
4. 空间分布：$1/r^3$ 衰减，极区径向、赤道切向；
5. IGRF 一阶近似：三个系数即可重现全球主磁场；
6. 局限性：无法描述局部异常、高阶结构、时变细节。

在接下来的章节中，我们将基于这些物理图像，详细推导 Basilisk 中居中偶极子磁场模型的数学实现。

---

1. 坐标系与符号约定

1.1 坐标系

• 行星惯性坐标系：$\mathcal{N}$（在 Basilisk 中为 J2000 惯性系），向量下标 $\_N$。
• 行星自转固定坐标系：$\mathcal{P}$（planet-fixed frame，在 Basilisk 中记为 P frame），向量下标 $\_P$。
  • ${\hat{\mathbf{p}}}_3$：自转轴方向（类似地球的北极方向）
  • ${\hat{\mathbf{p}}}_1$：赤道面内，指向本初子午面
  • ${\hat{\mathbf{p}}}_2$：补成右手系

在 Basilisk 中，行星状态消息 [SpicePlanetStateMsgPayload](file:///d:/small_paper/basilisk/src/architecture/msgPayloadDefC/SpicePlanetStateMsgPayload.h) 含有方向余弦矩阵
$$C_{P/N}=\text{planetState.J20002Pfix}$$
它将向量从 $\mathcal{N}$ 变换到 $\mathcal{P}$：
$${\mathbf{r}}_{B/P,P}=C_{P/N}{\mathbf{r}}_{B/P,N}$$

1.2 位置向量与半径

• ${\mathbf{r}}_{B/N,N}$：航天器相对于行星质心在惯性系中的位置。
• ${\mathbf{r}}_{P/N,N}$：行星质心在惯性系中的位置（通常为零向量，若原点在行星中心）。
• 相对位置向量：
  $${\mathbf{r}}_{B/P,N}={\mathbf{r}}_{B/N,N}-{\mathbf{r}}_{P/N,N}$$
• 在 Basilisk 基类 [MagneticFieldBase](file:///d:/small_paper/basilisk/src/simulation/environment/_GeneralModuleFiles/magneticFieldBase.cpp) 中：
  • 通过 updateRelativePos 计算
    $${\mathbf{r}}_{B/P,N}\to {\mathbf{r}}_{B/P,P}=C_{P/N}{\mathbf{r}}_{B/P,N}$$
  • 轨道半径（标量）
    $$r=\Vert {\mathbf{r}}_{B/P}\Vert$$
    在代码中记录为 this->orbitRadius。

1.3 IGRF 一阶高斯系数

在 IGRF 模型中，内源磁场的标量势在球坐标 $(r,\theta ,\varphi )$ 下可写为
$$V(r,\theta ,\varphi )=a\sum _{n=1}^N{\left(\frac{a}{r}\right)}^{n+1}\sum _{m=0}^n(g_n^m\cos m\varphi +h_n^m\sin m\varphi )P_n^m(\cos \theta )$$
其中：

• $a$：参考半径（地球通常取 $a\approx 6371.2\text{km}$），在 Basilisk 中为 planetRadius。
• $P_n^m$：缔合勒让德函数。
• $g_n^m,h_n^m$：高斯球谐系数，单位 nT。

当只保留 一阶项（n = 1） 时：
$$V_1(r,\theta ,\varphi )=a{\left(\frac{a}{r}\right)}^2\left[g_1^0\cos \theta +g_1^1\sin \theta \cos \varphi +h_1^1\sin \theta \sin \varphi \right]$$
这对应一个居中磁偶极子，其偶极矩与 $g_1^0,g_1^1,h_1^1$ 线性相关。

在 Basilisk 中，这三个系数（已从 nT 转为 T）存储在 MagneticFieldCenteredDipole 内：

// magneticFieldCenteredDipole.h
// IGRF 一阶系数，单位 [T]
double g10; //!< IGRF coefficient g_1^0
double g11; //!< IGRF coefficient g_1^1
double h11; //!< IGRF coefficient h_1^1

对于地球，这些系数在 Python 辅助函数
[simSetPlanetEnvironment.centeredDipoleMagField](file:///d:/small_paper/basilisk/src/utilities/simSetPlanetEnvironment.py#L54-L60) 中设置：

# 2020 IGRF 模型参数（单位从 nT 转成 T）
magFieldModule.g10 = -30926.00/1e9  # T
magFieldModule.g11 =  -2318.00/1e9  # T
magFieldModule.h11 =   5817.00/1e9  # T
magFieldModule.planetRadius = 6371.2*1000  # m

2. 从 IGRF 系数到偶极矩向量

2.1 一阶 IGRF 与偶极子标量势

对一阶项 $n=1$，标量势为
$$V_1(r,\theta ,\varphi )=a^3\frac{1}{r^2}\left[g_1^0\cos \theta +g_1^1\sin \theta \cos \varphi +h_1^1\sin \theta \sin \varphi \right]$$

另一方面，理想磁偶极子在球坐标下的标量势可写为（以参考半径 $a$ 归一化）
$$V_{\text{dip}}(r,\theta ,\varphi )={\left(\frac{a}{r}\right)}^{2}M(\cos \alpha \cos \theta +\sin \alpha \sin \theta \cos (\varphi -{\varphi }_m))$$
其中：

• $M$：偶极矩大小（常数）；
• $\alpha$：偶极倾角；
• ${\varphi }_m$：偶极在赤道面上的经度。

将两者对比，可以构造一个在行星固连坐标系 $\mathcal{P}$ 中的等效偶极矩向量 $\mathbf{m}\_P$：
$${\mathbf{m}}_P=\left[\begin{array}{c}{m}_x\\ m_y\\ m_z\end{array}\right]\propto \left[\begin{array}{c}{g}_1^1\\ h_1^1\\ g_1^0\end{array}\right]$$
差别只在于一个常数因子（单位和参考半径的处理）。Basilisk 将这些系数直接作为“偶极矩分量”的物理量使用，因此在代码中有：

// magneticFieldCenteredDipole.cpp
Eigen::Vector3d dipoleCoefficients; // first 3 IGRF coefficients
dipoleCoefficients << this->g11, this->h11, this->g10;

因此，在 Basilisk 中可以直接认为：
$${\mathbf{m}}_P=\left[\begin{array}{c}{g}_1^1\\ h_1^1\\ g_1^0\end{array}\right]\text{(单位 T)}$$
这里每个分量已经包含了“${\mu }_0/4\pi$”等物理常数以及从 nT 到 T 的单位转换。

3. 理想磁偶极子的磁场公式

3.1 磁偶极子场的矢量形式

在经典电磁学中，一个位于原点的磁偶极子 $\mathbf{m}$ 所产生的磁感应强度（忽略环境介质）可写为：
$$\mathbf{B}(\mathbf{r})=\frac{{\mu }_0}{4\pi }\frac{1}{r^3}(3(\mathbf{m}\cdot \hat{\mathbf{r}})\hat{\mathbf{r}}-\mathbf{m})$$
其中：

• $\mathbf{r}$：场点位置；$r=\Vert \mathbf{r}\Vert$；
• $\hat{\mathbf{r}}=\mathbf{r}/r$：径向单位向量；
• $\mathbf{m}$：偶极矩向量。

IGRF/WMM 模型中，球谐系数已经在参考半径 $a$ 处定义了磁标势，因此我们可以将常数和 ${\mu }_0/4\pi$ 吸收进 $\mathbf{m}$ 的定义，得到更工程化的形式：
$$\mathbf{B}(\mathbf{r})={\left(\frac{a}{r}\right)}^3(3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}-\mathbf{m})$$
其中 $\mathbf{m}$ 直接用 IGRF 系数构造（单位 T）。

3.2 Basilisk 采用的形式

Basilisk 里 [MagneticFieldCenteredDipole::evaluateMagneticFieldModel](file:///d:/small_paper/basilisk/src/simulation/environment/magneticFieldCenteredDipole/magneticFieldCenteredDipole.cpp) 的关键代码如下：

Eigen::Vector3d magField_P;         // [T] magnetic field in planet-fixed frame
Eigen::Vector3d rHat_P;             // [] normalized position vector in P frame
Eigen::Vector3d dipoleCoefficients; // [] first 3 IGRF coefficients

rHat_P = this->r_BP_P.normalized();

dipoleCoefficients << this->g11, this->h11, this->g10;
magField_P = pow(this->planetRadius/this->orbitRadius, 3)
             * (3 * rHat_P * rHat_P.dot(dipoleCoefficients)
                - dipoleCoefficients);

翻译成数学表达：

1. $\mathcal{P}$ 系下航天器位置单位向量：
   $${\hat{\mathbf{r}}}_P=\frac{{\mathbf{r}}_{B/P,P}}{r},r=\Vert {\mathbf{r}}_{B/P}\Vert$$
2. 偶极矩向量：
   $${\mathbf{m}}_P=\left[\begin{array}{c}{g}_1^1\\ h_1^1\\ g_1^0\end{array}\right]$$
3. 参考半径 $R=\text{planetRadius}$：地球为 $R=6371.2\text{km}$ 对应的米值。
4. 磁场（在行星固定系 $\mathcal{P}$ 分量）为：
   $$\boxed{{\mathbf{B}}_P({\mathbf{r}}_{B/P})={\left(\frac{R}{r}\right)}^3\left[3\right({\hat{\mathbf{r}}}_P\cdot {\mathbf{m}}_P){\hat{\mathbf{r}}}_P-{\mathbf{m}}_P]}$$

��以看出这正是理想磁偶极子公式，只是将常数合并进了 ${\mathbf{m}}_P$，并使用 IGRF 提供的一阶系数来构造该偶极矩。

3.3 磁偶极子公式的详细数学推导

下面从电磁学基本原理出发，详细推导上述公式。

3.3.1 磁偶极子的磁标势

在静磁学中，当磁场源远离观测点时，可以用磁标势 ${\varphi }_m$ 描述磁场（类比于静电学的电势）。对于一个位于原点、磁矩为 $\mathbf{m}$ 的点磁偶极子，其在空间任意点 $\mathbf{r}$ 处产生的磁标势为：

$${\varphi }_m(\mathbf{r})=\frac{{\mu }_0}{4\pi }\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}$$

物理意义：

• 磁标势 ${\varphi }_m$ 类似于电势，表示磁场的势能密度。
• $\mathbf{m}\cdot \mathbf{r}=mr\cos \theta$（其中 $\theta$ 是 $\mathbf{m}$ 与 $\mathbf{r}$ 的夹角），体现了磁偶极子的方向性：沿着磁矩方向势最大，垂直方向势为零。
• $1/r^3$ 的衰减反映了偶极场比单极场（$1/r^2$）衰减更快。

为了与 IGRF/WMM 的工程表达一致，我们引入参考半径 $R$（行星半径），并将常数 ${\mu }_0/(4\pi )$ 吸收进偶极矩 $\mathbf{m}$ 的定义中，得到归一化形式：

$$\boxed{{\varphi }_m(\mathbf{r})=\frac{R^3}{r^2}\frac{\mathbf{m}\cdot \mathbf{r}}{r}}$$

改写为：

$${\varphi }_m(\mathbf{r})=R^3\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}$$

其中 $\mathbf{m}$ 现在包含了所有物理常数，单位为 Tesla（在 Basilisk 中直接用 IGRF 系数构造）。

3.3.2 从磁标势到磁场：梯度运算

磁场强度 $\mathbf{B}$ 是磁标势的负梯度：

$$\mathbf{B}(\mathbf{r})=-\nabla {\varphi }_m(\mathbf{r})$$

代入上面的势函数：

$$\mathbf{B}=-\nabla \left(R^3\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}\right)$$

提取常数 $R^3$：

$$\mathbf{B}=-R^3\nabla \left(\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}\right)$$

3.3.3 计算梯度：分步推导

设 $\mathbf{m}$ 是常向量（在行星固定系中固定），$\mathbf{r}$ 是位置向量。我们需要计算：

$$\nabla \left(\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}\right)$$

使用乘积法则，设 $u=\mathbf{m}\cdot \mathbf{r}$，$v=1/r^3$：

$$\nabla (uv)=u\nabla v+v\nabla u$$

第一项：$\nabla (\mathbf{m}\cdot \mathbf{r})$

因为 $\mathbf{m}$ 是常向量：

$$\nabla (\mathbf{m}\cdot \mathbf{r})=\mathbf{m}$$

（这是因为 $\mathbf{m}\cdot \mathbf{r}=m_{x}x+m_{y}y+m_{z}z$，对 $\mathbf{r}$ 求梯度得 $\mathbf{m}$）

第二项：$\nabla (1/r^3)$

利用链式法则，$r=|\mathbf{r}|=\sqrt{x^2+y^2+z^2}$：

$$\frac{\partial r}{\partial x}=\frac{x}{r},\nabla r=\frac{\mathbf{r}}{r}=\hat{\mathbf{r}}$$

因此：

$$\nabla (r^{-3})=-3r^{-4}\nabla r=-3r^{-4}\hat{\mathbf{r}}=-\frac{3\hat{\mathbf{r}}}{r^4}$$

合并两项：

$$\nabla \left(\frac{\mathbf{m}\cdot \mathbf{r}}{r^3}\right)=(\mathbf{m}\cdot \mathbf{r})\nabla (r^{-3})+r^{-3}\nabla (\mathbf{m}\cdot \mathbf{r})$$

$$=(\mathbf{m}\cdot \mathbf{r})\left(-\frac{3\hat{\mathbf{r}}}{r^4}\right)+\frac{1}{r^3}\mathbf{m}$$

$$=-\frac{3(\mathbf{m}\cdot \mathbf{r})\hat{\mathbf{r}}}{r^4}+\frac{\mathbf{m}}{r^3}$$

3.3.4 代回磁场表达式

$$\mathbf{B}=-R^3\left[-\frac{3(\mathbf{m}\cdot \mathbf{r})\hat{\mathbf{r}}}{r^4}+\frac{\mathbf{m}}{r^3}\right]$$

$$=R^3\left[\frac{3(\mathbf{m}\cdot \mathbf{r})\hat{\mathbf{r}}}{r^4}-\frac{\mathbf{m}}{r^3}\right]$$

提取公因子 $1/r^3$：

$$\mathbf{B}=\frac{R^3}{r^3}\left[3\frac{\mathbf{m}\cdot \mathbf{r}}{r}\hat{\mathbf{r}}-\mathbf{m}\right]$$

注意到 $\mathbf{r}=r\hat{\mathbf{r}}$，所以：

$$\mathbf{m}\cdot \mathbf{r}=\mathbf{m}\cdot (r\hat{\mathbf{r}})=r(\mathbf{m}\cdot \hat{\mathbf{r}})$$

代入：

$$\mathbf{B}=\frac{R^3}{r^3}\left[3\frac{r(\mathbf{m}\cdot \hat{\mathbf{r}})}{r}\hat{\mathbf{r}}-\mathbf{m}\right]$$

$$=\frac{R^3}{r^3}\left[3(\mathbf{m}\cdot \hat{\mathbf{r}})\hat{\mathbf{r}}-\mathbf{m}\right]$$

即：

$$\boxed{\mathbf{B}(\mathbf{r})={\left(\frac{R}{r}\right)}^3\left[3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}-\mathbf{m}\right]}$$

这正是 Basilisk 代码中使用的公式！在行星固定系 $\mathcal{P}$ 中记为：

$${\mathbf{B}}_P({\mathbf{r}}_{B/P})={\left(\frac{R}{r}\right)}^3\left[3({\hat{\mathbf{r}}}_P\cdot {\mathbf{m}}_P){\hat{\mathbf{r}}}_P-{\mathbf{m}}_P\right]$$

3.3.5 公式各项的物理含义

将公式拆解为两部分：

$$\mathbf{B}={\left(\frac{R}{r}\right)}^3[\underset{\text{径向分量}}{\underbrace{3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}}}-\underset{\text{偶极矩方向}}{\underbrace{\mathbf{m}}}]$$

1. 衰减因子 $(R/r{)}^3$：

   • 磁偶极子场以 $1/r^3$ 速度衰减（比库仑场的 $1/r^2$ 快）。
   • $R$ 是参考半径，用于归一化，使得在 $r=R$ 处系数为 1。

2. 径向增强项 $3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}$：

   • $(\hat{\mathbf{r}}\cdot \mathbf{m})$ 是偶极矩在径向方向的投影，表征"沿着观测方向的偶极强度"。
   • 乘以 $3\hat{\mathbf{r}}$ 表示沿径向的磁场分量被放大 3 倍。
   • 物理上：在偶极轴线上（$\hat{\mathbf{r}}\parallel \mathbf{m}$），磁场沿轴向最强。

3. 偶极矩抵消项 $-\mathbf{m}$：

   • 减去偶极矩本身，使得在垂直于偶极轴的赤道面上，径向分量为零，磁场沿切向。
   • 当 $\hat{\mathbf{r}}\perp \mathbf{m}$ 时，$(\hat{\mathbf{r}}\cdot \mathbf{m})=0$，$\mathbf{B}=-(R/r{)}^3\mathbf{m}$，即磁场反向于偶极矩，指向南极。

4. 合成效果：

   • 极区（$\hat{\mathbf{r}}\parallel \mathbf{m}$）：$\mathbf{B}\approx (R/r{)}^3\cdot 2\mathbf{m}$，指向外（北极）或内（南极）。
   • 赤道（$\hat{\mathbf{r}}\perp \mathbf{m}$）：$\mathbf{B}\approx -(R/r{)}^3\mathbf{m}$，沿纬向环绕。

3.3.6 与 Basilisk 代码的对应

回顾代码：

rHat_P = this->r_BP_P.normalized();  // 单位径向向量
dipoleCoefficients << this->g11, this->h11, this->g10;  // 偶极矩
magField_P = pow(this->planetRadius/this->orbitRadius, 3)
             * (3 * rHat_P * rHat_P.dot(dipoleCoefficients)
                - dipoleCoefficients);

逐项对应：

• pow(planetRadius/orbitRadius, 3) → $(R/r{)}^3$
• rHat_P.dot(dipoleCoefficients) → ${\hat{\mathbf{r}}}_P\cdot {\mathbf{m}}_P$
• 3 * rHat_P * (...) → $3({\hat{\mathbf{r}}}_P\cdot {\mathbf{m}}_P){\hat{\mathbf{r}}}_P$
• - dipoleCoefficients → $-{\mathbf{m}}_P$

完全一致！

3.3.7 小结

从电磁学基本原理出发，我们推导了居中磁偶极子的磁场公式：

$$\mathbf{B}(\mathbf{r})={\left(\frac{R}{r}\right)}^3\left[3(\hat{\mathbf{r}}\cdot \mathbf{m})\hat{\mathbf{r}}-\mathbf{m}\right]$$

关键步骤：

1. 写出磁标势 ${\varphi }_m=R^3(\mathbf{m}\cdot \mathbf{r})/r^3$。
2. 通过 $\mathbf{B}=-\nabla {\varphi }_m$ 计算磁场。
3. 利用乘积法则和链式法则展开梯度。
4. 化简得到最终公式。

物理图像：

• 磁场由"径向增强"和"偶极抵消"两部分组成。
• 极区磁场指向外/内，赤道磁场沿纬向，典型的偶极场结构。
• Basilisk 通过 IGRF 一阶系数 $(g_1^1,h_1^1,g_1^0)$ 构造 ${\mathbf{m}}_P$，实现了地球磁场的偶极近似。

4. 坐标系变换：从行星固定系到惯性系

Basilisk 中各种环境模块输出的磁场都采用惯性系 $\mathcal{N}$（或仿真统一 N 系）的分量。对于居中偶极子模型，转换步骤如下：

1. 已有 $\mathcal{P}$ 系磁场 ${\mathbf{B}}_P$，如上公式所示。
2. 惯性系与行星固定系的方向余弦矩阵为 $C_{P/N}$，满足
   $${\mathbf{v}}_P=C_{P/N}{\mathbf{v}}_N$$
3. 反向变换为
   $${\mathbf{v}}_N=C_{N/P}{\mathbf{v}}_P=C_{P/N}^T{\mathbf{v}}_P$$

在代码中对应于：

// 将磁场从 P 系转换到 N 系，并写入消息
m33tMultV3(this->planetState.J20002Pfix, magField_P.data(), msg->magField_N);

m33tMultV3(A, x, y) 计算的是 $y=A^{T}x$，因此上述调用实现了：
$${\mathbf{B}}_N=C_{P/N}^T{\mathbf{B}}_P$$

最终，磁场消息 MagneticFieldMsgPayload 中的 magField_N 即为惯性系下的磁场分量。

5. 从 IGRF 系数到磁场矢量的完整流程总结

综上，Basilisk 中居中偶极子磁场模型采用 IGRF 一阶系数的完整数学关系可按以下步骤梳理：

1. 读取 IGRF 一阶系数与行星半径

   • 以地球为例：
     $$g_1^0,g_1^1,h_1^1\text{(单位 T)},R=R_{\text{planet}}$$
   • Python 中从 nT 转为 T，并设置 magFieldModule.g10/g11/h11 以及 planetRadius。

2. 构造偶极矩向量（行星固定系）

   $${\mathbf{m}}_P=\left[\begin{array}{c}{g}_1^1\\ h_1^1\\ g_1^0\end{array}\right]$$

3. 计算航天器相对行星的位置与单位径向向量

   • 惯性系相对位置：
     $${\mathbf{r}}_{B/P,N}={\mathbf{r}}_{B/N,N}-{\mathbf{r}}_{P/N,N}$$
   • 变换到行星固定系：
     $${\mathbf{r}}_{B/P,P}=C_{P/N}{\mathbf{r}}_{B/P,N}$$
   • 计算半径与单位向量：
     $$r=\Vert {\mathbf{r}}_{B/P}\Vert ,{\hat{\mathbf{r}}}_P=\frac{{\mathbf{r}}_{B/P,P}}{r}$$

4. 居中偶极子磁场（行星固定系）

   $${\mathbf{B}}_P({\mathbf{r}}_{B/P})={\left(\frac{R}{r}\right)}^3\left[3\right({\hat{\mathbf{r}}}_P\cdot {\mathbf{m}}_P){\hat{\mathbf{r}}}_P-{\mathbf{m}}_P]$$

5. 坐标变换到惯性系

   $${\mathbf{B}}_N=C_{P/N}^T{\mathbf{B}}_P$$

6. 磁场强度与方向

   • 磁场强度（标量）：
     $$\Vert \mathbf{B}\Vert =\sqrt{B_x^2+B_y^2+B_z^2}$$
   • 磁场方向单位向量：
     $$\hat{\mathbf{B}}=\frac{\mathbf{B}}{\Vert \mathbf{B}\Vert }$$

其中 $\mathbf{B}$ 可以取 ${\mathbf{B}}_P$（行星固定系）或 ${\mathbf{B}}_N$（惯性系），视具体使用需求而定。

6. 与 IGRF/WMM 完整模型的关系

• IGRF 与 WMM 完整模型都包含高阶 $n>1$ 的球谐系数，会提供更精细的空间结构。
• Basilisk 中的 MagneticFieldCenteredDipole 只使用 IGRF 的一阶项，构造一个居中偶极子模型：
  • 优点：实现简单、适合多行星（只需提供等效偶极的强度和方向）。
  • 缺点：忽略高阶项导致精度有限，无法描述局地磁异常等。
• 若需要高精度地球磁场，Basilisk 提供了基于 WMM 系数文件 WMM.COF 的 MagneticFieldWMM 模块，可以视为“完整球谐版”，而 MagneticFieldCenteredDipole 则是“IGRF 一阶近似版”。