satellogic / orbit-predictor

Python library to propagate satellite orbits.
MIT License
140 stars 43 forks source link

SmartLocationPredictors sometimes gets wrong result #113

Closed astrojuanlu closed 3 years ago

astrojuanlu commented 3 years ago
from datetime import datetime, timedelta

from orbit_predictor.locations import Location
from orbit_predictor.predictors.numerical import J2Predictor
from orbit_predictor.predictors.pass_iterators import (
    LocationPredictor,
    SmartLocationPredictor,
)

start = datetime(2020, 5, 9)
end = start + timedelta(days=1)
location = Location(name="loc", latitude_deg=11, longitude_deg=0, elevation_m=0)
predictor = J2Predictor(
    sma=475 + 6371, ecc=1.65e-3, inc=53.0, argp=90, raan=0, ta=300, epoch=start,
)

orbit_predictor_passes_new = list(
    predictor.passes_over(
        location,
        start,
        aos_at_dg=0,
        max_elevation_gt=0,
        limit_date=end,
        location_predictor_class=SmartLocationPredictor,
        tolerance_s=1e-3,
    )
)

orbit_predictor_passes_old = list(
    predictor.passes_over(
        location,
        start,
        aos_at_dg=0,
        max_elevation_gt=0,
        limit_date=end,
        location_predictor_class=LocationPredictor,
        tolerance_s=1e-3,
    )
)
print(len(orbit_predictor_passes_old))
op_pass = orbit_predictor_passes_old[1]
print(
    f"{op_pass.aos}, {op_pass.los}, "
    f"{op_pass.max_elevation_date}, "
    f"{op_pass.off_nadir_deg},  "
    f"{op_pass.max_elevation_deg}"
)

print("=============")

print(len(orbit_predictor_passes_new))
op_pass = orbit_predictor_passes_new[1]
print(
    f"{op_pass.aos}, {op_pass.los}, "
    f"{op_pass.max_elevation_date}, "
    f"{op_pass.off_nadir_deg},  "
    f"{op_pass.max_elevation_deg}"
)

Result:

$ python passes_op.py 
5
2020-05-09 09:13:29.377683, 2020-05-09 09:25:01.108565, 2020-05-09 09:19:15.942740, -9.67716027862082,  79.64934166047799
=============
5
2020-05-09 09:13:29.378581, 2020-05-09 09:25:01.109094, 2020-05-09 09:18:30.157464, -34.19785213567671,  52.85510803065424

(via @eguaio)

astrojuanlu commented 3 years ago

SmartLocationPredictor clearly gets a wrong result:

elev

``` import numpy as np import matplotlib.pyplot as plt import pandas as pd problematic_pass_ok = orbit_predictor_passes_old[1] at = problematic_pass_ok.aos dt = timedelta(seconds=1) elevations = {} while at < problematic_pass_ok.los: elevations[at] = location.elevation_for(predictor.get_only_position(at)) at += dt elevs_ser = np.degrees(pd.Series(elevations)) ax = elevs_ser.plot(label="Elevation") ax.plot( problematic_pass_ok.max_elevation_date, np.degrees( location.elevation_for( predictor.get_only_position(problematic_pass_ok.max_elevation_date) ) ), color="g", marker="x", mew=2, label="Old LocationPredictor", ) ax.plot( orbit_predictor_passes_new[1].max_elevation_date, np.degrees( location.elevation_for( predictor.get_only_position( orbit_predictor_passes_new[1].max_elevation_date ) ) ), color="r", marker="x", mew=2, label="SmartLocationPredictor", ) ax.legend() plt.show() ```
astrojuanlu commented 3 years ago

The problem is that this assumption is not necessarily true:

https://github.com/satellogic/orbit-predictor/blob/3d25be6b1de1ecf6f37586e9fd9531e01a47216e/orbit_predictor/predictors/pass_iterators.py#L235-L241

in other words: the next local maximum might not be between the current moment (last LOS) and the next orbital period.

In this particular case, if one reruns minimize_scalar between the computed AOS and LOS, the correct TCA is retrieved:

start_date + dt.timedelta(seconds=minimize_scalar(
    lambda t: -elevation(t),
    bounds=(t_aos, t_los),
    method=minimize_scalar_bounded_alt,
    options=dict(xatol=self.tolerance_s),
).x)

Pros:

Cons: