jupiters 0.0.3

Jupiter celestial simulation crate for the MilkyWay SolarSystem workspace
Documentation
use sciforge::hub::domain::common::constants::DEGREE_TO_RAD;

pub const SIDEREALDAYS: f64 = 35730.0;
pub const AXIALTILTDEG: f64 = 3.13;
pub const AXIALTILTRAD: f64 = AXIALTILTDEG * DEGREE_TO_RAD;
pub const PRECESSIONPERIODYEARS: f64 = 500000.0;

pub struct JupiterRotation {
    pub angularvelocityrads: f64,
    pub axialtiltrad: f64,
}

impl Default for JupiterRotation {
    fn default() -> Self {
        Self::new()
    }
}

impl JupiterRotation {
    pub fn new() -> Self {
        Self {
            angularvelocityrads: 2.0 * std::f64::consts::PI / SIDEREALDAYS,
            axialtiltrad: AXIALTILTRAD,
        }
    }

    pub fn surfacevelocityatlatitude(&self, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        self.angularvelocityrads * crate::JUPITEREQUATORIALRADIUS * latrad.cos()
    }

    pub fn centripetalaccelerationatlatitude(&self, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        let omega = self.angularvelocityrads;
        omega * omega * crate::JUPITEREQUATORIALRADIUS * latrad.cos()
    }

    pub fn coriolisparameter(&self, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        2.0 * self.angularvelocityrads * latrad.sin()
    }

    pub fn momentofinertia(&self) -> f64 {
        0.254 * crate::JUPITERMASS * crate::JUPITEREQUATORIALRADIUS * crate::JUPITEREQUATORIALRADIUS
    }

    pub fn rotationalkineticenergy(&self) -> f64 {
        0.5 * self.momentofinertia() * self.angularvelocityrads.powi(2)
    }

    pub fn angular_momentum(&self) -> f64 {
        self.momentofinertia() * self.angularvelocityrads
    }

    pub fn precession_rate_rad_per_year(&self) -> f64 {
        2.0 * std::f64::consts::PI / PRECESSIONPERIODYEARS
    }

    pub fn effective_tilt_rad(&self, juliancenturiessincej2000: f64) -> f64 {
        let rate = self.precession_rate_rad_per_year() / 100.0;
        self.axialtiltrad + rate * juliancenturiessincej2000 * 0.001
    }

    pub fn day_length_variation_due_to_tilt(&self, dayofyear: u32, latitudedeg: f64) -> f64 {
        let latrad = latitudedeg.to_radians();
        let orbitalperioddays = 4332.59;
        let declination = self.axialtiltrad
            * (2.0 * std::f64::consts::PI * (dayofyear as f64 - 81.0) / orbitalperioddays).sin();
        let coshourangle = -(latrad.tan() * declination.tan());
        if coshourangle > 1.0 {
            return 0.0;
        }
        if coshourangle < -1.0 {
            return SIDEREALDAYS / 3600.0;
        }
        2.0 * coshourangle.acos() * (SIDEREALDAYS / 3600.0) / (2.0 * std::f64::consts::PI)
    }
}