solarsystems 0.0.1

N-body solar system engine — gravitational dynamics, orbital mechanics, perturbations, event detection, and full celestial orchestration
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
89
90
91
92
93
94
95
96
97

pub const SECONDS_PER_MINUTE: f64 = 60.0;
pub const SECONDS_PER_HOUR: f64 = 3_600.0;
pub const SECONDS_PER_DAY: f64 = 86_400.0;
pub const SECONDS_PER_JULIAN_YEAR: f64 = 365.25 * SECONDS_PER_DAY;
pub const SECONDS_PER_SIDEREAL_YEAR: f64 = 365.256_363_004 * SECONDS_PER_DAY;
pub const DAYS_PER_JULIAN_CENTURY: f64 = 36_525.0;
pub const J2000_EPOCH_JD: f64 = 2_451_545.0;

#[derive(Clone, Debug)]
pub struct EphemerisTime {
    pub epoch_jd: f64,
    pub elapsed_s: f64,
}

impl EphemerisTime {
    pub fn j2000() -> Self {
        Self {
            epoch_jd: J2000_EPOCH_JD,
            elapsed_s: 0.0,
        }
    }

    pub fn from_jd(jd: f64) -> Self {
        Self {
            epoch_jd: J2000_EPOCH_JD,
            elapsed_s: (jd - J2000_EPOCH_JD) * SECONDS_PER_DAY,
        }
    }

    pub fn advance(&mut self, dt: f64) {
        self.elapsed_s += dt;
    }

    pub fn current_jd(&self) -> f64 {
        self.epoch_jd + self.elapsed_s / SECONDS_PER_DAY
    }

    pub fn centuries_since_j2000(&self) -> f64 {
        self.elapsed_s / (DAYS_PER_JULIAN_CENTURY * SECONDS_PER_DAY)
    }

    pub fn years_elapsed(&self) -> f64 {
        self.elapsed_s / SECONDS_PER_JULIAN_YEAR
    }

    pub fn days_elapsed(&self) -> f64 {
        self.elapsed_s / SECONDS_PER_DAY
    }
}

pub fn julian_date_to_calendar(jd: f64) -> (i32, u32, u32, u32, u32, f64) {
    let jd = jd + 0.5;
    let z = jd.floor() as i64;
    let f = jd - z as f64;
    let a = if z < 2_299_161 {
        z
    } else {
        let alpha = ((z as f64 - 1_867_216.25) / 36_524.25).floor() as i64;
        z + 1 + alpha - alpha / 4
    };
    let b = a + 1524;
    let c = ((b as f64 - 122.1) / 365.25).floor() as i64;
    let d = (365.25 * c as f64).floor() as i64;
    let e = ((b - d) as f64 / 30.6001).floor() as i64;
    let day = (b - d - (30.6001 * e as f64).floor() as i64) as u32;
    let month = if e < 14 { e - 1 } else { e - 13 } as u32;
    let year = if month > 2 { c - 4716 } else { c - 4715 } as i32;
    let hours_total = f * 24.0;
    let hour = hours_total.floor() as u32;
    let minutes_total = (hours_total - hour as f64) * 60.0;
    let minute = minutes_total.floor() as u32;
    let second = (minutes_total - minute as f64) * 60.0;
    (year, month, day, hour, minute, second)
}

pub fn calendar_to_julian_date(
    year: i32,
    month: u32,
    day: u32,
    hour: u32,
    minute: u32,
    second: f64,
) -> f64 {
    let (y, m) = if month <= 2 {
        (year - 1, month + 12)
    } else {
        (year, month)
    };
    let a = y / 100;
    let b = 2 - a + a / 4;
    let jd = (365.25 * (y + 4716) as f64).floor()
        + (30.6001 * (m + 1) as f64).floor()
        + day as f64
        + b as f64
        - 1524.5;
    jd + (hour as f64 + minute as f64 / 60.0 + second / 3600.0) / 24.0
}