diff options
Diffstat (limited to 'src/solar.rs')
| -rw-r--r-- | src/solar.rs | 314 |
1 files changed, 314 insertions, 0 deletions
diff --git a/src/solar.rs b/src/solar.rs new file mode 100644 index 0000000..af52afc --- /dev/null +++ b/src/solar.rs @@ -0,0 +1,314 @@ +#![allow(non_upper_case_globals)] +#![allow(invalid_value)] +#![allow(non_snake_case)] +#![allow(dead_code)] + +// based on https://prospertimes.neocities.org/solarterms.js +use lazy_static::lazy_static; +use num::{cast, traits::AsPrimitive}; +use std::num::Wrapping as W; +use std::{ + f64::consts::PI, + mem::MaybeUninit, + ops::{Div, Rem}, +}; + +// this is so stupid. i should not have to put pub every single line +pub struct Term { + pub solar_term: &'static str, + pub ang: usize, + pub month: usize, + pub earliest_day: usize, + pub latest_day: usize, +} + +#[cfg_attr(not(target_family = "wasm"), derive(Debug))] +pub struct SexagenaryDate { + pub year: &'static str, + pub month: &'static str, + pub day: &'static str, + pub term: Option<(&'static str, usize)>, +} + +const UNIT_YR: usize = 1984; + +// cant use [].map(|(...)| Term {}) on consts +// this sucks, i should not have2 write design8ed initializers each&everytime +macro_rules! terms { + ($({ $s: expr, $a: expr, $m: expr, $e: expr, $l: expr }),* $(,)?) => { + [ + $(Term { + solar_term: $s, + ang: $a, + month: $m, + earliest_day: $e, + latest_day: $l, + }),* + ] + }; +} + +pub const TERMS: [Term; 24] = terms![ + {"小寒", 285, 1, 3, 7}, + {"大寒", 300, 1, 18, 22}, + {"立春", 315, 2, 2, 6}, + {"雨水", 330, 2, 16, 20}, + {"驚蟄", 345, 3, 3, 7}, + {"春分", 360, 3, 18, 22}, + {"淸明", 15, 4, 2, 6}, + {"穀雨", 30, 4, 18, 22}, + {"立夏", 45, 5, 3, 7}, + {"小滿", 60, 5, 19, 23}, + {"芒種", 75, 6, 3, 7}, + {"夏至", 90, 6, 19, 23}, + {"小暑", 105, 7, 5, 9}, + {"大暑", 120, 7, 20, 24}, + {"立秋", 135, 8, 5, 9}, + {"處暑", 150, 8, 21, 25}, + {"白露", 165, 9, 5, 9}, + {"秋分", 180, 9, 21, 25}, + {"寒露", 195, 10, 6, 10}, + {"霜降", 210, 10, 21, 25}, + {"立冬", 225, 11, 5, 9}, + {"小雪", 240, 11, 20, 24}, + {"大雪", 255, 12, 5, 9}, + {"冬至", 270, 12, 19, 23}, +]; + +fn compute_JD(mut y: f64, mut m: f64, d: f64) -> f64 { + if m <= 2. { + y -= 1.; + m += 12.; + } + (365.25 * (y + 4716.)).floor() + (30.6001 * (m + 1.)).floor() + d + - 13. + - 1524.5 +} + +macro_rules! JD { + ($y: expr, $m: expr, $d: expr) => { + compute_JD($y as f64, $m as f64, $d as f64) + }; +} + +// dont ask me how this works +fn compute_ang(jd: f64) -> f64 { + let T = (jd - 2451545.) / 36525.; + let Lmean = (280.46646 + (36000.76983 * T) + (0.0003032 * T * T)) % 360.; + let M = (357.52911 + (35999.05029 * T) - (0.0001537 * T * T)) % 360.; + let C = ((1.914602 - (0.004817 * T) - (0.000014 * T * T)) + * (M * PI / 180.).sin()) + + ((0.019993 - (0.000101 * T)) * (2. * M * PI / 180.).sin()) + + (0.000289 * (3. * M * PI / 180.).sin()); + let Ltrue = Lmean + C; + let Lapp = Ltrue + - 0.00569 + - (0.00478 * ((125.04 - 1934.136 * T) * PI / 180.).sin()); + Lapp +} + +fn bisect<T>( + lo: &mut f64, + hi: &mut f64, + base: T, + scale: T, + target: f64, + y: usize, + m: f64, +) where + T: Into<f64> + Copy, +{ + while lo <= hi { + let mid = ((*lo + *hi) / 2.).floor(); + if compute_ang(JD!(y, m, base.into() + mid * scale.into())) < target { + *lo = mid + 1.; + } else { + *hi = mid - 1.; + } + } +} + +fn compute_solarterm_day(i: usize, y: usize) -> f64 { + //double solarterm_day, test_day, test_hour, test_minute, test_long; + let mut solarterm_day; + let m: f64 = TERMS[i].month as f64; + let mut epd: f64 = TERMS[i].earliest_day as f64; // earliest possible day + let mut lpd: f64 = TERMS[i].latest_day as f64; // latest possible day + let mut eph: f64 = 0.; // earliest possible hour + let mut lph: f64 = 23.; // latest possible hour + let mut epm: f64 = 0.; // earliest possible minute + let mut lpm: f64 = 59.; // latest possible minute + let target_long: f64 = TERMS[i].ang as f64; + bisect(&mut epd, &mut lpd, 0, 1, target_long, y, m); + solarterm_day = lpd; + bisect( + &mut eph, + &mut lph, + solarterm_day, + 1. / 24., + target_long, + y, + m, + ); + solarterm_day += lph / 24.; + bisect( + &mut epm, + &mut lpm, + solarterm_day, + 1. / 1440., + target_long, + y, + m, + ); + solarterm_day + (lpm + 1.) / 1440. +} + +#[inline(always)] +fn align<T>(a: T, to: T, zero: T) -> usize +where + T: cast::AsPrimitive<isize>, +{ + let (a, to, zero) = (a.as_(), to.as_(), zero.as_()); + let tmp = a + (to - zero); + (tmp - to * (tmp >= to) as isize) as usize +} + +#[inline(always)] +fn mod2ganzhi<T>(g: T, z: T) -> usize +where + T: cast::AsPrimitive<isize>, +{ + let (g, z) = (g.as_(), z.as_()); + // chinese remainder theorem + // x ≡ g (mod 10) + // x ≡ z (mod 12) + // 10x + 12y = gcd(10,12) = 2 + // x=-1, y=1 + // X=(g-z)/gcd(10,12)=(g-z)/2 + // 10xX + 12yX = 2X + // => -10X + 12X = 2X + // => -5(g-z) + 6(g-z) = g-z + // x = g-(-5(g-z)) (case 1) = 6(g-z)+z (case 2) + // x = g+5(g-z) = 6g-5z + // dont ask me how the math works idk + align(6 * g - 5 * z, 60, 0) +} + +fn stday(i: usize, y: usize) -> f64 { + const YEARS: usize = 200; + // rust doesnt have (*a)[n], | @ least ill have2 use a cr8 4 it + // also lazy_static doesnt work with mutables (im not using mutex) + assert!(y < UNIT_YR + YEARS && UNIT_YR <= y); + static mut STDAYS: Vec<f64> = vec![]; + static mut INIT: bool = false; + unsafe { + if !INIT { + STDAYS = vec![0.; YEARS * TERMS.len()]; + INIT = true; + } + let idx = y - UNIT_YR; + let ret = &mut STDAYS[idx * TERMS.len() + i]; + (if int(*ret) != 0 { + *ret + } else { + let d = compute_solarterm_day(i, y); + *ret = d; + d + }) + .floor() + } +} + +lazy_static! { + static ref GANZHIS: Vec<String> = { + let mut tmp: Vec<String> = vec![]; + tmp.reserve(60); + (0..60).for_each(|i| { + // looks so much worse than using format!() + // cant be bothered to benchmark + let mut s = gan[i % 10].to_string(); + s.push(zhi[i % 12]); + tmp.push(s); + }); + tmp + }; +} +const gan: [char; 10] = + ['甲', '乙', '丙', '丁', '戊', '己', '庚', '辛', '壬', '癸']; +const zhi: [char; 12] = [ + '子', '丑', '寅', '卯', '辰', '巳', '午', '未', '申', '酉', '戌', '亥', +]; +pub fn ganzhi(i: usize) -> &'static str { + &GANZHIS[i] +} + +#[inline(always)] +fn int(x: f64) -> usize { + x as usize +} + +lazy_static! { + static ref JIAZI: usize = int(JD!(UNIT_YR, 1, 31)); +} +pub fn solar(y: usize, m: usize, d: usize) -> SexagenaryDate { + let jdf = JD!(y, m, d); + let jd = int(jdf); + let a = compute_ang(jdf); + let ygz = (y + - UNIT_YR + - ((m <= 2 && a < 316.) && (m < 2 || d < int(stday(2, y)))) as usize) + % 60; + let rem = a % 15.; + let div = int(a.div(15.).floor()); + let mut dz = align((div + 1) / 2, 12, 9); + let mut termb = rem > 14.; + let mut term: usize = unsafe { MaybeUninit::uninit().assume_init() }; + if termb { + term = align(div, 24, 18); + let termday = stday(term, y); + termb = d == int(termday); + if (div & 1) == 0 { + dz += termb as usize; + dz = align(dz, 12, 12); + } + } + let mut tmp = ygz.rem(5); + tmp = tmp * 2 + 2; + tmp -= 10 * (tmp == 10) as usize; + let mgz = mod2ganzhi(tmp + align(dz, 12, 2).rem(10), dz); + let dgz = (jd - *JIAZI).rem(60); + SexagenaryDate { + year: ganzhi(ygz), + month: ganzhi(mgz), + day: ganzhi(dgz), + term: termb.then(|| (TERMS[term].solar_term, term)), + } +} + +// https://en.wikipedia.org/wiki/Determination_of_the_day_of_the_week#Disparate_variation +pub fn dow(y: usize, mut m: usize, d: usize) -> (&'static str, &'static str) { + let Y = y - (m <= 2) as usize; + let y2 = Y % 100; + let c = Y / 100; + m += 9; + m -= 12 * (m >= 12) as usize; + let a = (d + (26 * (m + 1) - 2) / 10 + y2 + y2 / 4 + c / 4 - 2 * c) % 7; + [ + ("日", "Sun"), + ("月", "Mon"), + ("火", "Tue"), + ("水", "Wed"), + ("木", "Thu"), + ("金", "Fri"), + ("土", "Sat"), + ][a] +} + +#[cfg(not(target_family = "wasm"))] +fn main() { + (1..=12).for_each(|m| { + (1..28).for_each(|d| { + println!("{:?}", solar(2026, m, d)); + }); + }); +} |