use std::fmt;

use crate::periodic_terms::*;

const SUN_RADIUS: f64 = 0.26667;

const L_COUNT: usize = 6;
const B_COUNT: usize = 2;
const R_COUNT: usize = 5;
const Y_COUNT: usize = 63;

trait AngleExt {
    fn limit_degrees(self) -> Self;
    fn limit_degrees180(self) -> Self;
    fn limit_degrees180pm(self) -> Self;
    fn limit_minutes(self) -> Self;
    fn limit_zero2one(self) -> Self;
}

impl AngleExt for f64 {
    fn limit_degrees(self) -> f64 {
        self.rem_euclid(360.0)
    }

    fn limit_degrees180(self) -> f64 {
        self.rem_euclid(180.0)
    }

    fn limit_degrees180pm(self) -> f64 {
        let limited = self.rem_euclid(360.0);
        if limited > 180.0 { limited - 360.0 } else { limited }
    }

    fn limit_minutes(self) -> f64 {
        if self < -20.0 {
            self + 1440.0
        } else if self > 20.0 {
            self - 1440.0
        } else {
            self
        }
    }

    fn limit_zero2one(self) -> f64 {
        self.rem_euclid(1.0)
    }
}

fn julian_day(
    year: i64,
    month: i64,
    day: i64,
    hour: i64,
    minute: i64,
    second: f64,
    dut1: f64,
    tz: f64,
) -> f64 {
    let mut year = year as f64;
    let mut month = month as f64;
    let day = day as f64;
    let hour = hour as f64;
    let minute = minute as f64;

    let day_decimal = day + (hour - tz + (minute + (second + dut1) / 60.0) / 60.0) / 24.0;

    if month < 3.0 {
        month += 12.0;
        year -= 1.0;
    }

    let mut jd = (365.25 * (year + 4716.0)).floor()
        + (30.6001 * (month + 1.0)).floor()
        + day_decimal
        - 1524.5;

    if jd > 2299160.0 {
        let a = (year / 100.0).floor();
        jd += 2.0 - a + (a / 4.0).floor();
    }

    jd
}

fn julian_century(jd: f64) -> f64 {
    (jd - 2451545.0) / 36525.0
}

fn julian_ephemeris_day(jd: f64, delta_t: f64) -> f64 {
    jd + delta_t / 86400.0
}

fn julian_ephemeris_century(jde: f64) -> f64 {
    (jde - 2451545.0) / 36525.0
}

fn julian_ephemeris_millennium(jce: f64) -> f64 {
    jce / 10.0
}

fn third_order_polynomial(a: f64, b: f64, c: f64, d: f64, x: f64) -> f64 {
    ((a * x + b) * x + c) * x + d
}

fn earth_periodic_term_summation(terms: &[[f64; 3]], jme: f64) -> f64 {
    terms.iter().map(|&[a, b, c]| a * (b + c * jme).cos()).sum()
}

fn earth_values(term_sum: &[f64], jme: f64) -> f64 {
    term_sum
        .iter()
        .enumerate()
        .map(|(i, &s)| s * jme.powi(i as i32))
        .sum::<f64>()
        / 1.0e8
}

fn earth_heliocentric_longitude(jme: f64) -> f64 {
    let mut sum = [0.0; L_COUNT];
    for i in 0..L_COUNT {
        sum[i] = earth_periodic_term_summation(L_TERMS[i], jme);
    }
    earth_values(&sum, jme).to_degrees().limit_degrees()
}

fn earth_heliocentric_latitude(jme: f64) -> f64 {
    let mut sum = [0.0; B_COUNT];
    for i in 0..B_COUNT {
        sum[i] = earth_periodic_term_summation(B_TERMS[i], jme);
    }
    earth_values(&sum, jme).to_degrees()
}

fn earth_radius_vector(jme: f64) -> f64 {
    let mut sum = [0.0; R_COUNT];
    for i in 0..R_COUNT {
        sum[i] = earth_periodic_term_summation(R_TERMS[i], jme);
    }
    earth_values(&sum, jme)
}

fn geocentric_longitude(l: f64) -> f64 {
    (l + 180.0).limit_degrees()
}

fn geocentric_latitude(b: f64) -> f64 {
    -b
}

fn mean_elongation_moon_sun(jce: f64) -> f64 {
    third_order_polynomial(1.0 / 189474.0, -0.0019142, 445267.11148, 297.85036, jce)
}

fn mean_anomaly_sun(jce: f64) -> f64 {
    third_order_polynomial(-1.0 / 300000.0, -0.0001603, 35999.05034, 357.52772, jce)
}

fn mean_anomaly_moon(jce: f64) -> f64 {
    third_order_polynomial(1.0 / 56250.0, 0.0086972, 477198.867398, 134.96298, jce)
}

fn argument_latitude_moon(jce: f64) -> f64 {
    third_order_polynomial(1.0 / 327270.0, -0.0036825, 483202.017538, 93.27191, jce)
}

fn ascending_longitude_moon(jce: f64) -> f64 {
    third_order_polynomial(1.0 / 450000.0, 0.0020708, -1934.136261, 125.04452, jce)
}

fn xy_term_summation(i: usize, x: &[f64]) -> f64 {
    x.iter()
        .zip(Y_TERMS[i].iter())
        .map(|(xi, &yi)| xi * yi as f64)
        .sum()
}

fn nutation_longitude_and_obliquity(jce: f64, x: &[f64]) -> (f64, f64) {
    let mut sum_psi = 0.0;
    let mut sum_epsilon = 0.0;

    for i in 0..Y_COUNT {
        let xy = xy_term_summation(i, x).to_radians();
        let t = &PE_TERMS[i];
        sum_psi += (t[0] + jce * t[1]) * xy.sin();
        sum_epsilon += (t[2] + jce * t[3]) * xy.cos();
    }

    (sum_psi / 36_000_000.0, sum_epsilon / 36_000_000.0)
}

fn ecliptic_mean_obliquity(jme: f64) -> f64 {
    let u = jme / 10.0;
    84381.448
        + u * (-4680.93
            + u * (-1.55
                + u * (1999.25
                    + u * (-51.38
                        + u * (-249.67
                            + u * (-39.05
                                + u * (7.12 + u * (27.87 + u * (5.79 + u * 2.45)))))))))
}

fn ecliptic_true_obliquity(delta_epsilon: f64, epsilon0: f64) -> f64 {
    delta_epsilon + epsilon0 / 3600.0
}

fn aberration_correction(r: f64) -> f64 {
    -20.4898 / (3600.0 * r)
}

fn apparent_sun_longitude(theta: f64, delta_psi: f64, delta_tau: f64) -> f64 {
    theta + delta_psi + delta_tau
}

fn greenwich_mean_sidereal_time(jd: f64, jc: f64) -> f64 {
    (280.46061837
        + 360.98564736629 * (jd - 2451545.0)
        + jc * jc * (0.000387933 - jc / 38710000.0))
        .limit_degrees()
}

fn greenwich_sidereal_time(nu0: f64, delta_psi: f64, epsilon: f64) -> f64 {
    nu0 + delta_psi * epsilon.to_radians().cos()
}

fn geocentric_right_ascension(lamda: f64, epsilon: f64, beta: f64) -> f64 {
    let sin_lambda = lamda.to_radians().sin();
    let cos_lambda = lamda.to_radians().cos();
    let cos_epsilon = epsilon.to_radians().cos();
    let sin_epsilon = epsilon.to_radians().sin();
    let tan_beta = beta.to_radians().tan();

    (sin_lambda * cos_epsilon - tan_beta * sin_epsilon)
        .atan2(cos_lambda)
        .to_degrees()
        .limit_degrees()
}

fn geocentric_declination(beta: f64, epsilon: f64, lamda: f64) -> f64 {
    let sin_beta = beta.to_radians().sin();
    let cos_beta = beta.to_radians().cos();
    let cos_epsilon = epsilon.to_radians().cos();
    let sin_epsilon = epsilon.to_radians().sin();
    let sin_lambda = lamda.to_radians().sin();

    (sin_beta * cos_epsilon + cos_beta * sin_epsilon * sin_lambda)
        .asin()
        .to_degrees()
}

fn observer_hour_angle(nu: f64, longitude: f64, alpha_deg: f64) -> f64 {
    (nu + longitude - alpha_deg).limit_degrees()
}

fn sun_equatorial_horizontal_parallax(r: f64) -> f64 {
    8.794 / (3600.0 * r)
}

fn right_ascension_parallax_and_topocentric_dec(
    latitude: f64,
    elevation: f64,
    xi: f64,
    h: f64,
    delta: f64,
) -> (f64, f64) {
    let lat_rad = latitude.to_radians();
    let xi_rad = xi.to_radians();
    let h_rad = h.to_radians();
    let delta_rad = delta.to_radians();

    let u = (0.99664719 * lat_rad.tan()).atan();
    let y = 0.99664719 * u.sin() + elevation * lat_rad.sin() / 6378140.0;
    let x = u.cos() + elevation * lat_rad.cos() / 6378140.0;

    let delta_alpha_rad =
        (-x * xi_rad.sin() * h_rad.sin()).atan2(delta_rad.cos() - x * xi_rad.sin() * h_rad.cos());

    let delta_prime_rad = (delta_rad.sin() - y * xi_rad.sin() * delta_alpha_rad.cos())
        .atan2(delta_rad.cos() - x * xi_rad.sin() * h_rad.cos());

    (delta_alpha_rad.to_degrees(), delta_prime_rad.to_degrees())
}

fn topocentric_local_hour_angle(h: f64, delta_alpha: f64) -> f64 {
    h - delta_alpha
}

fn topocentric_elevation_angle(latitude: f64, delta_prime: f64, h_prime: f64) -> f64 {
    let lat_rad = latitude.to_radians();
    let delta_prime_rad = delta_prime.to_radians();
    let h_prime_rad = h_prime.to_radians();
    (lat_rad.sin() * delta_prime_rad.sin()
        + lat_rad.cos() * delta_prime_rad.cos() * h_prime_rad.cos())
    .asin()
    .to_degrees()
}

fn atmospheric_refraction_correction(
    pressure: f64,
    temperature: f64,
    atmos_refract: f64,
    e0: f64,
) -> f64 {
    if e0 < -(SUN_RADIUS + atmos_refract) {
        return 0.0;
    }
    let a = (pressure / 1010.0) * (283.0 / (273.0 + temperature));
    let e = (e0 + 10.3 / (e0 + 5.11)).to_radians();
    a * 1.02 / (60.0 * e.tan())
}

fn topocentric_elevation_angle_corrected(e0: f64, delta_e: f64) -> f64 {
    e0 + delta_e
}

fn topocentric_zenith_angle(e: f64) -> f64 {
    90.0 - e
}

fn topocentric_azimuth_angle_astro(h_prime: f64, latitude: f64, delta_prime: f64) -> f64 {
    let h_prime_rad = h_prime.to_radians();
    let lat_rad = latitude.to_radians();
    let delta_prime_rad = delta_prime.to_radians();
    let a = h_prime_rad.cos() * lat_rad.sin() - delta_prime_rad.tan() * lat_rad.cos();
    h_prime_rad.sin().atan2(a).to_degrees().limit_degrees()
}

fn surface_incidence_angle(zenith: f64, azimuth_astro: f64, azm_rotation: f64, slope: f64) -> f64 {
    let zenith_rad = zenith.to_radians();
    let slope_rad = slope.to_radians();
    let azimuth_astro_rad = azimuth_astro.to_radians();
    let azm_rotation_rad = azm_rotation.to_radians();

    (zenith_rad.cos() * slope_rad.cos()
        + slope_rad.sin() * zenith_rad.sin() * (azimuth_astro_rad - azm_rotation_rad).cos())
    .acos()
    .to_degrees()
}

fn sum_mean_longitude(jme: f64) -> f64 {
    (280.4664567
        + jme * (360007.6982779
            + jme * (0.03032028
                + jme * (1.0 / 49931.0
                    + jme * (-1.0 / 15300.0 + jme * (-1.0 / 2000000.0))))))
    .limit_degrees()
}

fn eot(m: f64, alpha: f64, del_psi: f64, epsilon: f64) -> f64 {
    (4.0 * (m - 0.0057183 - alpha + del_psi * epsilon.to_radians().cos())).limit_minutes()
}

fn dayfrac_to_local_hr(dayfrac: f64, timezone: f64) -> f64 {
    24.0 * (dayfrac + timezone / 24.0).limit_zero2one()
}

fn approx_sun_transit_time(alpha_zero: f64, longitude: f64, nu: f64) -> f64 {
    (alpha_zero - longitude - nu) / 360.0
}

fn sun_hour_angle_at_rise_set(latitude: f64, delta_zero: f64, h0_prime: f64) -> f64 {
    let argument = (h0_prime.to_radians().sin()
        - latitude.to_radians().sin() * delta_zero.to_radians().sin())
        / (latitude.to_radians().cos() * delta_zero.to_radians().cos());

    if argument.abs() <= 1.0 {
        argument.acos().to_degrees().limit_degrees180()
    } else {
        -99999.0
    }
}

fn approx_sun_rise_and_set(sun_transit: f64, h0: f64) -> (f64, f64, f64) {
    let h0_dfrac = h0 / 360.0;
    (
        (sun_transit - h0_dfrac).limit_zero2one(),
        (sun_transit + h0_dfrac).limit_zero2one(),
        sun_transit.limit_zero2one(),
    )
}

fn rts_alpha_delta_prime(jd_zero: f64, jd_minus: f64, jd_plus: f64, n: f64) -> f64 {
    let mut a = jd_zero - jd_minus;
    let mut b = jd_plus - jd_zero;

    if a.abs() >= 2.0 {
        a = a.limit_zero2one();
    }
    if b.abs() >= 2.0 {
        b = b.limit_zero2one();
    }
    jd_zero + n * (a + b + (b - a) * n) / 2.0
}

fn rts_sun_altitude(latitude: f64, delta_prime: f64, h_prime: f64) -> f64 {
    let a = latitude.to_radians().sin() * delta_prime.to_radians().sin();
    let b = latitude.to_radians().cos() * delta_prime.to_radians().cos() * h_prime.to_radians().cos();
    (a + b).asin().to_degrees()
}

fn sun_rise_and_set(
    sun_transit: f64,
    h_rts: f64,
    delta_prime: f64,
    latitude: f64,
    h_prime: f64,
    h0_prime: f64,
) -> f64 {
    sun_transit + (h_rts - h0_prime)
        / (360.0
            * delta_prime.to_radians().cos()
            * latitude.to_radians().cos()
            * h_prime.to_radians().sin())
}

fn topocentric_azimuth_angle(azimuth_astro: f64) -> f64 {
    (azimuth_astro + 180.0).limit_degrees()
}

/// Estimate atmospheric pressure (mbar) from elevation (metres).
///
/// Uses the standard atmosphere formula. Useful when measured pressure is
/// unavailable — pass the result as [`SpaInput::pressure`].
#[must_use]
pub fn pressure_from_elevation(elevation_m: f64) -> f64 {
    1013.25 * (1.0 - 2.25577e-5 * elevation_m).powf(5.25588)
}


/// Validation error returned by [`SpaInput::calculate`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum SpaError {
    /// Year outside −2000 ..= 6000.
    Year,
    /// Month outside 1 ..= 12.
    Month,
    /// Day outside 1 ..= 31.
    Day,
    /// Hour outside 0 ..= 24.
    Hour,
    /// Minute outside 0 ..= 59 (or > 0 when hour = 24).
    Minute,
    /// Second outside 0.0 ..< 60.0 (or > 0 when hour = 24).
    Second,
    /// |ΔT| > 8000 seconds.
    DeltaT,
    /// |timezone| > 18 hours.
    Timezone,
    /// |longitude| > 180°.
    Longitude,
    /// |latitude| > 90°.
    Latitude,
    /// Elevation below −6 500 000 m.
    Elevation,
    /// Pressure outside 0 ..= 5000 mbar.
    Pressure,
    /// Temperature outside −273 ..= 6000 °C.
    Temperature,
    /// |slope| > 360°.
    Slope,
    /// |azm_rotation| > 360°.
    AzmRotation,
    /// |atmos_refract| > 5°.
    AtmosRefract,
    /// |delta_ut1| >= 1 second.
    DeltaUt1,
}

impl fmt::Display for SpaError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            SpaError::Year        => write!(f, "year out of range (-2000..=6000)"),
            SpaError::Month       => write!(f, "month out of range (1..=12)"),
            SpaError::Day         => write!(f, "day out of range (1..=31)"),
            SpaError::Hour        => write!(f, "hour out of range (0..=24)"),
            SpaError::Minute      => write!(f, "minute out of range (0..=59)"),
            SpaError::Second      => write!(f, "second out of range (0.0..<60.0)"),
            SpaError::DeltaT      => write!(f, "delta_t magnitude exceeds 8000 s"),
            SpaError::Timezone    => write!(f, "timezone out of range (-18..=18)"),
            SpaError::Longitude   => write!(f, "longitude out of range (-180..=180)"),
            SpaError::Latitude    => write!(f, "latitude out of range (-90..=90)"),
            SpaError::Elevation   => write!(f, "elevation below -6 500 000 m"),
            SpaError::Pressure    => write!(f, "pressure out of range (0..=5000 mbar)"),
            SpaError::Temperature => write!(f, "temperature out of range (-273..=6000 °C)"),
            SpaError::Slope       => write!(f, "slope magnitude exceeds 360°"),
            SpaError::AzmRotation => write!(f, "azm_rotation magnitude exceeds 360°"),
            SpaError::AtmosRefract => write!(f, "atmos_refract magnitude exceeds 5°"),
            SpaError::DeltaUt1    => write!(f, "delta_ut1 magnitude must be < 1 s"),
        }
    }
}

impl std::error::Error for SpaError {}

/// All inputs to the solar position algorithm.
///
/// Construct with a struct literal and use `..Default::default()` to fill
/// fields you do not need to customise:
///
/// ```rust
/// # use from_rust_till_dawn::SpaInput;
/// let out = SpaInput {
///     year: 2024, month: 6, day: 21,
///     hour: 12, latitude: 48.8566, longitude: 2.3522,
///     ..Default::default()
/// }.calculate().unwrap();
/// ```
#[derive(Debug, Clone)]
pub struct SpaInput {
    /// Calendar year (−2000 ..= 6000).
    pub year: i32,
    /// Calendar month (1 ..= 12).
    pub month: i32,
    /// Calendar day (1 ..= 31).
    pub day: i32,
    /// Hour of day in local time (0 ..= 24).
    pub hour: i32,
    /// Minute of hour (0 ..= 59).
    pub minute: i32,
    /// Second of minute (0.0 ..< 60.0).
    pub second: f64,
    /// Fractional UT1 − UTC correction in seconds (|delta_ut1| < 1).
    pub delta_ut1: f64,
    /// Difference between Terrestrial Time and UT1 (ΔT, seconds).
    pub delta_t: f64,
    /// Observer's UTC offset in hours (−18 ..= 18).
    pub timezone: f64,
    /// Observer longitude in degrees, east positive (−180 ..= 180).
    pub longitude: f64,
    /// Observer latitude in degrees, north positive (−90 ..= 90).
    pub latitude: f64,
    /// Observer elevation above sea level in metres.
    pub elevation: f64,
    /// Atmospheric pressure in mbar. `None` estimates it from `elevation`.
    pub pressure: Option<f64>,
    /// Air temperature in °C (affects refraction correction).
    pub temperature: f64,
    /// Slope of the surface for incidence-angle calculation (degrees).
    pub slope: f64,
    /// Surface azimuth rotation from south, east positive (degrees).
    pub azm_rotation: f64,
    /// Atmospheric refraction at sunrise/sunset (degrees, typically 0.5667).
    pub atmos_refract: f64,
}

impl Default for SpaInput {
    fn default() -> Self {
        Self {
            year: 2000,
            month: 1,
            day: 1,
            hour: 0,
            minute: 0,
            second: 0.0,
            delta_ut1: 0.0,
            delta_t: 0.0,
            timezone: 0.0,
            longitude: 0.0,
            latitude: 0.0,
            elevation: 0.0,
            pressure: None,
            temperature: 15.0,
            slope: 0.0,
            azm_rotation: 0.0,
            atmos_refract: 0.5667,
        }
    }
}

/// Computed outputs from the solar position algorithm.
#[derive(Debug, Clone)]
pub struct SpaOutput {
    /// Topocentric zenith angle (degrees).
    pub zenith: f64,
    /// Topocentric azimuth angle, astronomers' convention N=0° clockwise (degrees).
    pub azimuth_astro: f64,
    /// Topocentric azimuth angle, navigators' convention N=0° clockwise (degrees).
    pub azimuth: f64,
    /// Refraction-corrected topocentric elevation angle (degrees above horizon).
    pub elevation: f64,
    /// Surface incidence angle (degrees). Only computed for `ZenithAzimuthIncidence` / `All`.
    pub incidence: f64,
    /// Equation of time (minutes). Only computed for `ZenithAzimuthRts` / `All`.
    pub eot: f64,
    /// Solar noon in local fractional hours. Only computed for `ZenithAzimuthRts` / `All`.
    pub suntransit: f64,
    /// Sunrise in local fractional hours. Only computed for `ZenithAzimuthRts` / `All`.
    pub sunrise: f64,
    /// Sunset in local fractional hours. Only computed for `ZenithAzimuthRts` / `All`.
    pub sunset: f64,
}

impl SpaInput {
    /// Run the solar position algorithm and return the computed outputs.
    #[must_use = "the computed solar position is unused"]
    pub fn calculate(self) -> Result<SpaOutput, SpaError> {
        validate_inputs(&self)?;
        let pressure = self.pressure.unwrap_or_else(|| pressure_from_elevation(self.elevation));
        Ok(run_calculation(&self, pressure))
    }
}

// Calculation 

fn validate_inputs(i: &SpaInput) -> Result<(), SpaError> {
    if i.year < -2000 || i.year > 6000 { return Err(SpaError::Year); }
    if i.month < 1 || i.month > 12 { return Err(SpaError::Month); }
    if i.day < 1 || i.day > 31 { return Err(SpaError::Day); }
    if i.hour < 0 || i.hour > 24 { return Err(SpaError::Hour); }
    if i.minute < 0 || i.minute > 59 { return Err(SpaError::Minute); }
    if i.second < 0.0 || i.second >= 60.0 { return Err(SpaError::Second); }
    if let Some(p) = i.pressure {
        if p < 0.0 || p > 5000.0 { return Err(SpaError::Pressure); }
    }
    if i.temperature <= -273.0 || i.temperature > 6000.0 { return Err(SpaError::Temperature); }
    if i.delta_ut1 <= -1.0 || i.delta_ut1 >= 1.0 { return Err(SpaError::DeltaUt1); }
    if i.hour == 24 && i.minute > 0 { return Err(SpaError::Minute); }
    if i.hour == 24 && i.second > 0.0 { return Err(SpaError::Second); }
    if i.delta_t.abs() > 8000.0 { return Err(SpaError::DeltaT); }
    if i.timezone.abs() > 18.0 { return Err(SpaError::Timezone); }
    if i.longitude.abs() > 180.0 { return Err(SpaError::Longitude); }
    if i.latitude.abs() > 90.0 { return Err(SpaError::Latitude); }
    if i.atmos_refract.abs() > 5.0 { return Err(SpaError::AtmosRefract); }
    if i.elevation < -6_500_000.0 { return Err(SpaError::Elevation); }
    if i.slope.abs() > 360.0 { return Err(SpaError::Slope); }
    if i.azm_rotation.abs() > 360.0 { return Err(SpaError::AzmRotation); }
    Ok(())
}

struct GeoData {
    jme: f64,
    r: f64,
    nu: f64,
    alpha: f64,
    delta: f64,
    del_psi: f64,
    epsilon: f64,
}

fn geocentric_sun_ra_and_dec(jd: f64, delta_t: f64) -> GeoData {
    let jc = julian_century(jd);
    let jde = julian_ephemeris_day(jd, delta_t);
    let jce = julian_ephemeris_century(jde);
    let jme = julian_ephemeris_millennium(jce);

    let l = earth_heliocentric_longitude(jme);
    let b = earth_heliocentric_latitude(jme);
    let r = earth_radius_vector(jme);

    let theta = geocentric_longitude(l);
    let beta = geocentric_latitude(b);

    let x = [
        mean_elongation_moon_sun(jce),
        mean_anomaly_sun(jce),
        mean_anomaly_moon(jce),
        argument_latitude_moon(jce),
        ascending_longitude_moon(jce),
    ];
    let (del_psi, del_epsilon) = nutation_longitude_and_obliquity(jce, &x);

    let epsilon0 = ecliptic_mean_obliquity(jme);
    let epsilon = ecliptic_true_obliquity(del_epsilon, epsilon0);

    let del_tau = aberration_correction(r);
    let lamda = apparent_sun_longitude(theta, del_psi, del_tau);
    let nu0 = greenwich_mean_sidereal_time(jd, jc);
    let nu = greenwich_sidereal_time(nu0, del_psi, epsilon);

    let alpha = geocentric_right_ascension(lamda, epsilon, beta);
    let delta = geocentric_declination(beta, epsilon, lamda);

    GeoData { jme, r, nu, alpha, delta, del_psi, epsilon }
}

struct RtsOutput {
    eot: f64,
    suntransit: f64,
    sunrise: f64,
    sunset: f64,
}

fn eot_and_rts(
    geo: &GeoData,
    year: i32, month: i32, day: i32,
    delta_t: f64, longitude: f64, latitude: f64,
    timezone: f64, atmos_refract: f64,
) -> RtsOutput {
    const SUN_TRANSIT: usize = 0;
    const SUN_RISE: usize = 1;
    const SUN_SET: usize = 2;
    const JD_MINUS: usize = 0;
    const JD_ZERO: usize = 1;
    const JD_PLUS: usize = 2;

    let m = sum_mean_longitude(geo.jme);
    let eot_val = eot(m, geo.alpha, geo.del_psi, geo.epsilon);

    let jd_midnight = julian_day(year as i64, month as i64, day as i64, 0, 0, 0.0, 0.0, 0.0);
    let nu = geocentric_sun_ra_and_dec(jd_midnight, delta_t).nu;

    let geo_days = [
        geocentric_sun_ra_and_dec(jd_midnight - 1.0, 0.0),
        geocentric_sun_ra_and_dec(jd_midnight,       0.0),
        geocentric_sun_ra_and_dec(jd_midnight + 1.0, 0.0),
    ];
    let alpha = [geo_days[JD_MINUS].alpha, geo_days[JD_ZERO].alpha, geo_days[JD_PLUS].alpha];
    let delta = [geo_days[JD_MINUS].delta, geo_days[JD_ZERO].delta, geo_days[JD_PLUS].delta];

    let h0_prime = -(SUN_RADIUS + atmos_refract);
    let mut m_rts = [0.0f64; 3];
    m_rts[SUN_TRANSIT] = approx_sun_transit_time(alpha[JD_ZERO], longitude, nu);
    let h0 = sun_hour_angle_at_rise_set(latitude, delta[JD_ZERO], h0_prime);

    if h0 < 0.0 {
        return RtsOutput { eot: eot_val, suntransit: -99999.0, sunrise: -99999.0, sunset: -99999.0 };
    }

    let (sun_rise, sun_set, sun_transit) = approx_sun_rise_and_set(m_rts[SUN_TRANSIT], h0);
    m_rts[SUN_TRANSIT] = sun_transit;
    m_rts[SUN_RISE] = sun_rise;
    m_rts[SUN_SET] = sun_set;

    let mut nu_rts = [0.0f64; 3];
    let mut h_rts = [0.0f64; 3];
    let mut alpha_prime = [0.0f64; 3];
    let mut delta_prime = [0.0f64; 3];
    let mut h_prime = [0.0f64; 3];

    for i in 0..3 {
        nu_rts[i] = nu + 360.985647 * m_rts[i];
        let n = m_rts[i] + delta_t / 86400.0;
        alpha_prime[i] = rts_alpha_delta_prime(alpha[JD_ZERO], alpha[JD_MINUS], alpha[JD_PLUS], n);
        delta_prime[i] = rts_alpha_delta_prime(delta[JD_ZERO], delta[JD_MINUS], delta[JD_PLUS], n);
        h_prime[i] = (nu_rts[i] + longitude - alpha_prime[i]).limit_degrees180pm();
        h_rts[i] = rts_sun_altitude(latitude, delta_prime[i], h_prime[i]);
    }

    RtsOutput {
        eot: eot_val,
        suntransit: dayfrac_to_local_hr(
            m_rts[SUN_TRANSIT] - h_prime[SUN_TRANSIT] / 360.0,
            timezone,
        ),
        sunrise: dayfrac_to_local_hr(
            sun_rise_and_set(
                m_rts[SUN_RISE], h_rts[SUN_RISE], delta_prime[SUN_RISE],
                latitude, h_prime[SUN_RISE], h0_prime,
            ),
            timezone,
        ),
        sunset: dayfrac_to_local_hr(
            sun_rise_and_set(
                m_rts[SUN_SET], h_rts[SUN_SET], delta_prime[SUN_SET],
                latitude, h_prime[SUN_SET], h0_prime,
            ),
            timezone,
        ),
    }
}

fn run_calculation(input: &SpaInput, pressure: f64) -> SpaOutput {
    let jd = julian_day(
        input.year as i64, input.month as i64, input.day as i64,
        input.hour as i64, input.minute as i64, input.second,
        input.delta_ut1, input.timezone,
    );

    let geo = geocentric_sun_ra_and_dec(jd, input.delta_t);

    let h = observer_hour_angle(geo.nu, input.longitude, geo.alpha);
    let xi = sun_equatorial_horizontal_parallax(geo.r);
    let (del_alpha, delta_prime) = right_ascension_parallax_and_topocentric_dec(
        input.latitude, input.elevation, xi, h, geo.delta,
    );
    let h_prime = topocentric_local_hour_angle(h, del_alpha);
    let e0 = topocentric_elevation_angle(input.latitude, delta_prime, h_prime);
    let del_e = atmospheric_refraction_correction(pressure, input.temperature, input.atmos_refract, e0);
    let e = topocentric_elevation_angle_corrected(e0, del_e);
    let zenith = topocentric_zenith_angle(e);
    let azimuth_astro = topocentric_azimuth_angle_astro(h_prime, input.latitude, delta_prime);
    let azimuth = topocentric_azimuth_angle(azimuth_astro);

    let incidence = surface_incidence_angle(zenith, azimuth_astro, input.azm_rotation, input.slope);

    let rts = eot_and_rts(
        &geo,
        input.year, input.month, input.day,
        input.delta_t, input.longitude, input.latitude,
        input.timezone, input.atmos_refract,
    );

    SpaOutput {
        zenith,
        azimuth_astro,
        azimuth,
        elevation: e,
        incidence,
        eot: rts.eot,
        suntransit: rts.suntransit,
        sunrise: rts.sunrise,
        sunset: rts.sunset,
    }
}

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

    // Reference values from spa_tester.c (NREL SPA, August 12 2004).
    // Input: 2003-10-17 12:30:30 UTC-7, Golden CO, pressure=820 mbar, temp=11°C.
    #[test]
    fn spa_reference_2003_10_17() {
        let input = SpaInput {
            year: 2003,
            month: 10,
            day: 17,
            hour: 12,
            minute: 30,
            second: 30.0,
            timezone: -7.0,
            delta_ut1: 0.0,
            delta_t: 67.0,
            longitude: -105.1786,
            latitude: 39.742476,
            elevation: 1830.14,
            pressure: Some(820.0),
            temperature: 11.0,
            slope: 30.0,
            azm_rotation: -10.0,
            atmos_refract: 0.5667,
        };


        let jd = julian_day(2003, 10, 17, 12, 30, 30.0, 0.0, -7.0);
        assert!((jd - 2452930.312847).abs() < 1e-6, "jd={jd}");

        let geo = geocentric_sun_ra_and_dec(jd, 67.0);

        let l = earth_heliocentric_longitude(geo.jme);
        assert!((l - 24.01826).abs() < 1e-4, "l={l}");

        let b = earth_heliocentric_latitude(geo.jme);
        assert!((b - (-1.011219e-4)).abs() < 1e-9, "b={b}");

        assert!((geo.r - 0.996542).abs() < 1e-6, "r={}", geo.r);

        let h = observer_hour_angle(geo.nu, -105.1786, geo.alpha);
        assert!((h - 11.105902).abs() < 1e-5, "h={h}");

        assert!((geo.del_psi - (-3.998404e-3)).abs() < 1e-8, "del_psi={}", geo.del_psi);
        assert!((geo.epsilon - 23.440465).abs() < 1e-5, "epsilon={}", geo.epsilon);


        let out = input.calculate().expect("SPA calculation failed");

        assert!((out.zenith - 50.111622).abs() < 1e-5, "zenith={}", out.zenith);
        assert!((out.azimuth - 194.340241).abs() < 1e-5, "azimuth={}", out.azimuth);
        assert!((out.incidence - 25.187000).abs() < 1e-4, "incidence={}", out.incidence);

        // Sunrise 06:12:43 and sunset 17:20:19 — tolerance: 1 second
        let sunrise_expected = 6.0 + 12.0 / 60.0 + 43.0 / 3600.0;
        let sunset_expected = 17.0 + 20.0 / 60.0 + 19.0 / 3600.0;
        assert!((out.sunrise - sunrise_expected).abs() < 1.0 / 3600.0, "sunrise={}", out.sunrise);
        assert!((out.sunset - sunset_expected).abs() < 1.0 / 3600.0, "sunset={}", out.sunset);
    }
}
