use crate::temporal::epoch::J2000_EPOCH;
use sciforge::hub::prelude::constants::AU;
pub const EARTH_AXIAL_TILT_DEG: f64 = 23.4393;
pub struct SolarPosition {
pub azimuth_deg: f64,
pub elevation_deg: f64,
pub distance_au: f64,
pub direction: [f64; 3],
}
impl SolarPosition {
pub fn compute(julian_date: f64, lat_deg: f64, lon_deg: f64) -> Self {
let d = julian_date - J2000_EPOCH;
let mean_anomaly = (357.5291 + 0.98560028 * d) % 360.0;
let m_rad = mean_anomaly.to_radians();
let equation_of_center =
1.9148 * m_rad.sin() + 0.0200 * (2.0 * m_rad).sin() + 0.0003 * (3.0 * m_rad).sin();
let l_venus = (3.176_146_7 + 1_021.328_554_4 * d / 36525.0) % (2.0 * std::f64::consts::PI);
let l_jupiter = (0.599_546 + 145.270_285_3 * d / 36525.0) % (2.0 * std::f64::consts::PI);
let perturbation = 0.00313 * (2.0 * l_venus - m_rad).sin()
+ 0.00178 * (3.0 * m_rad).sin()
+ 0.00077 * (l_jupiter - m_rad).sin();
let ecliptic_lon =
(mean_anomaly + equation_of_center + perturbation + 180.0 + 102.9372) % 360.0;
let ecl_rad = ecliptic_lon.to_radians();
let obliquity_rad = EARTH_AXIAL_TILT_DEG.to_radians();
let declination = (obliquity_rad.sin() * ecl_rad.sin()).asin();
let right_ascension = (obliquity_rad.cos() * ecl_rad.sin()).atan2(ecl_rad.cos());
let sidereal_time = (280.16 + 360.9856235 * d + lon_deg) % 360.0;
let hour_angle = sidereal_time.to_radians() - right_ascension;
let lat_rad = lat_deg.to_radians();
let altitude = (lat_rad.sin() * declination.sin()
+ lat_rad.cos() * declination.cos() * hour_angle.cos())
.asin();
let azimuth = (declination.sin() - altitude.sin() * lat_rad.sin())
.atan2(altitude.cos() * lat_rad.cos());
let dir_x = altitude.cos() * azimuth.sin();
let dir_y = altitude.sin();
let dir_z = altitude.cos() * azimuth.cos();
Self {
azimuth_deg: azimuth.to_degrees(),
elevation_deg: altitude.to_degrees(),
distance_au: 1.0 - 0.0167 * m_rad.cos(),
direction: [dir_x, dir_y, dir_z],
}
}
pub fn distance_m(&self) -> f64 {
self.distance_au * AU
}
pub fn is_above_horizon(&self) -> bool {
self.elevation_deg > 0.0
}
}