mingli-astro 1.1.0

Computational astronomy and calendrics: Julian day, apparent solar longitude and the 24 solar terms, lunar phase and lunisolar intercalation, and the sexagenary cycle.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

//! 朔（新月）时刻：Meeus ch.49 截断模型，精度约数分钟，足以判定朔落在哪一民用日。

use crate::{delta_t_seconds, year_of_jd};

fn mean_phase_jde(k: f64, t: f64) -> f64 {
    2451550.09766 + 29.530588861 * k + 0.00015437 * t * t - 0.000000150 * t * t * t
        + 0.00000000073 * t * t * t * t
}

/// 第 `k` 个朔（phase=0，新月）的 JD(UT)。`k=0` 对应 2000-01-06 前后的朔。
#[must_use]
pub fn new_moon_jd_ut(k: i64) -> f64 {
    let k = k as f64;
    let t = k / 1236.85;
    let mut jde = mean_phase_jde(k, t);
    let e = 1.0 - 0.002516 * t - 0.0000074 * t * t;
    let m = (2.5534 + 29.10535670 * k - 0.0000014 * t * t - 0.00000011 * t * t * t).to_radians();
    let mp = (201.5643 + 385.81693528 * k + 0.0107582 * t * t + 0.00001238 * t * t * t
        - 0.000000058 * t * t * t * t)
        .to_radians();
    let f = (160.7108 + 390.67050284 * k - 0.0016118 * t * t - 0.00000227 * t * t * t
        + 0.000000011 * t * t * t * t)
        .to_radians();
    let omega = (124.7746 - 1.56375588 * k + 0.0020672 * t * t + 0.00000215 * t * t * t).to_radians();

    let corr = -0.40720 * mp.sin()
        + 0.17241 * e * m.sin()
        + 0.01608 * (2.0 * mp).sin()
        + 0.01039 * (2.0 * f).sin()
        + 0.00739 * e * (mp - m).sin()
        - 0.00514 * e * (mp + m).sin()
        + 0.00208 * e * e * (2.0 * m).sin()
        - 0.00111 * (mp - 2.0 * f).sin()
        - 0.00057 * (mp + 2.0 * f).sin()
        + 0.00056 * e * (2.0 * mp + m).sin()
        - 0.00042 * (3.0 * mp).sin()
        + 0.00042 * e * (m + 2.0 * f).sin()
        + 0.00038 * e * (m - 2.0 * f).sin()
        - 0.00024 * e * (2.0 * mp - m).sin()
        - 0.00017 * omega.sin()
        - 0.00007 * (mp + 2.0 * m).sin()
        + 0.00004 * (2.0 * mp - 2.0 * f).sin()
        + 0.00004 * (3.0 * m).sin()
        + 0.00003 * (mp + m - 2.0 * f).sin()
        + 0.00003 * (2.0 * mp + 2.0 * f).sin()
        - 0.00003 * (mp + m + 2.0 * f).sin()
        + 0.00003 * (mp - m + 2.0 * f).sin()
        - 0.00002 * (mp - m - 2.0 * f).sin()
        - 0.00002 * (3.0 * mp + m).sin()
        + 0.00002 * (4.0 * mp).sin();
    jde += corr;

    // 附加行星摄动项（A1..A14）
    let angles = [
        299.77 + 0.107408 * k - 0.009173 * t * t,
        251.88 + 0.016321 * k,
        251.83 + 26.651886 * k,
        349.42 + 36.412478 * k,
        84.66 + 18.206239 * k,
        141.74 + 53.303771 * k,
        207.14 + 2.453732 * k,
        154.84 + 7.306860 * k,
        34.52 + 27.261239 * k,
        207.19 + 0.121824 * k,
        291.34 + 1.844379 * k,
        161.72 + 24.198154 * k,
        239.56 + 25.513099 * k,
        331.55 + 3.592518 * k,
    ];
    let coef = [
        0.000325, 0.000165, 0.000164, 0.000126, 0.000110, 0.000062, 0.000060, 0.000056, 0.000047,
        0.000042, 0.000040, 0.000037, 0.000035, 0.000023,
    ];
    let add: f64 = coef
        .iter()
        .zip(angles.iter())
        .map(|(c, ang)| c * ang.to_radians().sin())
        .sum();
    jde += add;

    // JDE（力学时） → JD(UT)
    jde - delta_t_seconds(year_of_jd(jde)) / 86400.0
}

/// 估计某 JD 附近朔的整数序号 `k`（仅 crate 内使用）。
pub(crate) fn new_moon_k_near(jd: f64) -> i64 {
    ((jd - 2451550.09766) / 29.530588861).round() as i64
}