mawaqit 0.2.4

use chrono::{DateTime, Datelike, TimeZone, Utc};

use crate::astronomy::ops;
use crate::astronomy::unit::{Angle, Coordinates, Stride};

/// Geocentric solar coordinates (declination, right ascension, apparent
/// sidereal time) for a given Julian day.
#[derive(PartialEq, Debug, Copy, Clone)]
pub struct SolarCoordinates {
    // The declination of the sun, the angle between
    // the rays of the Sun and the plane of the Earth's equator.
    declination: Angle,

    // Right ascension of the Sun, the angular distance on the
    // celestial equator from the vernal equinox to the hour circle.
    right_ascension: Angle,

    // Apparent sidereal time, the hour angle of the vernal equinox.
    apparent_sidereal_time: Angle,
}

impl SolarCoordinates {
    fn new(julian_day: f64) -> Self {
        let julian_century = ops::julian_century(julian_day);
        let mean_solar_longitude = ops::mean_solar_longitude(julian_century);
        let mean_lunar_longitude = ops::mean_lunar_longitude(julian_century);
        let ascending_lunar_node = ops::ascending_lunar_node_longitude(julian_century);
        let apparent_solar_longitude =
            ops::apparent_solar_longitude(julian_century, mean_solar_longitude).radians();

        let mean_sidereal_time = ops::mean_sidereal_time(julian_century);
        let nutation_longitude = ops::nutation_in_longitude(
            mean_solar_longitude,
            mean_lunar_longitude,
            ascending_lunar_node,
        );
        let nutation_obliq = ops::nutation_in_obliquity(
            mean_solar_longitude,
            mean_lunar_longitude,
            ascending_lunar_node,
        );

        let mean_obliq_ecliptic = ops::mean_obliquity_of_the_ecliptic(julian_century);
        let apparent_obliq_ecliptic =
            ops::apparent_obliquity_of_the_ecliptic(julian_century, mean_obliq_ecliptic).radians();

        // Equation from Astronomical Algorithms page 165
        let declination = Angle::from_radians(
            (apparent_obliq_ecliptic.sin() * apparent_solar_longitude.sin()).asin(),
        );

        // Equation from Astronomical Algorithms page 165
        let right_ascension = Angle::from_radians(
            (apparent_obliq_ecliptic.cos() * apparent_solar_longitude.sin())
                .atan2(apparent_solar_longitude.cos()),
        )
        .unwound();

        // Equation from Astronomical Algorithms page 88
        let apparent_sidereal_time = Angle::new(
            mean_sidereal_time.degrees
                + ((nutation_longitude * 3600.0)
                    * Angle::new(mean_obliq_ecliptic.degrees + nutation_obliq)
                        .radians()
                        .cos())
                    / 3600.0,
        );

        SolarCoordinates {
            declination,
            right_ascension,
            apparent_sidereal_time,
        }
    }
}

/// Solar time data and calculations for a given date and location.
#[derive(Debug, Copy, Clone)]
pub struct SolarTime {
    date: DateTime<Utc>,
    observer: Coordinates,
    solar: SolarCoordinates,
    pub transit: DateTime<Utc>,
    pub sunrise: Option<DateTime<Utc>>,
    pub sunset: Option<DateTime<Utc>>,
    prev_solar: SolarCoordinates,
    next_solar: SolarCoordinates,
    approx_transit: f64,
}

impl SolarTime {
    /// Compute solar time data for the given date and coordinates.
    ///
    /// Returns `Err` when the sun never rises or sets (polar day/night),
    /// or when solar transit is at or below the geometric horizon
    /// (no real daylight).
    pub fn new(date: DateTime<Utc>, coordinates: Coordinates) -> Result<SolarTime, &'static str> {
        // All calculation need to occur at 0h0m UTC
        let today = Utc
            .with_ymd_and_hms(date.year(), date.month(), date.day(), 0, 0, 0)
            .single()
            .ok_or("invalid date for solar time")?;
        let tomorrow = today.tomorrow();
        let yesterday = today.yesterday();
        let prev_solar = SolarCoordinates::new(yesterday.julian_day());
        let solar = SolarCoordinates::new(today.julian_day());
        let next_solar = SolarCoordinates::new(tomorrow.julian_day());
        let solar_altitude = Angle::new(-50.0 / 60.0);
        let approx_transit = ops::approximate_transit(
            coordinates.longitude_angle(),
            solar.apparent_sidereal_time,
            solar.right_ascension,
        );
        let transit_time = ops::corrected_transit(
            approx_transit,
            coordinates.longitude_angle(),
            solar.apparent_sidereal_time,
            solar.right_ascension,
            prev_solar.right_ascension,
            next_solar.right_ascension,
        );

        let sunrise_time = ops::corrected_hour_angle(
            approx_transit,
            solar_altitude,
            coordinates,
            false,
            solar.apparent_sidereal_time,
            ops::InterpolatedAngle::new(
                solar.right_ascension,
                prev_solar.right_ascension,
                next_solar.right_ascension,
            ),
            ops::InterpolatedAngle::new(
                solar.declination,
                prev_solar.declination,
                next_solar.declination,
            ),
        );
        let sunset_time = ops::corrected_hour_angle(
            approx_transit,
            solar_altitude,
            coordinates,
            true,
            solar.apparent_sidereal_time,
            ops::InterpolatedAngle::new(
                solar.right_ascension,
                prev_solar.right_ascension,
                next_solar.right_ascension,
            ),
            ops::InterpolatedAngle::new(
                solar.declination,
                prev_solar.declination,
                next_solar.declination,
            ),
        );

        Ok(SolarTime {
            date,
            observer: coordinates,
            solar,
            transit: SolarTime::setting_hour(transit_time, &date)
                .ok_or("transit computation failed")?,
            sunrise: Some(
                SolarTime::setting_hour(sunrise_time, &date).ok_or("sunrise computation failed")?,
            ),
            sunset: Some(
                SolarTime::setting_hour(sunset_time, &date).ok_or("sunset computation failed")?,
            ),
            prev_solar,
            next_solar,
            approx_transit,
        })
    }

    /// Compute the time when the Sun reaches the given angle below the horizon.
    /// Set `after_transit` to `true` for times after solar transit (e.g. Isha, sunset),
    /// `false` for times before transit (e.g. Fajr, sunrise).
    pub fn time_for_solar_angle(&self, angle: Angle, after_transit: bool) -> Option<DateTime<Utc>> {
        let hours = ops::corrected_hour_angle(
            self.approx_transit,
            angle,
            self.observer,
            after_transit,
            self.solar.apparent_sidereal_time,
            ops::InterpolatedAngle::new(
                self.solar.right_ascension,
                self.prev_solar.right_ascension,
                self.next_solar.right_ascension,
            ),
            ops::InterpolatedAngle::new(
                self.solar.declination,
                self.prev_solar.declination,
                self.next_solar.declination,
            ),
        );

        SolarTime::setting_hour(hours, &self.date).filter(|result| {
            if after_transit {
                *result > self.transit
            } else {
                *result < self.transit
            }
        })
    }

    /// Return the Sun's declination at this date.
    #[must_use]
    pub fn declination(&self) -> Angle {
        self.solar.declination
    }

    /// Calculates the afternoon time when the shadow is `shadow_length` object-heights longer
    /// than its transit minimum. The two commonly-used values for `shadow_length` are `1.0`
    /// and `2.0`. Returns `None` if the required shadow is astronomically unattainable on this
    /// date.
    #[must_use]
    pub fn time_for_shadow(&self, shadow_length: f64) -> Option<DateTime<Utc>> {
        let absolute_degrees = (self.observer.latitude - self.declination().degrees).abs();
        let tangent = Angle::new(absolute_degrees);
        let inverse = shadow_length + tangent.radians().tan();
        let angle = Angle::from_radians((1.0 / inverse).atan());

        if angle.degrees < 0.0 {
            return None;
        }

        self.time_for_solar_angle(angle, true)
    }

    fn setting_hour(value: f64, date: &DateTime<Utc>) -> Option<DateTime<Utc>> {
        if value.is_normal() {
            let calculated_hours = value.floor();
            let calculated_minutes = ((value - calculated_hours) * 60.0).floor();
            let calculated_seconds =
                ((value - (calculated_hours + calculated_minutes / 60.0)) * 60.0 * 60.0).floor();

            let (adjusted_hour, adjusted_date) = SolarTime::hour_adjustment(calculated_hours, date);

            // Round to the nearest minute
            let adjusted_mins = (calculated_minutes + calculated_seconds / 60.0).round() as u32;
            let adjusted_secs: u32 = 0;

            let (hour, mins, secs, result_date) = if adjusted_mins > 59 {
                if adjusted_hour >= 23 {
                    // Minute-rounding can push an already-wrapped hour past
                    // 23:59 into invalid 24:00. Carry to next day.
                    (0, 0, adjusted_secs, adjusted_date.tomorrow())
                } else {
                    (adjusted_hour + 1, 0, adjusted_secs, adjusted_date)
                }
            } else {
                (adjusted_hour, adjusted_mins, adjusted_secs, adjusted_date)
            };

            Utc.with_ymd_and_hms(
                result_date.year(),
                result_date.month(),
                result_date.day(),
                hour,
                mins,
                secs,
            )
            .single()
        } else {
            None
        }
    }

    fn hour_adjustment(calculated_hours: f64, date: &DateTime<Utc>) -> (u32, DateTime<Utc>) {
        // Adjust the hour to be within 0..=23,
        // wrapping around as needed; otherwise
        // chrono method will panic.
        if calculated_hours < 0.0 {
            ((calculated_hours + 24.0) as u32, date.yesterday())
        } else if calculated_hours >= 24.0 {
            ((calculated_hours - 24.0) as u32, date.tomorrow())
        } else {
            (calculated_hours as u32, *date)
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use chrono::Local;

    #[test]
    fn solar_coordinates() {
        let julian_day = ops::julian_day(1992, 10, 13, 0.0);
        let solar = SolarCoordinates::new(julian_day);

        assert!((solar.declination.degrees - (-7.7850685152648795)).abs() < 1e-12);
        assert!((solar.right_ascension.degrees - 198.380_822_142_518_8).abs() < 1e-12);
        assert!((solar.right_ascension.unwound().degrees - 198.380_822_142_518_8).abs() < 1e-12);
    }

    #[test]
    fn zero_out_time_for_a_date() {
        // Local date below is 2019-01-11T04:41:19Z in UTC
        let utc_date = Utc
            .with_ymd_and_hms(2019, 1, 11, 23, 41, 19)
            .single()
            .expect("Invalid date and time provided");
        let updated = Utc
            .with_ymd_and_hms(utc_date.year(), utc_date.month(), utc_date.day(), 0, 0, 0)
            .single();

        assert_eq!(updated, Utc.with_ymd_and_hms(2019, 1, 11, 0, 0, 0).single());
    }

    #[test]
    fn calculate_date_for_tomorrow() {
        let date = Local
            .with_ymd_and_hms(2019, 1, 10, 0, 0, 0)
            .single()
            .expect("Invalid date and time provided");
        let tomorrow = date.tomorrow();

        assert_eq!(tomorrow.year(), 2019);
        assert_eq!(tomorrow.month(), 1);
        assert_eq!(tomorrow.day(), 11);
    }

    #[test]
    fn calculate_julian_date() {
        let utc = Utc
            .with_ymd_and_hms(1992, 10, 13, 0, 0, 0)
            .single()
            .expect("Invalid date and time provided");
        let julian_day = ops::julian_day(1992, 10, 13, 0.0);

        assert_eq!(utc.julian_day(), julian_day);
    }

    #[test]
    fn calculate_solar_time() {
        let coordinates = Coordinates::new(35.0 + 47.0 / 60.0, -78.0 - 39.0 / 60.0);
        let date = Utc
            .with_ymd_and_hms(2015, 7, 12, 0, 0, 0)
            .single()
            .expect("Invalid date and time provided");
        let solar = SolarTime::new(date, coordinates).unwrap();
        let transit_date = Utc.with_ymd_and_hms(2015, 7, 12, 17, 20, 0).unwrap();
        let sunrise_date = Utc.with_ymd_and_hms(2015, 7, 12, 10, 8, 0).unwrap();
        let sunset_date = Utc.with_ymd_and_hms(2015, 7, 13, 0, 32, 0).unwrap();

        assert_eq!(solar.transit, transit_date);
        assert_eq!(solar.sunrise, Some(sunrise_date));
        assert_eq!(solar.sunset, Some(sunset_date));
    }

    #[test]
    fn calculate_time_for_solar_angle() {
        let coordinates = Coordinates::new(35.0 + 47.0 / 60.0, -78.0 - 39.0 / 60.0);
        let date = Utc
            .with_ymd_and_hms(2015, 7, 12, 0, 0, 0)
            .single()
            .expect("Invalid date and time provided");
        let solar = SolarTime::new(date, coordinates).unwrap();
        let angle = Angle::new(-6.0);
        let twilight_start = solar.time_for_solar_angle(angle, false).unwrap();
        let twilight_end = solar.time_for_solar_angle(angle, true).unwrap();

        assert_eq!(twilight_start.format("%-k:%M").to_string(), "9:38");
        assert_eq!(twilight_end.format("%-k:%M").to_string(), "1:02");
    }

    #[test]
    fn calculate_corrected_hour_angle() {
        let coordinates = Coordinates::new(35.0 + 47.0 / 60.0, -78.0 - 39.0 / 60.0);
        let date = Utc
            .with_ymd_and_hms(2015, 7, 12, 0, 0, 0)
            .single()
            .expect("Invalid date and time provided");
        let today = Utc
            .with_ymd_and_hms(date.year(), date.month(), date.day(), 0, 0, 0)
            .single()
            .expect("Invalid date and time provided.");
        let tomorrow = today.tomorrow();
        let yesterday = today.yesterday();
        let prev_solar = SolarCoordinates::new(yesterday.julian_day());
        let solar = SolarCoordinates::new(today.julian_day());
        let next_solar = SolarCoordinates::new(tomorrow.julian_day());
        let solar_altitude = Angle::new(-50.0 / 60.0);
        let approx_transit = ops::approximate_transit(
            coordinates.longitude_angle(),
            solar.apparent_sidereal_time,
            solar.right_ascension,
        );
        let sunrise_time = ops::corrected_hour_angle(
            approx_transit,
            solar_altitude,
            coordinates,
            false,
            solar.apparent_sidereal_time,
            ops::InterpolatedAngle::new(
                solar.right_ascension,
                prev_solar.right_ascension,
                next_solar.right_ascension,
            ),
            ops::InterpolatedAngle::new(
                solar.declination,
                prev_solar.declination,
                next_solar.declination,
            ),
        );

        assert!((sunrise_time - 10.131_800_480_632_85).abs() < 1e-12);
    }
}