predict_rs/
geodetic.rs

1use core::f64::consts::{FRAC_PI_2, TAU};
2
3use crate::consts::{EARTH_RADIUS_KM_WGS84, FLATTENING_FACTOR, JULIAN_TIME_DIFF};
4use crate::julian_date::theta_g_jd;
5use crate::math::{fmod2p, Vec3};
6
7#[derive(Debug)]
8pub struct Geodetic {
9    pub lat: f64,
10    pub lon: f64,
11    pub alt: f64,
12    pub theta: f64,
13}
14
15/// Calculate the Geodetic position of a satellite given a unix epoch time
16/// and it's ECI position
17#[must_use]
18pub fn calculate_lat_lon_alt(time: f64, pos: Vec3) -> Geodetic {
19    /* Procedure Calculate_LatLonAlt will calculate the geodetic  */
20    /* position of an object given its ECI position pos and time. */
21    /* It is intended to be used to determine the ground track of */
22    /* a satellite.  The calculations  assume the earth to be an  */
23    /* oblate spheroid as defined in WGS '72.                     */
24
25    /* Reference:  The 1992 Astronomical Almanac, page K12. */
26
27    //Convert to julian time:
28    let newtime = time + JULIAN_TIME_DIFF;
29    let theta = pos.1.atan2(pos.0); /* radians */
30    let lon = fmod2p(theta - theta_g_jd(newtime)); /* radians */
31    let r = (pos.0.powf(2.0) + pos.1.powf(2.0)).sqrt();
32    let e2 = FLATTENING_FACTOR * (2.0 - FLATTENING_FACTOR);
33    let mut lat = pos.2.atan2(r); /* radians */
34    let mut phi;
35    let mut c;
36    loop {
37        phi = lat;
38        c = 1.0 / (1.0 - e2 * phi.sin().powf(2.0)).sqrt();
39        lat = (pos.2 + EARTH_RADIUS_KM_WGS84 * c * e2 * phi.sin()).atan2(r);
40        // geodetic->lat=atan2(pos[2]+EARTH_RADIUS_KM_WGS84*c*e2*sin(phi),r);
41        if (lat - phi).abs() < 1E-10 {
42            break;
43        }
44    }
45
46    let alt = r / lat.cos() - EARTH_RADIUS_KM_WGS84 * c; /* kilometers */
47
48    if lat > FRAC_PI_2 {
49        lat -= TAU;
50    }
51    Geodetic {
52        lat,
53        lon,
54        alt,
55        theta,
56    }
57}
predict_rs/geodetic.rs

predict_rs/
geodetic.rs
