# 十一月新增功能说明

本期（2026 年 11 月）在 moongeodesy 里新增了两个包：

- `magnetic`：世界地磁模型 WMM2025；
- `conic`：兰勃特等角圆锥投影和阿尔伯斯等积圆锥投影。

另外，把 README 里的每一条说法都对上了测试。本文讲清楚这些代码在做什么、为什么这样写、正确性怎么验证。比赛规则要求参赛者能解释 AI 写的代码，常见问题放在每一节的最后。

为什么选这两个功能，见 [`十一月新增功能选型.md`](十一月新增功能选型.md)。

---

## 一、`magnetic`：WMM2025 地磁模型

### 解决什么问题

罗盘指向磁北，地图用的是真北，两者之间的夹角叫**磁偏角**。2026 年 11 月，北京的磁偏角约 −7.6°（磁北偏西），漠河约 −12.8°，喀什约 +4.1°。如果按地图上量出的方位角直接走罗盘，偏 7.5° 时每走 2 公里会偏出约 260 米（场景 s5 计算并检查了这个数）。无人机、手机指南针、野外测量都要做这个改正。

WMM 是美国国家地理空间情报局（NGA）和英国国防地理中心（DGC）的官方模型，由 NOAA NCEI 和英国地质调查局（BGS）制作，每 5 年更新一次。WMM2025 的适用时间是 2025.0 到 2030.0。

### 模型是什么

地球主磁场可以写成一个标量位势 V 的负梯度，V 用 12 阶球谐函数展开：

```
V(r, φ', λ, t) = a · Σ(n=1..12) (a/r)^(n+1) · Σ(m=0..n) [g(n,m,t)·cos(mλ) + h(n,m,t)·sin(mλ)] · P(n,m)(sin φ')
```

- a = 6371200 m 是地磁参考半径；
- (r, φ', λ) 是地心球坐标：半径、地心纬度、经度；
- g、h 是高斯系数，共 90 组，随时间线性变化：g(t) = g + (t − 2025.0)·ġ；
- P(n,m) 是施密特半归一化的连带勒让德函数。

### 计算步骤（`src/magnetic/magnetic.mbt`）

1. **大地坐标 → 地心球坐标**（`geocentric`）：用 WGS-84 椭球把纬度 φ、椭球高 h 换成地心纬度 φ' 和地心距 r。
2. **当日系数**（`coefficients_at`）：g(t) = g + (t − t0)·ġ。
3. **勒让德函数**（`legendre_poly`）：见下面“关键设计”。
4. **求和得到地心坐标系下的磁场**（`sum_field`）：
   - 北向分量 X' 对 P 求关于 φ' 的导数；
   - 东向分量 Y' 要用 P / cos φ'；
   - 垂直分量 Z' 带 (n+1) 因子。

   年变化率用同样的求和，只是把系数换成 ġ、ḣ。
5. **旋转到椭球坐标系**：地心纬度和大地纬度差 ψ = φ' − φ，X = X'cos ψ − Z'sin ψ，Z = X'sin ψ + Z'cos ψ。
6. **派生量**：
   - 水平强度 H = √(X² + Y²)，总强度 F = √(H² + Z²)；
   - 磁倾角 I = atan2(Z, H)，磁偏角 D = atan2(Y, X)；
   - 它们的年变化率按求导公式算。

### 关键设计：极点不需要特殊处理

东向分量 Y' 要用 P(n,m)/cos φ'。在地理极点 cos φ' = 0，直接除会得到 0/0。NOAA 官方 C 程序在极点走单独的分支。

我们的写法是把勒让德函数拆成两部分：

```
P(n,m)(x) = cos^m(φ') · T(n,m)(x)，  x = sin φ'
```

T(n,m) 是多项式，用和 P 相同的递推式计算：

- T(m,m) = √((2m−1)/(2m)) · T(m−1,m−1)；
- T(n,m) = [(2n−1)·x·T(n−1,m) − √((n−1)²−m²)·T(n−2,m)] / √(n²−m²)。

导数 T' 也一起递推。这样：

- P/cos φ' = cos^(m−1)·T，在极点有限；
- dP/dφ' = −m·x·cos^(m−1)·T + cos^(m+1)·T'。

数组 `cpow[k] = cos^k`，其中 `cpow[0] = 1`，所以极点上 m = 1 的项是 1 × T，自然有限，不需要特殊分支。测试“field is finite and continuous at the geographic poles”检查了这一点：极点与离极点 1e-7°（约 1 厘米）的点相差不到 1e-3 nT。

### 其他接口

- `grivation(decl, gamma)`：格网磁偏角，即从格网北到磁北的角度，等于磁偏角减子午线收敛角。UTM 的收敛角直接用本库 `utm` 包的 `from_latlon(...).gamma`。场景 s5 演示了两种算法：按真北算的罗盘方位，和从 UTM 地图格网方位算的罗盘方位，两者必须相等，否则程序中止。
- `MagneticElements::zone()`：按 WMM2025 的定义判断罗盘还能不能用。H < 2000 nT 为“禁用区”（Blackout），H < 6000 nT 为“警告区”（Caution）。
- `decimal_year(y, m, d)`：与 NOAA 软件一致，等于“年 + (年内第几天 − 1) / 当年天数”。
- 输入超出模型范围时抛出 `MagneticError`：纬度不在 ±90° 内、经度不是有限数、高度不在 −1 km 到 850 km、日期不在 2025.0 到 2030.0。这些限制都来自模型规范 MIL-PRF-89500B。

### 怎么验证

| 参考数据 | 来源 | 结果 |
| --- | --- | --- |
| 100 个双精度官方测试点（6 位小数） | NOAA 随系数一起发布的 `WMM2025_TestValues.txt` | Y 在文件舍入范围内（5e-7 nT）；Z < 2.5e-6 nT；X、H、F < 1e-3 nT（差异出在文件里，见下） |
| 技术报告表 6 的 12 个标准值（含格网磁偏角） | NOAA | 全部在舍入范围内 |
| 技术报告表 3b 的逐步算例（10 位小数） | NOAA | 最终结果每一项都在最后一位的半个单位内；中间量（地心纬度、地心半径、当日系数、地心分量及其变化率）在最后一位的一个单位内 |
| 技术报告表 3d 的 UPS/UTM 格网磁偏角 | NOAA；UTM 收敛角用本库 `utm` 算 | < 1e-9° |
| 119 个点的独立高精度计算 | `tools/gen_wmm_independent.py`（mpmath，40 位精度） | < 1e-10 nT |

独立计算刻意用了完全不同的算法，和 MoonBit 代码没有共享任何步骤：

- 勒让德函数用罗德里格斯公式的精确有理系数，而不是递推；
- 磁场用位势在地心直角坐标系里的数值梯度，而不是解析求导；
- 分量用椭球局部坐标系的单位向量投影，而不是按 ψ 旋转。

**一个发现：** NOAA 的 100 点文件里，X 分量和我们差最多 7.2e-4 nT（相对误差 1.3e-8），Z 差最多 2.2e-6 nT，Y 在文件的舍入范围内。独立计算给出的 X、Z 都和我们一样；在赤道上（ψ = 0）文件的 X 也和我们一致到舍入（有测试检查）。所以这些微小差异出在 NOAA 的文件里，不是我们的代码。7.2e-4 nT 远小于报告里说的单精度软件误差 0.1 nT，更远小于模型本身约 100 nT 的精度，对使用没有影响。测试按实际情况设限，并在注释和 README 里写明了原因。

### 常见问题

- **为什么不直接移植 NOAA 的 C 代码或 GeographicLib？** 按报告的方程自己写，NOAA 官方测试值（由 NOAA 自己的软件算出）就成了独立的参考；再加上一个算法完全不同的 40 位精度实现做第二个参考。如果是移植，就只能用同源的数据来测。
- **高度用的是海拔还是椭球高？** 椭球高。两者差一个大地水准面高。报告第 1.2 节说这对磁场的影响约 1 nT 或更小，所以没有内置大地水准面模型。
- **2030 年以后怎么办？** NOAA 会发布 WMM2030。到时用 `tools/gen_wmm.py` 换系数文件，脚本会校验文件的 SHA-256，再加一个模型构造函数即可。

---

## 二、`conic`：兰勃特等角圆锥和阿尔伯斯等积圆锥

### 解决什么问题

- **阿尔伯斯等积圆锥投影（Albers）**：投影后的面积和地面面积相等，用于面积统计。国内常用 CGCS2000 或克拉索夫斯基椭球，双标准纬线 25°N 和 47°N，中央经线 105°E。
- **兰勃特等角圆锥投影（LCC）**：保持角度不变，用于中小比例尺地图和很多国家的平面坐标系，例如美国州平面坐标、比利时 Lambert 72、法国 Lambert II。

场景 s6 用 0.5° × 0.5° 的格子对比了三种投影的面积误差，从三亚到漠河：

| 投影 | 面积误差 |
| --- | --- |
| Albers | 与椭球面积一致到 1e-9 |
| 同参数的兰勃特投影 | 最多差 6% |
| 网页地图常用的 Web Mercator | 在漠河差 176% |

### 实现（`src/conic/`）

- `lcc.mbt`、`albers.mbt`：移植 GeographicLib 的 `LambertConformalConic` 和 `AlbersEqualArea`（Karney，MIT/X11 许可证）。
- `grid.mbt`：EPSG 风格的格网，包括假原点纬度、中央经线、假东和假北，对应 EPSG 方法 9801、9802 和 9822。
- `common.mbt`：错误类型、结果结构，以及“均差”（divided differences）工具函数。

**为什么要用均差。** 教科书（Snyder，EPSG 指南 7-2）里，兰勃特投影的圆锥常数是：

```
n = (ln m1 − ln m2) / (ln t1 − ln t2)
```

两条标准纬线很接近时，分子分母都接近 0，双精度会丢掉很多位有效数字。Karney 的写法把 (f(x) − f(y)) / (x − y) 改写成不需要相减的形式，叫均差，例如：

- `dasinh`
- `dlog1p`
- `dsn`

这样在下面几种情况下都保持精度：

- 两条标准纬线很接近；
- 标准纬线在赤道附近，投影退化为墨卡托；
- 标准纬线在极点附近，投影退化为极球面。

测试里有一组标准纬线 40° 和 40.0001° 的例子。PROJ 在这组上差了 6 毫米，正是这个消去问题；我们和 50 位精度的计算一致到纳米级。另一个测试把两条纬线拉开 1e-9°，我们算出的几微米差别与 50 位精度的值一致到 2e-8 米。

**移植时的一处判断。** GeographicLib 2.5 的 `AlbersEqualArea::Forward` 对纬度乘了两次半球符号：源码第 499 行和第 501 行各乘一次。我们只乘一次。南半球的测试证实这样是对的：

- GIGS 5109 澳大利亚 Albers；
- PROJ 对照里的澳大利亚 Albers −18°/−36°。

两者都通过了。

**`ConicGrid` 的假原点。** EPSG 的双标准纬线方法规定了“假原点纬度”。我们先算出这个纬度在中央经线上的北坐标 y0，之后所有北坐标都减去 y0 再加假北。比利时 Lambert 72 的假原点在北极（纬度 90°），测试也覆盖到了。

### 怎么验证

我们移植的是 GeographicLib 的写法，下面这些参考数据都不用这种写法：

| 参考数据 | 内容 | 结果 |
| --- | --- | --- |
| EPSG 指南 7-2 / Snyder 闭式公式的 50 位精度计算（`tools/gen_conic.py`） | 18 套投影，1792 个点，覆盖中国 CGCS2000 和克拉索夫斯基、美国、欧洲、澳大利亚、南半球、极区、赤道附近、单标准纬线、球面 | 坐标 < 2.5e-8 m，收敛角 < 1e-12°，比例因子 < 1e-13 |
| PROJ 9.5（C 语言的 Snyder 公式实现） | 同样的点 | < 2.5e-7 m；标准纬线极近的那组 < 1 cm，原因见上 |
| IOGP GIGS 官方测试集 5102、5103、5109 | 法国 EuroLambert、法国 NTF Lambert II（百分度、巴黎本初子午线）、比利时 Lambert 72、美国犹他北区（英尺和美国测量英尺）、澳大利亚 Albers | 兰勃特在 1 毫米内，Albers 在 5 毫米内。比文件印刷的最后一位略大，但 PROJ 在同样的点上差得一样（例如 4.5 毫米），所以差异在 GIGS 的数值里；GIGS 自己的容差是 3–5 厘米 |
| EPSG 指南 7-2 算例 | 德州中南区（双标准纬线，美国测量英尺）、牙买加（单标准纬线） | 印刷精度内 |
| Snyder《Map Projections》附录 A 算例 | Albers 和兰勃特各一个，克拉克 1866 椭球 | 0.1 米内（Snyder 手算约 8 位有效数字） |

还有几项性质测试，都用本库的测地线来“量”：

- **等角**：沿北、东、东北三个方向量出的比例尺，都等于 `k`，误差 < 1e-9。
- **等积**：中国范围一个 30° × 50° 的大框，Albers 平面面积与椭球面积相对差 1.8e-9。
- **退化情况**：标准纬线取 0° 时等于墨卡托（误差 < 1e-7 m）；取 90° 时等于本库的 UPS 极球面投影（< 1e-8 m）。
- **长轴在极轴方向的椭球（扁率为负）**：PROJ 不支持，用本库的测地线检查等角、等积和往返。
- **极点**：Albers 在极点附近把经线压得极扁，离极点约 1e-6°（10 厘米）以内的点在地图上和极点本身的坐标在舍入意义下相同。反算遇到这种输入就直接返回极点，所以极点本身往返精确，离极点 1e-4° 以外的点仍精确到 1e-7°。兰勃特投影的两个极点也都能往返。

### 常见问题

- **中国有官方的 Albers 参数吗？** 没有单一的官方参数。常见做法是 25°N/47°N、中央经线 105°E 或 110°E、原点在赤道，椭球用 CGCS2000 或克拉索夫斯基，所以库里没有内置“官方预设”，参数由用户传入。
- **GIGS 数据的许可？** GIGS 版权属于 IOGP。使用条款允许在注明 IOGP 所有权的前提下使用、复制和分发，禁止把 GIGS 本身拿去营利，商业软件可以包含它，只要价值在于软件本身。数据只放在测试文件 `src/conic/gigs_data_test.mbt` 里，文件头、README 和 `THIRD_PARTY_NOTICES.md` 都写明了来源和条款，并说明这个文件不适用本项目的 MIT 许可。

---

## 三、本期对十月验收问题的整改和 README 核对

十月验收没过的三条已在 11/08 修复（见 PR #32）：

1. 国外坐标被 GCJ 偏移（首尔、河内）；
2. 矩形边缘的反算误差；
3. ±90° 处 UPS 计算出 NaN。

教训是：测试全过，但 README 里有一条承诺没有测试，验收照样不过。所以本期把 README 和各包说明逐条核对了两遍，清单见 [`十一月验收清单.md`](十一月验收清单.md)。主要改动：

- **国外沿海城镇和小岛也不再偏移。** 第二遍核对发现，十月的修复只覆盖了内陆：远离中国的外国海岸线被简化到约 5.5 公里，20 平方公里以下的小岛被丢掉，所以那霸、槟城、亚庇、下龙、郁陵岛、南浦等沿海地点仍被当成海面加了偏移。现在 `tools/gen_china_boundary.py` 用 Natural Earth 的国家和小岛两个图层，外国陆地向海外扩约 9 公里，再精确扣掉中国陆地；查询改成按 1° 网格索引，速度也更快（约 100 纳秒）。检验不再只靠手选城市，而是用与 Natural Earth 无关的 GeoNames 全部 5000 人以上城镇：12115 个外国城镇不偏移，3234 个中国（含港澳台）城镇偏移；只有 3 个位于 Natural Earth 未收录小岛上的外国城镇例外，README 写明了它们的名字。
- **负扁率椭球上的反算收敛。** 从 GeographicLib 移植的 `tauf` 用 1 − es² 代替 1 − e²，扁率为负时牛顿迭代只线性收敛、提前停止，兰勃特往返差 4.6e-8°。改用带符号的 e² 后精确；扁率为正的常用椭球不受影响。

- **ECEF / 局部 ENU 的参考数据换成独立来源。** 原来的参考数据是同一份 GeographicLib 代码的 Python 转写，和被测代码同源。现在改为 50 位精度的教科书公式加 PROJ，结果一致到 7.5e-9 米。
- **测试断言收紧到 README 公布的数值。** 原来有的断言比 README 宽 2 到 270 倍。现在 README 写的就是测试断言的界限，括号里注明实测值。
- **场景程序结果不对就中止。** 原来只打印结果，现在 s1–s6 都会检查自己打印的关键数字。CI 在 4 个后端都跑这些场景。
- **README 的代码块自动生成测试。** `tools/sync_readme_tests.py` 把 README 的代码块变成测试；如果两者不一致，CI 失败。
- **补上没有测试的 API 承诺。** 例如：MGRS 解码默认返回 1 米格子中心、`Geodesic::line` 在设定距离前距离是 NaN、UTM 各种越界报错。
- **改正写错或过时的说法。** 例如：
  - GCJ 单步近似反算的误差实际可达 3.04 米，原来写的是 1–2 米；
  - UPS 极点为 NaN 的旧说明，已经不成立；
  - 测试数量；
  - 基准测试数据（加入国界判断后，GCJ 正算变慢了）；
  - 代码行数。

## 四、复现

```bash
bash tools/acceptance.sh     # 跑清单里全部自动检查（首次编译需几分钟）
moon test --target all       # 只跑测试
moon run src/examples/s5_compass
moon run src/examples/s6_conic_area
```

重新生成参考数据需要 Python 和 `pip install mpmath pyproj geographiclib shapely`。每个脚本开头都写明了数据来源和下载地址。