use core::f64::consts::{FRAC_PI_2, TAU};
use crate::consts::{EARTH_RADIUS_KM_WGS84, FLATTENING_FACTOR, JULIAN_TIME_DIFF};
use crate::julian_date::theta_g_jd;
use crate::math::{fmod2p, Vec3};
#[derive(Debug)]
pub struct Geodetic {
pub lat: f64,
pub lon: f64,
pub alt: f64,
pub theta: f64,
}
#[must_use]
pub fn calculate_lat_lon_alt(time: f64, pos: Vec3) -> Geodetic {
let newtime = time + JULIAN_TIME_DIFF;
let theta = pos.1.atan2(pos.0);
let lon = fmod2p(theta - theta_g_jd(newtime));
let r = (pos.0.powf(2.0) + pos.1.powf(2.0)).sqrt();
let e2 = FLATTENING_FACTOR * (2.0 - FLATTENING_FACTOR);
let mut lat = pos.2.atan2(r);
let mut phi;
let mut c;
loop {
phi = lat;
c = 1.0 / (1.0 - e2 * phi.sin().powf(2.0)).sqrt();
lat = (pos.2 + EARTH_RADIUS_KM_WGS84 * c * e2 * phi.sin()).atan2(r);
if (lat - phi).abs() < 1E-10 {
break;
}
}
let alt = r / lat.cos() - EARTH_RADIUS_KM_WGS84 * c;
if lat > FRAC_PI_2 {
lat -= TAU;
}
Geodetic {
lat,
lon,
alt,
theta,
}
}