use crate::AXIAL_TILT_RAD;
use crate::ECCENTRICITY;
use crate::SEMI_MAJOR_AXIS_M;
use sciforge::hub::domain::meteorology::radiation::solar_zenith_angle;

pub struct SolarPosition {
    pub elevation_rad: f64,
    pub azimuth_rad: f64,
    pub distance_au: f64,
}

const AU: f64 = 1.496e11;

impl SolarPosition {
    pub fn compute(
        lat_deg: f64,
        lon_deg: f64,
        true_anomaly_rad: f64,
        subsolar_lon_deg: f64,
    ) -> Self {
        let r = SEMI_MAJOR_AXIS_M * (1.0 - ECCENTRICITY * ECCENTRICITY)
            / (1.0 + ECCENTRICITY * true_anomaly_rad.cos());
        let distance_au = r / AU;

        let solar_declination = AXIAL_TILT_RAD.sin() * true_anomaly_rad.sin();
        let dec_rad = solar_declination.asin();

        let hour_angle_deg = lon_deg - subsolar_lon_deg;
        let hour_angle_rad = hour_angle_deg.to_radians();

        let lat_rad = lat_deg.to_radians();
        let sin_elev =
            lat_rad.sin() * dec_rad.sin() + lat_rad.cos() * dec_rad.cos() * hour_angle_rad.cos();
        let elevation_rad = sin_elev.asin();

        let cos_az = (dec_rad.sin() - lat_rad.sin() * sin_elev)
            / (lat_rad.cos() * elevation_rad.cos() + 1e-30);
        let mut azimuth_rad = cos_az.clamp(-1.0, 1.0).acos();
        if hour_angle_rad.sin() > 0.0 {
            azimuth_rad = 2.0 * std::f64::consts::PI - azimuth_rad;
        }

        Self {
            elevation_rad,
            azimuth_rad,
            distance_au,
        }
    }

    pub fn elevation_deg(&self) -> f64 {
        self.elevation_rad.to_degrees()
    }

    pub fn is_above_horizon(&self) -> bool {
        self.elevation_rad > 0.0
    }

    pub fn irradiance(&self) -> f64 {
        let base = sciforge::hub::domain::meteorology::radiation::solar_constant();
        base / (self.distance_au * self.distance_au)
    }
}

pub fn zenith_angle(lat_deg: f64, dec_deg: f64, hour_angle_deg: f64) -> f64 {
    solar_zenith_angle(lat_deg, dec_deg, hour_angle_deg)
}
