Space Math in Rust

math.rs

#![allow(unused)]

use chrono::prelude::*;
use std::f64::consts::PI;

type Float = f64;

#[derive(Copy, Clone)]
pub struct GpsCoordinate {
    pub latitude: Angle,
    pub longitude: Angle,
}

#[derive(Copy, Clone)]
pub struct Dms {
    pub degrees: Float,
    pub minutes: Float,
    pub seconds: Float,
}

#[derive(Copy, Clone)]
pub struct Hms {
    pub hours: Float,
    pub minutes: Float,
    pub seconds: Float,
}

#[derive(Copy, Clone)]
pub struct Angle {
    pub rads: Float,
}

#[derive(Copy, Clone)]
pub struct SphericalCoordinate {
    pub r: Float,
    pub ra: Angle,
    pub decl: Angle,
}

#[derive(Copy, Clone)]
pub struct RectangularCoordinate {
    pub x: Angle,
    pub y: Angle,
    pub z: Angle,
}

#[derive(Copy, Clone)]
pub struct HorizontalCoordinate {
    pub altitude: Angle,
    pub azimuth: Angle,
}

impl Angle {
    pub fn from_rads(rads: Float) -> Self {
        return Angle { rads: rads };
    }

    pub fn from_deg(degs: Float) -> Self {
        return Angle {
            rads: (degs / 180.0) * PI,
        };
    }

    pub fn as_rads(&self) -> Float {
        return self.rads;
    }

    pub fn as_deg(&self) -> Float {
        return (self.rads * 180.0) / PI;
    }

    pub fn as_dms(&self) -> Dms {
        let deg = self.as_deg();
        let min = (deg - deg.floor()) * 60.0;
        let sec = (min - min.floor()) * 60.0;

        return Dms {
            degrees: self.as_deg(),
            minutes: min,
            seconds: sec,
        };
    }

    pub fn as_hms(&self) -> Hms {
        let degrees = self.as_deg();
        let hours = degrees / 15.0;
        let minutes = (hours - hours.floor()) * 60.0;
        let seconds = (minutes - minutes.floor()) * 60.0;

        return Hms {
            hours: hours,
            minutes: minutes,
            seconds: seconds,
        };
    }

    pub fn rev(&self) -> Self {
        return Angle::from_rads(self.rads - (self.rads / (2.0 * PI)).floor() * (2.0 * PI));
    }
}

pub struct SpaceMath {}

impl SpaceMath {
    pub fn sum_and_exponentially_multiply(values: Vec<Float>, multiplier: Float) -> Float {
        let mut result = 0.0;
        for i in 0..values.len() {
            if i == 0 {
                result += values.get(i).unwrap();
            } else {
                result += values.get(i).unwrap() * multiplier.powf(i as Float);
            }
        }
        return result;
    }

    pub fn rads(deg: Float) -> Float {
        return (deg / 180.0) * PI;
    }

    pub fn time_to_dec(hours: Float, minutes: Float, seconds: Float) -> Float {
        return hours + (minutes / 60.0) + (seconds / 3600.0);
    }

    pub fn time_to_angle(hours: Float, minutes: Float, seconds: Float) -> Angle {
        // NOTE: according this website, we should consider
        // 15.04107 instead of 15.0 http://www.stjarnhimlen.se/comp/riset.html
        return Angle::from_deg(SpaceMath::time_to_dec(hours, minutes, seconds) * 15.0);
    }

    pub fn instant_to_days_julian(instant: NaiveDateTime) -> Float {
        // Add the hour decimal.
        let day = instant.day() as Float
            + SpaceMath::time_to_dec(
                instant.hour() as Float,
                instant.minute() as Float,
                instant.second() as Float,
            ) / 24.0;

        let mut month = instant.month() as Float;
        let mut year = instant.year() as Float;

        // If we are calculating for january or february, consider it the "13th and 14th"
        // months.
        if month == 1.0 || month == 2.0 {
            year = year - 1.0;
            month = month + 12.0;
        }

        let a = (year / 100.0).floor();
        let b = 2.0 - a + (a / 4.0).floor();
        return (365.25 * (year + 4716.0)).floor() + (30.6001 * (month + 1.0)).floor() + day + b
            - 1524.5;
    }

    pub fn julian_to_day_number(days_julian: Float) -> Float {
        return days_julian - 2451543.5;
    }

    pub fn julian_to_century_time(julian_date: Float) -> Float {
        return (julian_date - 2451545.0) / 36525.0;
    }

    pub fn julian_millennia_from_epoch_2000(julian_days: Float) -> Float {
        return (julian_days - 2451545.0) / 365250.0;
    }

    pub fn julian_centuries_from_epoch_2000(julian_days: Float) -> Float {
        return 10.0 * SpaceMath::julian_millennia_from_epoch_2000(julian_days);
    }

    // This gets the mean longitude of the sun at an instant using keplerian formulas
    pub fn sol_l(d: Float) -> Angle {
        let w = Angle::from_deg(SpaceMath::sum_and_exponentially_multiply(
            vec![282.9404, 4.70935E-5],
            d,
        ))
        .as_deg();

        let m = Angle::from_deg(SpaceMath::sum_and_exponentially_multiply(
            vec![356.0470, 0.9856002585],
            d,
        ))
        .rev()
        .as_deg();

        return Angle::from_deg(w + m).rev();
    }

    pub fn calculate_hour_angle<T: TimeZone>(
        instant: DateTime<T>,
        longitude: Angle,
        ra: Angle,
    ) -> Angle {
        // Get the current mean longitude of the sun
        let instant_utc = instant.naive_utc();
        let jd = SpaceMath::instant_to_days_julian(instant_utc);
        let d = SpaceMath::julian_to_day_number(jd);

        let gmst0 = Angle::from_deg(SpaceMath::sol_l(d).as_deg() + 180.0)
            .rev()
            .as_deg()
            / 15.0; // Hours
        let sidtime = gmst0
            + SpaceMath::time_to_dec(
                instant.hour() as Float,
                instant.minute() as Float,
                instant.second() as Float,
            )
            + (longitude.as_deg() / 15.0);
        let ha = sidtime - (ra.as_deg() / 15.0);
        return Angle::from_deg(ha * 15.0);
    }

    pub fn sphere_to_rect(sphere: SphericalCoordinate) -> RectangularCoordinate {
        return RectangularCoordinate {
            x: Angle::from_rads(sphere.r * sphere.ra.as_rads().cos() * sphere.decl.as_rads().cos()),
            y: Angle::from_rads(sphere.r * sphere.ra.as_rads().sin() * sphere.decl.as_rads().cos()),
            z: Angle::from_rads(sphere.r * sphere.decl.as_rads().sin()),
        };
    }

    pub fn rect_to_sphere(rect: RectangularCoordinate) -> SphericalCoordinate {
        return SphericalCoordinate {
            r: (rect.x.as_rads().powf(2.0)
                + rect.y.as_rads().powf(2.0)
                + rect.z.as_rads().powf(2.0))
            .sqrt(),
            ra: Angle::from_rads(rect.y.as_rads().atan2(rect.x.as_rads())).rev(),
            decl: Angle::from_rads(
                rect.z
                    .as_rads()
                    .atan2((rect.x.as_rads().powf(2.0) + rect.y.as_rads().powf(2.0)).sqrt()),
            ),
        };
    }

    pub fn rotate_around_x(rect: RectangularCoordinate, oblecl: Angle) -> RectangularCoordinate {
        return RectangularCoordinate {
            x: rect.x,
            y: Angle::from_rads(
                rect.y.as_rads() * oblecl.as_rads().cos()
                    - rect.z.as_rads() * oblecl.as_rads().sin(),
            ),
            z: Angle::from_rads(
                rect.y.as_rads() * oblecl.as_rads().sin()
                    + rect.z.as_rads() * oblecl.as_rads().cos(),
            ),
        };
    }

    pub fn sphere_to_horizontal<T: TimeZone>(
        gps: GpsCoordinate,
        target: SphericalCoordinate,
        instant: DateTime<T>,
    ) -> HorizontalCoordinate {
        let ha = SpaceMath::calculate_hour_angle(instant, gps.longitude, target.ra);
        let x = ha.as_rads().cos() * target.decl.as_rads().cos();
        let y = ha.as_rads().sin() * target.decl.as_rads().cos();
        let z = target.decl.as_rads().sin();
        let d90 = SpaceMath::rads(90.0);
        let xhor =
            x * (d90 - gps.latitude.as_rads()).cos() - z * (d90 - gps.latitude.as_rads()).sin();
        let yhor = y;
        let zhor =
            x * (d90 - gps.latitude.as_rads()).sin() + z * (d90 - gps.latitude.as_rads()).cos();
        let azimuth = yhor.atan2(xhor) + PI;
        let altitude = zhor.atan2((xhor.powf(2.0) + yhor.powf(2.0)).sqrt());

        return HorizontalCoordinate {
            altitude: Angle::from_rads(altitude),
            azimuth: Angle::from_rads(azimuth),
        };
    }
}

#[cfg(test)]
mod test_coordinates {
    use super::*;

    fn round(x: Float, decimals: u32) -> Float {
        let y = 10i32.pow(decimals) as Float;
        (x * y).round() / y
    }

    #[test]
    pub fn test_spherical_to_ecliptic() {
        let sol = SphericalCoordinate {
            ra: Angle::from_deg(90.0),
            decl: Angle { rads: 0.0 },
            r: 1.0,
        };

        // Verify equatorial conversion works
        let xyz = SpaceMath::sphere_to_rect(sol);
        assert_eq!(round(xyz.x.as_rads(), 8), 0.0);
        assert_eq!(round(xyz.y.as_rads(), 8), 1.0);
        assert_eq!(xyz.z.as_rads(), (0.0 as Float).sin());

        // Now verify ecliptic rotation works
        let xyz_ecliptic = SpaceMath::rotate_around_x(xyz, Angle::from_deg(23.4));

        assert_eq!(round(xyz_ecliptic.x.as_rads(), 8), 0.0);
        assert_eq!(round(xyz_ecliptic.y.as_rads(), 4), 0.9178);
        assert_eq!(round(xyz_ecliptic.z.as_rads(), 4), 0.3971);

        // Convert ecliptic to spherical
        let spherical = SpaceMath::rect_to_sphere(xyz_ecliptic);
        assert_eq!(spherical.r.round(), 1.0);
        assert_eq!(spherical.ra.as_deg().round(), 90.0);
        assert_eq!(round(spherical.decl.as_deg(), 2), 23.4);
    }
}