ut-astria / orbdetpy

Orbit Determination with Python
https://ut-astria.github.io/orbdetpy
GNU General Public License v3.0
115 stars 38 forks source link

Failing to interpolate known GPS data #15

Closed edvardm closed 4 years ago

edvardm commented 4 years ago

Hi,

First off, thanks for a really good package! I hope I could contribute something to the project. However, there's currently an issue I'm having: for only few lines of data from given file I'm able to get nice interpolated data, but this fails with nans as values, based on sample test_ephemeris.py that is included in the distribution, but instead of propagating an orbit I use existing GPS data:

import csv
import sys
from typing import Union

try:
    import pendulum
    from loguru import logger
except ImportError:
    sys.exit("Please run 'pip install -U pendulum loguru' to run this")

from orbdetpy import Frame
from orbdetpy.conversion import get_J2000_epoch_offset
from orbdetpy.utilities import interpolate_ephemeris

STEP_SIZE = 60.0

Numeric = Union[int, float]

def main(fname):
    epoch = pendulum.parse("31 Dec 1999 23:59:28", strict=False).timestamp()
    logger.debug("(Unix) epoch set to {}, using file {}", epoch, fname)
    svs = list(parse_file(epoch, fname))  # state vectors; time, pos<3d>, vel<3d>

    times = [j2k_epoch_offset(sv[0]) for sv in svs]
    start_t, end_t = times[0], times[-1]
    assert start_t < end_t

    states = [e[1:] for e in svs]  # extract only position and velocity data

    logger.debug(
        "Interpolating {:d} values using step size {:.3f}s between {} .. {}",
        len(times),
        STEP_SIZE,
        start_t,
        end_t,
    )

    for row in interpolate_ephemeris(
        source_frame=Frame.GCRF,
        times=times,
        states=states,
        num_points=len(times),
        dest_frame=Frame.GCRF,
        interp_start=start_t,
        interp_end=end_t,
        step_size=STEP_SIZE,
    ):
        print(row.time, row.true_state)

def j2k_epoch_offset(stk_epoch: Numeric):
    time_str = pendulum.from_timestamp(stk_epoch).isoformat()
    return get_J2000_epoch_offset(time_str)

def parse_file(epoch, fpath):
    with open(fpath) as fh:
        reader = csv.reader(fh, delimiter=";")
        for row in reader:
            offset, *rest = [float(e) for e in row]
            yield [offset + epoch, *rest]  # we want first entry as time, which is epoch + offset from that

if __name__ == "__main__":
    import sys

    main(sys.argv[1])

Data file (ephemerides) attached, can be run with

pip install -U loguru pendulum  # needed only if missing those packages
python interpolate.py raw.txt

So the output I get is

orbdetpy-server version 2.0.2
2020-09-02 12:43:03.585 | DEBUG    | __main__:main:22 - (Unix) epoch set to 946684768.0, using file raw.csv
2020-09-02 12:43:03.752 | DEBUG    | __main__:main:40 - Interpolating 48 values using step size 60.000s between 631108873.184 .. 631109286.184
...
631109233.184 [1.5458659847874403e+47, -8.419914678062475e+46, 3.36869557803364e+47, 1.2458308352515215e+47, 8.523956895271134e+46, 1.7033524692088471e+47]
631109286.184 [nan, nan, nan, nan, nan, nan]
Shutting down server...

I'm pretty certain that the coordinate system is wrong (likely source is ECEF(?)), but shoudn't interpolation still work? I'm assuming it's simply using Hermite or Lagrange algorithm to interpolate points, which means that it just treats those points as plain numbers?

I'm assuming I'm doing something very silly here, so any help is appreciated.

raw.txt

edvardm commented 4 years ago

As my colleague pointed out, only now realized that interpolated values are way too large. So I guess there's still kind of a bug -- instead of nans I think I'd rather have an exception telling user is a doofus, but noticed that this is wrong tool to use in this situation. What I want is actually orbit propagation, not interpolation.