astroceleste_engine/chinese/
calendar.rs

1//! The Chinese calendar around a moment: the solar terms, the lunar months and the
2//! equation of time, which Ba Zi and Zi Wei Dou Shu are cast from.
3//!
4//! The rules are those of the modern calendar: true solar terms and true new moons (in
5//! force since 1645), dated in China Standard Time (UTC+8) from 1929 and on the Beijing
6//! meridian (116°25′ E) before. They are applied to every date, so dates before 1645
7//! follow them proleptically rather than the mean-term calendars then in use.
8//!
9//! A lunar month starts on the day (China time) of a new moon. Month 11 is the month that
10//! holds the winter solstice. When 13 months run from one month 11 to the next, the first
11//! of them without a major term (中氣) is intercalary and repeats the number before it.
12
13use serde::{Deserialize, Serialize};
14
15use crate::constants::TAU;
16use crate::ephemeris::observe::{apparent, observe};
17use crate::ephemeris::{KernelSet, SpkError};
18use crate::error::EngineError;
19use crate::frames::nutation::iau2000b_radians;
20use crate::frames::{mxv, Orientation};
21use crate::instant::UtcInstant;
22use crate::planets::require_kernel;
23use crate::search::{instant_of, Sky, SUN, SYNODIC_MONTH, TROPICAL_YEAR};
24use crate::time::Time;
25
26/// Julian day number of 1970-01-01.
27const UNIX_EPOCH_JDN: i64 = 2_440_588;
28const MICROS_PER_DAY: i64 = 86_400_000_000;
29/// Births this close (days) to a lunar month boundary may fall either side of it once
30/// converted to the local civil or solar date, so the months around are kept too.
31const DATE_MARGIN: i64 = 3;
32
33/// The 24 solar terms by the Sun's longitude / 15°, from the spring equinox.
34const TERM_NAMES: [&str; 24] = [
35    "chunfen",
36    "qingming",
37    "guyu",
38    "lixia",
39    "xiaoman",
40    "mangzhong",
41    "xiazhi",
42    "xiaoshu",
43    "dashu",
44    "liqiu",
45    "chushu",
46    "bailu",
47    "qiufen",
48    "hanlu",
49    "shuangjiang",
50    "lidong",
51    "xiaoxue",
52    "daxue",
53    "dongzhi",
54    "xiaohan",
55    "dahan",
56    "lichun",
57    "yushui",
58    "jingzhe",
59];
60
61/// The longitude of 立春, the Beginning of Spring, which opens the year of the pillars.
62pub(crate) const LICHUN: u16 = 315;
63
64/// One of the 24 solar terms (節氣): the moment the Sun reaches a multiple of 15°.
65#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
66pub struct SolarTerm {
67    /// Pinyin code, e.g. "lichun".
68    pub name: String,
69    /// The Sun's apparent ecliptic longitude, degrees (a multiple of 15).
70    pub longitude: u16,
71    /// A "jie" (節) term, which opens a solar month (longitudes 15° + 30°·k); otherwise a
72    /// major term (中氣, multiples of 30°).
73    pub jie: bool,
74    /// The moment, ISO 8601 UTC to the second.
75    pub utc: String,
76}
77
78/// A month of the lunisolar calendar.
79#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
80pub struct LunarMonth {
81    /// The lunar year: the Gregorian year its New Year falls in.
82    pub year: i32,
83    /// Month number, 1-12.
84    pub month: u8,
85    /// Whether it is the intercalary month (閏月) repeating `month`.
86    pub leap: bool,
87    /// First day, `YYYY-MM-DD` (China time).
88    pub start: String,
89    /// Length in days: 29 or 30.
90    pub days: u8,
91}
92
93/// The Chinese calendar around a moment: what the Four Pillars and the Purple Star chart
94/// need from the ephemeris.
95#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
96pub struct ChineseCalendar {
97    /// Apparent minus mean solar time at the moment, minutes.
98    pub equation_of_time: f64,
99    /// Every solar term from the last 立春 at or before the moment to the first jie term
100    /// after it.
101    pub solar_terms: Vec<SolarTerm>,
102    /// The lunar months around the date of the moment (a few days either side).
103    pub lunar_months: Vec<LunarMonth>,
104}
105
106impl SolarTerm {
107    /// The moment of the term.
108    pub(crate) fn instant(&self) -> Result<UtcInstant, EngineError> {
109        UtcInstant::parse(&self.utc).map_err(|e| EngineError::InvalidInput(e.to_string()))
110    }
111}
112
113impl LunarMonth {
114    /// Julian day number of the first day.
115    pub(crate) fn start_day(&self) -> Result<i64, EngineError> {
116        let start = UtcInstant::parse(&format!("{}T00:00:00Z", self.start))
117            .map_err(|e| EngineError::InvalidInput(e.to_string()))?;
118        Ok(day_number(start))
119    }
120}
121
122/// Julian day number of the UTC date of `instant`.
123pub(crate) fn day_number(instant: UtcInstant) -> i64 {
124    instant.micros().div_euclid(MICROS_PER_DAY) + UNIX_EPOCH_JDN
125}
126
127/// `YYYY-MM-DD` of a Julian day number.
128fn date_of(day: i64) -> String {
129    let (y, m, d, ..) = UtcInstant::from_micros((day - UNIX_EPOCH_JDN) * MICROS_PER_DAY).civil();
130    format!("{y:04}-{m:02}-{d:02}")
131}
132
133/// Offset of China time from UT, days: UTC+8 from 1929, the Beijing meridian before.
134fn china_offset(jd: f64) -> f64 {
135    // 1929-01-01T00:00 at UTC+8.
136    const CST_SINCE_JD: f64 = 2_425_612.5 - 8.0 / 24.0;
137    if jd < CST_SINCE_JD {
138        (116.0 + 25.0 / 60.0) / 360.0
139    } else {
140        8.0 / 24.0
141    }
142}
143
144/// Julian day number of the date in China at UT Julian date `jd`.
145fn china_day(jd: f64) -> i64 {
146    (jd + 0.5 + china_offset(jd)).floor() as i64
147}
148
149impl Sky<'_> {
150    /// The new moon that starts the month holding the China-time day of the winter
151    /// solstice `solstice`: month 11.
152    fn month_eleven(&self, solstice: f64) -> Result<f64, EngineError> {
153        let solstice_day = china_day(solstice);
154        let moon = self.new_moon_before(solstice + 1.0)?;
155        if china_day(moon) > solstice_day {
156            self.new_moon_before(moon - 1.0)
157        } else {
158            Ok(moon)
159        }
160    }
161
162    /// The months of the year (歲) from the month 11 holding `solstice` to the next
163    /// month 11, with the first day of that next month 11 and the next solstice.
164    fn sui(&self, solstice: f64) -> Result<(Vec<Month>, i64, f64), EngineError> {
165        let next_solstice = self.sun_at(270.0, solstice + TROPICAL_YEAR)?;
166        let end = china_day(self.month_eleven(next_solstice)?);
167
168        let mut moon = self.month_eleven(solstice)?;
169        let mut starts = vec![china_day(moon)];
170        loop {
171            moon = self.new_moon_near(moon + SYNODIC_MONTH)?;
172            let day = china_day(moon);
173            if day >= end {
174                break;
175            }
176            starts.push(day);
177            if starts.len() > 13 {
178                return Err(diverged("a year of more than 13 months"));
179            }
180        }
181
182        // The major terms: the solstice and the eleven after it (the next solstice
183        // already starts the next year).
184        let mut major_days = vec![china_day(solstice)];
185        let mut term = solstice;
186        for k in 1..12 {
187            term = self.sun_at(f64::from((270 + 30 * k) % 360), term + 30.44)?;
188            major_days.push(china_day(term));
189        }
190
191        let (first_year, ..) =
192            UtcInstant::from_micros((starts[0] - UNIX_EPOCH_JDN) * MICROS_PER_DAY).civil();
193        let mut ends: Vec<i64> = starts[1..].to_vec();
194        ends.push(end);
195        let mut leap_pending = starts.len() == 13;
196        let mut number = 11u8;
197        let mut year = first_year;
198        let mut months = Vec::with_capacity(starts.len());
199        for (i, (&start, &month_end)) in starts.iter().zip(&ends).enumerate() {
200            let has_major = major_days.iter().any(|&d| start <= d && d < month_end);
201            let leap = i > 0 && leap_pending && !has_major;
202            if leap {
203                leap_pending = false;
204            } else if i > 0 {
205                number = number % 12 + 1;
206                if number == 1 {
207                    year = first_year + 1;
208                }
209            }
210            months.push(Month {
211                start,
212                end: month_end,
213                number,
214                leap,
215                year,
216            });
217        }
218        Ok((months, end, next_solstice))
219    }
220
221    /// Equation of time at `jd`, minutes: apparent minus mean solar time.
222    fn equation_of_time(&self, jd: f64) -> Result<f64, EngineError> {
223        let kernel = require_kernel(self.kernels, jd)?;
224        let t = Time::from_ut1(jd);
225        let orientation = Orientation::with_nutation(&t, iau2000b_radians(t.tt()));
226        let position = apparent(kernel, &observe(kernel, SUN, &t)?, &t)?;
227        let equatorial = mxv(&orientation.m, &position);
228        let ra_hours = equatorial[1].atan2(equatorial[0]).rem_euclid(TAU) * 24.0 / TAU;
229        let apparent_hours = orientation.gast_hours - ra_hours + 12.0;
230        let mean_hours = (jd + 0.5).rem_euclid(1.0) * 24.0;
231        Ok(((apparent_hours - mean_hours + 12.0).rem_euclid(24.0) - 12.0) * 60.0)
232    }
233}
234
235/// A lunar month, by Julian day numbers.
236struct Month {
237    start: i64,
238    /// First day of the next month.
239    end: i64,
240    number: u8,
241    leap: bool,
242    year: i32,
243}
244
245/// A search the ephemeris cannot satisfy: only a corrupt kernel gets here.
246fn diverged(what: &str) -> EngineError {
247    EngineError::Ephemeris(SpkError::Format(format!("Chinese calendar search: {what}")))
248}
249
250/// The Chinese calendar around `instant`, or `None` when the loaded kernels do not cover
251/// the year or so of ephemeris it needs (about 14 months before the moment and 2 after).
252pub fn chinese_calendar(
253    kernels: &KernelSet,
254    instant: UtcInstant,
255) -> Result<Option<ChineseCalendar>, EngineError> {
256    match calendar(kernels, instant) {
257        Err(EngineError::OutOfRange { .. }) => Ok(None),
258        other => other.map(Some),
259    }
260}
261
262fn calendar(kernels: &KernelSet, instant: UtcInstant) -> Result<ChineseCalendar, EngineError> {
263    let sky = Sky { kernels };
264    let jd = instant.julian_day();
265
266    let mut solar_terms = Vec::new();
267    let mut term = sky.sun_before(f64::from(LICHUN), jd)?;
268    let mut longitude = LICHUN;
269    loop {
270        let jie = longitude % 30 == 15;
271        solar_terms.push(SolarTerm {
272            name: TERM_NAMES[usize::from(longitude / 15)].to_string(),
273            longitude,
274            jie,
275            utc: instant_of(term).isoformat(),
276        });
277        if jie && term > jd {
278            break;
279        }
280        if solar_terms.len() > 48 {
281            return Err(diverged("no solar term after the moment"));
282        }
283        longitude = (longitude + 15) % 360;
284        term = sky.sun_at(f64::from(longitude), term + 15.2)?;
285    }
286
287    let day = china_day(jd);
288    let solstice = sky.sun_before(270.0, jd)?;
289    let (mut months, mut end, mut next_solstice) = sky.sui(solstice)?;
290    if day - DATE_MARGIN < months[0].start {
291        let (mut earlier, ..) = sky.sui(sky.sun_at(270.0, solstice - TROPICAL_YEAR)?)?;
292        earlier.append(&mut months);
293        months = earlier;
294    }
295    while day + DATE_MARGIN >= end {
296        if months.len() > 40 {
297            return Err(diverged("no lunar year after the moment"));
298        }
299        let (mut later, later_end, later_solstice) = sky.sui(next_solstice)?;
300        months.append(&mut later);
301        (end, next_solstice) = (later_end, later_solstice);
302    }
303    let lunar_months = months
304        .into_iter()
305        .filter(|m| m.start <= day + DATE_MARGIN && m.end > day - DATE_MARGIN)
306        .map(|m| LunarMonth {
307            year: m.year,
308            month: m.number,
309            leap: m.leap,
310            start: date_of(m.start),
311            days: (m.end - m.start) as u8,
312        })
313        .collect();
314
315    Ok(ChineseCalendar {
316        equation_of_time: sky.equation_of_time(jd)?,
317        solar_terms,
318        lunar_months,
319    })
320}
321
322/// A lunar year of the calendar: its New Year and its months in order.
323#[derive(Debug, Clone, PartialEq, Eq)]
324pub(crate) struct LunarYear {
325    /// The Gregorian year its New Year falls in.
326    pub year: i32,
327    /// Julian day number of the New Year.
328    pub start: i64,
329    /// Lengths of its months in order (29 or 30), the leap month in its place.
330    pub days: Vec<u8>,
331    /// Number of the month the leap month repeats, 0 without one.
332    pub leap: u8,
333}
334
335/// The lunar years `first..=last`, from the ephemeris (for the embedded table and its test).
336#[cfg(test)]
337pub(crate) fn lunar_years(
338    kernels: &KernelSet,
339    first: i32,
340    last: i32,
341) -> Result<Vec<LunarYear>, EngineError> {
342    let sky = Sky { kernels };
343    let december = |year: i32| UtcInstant::from_civil(year, 12, 10, 0, 0, 0, 0).julian_day();
344    let mut solstice = sky.sun_at(270.0, december(first - 1))?;
345    let mut months = Vec::new();
346    loop {
347        let (mut sui, _, next) = sky.sui(solstice)?;
348        let done = sui.last().is_some_and(|m| m.year > last);
349        months.append(&mut sui);
350        if done {
351            break;
352        }
353        solstice = next;
354    }
355    Ok((first..=last)
356        .map(|year| {
357            let of_year: Vec<&Month> = months.iter().filter(|m| m.year == year).collect();
358            LunarYear {
359                year,
360                start: of_year[0].start,
361                days: of_year.iter().map(|m| (m.end - m.start) as u8).collect(),
362                leap: of_year.iter().find(|m| m.leap).map_or(0, |m| m.number),
363            }
364        })
365        .collect())
366}
367
368/// The moments (UT Julian dates) of the 12 jie of a Gregorian year, from 小寒 (285°) to
369/// 大雪 (255°), from the ephemeris (for the embedded table and its test).
370#[cfg(test)]
371pub(crate) fn jie_of_year(kernels: &KernelSet, year: i32) -> Result<[f64; 12], EngineError> {
372    let sky = Sky { kernels };
373    let january = UtcInstant::from_civil(year, 1, 6, 0, 0, 0, 0).julian_day();
374    let mut moments = [0.0; 12];
375    for (i, moment) in moments.iter_mut().enumerate() {
376        let longitude = f64::from((285 + 30 * i as u32) % 360);
377        *moment = sky.sun_at(longitude, january + 30.44 * i as f64)?;
378    }
379    Ok(moments)
380}
astroceleste_engine/chinese/calendar.rs

astroceleste_engine/chinese/
calendar.rs
