Skip to content

use subsecond precision when propagating #34

Description

@matiasg

I noticed that there are a few kms differences for propagate_to() with respect to other implementations of SGP4

let lines = "1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992
             2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767";
let tle = TwoLineElement::from_lines(lines).unwrap();
let t = NaiveDate::from_ymd_opt(2020, 5, 29)
    .unwrap()
    .and_hms_micro_opt(1, 2, 3, 0)
    .unwrap()
    .and_utc();
let position = tle.propagate_to(t).unwrap().position;
println!("{:?}", position);

yields

[4262.419202151204, 406.5812008219196, -5284.072095504458]

while for instance by using Brandon Rhodes' python implementation

from sgp4.api import Satrec, jday
from sgp4.model import WGS84

TLE_LINES = ("1 25544U 98067A   20148.21301450  .00001715  00000-0  38778-4 0  9992",
             "2 25544  51.6435  92.2789 0002570 358.0648 144.9972 15.49396855228767")

iss = Satrec.twoline2rv(*TLE_LINES, WGS84)
jd = jday(2020, 5, 29, 1, 2, 3)
e, r, v = iss.sgp4(jd[0], jd[1])
r

yields

(4262.298413708936, 403.13838162020915, -5284.434613388618)

which is 3.465 km appart from the previous one.
This comes from not using subsecond precision on julian_day_to_datetime fn.

If changed to this

pub(crate) fn julian_day_to_datetime(jd: c_double) -> DateTime<Utc> {
    let mut year = c_int::default();
    let mut month = c_int::default();
    let mut day = c_int::default();
    let mut hour = c_int::default();
    let mut minute = c_int::default();
    let mut second = c_double::default();

    unsafe {
        invjday(
            jd,
            &mut year,
            &mut month,
            &mut day,
            &mut hour,
            &mut minute,
            &mut second,
        );
    }

    let usec = (second.fract() * 1e6).floor();
    NaiveDate::from_ymd_opt(year, month as u32, day as u32)
        .unwrap()
        .and_hms_micro_opt(hour as u32, minute as u32, second as u32, usec as u32)
        .unwrap()
        .and_utc()
}

it now yields

[4262.298360113213, 403.13686091734655, -5284.434773208045]

This is much closer to the python implementation. And, in fact, if changing the jd definition by 0.2 msecs

jd = jday(2020, 5, 29, 1, 2, 2.9998)

it gives the exact same result.

This last offset of 0.2 msecs seems to come from the use of millisecond precision in propagate_to:

        let min_since_epoch = (t - tle_epoch).num_milliseconds() as f64 / 60_000.;

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions