suncalc 0.4.0

//! The `sun` crate is a library for calculating the position of the sun.
//! It is a port of the `JavaScript` library
//! [suncalc](https://github.com/mourner/suncalc).
//!
//! # Example
//!
//! ```
//! let unixtime = suncalc::Timestamp(1362441600000);
//! let lat = 48.0;
//! let lon = 9.0;
//! let moon_pos = suncalc::get_position(unixtime,lat,lon);
//! let az  = moon_pos.azimuth.to_degrees();
//! let alt = moon_pos.altitude.to_degrees();
//! println!("The position of the sun is {}/{}", az, alt);
//! ```

use std::f64::consts::PI;

// base definitions

/// Holds the [azimuth](https://en.wikipedia.org/wiki/Azimuth)
/// and [altitude](https://en.wikipedia.org/wiki/Horizontal_coordinate_system)
/// angles of the sun position.
#[derive(Debug, Clone, Copy)]
pub struct Position {
    pub azimuth: f64,
    pub altitude: f64,
    pub distance: Option<f64>,
    pub parallactic_angle: Option<f64>,
}

#[derive(Debug, Clone, Copy)]
struct Coords {
    pub right_ascension: f64,
    pub declination: f64,
    pub distance: Option<f64>,
}

/// Holds illuminated _fraction_ of the moon, the _phase_, and the _angle_.
#[derive(Debug, Clone, Copy)]
pub struct Illumination {
    pub fraction: f64,
    pub phase: f64,
    pub angle: f64,
}

/// UNIX timestamp in milliseconds
#[derive(Debug, PartialEq, Eq, Clone, Copy)]
pub struct Timestamp(pub i64);

// date/time constants and conversions
const MILLISECONDS_PER_DAY: u32 = 1000 * 60 * 60 * 24;
const J1970: u32 = 2_440_588;
const J2000: u32 = 2_451_545;
const TO_RAD: f64 = PI / 180.0;
const OBLIQUITY_OF_EARTH: f64 = 23.4397 * TO_RAD;
const PERIHELION_OF_EARTH: f64 = 102.9372 * TO_RAD;

fn to_julian(timestamp: Timestamp) -> f64 {
    timestamp.0 as f64 / (MILLISECONDS_PER_DAY as f64) - 0.5 + J1970 as f64
}

fn from_julian(j: f64) -> Timestamp {
    Timestamp(((j + 0.5 - J1970 as f64) * MILLISECONDS_PER_DAY as f64).round() as i64)
}

fn to_days(timestamp: Timestamp) -> f64 {
    to_julian(timestamp) - J2000 as f64
}

// general calculations for position

fn right_ascension(l: f64, b: f64) -> f64 {
    (l.sin() * OBLIQUITY_OF_EARTH.cos() - b.tan() * OBLIQUITY_OF_EARTH.sin()).atan2(l.cos())
}

fn declination(l: f64, b: f64) -> f64 {
    (b.sin() * OBLIQUITY_OF_EARTH.cos() + b.cos() * OBLIQUITY_OF_EARTH.sin() * l.sin()).asin()
}

fn azimuth(h: f64, phi: f64, dec: f64) -> f64 {
    h.sin().atan2(h.cos() * phi.sin() - dec.tan() * phi.cos())
}

fn altitude(h: f64, phi: f64, dec: f64) -> f64 {
    (phi.sin() * dec.sin() + phi.cos() * dec.cos() * h.cos()).asin()
}

fn sidereal_time(d: f64, lw: f64) -> f64 {
    (280.16 + 360.985_623_5 * d).to_radians() - lw
}

fn astro_refraction(h: f64) -> f64 {
    let hh = if h < 0.0 { 0.0 } else { h };

    0.0002967 / (hh + 0.00312536 / (hh + 0.08901179)).tan()
}

// general sun calculations

fn solar_mean_anomaly(d: f64) -> f64 {
    (357.5291 + 0.985_600_28 * d).to_radians()
}

fn equation_of_center(m: f64) -> f64 {
    (1.9148 * (1.0 * m).sin() + 0.02 * (2.0 * m).sin() + 0.0003 * (3.0 * m).sin()).to_radians()
}

fn ecliptic_longitude(m: f64) -> f64 {
    m + equation_of_center(m) + PERIHELION_OF_EARTH + PI
}

// calculations for sun times

const J0: f64 = 0.0009;

fn julian_cycle(d: f64, lw: f64) -> f64 {
    (d - J0 - lw / (2.0 * PI)).round()
}

fn approx_transit(h_t: f64, lw: f64, n: f64) -> f64 {
    J0 + (h_t + lw) / (2.0 * PI) + n
}

fn solar_transit_j(ds: f64, m: f64, l: f64) -> f64 {
    J2000 as f64 + ds + 0.0053 * m.sin() - 0.0069 * (2.0 * l).sin()
}

fn hour_angle(h: f64, phi: f64, dec: f64) -> f64 {
    ((h.sin() - phi.sin() * dec.sin()) / (phi.cos() * dec.cos())).acos()
}

fn observer_angle(height: f64) -> f64 {
    -2.076 * height.sqrt() / 60.0
}

fn get_set_j(h: f64, lw: f64, phi: f64, dec: f64, n: f64, m: f64, l: f64) -> f64 {
    let w = hour_angle(h, phi, dec);
    let a = approx_transit(w, lw, n);

    solar_transit_j(a, m, l)
}

fn sun_coords(d: f64) -> Coords {
    let m = solar_mean_anomaly(d);
    let l = ecliptic_longitude(m);

    Coords {
        right_ascension: right_ascension(l, 0.0),
        declination: declination(l, 0.0),
        distance: None,
    }
}

/// Calculates the sun position for a given date and latitude/longitude.
/// The angles are calculated as [radians](https://en.wikipedia.org/wiki/Radian).
///
/// * `unixtime`  - [unix time](https://en.wikipedia.org/wiki/Unix_time) in milliseconds.
/// * `lat`       - [latitude](https://en.wikipedia.org/wiki/Latitude) in degrees.
/// * `lon`       - [longitude](https://en.wikipedia.org/wiki/Longitude) in degrees.
/// calculates the sun position for a given date and latitude/longitude
pub fn get_position(timestamp: Timestamp, lat: f64, lon: f64) -> Position {
    let lw = -lon.to_radians();
    let phi = lat.to_radians();
    let d = to_days(timestamp);
    let m = solar_mean_anomaly(d);
    let l = ecliptic_longitude(m);
    let dec = declination(l, 0.0);
    let ra = right_ascension(l, 0.0);
    let h = sidereal_time(d, lw) - ra;

    Position {
        azimuth: azimuth(h, phi, dec),
        altitude: altitude(h, phi, dec),
        distance: None,
        parallactic_angle: None,
    }
}

#[derive(Debug, Clone, Copy)]
pub struct SunTimes {
    pub solar_noon: Timestamp,
    pub nadir: Timestamp,

    pub sunrise: Timestamp,
    pub sunset: Timestamp,

    pub sunrise_end: Timestamp,
    pub sunset_start: Timestamp,

    pub dawn: Timestamp,
    pub dusk: Timestamp,

    pub nautical_dawn: Timestamp,
    pub nautical_dusk: Timestamp,

    pub night_end: Timestamp,
    pub night: Timestamp,

    pub golden_hour_end: Timestamp,
    pub golden_hour: Timestamp,
}

struct SunTime(f64);

impl SunTime {
    pub const SUNRISE_SUNSET: SunTime = SunTime(-0.833);
    pub const SUNRISE_END_SUNSET_START: SunTime = SunTime(-0.3);
    pub const DAWN_DUSK: SunTime = SunTime(-6.0);
    pub const NAUTICAL_DAWN_NAUTICAL_DUSK: SunTime = SunTime(-12.0);
    pub const NIGHT_END_NIGHT: SunTime = SunTime(-18.0);
    pub const GOLDEN_HOUR_END_GOLDEN_HOUR: SunTime = SunTime(6.0);

    pub fn calc(
        &self,
        j_noon: f64,
        dh: f64,
        lw: f64,
        phi: f64,
        dec: f64,
        n: f64,
        m: f64,
        l: f64,
    ) -> (Timestamp, Timestamp) {
        let h0 = (self.0 + dh).to_radians();
        let j_set = get_set_j(h0, lw, phi, dec, n, m, l);
        let j_rise = j_noon - (j_set - j_noon);
        (from_julian(j_rise), from_julian(j_set))
    }
}

pub fn get_times(timestamp: Timestamp, lat: f64, lon: f64, height: Option<f64>) -> SunTimes {
    let height = height.unwrap_or(0.0);

    let lw = (-lon).to_radians();
    let phi = lat.to_radians();

    let dh = observer_angle(height);

    let d = to_days(timestamp);
    let n = julian_cycle(d, lw);
    let ds = approx_transit(0.0, lw, n);

    let m = solar_mean_anomaly(ds);
    let l = ecliptic_longitude(m);
    let dec = declination(l, 0.0);

    let j_noon = solar_transit_j(ds, m, l);

    let (sunrise, sunset) = SunTime::SUNRISE_SUNSET.calc(j_noon, dh, lw, phi, dec, n, m, l);
    let (sunrise_end, sunset_start) =
        SunTime::SUNRISE_END_SUNSET_START.calc(j_noon, dh, lw, phi, dec, n, m, l);
    let (dawn, dusk) = SunTime::DAWN_DUSK.calc(j_noon, dh, lw, phi, dec, n, m, l);
    let (nautical_dawn, nautical_dusk) =
        SunTime::NAUTICAL_DAWN_NAUTICAL_DUSK.calc(j_noon, dh, lw, phi, dec, n, m, l);
    let (night_end, night) = SunTime::NIGHT_END_NIGHT.calc(j_noon, dh, lw, phi, dec, n, m, l);
    let (golden_hour_end, golden_hour) =
        SunTime::GOLDEN_HOUR_END_GOLDEN_HOUR.calc(j_noon, dh, lw, phi, dec, n, m, l);

    SunTimes {
        solar_noon: from_julian(j_noon),
        nadir: from_julian(j_noon - 0.5),
        sunrise,
        sunset,
        sunrise_end,
        sunset_start,
        dawn,
        dusk,
        nautical_dawn,
        nautical_dusk,
        night_end,
        night,
        golden_hour_end,
        golden_hour,
    }
}

// general moon calculations, based on http://aa.quae.nl/en/reken/hemelpositie.html formulas

fn lunar_mean_anomaly(d: f64) -> f64 {
    (134.963 + 13.064993 * d).to_radians()
}

fn lunar_ecliptic_longitude(d: f64) -> f64 {
    (218.316 + 13.176396 * d).to_radians()
}

fn lunar_mean_distance(d: f64) -> f64 {
    (93.272 + 13.229350 * d).to_radians()
}

fn moon_coords(d: f64) -> Coords {
    let l = lunar_ecliptic_longitude(d);
    let m = lunar_mean_anomaly(d);
    let f = lunar_mean_distance(d);

    let lng = l + TO_RAD * 6.289 * m.sin();
    let lat = TO_RAD * 5.128 * f.sin();
    let distance = 385001.0 - 20905.0 * m.cos(); // in km

    Coords {
        right_ascension: right_ascension(lng, lat),
        declination: declination(lng, lat),
        distance: Some(distance),
    }
}

/// Calculates the moon position for a given date and latitude/longitude
pub fn moon_pos(timestamp: Timestamp, lat: f64, lon: f64) -> Position {
    let lw = TO_RAD * -lon;
    let phi = TO_RAD * lat;
    let d = to_days(timestamp);

    let c = moon_coords(d);

    let h = sidereal_time(d, lw) - c.right_ascension;
    let mut alt = altitude(h, phi, c.declination);
    alt = alt + astro_refraction(alt);

    // formula 14.1 of "Astronomical Algorithms" 2nd edition by Jean Meeus (Willmann-Bell, Richmond) 1998.
    let pa = h
        .sin()
        .atan2(phi.tan() * c.declination.cos() - c.declination.sin() * h.cos());

    Position {
        azimuth: azimuth(h, phi, c.declination),
        altitude: alt,
        distance: c.distance,
        parallactic_angle: Some(pa),
    }
}

/// Calculates the moon illumination, phase, and angle for a given date
pub fn moon_illumination(timestamp: Timestamp) -> Illumination {
    let d = to_days(timestamp);
    let s = sun_coords(d);
    let m = moon_coords(d);
    let a = lunar_mean_anomaly(d);

    let distance = 385001.0 - 20905.0 * a.cos(); // distance to the moon in km

    let sdist = 149598000 as f64;

    let phi = (s.declination.sin() * m.declination.sin()
        + s.declination.cos()
            * m.declination.cos()
            * (s.right_ascension - m.right_ascension).cos())
    .acos();

    let inc = (sdist * phi.sin()).atan2(distance - sdist * phi.cos());
    let angle = (s.declination.cos() * (s.right_ascension - m.right_ascension).sin()).atan2(
        s.declination.sin() * m.declination.cos()
            - s.declination.cos()
                * (m.declination).sin()
                * (s.right_ascension - m.right_ascension).cos(),
    );

    let sign = if angle < 0.0 { -1.0 } else { 1.0 };

    Illumination {
        fraction: (1.0 + inc.cos()) / 2.0,
        phase: 0.5 + 0.5 * inc * sign / PI,
        angle: angle,
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_abs_diff_eq;

    // 2013-03-05 UTC
    const DATE: Timestamp = Timestamp(1362441600000);
    const LAT: f64 = 50.5;
    const LON: f64 = 30.5;

    #[test]
    fn test_pos() {
        let moon_pos = get_position(DATE, LAT, LON);

        assert_abs_diff_eq!(-2.5003175907168385, moon_pos.azimuth);
        assert_abs_diff_eq!(-0.7000406838781611, moon_pos.altitude);
    }

    macro_rules! assert_eq_secs {
        ($a:expr, $b:expr) => {
            assert_eq!($a.0 / 1000, $b.0 / 1000);
        };
    }

    #[test]
    fn test_times() {
        // 2013-03-05T10:10:57Z
        let solar_noon = Timestamp(1362478257000);
        // 2013-03-04T22:10:57Z
        let nadir = Timestamp(1362435057000);
        // 2013-03-05T04:34:56Z
        let sunrise = Timestamp(1362458096000);
        // 2013-03-05T15:46:57Z
        let sunset = Timestamp(1362498417000);
        // 2013-03-05T04:38:19Z
        let sunrise_end = Timestamp(1362458299000);
        // 2013-03-05T15:43:34Z
        let sunset_start = Timestamp(1362498214000);
        // 2013-03-05T04:02:17Z
        let dawn = Timestamp(1362456137000);
        // 2013-03-05T16:19:36Z
        let dusk = Timestamp(1362500376000);
        // 2013-03-05T03:24:31Z
        let nautical_dawn = Timestamp(1362453871000);
        // 2013-03-05T16:57:22Z
        let nautical_dusk = Timestamp(1362502642000);
        // 2013-03-05T02:46:17Z
        let night_end = Timestamp(1362451577000);
        // 2013-03-05T17:35:36Z
        let night = Timestamp(1362504936000);
        // 2013-03-05T05:19:01Z
        let golden_hour_end = Timestamp(1362460741000);
        // 2013-03-05T15:02:52Z
        let golden_hour = Timestamp(1362495772000);

        let times = get_times(DATE, LAT, LON, None);

        assert_eq_secs!(solar_noon, times.solar_noon);
        assert_eq_secs!(nadir, times.nadir);
        assert_eq_secs!(sunrise, times.sunrise);
        assert_eq_secs!(sunset, times.sunset);
        assert_eq_secs!(sunrise_end, times.sunrise_end);
        assert_eq_secs!(sunset_start, times.sunset_start);
        assert_eq_secs!(dawn, times.dawn);
        assert_eq_secs!(dusk, times.dusk);
        assert_eq_secs!(nautical_dawn, times.nautical_dawn);
        assert_eq_secs!(nautical_dusk, times.nautical_dusk);
        assert_eq_secs!(night_end, times.night_end);
        assert_eq_secs!(night, times.night);
        assert_eq_secs!(golden_hour_end, times.golden_hour_end);
        assert_eq_secs!(golden_hour, times.golden_hour)
    }

    #[test]
    fn test_times_height() {
        let height = 2000.0;

        // 2013-03-05T10:10:57Z
        let solar_noon = Timestamp(1362478257000);
        // 2013-03-04T22:10:57Z
        let nadir = Timestamp(1362435057000);
        // 2013-03-05T04:25:07Z
        let sunrise = Timestamp(1362457507000);
        // 2013-03-05T15:56:46Z
        let sunset = Timestamp(1362499006000);

        let times = get_times(DATE, LAT, LON, Some(height));

        assert_eq_secs!(solar_noon, times.solar_noon);
        assert_eq_secs!(nadir, times.nadir);
        assert_eq_secs!(sunrise, times.sunrise);
        assert_eq_secs!(sunset, times.sunset);
    }

    #[test]
    fn test_julian() {
        // 1. Jan. 2015
        assert_abs_diff_eq!(2457054.5, to_julian(Timestamp(1422748800000)));
        assert_eq!(Timestamp(1422748800000), from_julian(2457054.5));
    }

    #[test]
    fn test_to_days() {
        // 1. Jan. 2015
        assert_abs_diff_eq!(5509.5, to_days(Timestamp(1422748800000)));
    }

    #[test]
    fn test_moon_pos() {
        let moon_pos = moon_pos(DATE, LAT, LON);

        assert_abs_diff_eq!(moon_pos.azimuth, -0.9783999522438226);
        assert_abs_diff_eq!(moon_pos.altitude, 0.014551482243892251);
        assert_abs_diff_eq!(moon_pos.distance.unwrap(), 364121.37256256194);
    }

    #[test]
    fn test_moon_illumination() {
        let moon_illum = moon_illumination(DATE);

        assert_abs_diff_eq!(moon_illum.fraction, 0.4848068202456373);
        assert_abs_diff_eq!(moon_illum.phase, 0.7548368838538762);
        assert_abs_diff_eq!(moon_illum.angle, 1.6732942678578346);
    }
}