use crate::temporal::epoch::J2000EPOCH;
use sciforge::hub::prelude::constants::AU;
pub const JUPITERAXIALTILTDEG: f64 = 3.13;
pub struct SolarPosition {
pub azimuthdeg: f64,
pub elevationdeg: f64,
pub distanceau: f64,
pub direction: [f64; 3],
}
impl SolarPosition {
pub fn compute(juliandate: f64, latdeg: f64, londeg: f64) -> Self {
let d = juliandate - J2000EPOCH;
let jupiterorbitalperioddays = 4332.59;
let meananomaly = (20.020 + 360.0 * d / jupiterorbitalperioddays) % 360.0;
let mrad = meananomaly.to_radians();
let equationofcenter =
5.5549 * mrad.sin() + 0.1726 * (2.0 * mrad).sin() + 0.004 * (3.0 * mrad).sin();
let eclipticlon = (meananomaly + equationofcenter + 180.0 + 14.331) % 360.0;
let eclrad = eclipticlon.to_radians();
let obliquityrad = JUPITERAXIALTILTDEG.to_radians();
let declination = (obliquityrad.sin() * eclrad.sin()).asin();
let rightascension = (obliquityrad.cos() * eclrad.sin()).atan2(eclrad.cos());
let jupitersiderealrotperiodh = 9.925;
let siderealtime =
(280.16 + 360.0 / (jupitersiderealrotperiodh / 24.0) * d + londeg) % 360.0;
let hourangle = siderealtime.to_radians() - rightascension;
let latrad = latdeg.to_radians();
let altitude = (latrad.sin() * declination.sin()
+ latrad.cos() * declination.cos() * hourangle.cos())
.asin();
let azimuth = (declination.sin() - altitude.sin() * latrad.sin())
.atan2(altitude.cos() * latrad.cos());
let dirx = altitude.cos() * azimuth.sin();
let diry = altitude.sin();
let dirz = altitude.cos() * azimuth.cos();
let distanceau = 5.2044 - 0.2543 * mrad.cos();
Self {
azimuthdeg: azimuth.to_degrees(),
elevationdeg: altitude.to_degrees(),
distanceau,
direction: [dirx, diry, dirz],
}
}
pub fn distancem(&self) -> f64 {
self.distanceau * AU
}
pub fn isabovehorizon(&self) -> bool {
self.elevationdeg > 0.0
}
pub fn solarfluxwm2(&self) -> f64 {
1361.0 / (self.distanceau * self.distanceau)
}
}