TLE tracking functions

[varenius]

To track GNSS satellites we need TLE tracking functions. Old stuff is in https://github.com/varenius/salsa/blob/main/Control_program/satellites.py and dependencies.

[varenius]

Links to TLE data in https://github.com/varenius/salsa/blob/main/Control_program/TLElinks.txt

[varenius]

The old code relies on pyephem. I didn't find any TLE->Horizontal code for rust, yet. Maybe we need to re-implement it. A general description with some code-example links (including pyephem, but also PHP etc.) is ar https://en.wikipedia.org/wiki/Simplified_perturbations_models

[varenius]

This PHP function claims to obtain az/el given TLE data, observer location and time: https://github.com/shupp/Predict/blob/8cdba40109f40c90217965b5b9463a0c5129636c/Predict.php#L501. Maybe worth porting to rust and see if it agrees with the current (pyephem-based) software?

[varenius]

Hm, this is complicated. But maybe this will solve most (all?) of our problems: https://docs.rs/sgp4/latest/sgp4/ ?

[varenius]

Seems the sgp4 crate can nicely give x,y,z position (see https://github.com/neuromorphicsystems/sgp4#propagation-variables). Question is how to convert these to something we can handle, for example equatorial R.A./Dec.?

[varenius]

Actually, we probably can use https://github.com/MNahad/star-trak which in turn uses the sgp4 crate.

[varenius]

I don't see the star-track software released as an "official" crate. Could still use it, or possibly the parts (e.g. https://github.com/MNahad/star-trak/blob/master/src/engine/transforms.rs) which we need to get AER (Azimuth, Elevation, Range) coordinates from the sgp4-crate results?

[varenius]

If we assume the SGP4 crate works and can be covinced to give Earth Centered Interial carteesian x,y,z coordinates for a satellite at a given time, then we should be able to use existing gmst calculation + the description at https://celestrak.org/columns/v02n02/ to get az/el.

[varenius]

This code can be used with sgp4 to obtain satellite positions for GALILEO (and others) GNSS as az,el (needs ureq, sgp4):

fn main() -> anyhow::Result<()> {
    let response = ureq::get("https://celestrak.com/NORAD/elements/gp.php")
        .query("GROUP", "galileo")
        .query("FORMAT", "json")
        .call()?;
    let elements_vec: Vec<sgp4::Elements> = response.into_json()?;
    for elements in &elements_vec {
        let satn = elements.object_name.as_ref().unwrap();
        let constants = sgp4::Constants::from_elements(elements)?;
        let timestamp = Utc::now().naive_utc();
        let elapsed_ms = timestamp
            .signed_duration_since(elements.datetime)
            .num_milliseconds();
        if let Ok(prediction) = constants.propagate((elapsed_ms as f64) / (1000.0 * 60.0)) {
            let eci = prediction.position;
            let lon: f64 = 0.20802143022;
            let lat: f64 = 1.00170457462;
            let hor = horizontal_from_eci(eci[0], eci[1], eci[2], lat, lon, 0.0);
            println!("{} {:?}", satn, hor);
        }
    }
    Ok(())
}

To convert from ECI to horizontal we can use a function from https://celestrak.org/columns/v02n02/ as:

// Function to calculate satellite position for observe at given time
// from https://celestrak.org/columns/v02n02/
// satellite ECI x,y,z in km
// observer lat,long in radians, alt in km
fn horizontal_from_eci(xs: f64, ys: f64, zs: f64, lat: f64, lon: f64, alt: f64) -> (f64, f64) {
   let re = 6378.135; // Earth radius in km
   let full_circle = 2.0*PI;
   let theta = (gmst_now() + lon ) % full_circle;
   let r = (re + alt)*lat.cos();
   let xo = r*theta.cos();
   let yo = r*theta.sin();
   let zo = (re + alt)*lat.sin();
   let rx = xs - xo;
   let ry = ys - yo;
   let rz = zs - zo;
   let top_s = lat.sin()* theta.cos()*rx + lat.sin() * theta.sin()*ry - lat.cos()*rz;
   let top_e = -theta.sin()*rx + theta.cos()*ry;
   let top_z = lat.cos()* theta.cos()*rx + lat.cos() * theta.sin()*ry + lat.sin()*rz;

   let az = (-top_e/top_s).atan();
   if top_s > 0.0 {
       let az = az + PI;
   }
   if az < 0.0 {
       let az = az + full_circle;
   }
   let rg = (rx*rx + ry*ry + rz*rz).sqrt();
   let el = (top_z/rg).asin();

   (az.to_degrees(),el.to_degrees())
}

Note that this refers to gmst_now() which is a local test function I have, will do a proper commit where the above has been incorporated in to coords.rs with gmst(when) etc. later. The above has been verified to work against current SALSA software and e.g. https://www.n2yo.com/?s=37846.

[varenius]

eci-->hor function added now in https://github.com/centowen/salsa/pull/45/commits/64fb6b720a1af1ec836e3cd67fef045e6053eca0.