From 88a2c7ea80c29c6ca7f39c7c84923934922d5187 Mon Sep 17 00:00:00 2001 From: Jeff Klenzing Date: Thu, 23 Dec 2021 13:21:42 -0500 Subject: [PATCH 1/4] BUG: use consistent epoch --- pysatMissions/instruments/missions_sgp4.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/pysatMissions/instruments/missions_sgp4.py b/pysatMissions/instruments/missions_sgp4.py index 454f7cbb..cb2973cc 100644 --- a/pysatMissions/instruments/missions_sgp4.py +++ b/pysatMissions/instruments/missions_sgp4.py @@ -61,6 +61,9 @@ def init(self): 'August 21–24.')) logger.info(self.acknowledgements) + if 'epoch' not in self.kwargs['load'].keys(): + self.kwargs['load']['epoch'] = self.files.files.index[0] + return @@ -71,7 +74,7 @@ def init(self): def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None, alt_periapsis=None, alt_apoapsis=None, inclination=None, raan=0., arg_periapsis=0., mean_anomaly=0., - bstar=0., one_orbit=False, num_samples=None, cadence='1S'): + epoch=None, bstar=0., one_orbit=False, num_samples=None, cadence='1S'): """Generate position of satellite in ECI co-ordinates. Note @@ -125,6 +128,10 @@ def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None, orbital plane based on the orbital period. Optional when instantiating via orbital elements. (default=0.) + epoch : dt.datetime or NoneType + The epoch used for calculating orbital propagation from Keplerians. + If None, then use the first date in the file list for consistency across + multiple days. Note that this will be set at init. (default=None) bstar : float Inverse of the ballistic coefficient. Used to model satellite drag. Measured in inverse distance (1 / earth radius). Optional when @@ -176,7 +183,7 @@ def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None, times, index, dates = ps_meth.generate_times(fnames, num_samples, freq=cadence) # Calculate epoch for orbital propagator - epoch = (dates[0] - dt.datetime(1949, 12, 31)).days + epoch_days = (epoch - dt.datetime(1949, 12, 31)).days jd, _ = sapi.jday(dates[0].year, dates[0].month, dates[0].day, 0, 0, 0) if inclination is not None: @@ -185,7 +192,7 @@ def load(fnames, tag=None, inst_id=None, TLE1=None, TLE2=None, alt_apoapsis) satellite = sapi.Satrec() # according to module webpage, wgs72 is common - satellite.sgp4init(sapi.WGS72, 'i', 0, epoch, bstar, 0, 0, + satellite.sgp4init(sapi.WGS72, 'i', 0, epoch_days, bstar, 0, 0, eccentricity, np.radians(arg_periapsis), np.radians(inclination), np.radians(mean_anomaly), mean_motion, np.radians(raan)) From 6655f4257186a86161fe8142079f7cd4fe3bf513 Mon Sep 17 00:00:00 2001 From: Jeff Klenzing Date: Thu, 23 Dec 2021 13:21:59 -0500 Subject: [PATCH 2/4] TST: test for consistency between days --- pysatMissions/tests/test_instruments.py | 33 +++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/pysatMissions/tests/test_instruments.py b/pysatMissions/tests/test_instruments.py index 1a023207..3b5fd3d2 100644 --- a/pysatMissions/tests/test_instruments.py +++ b/pysatMissions/tests/test_instruments.py @@ -90,6 +90,39 @@ def teardown_class(self): # Custom package unit tests can be added here + @pytest.mark.parametrize("kwargs", [{}, {'epoch': dt.datetime(2019, 1, 1)}]) + def test_sgp4_data_continuity(self, kwargs): + """Test that data is continuous for sequential days. + + Parameters + ---------- + kwargs : dict + Optional kwargs to pass through. + """ + + # Define sat with custom Keplerian inputs + sat = pysat.Instrument('missions', 'sgp4', inclination=20, + alt_periapsis=400, **kwargs) + + # Get last 10 points of day 1 + sat.load(2018, 1) + day1 = sat.data[-10:] + + # Get first 10 points of day 2 + sat.load(2018, 2) + day2 = sat.data[:10] + + average_gradient = day1.diff().mean() + std_gradient = day1.diff().std() + gradient_between_days = day2.iloc[0] - day1.iloc[-1] + + # Check that the jump between days is within 3 sigma of average gradient + for key in average_gradient.keys(): + del_g = np.abs(average_gradient[key] - gradient_between_days[key]) + assert del_g < 3. * std_gradient[key] + + return + @pytest.mark.parametrize("inst_dict", [x for x in instruments['download']]) @pytest.mark.parametrize("kwarg,output", [(None, 1), ('10s', 10)]) def test_inst_cadence(self, inst_dict, kwarg, output): From 6b62c1bc83fdf0094d648e386af169e9ff70a53a Mon Sep 17 00:00:00 2001 From: Jeff Klenzing Date: Thu, 23 Dec 2021 14:20:54 -0500 Subject: [PATCH 3/4] TST: add parametrize --- pysatMissions/tests/test_instruments.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/pysatMissions/tests/test_instruments.py b/pysatMissions/tests/test_instruments.py index 3b5fd3d2..d5186a1b 100644 --- a/pysatMissions/tests/test_instruments.py +++ b/pysatMissions/tests/test_instruments.py @@ -90,7 +90,12 @@ def teardown_class(self): # Custom package unit tests can be added here - @pytest.mark.parametrize("kwargs", [{}, {'epoch': dt.datetime(2019, 1, 1)}]) + @pytest.mark.parametrize("kwargs", + [{}, + {'inclination': 20, 'alt_periapsis': 400}, + {'inclination': 80, 'alt_periapsis': 500, + 'alt_apoapsis': 600, + 'epoch': dt.datetime(2019, 1, 1)}]) def test_sgp4_data_continuity(self, kwargs): """Test that data is continuous for sequential days. @@ -101,8 +106,9 @@ def test_sgp4_data_continuity(self, kwargs): """ # Define sat with custom Keplerian inputs - sat = pysat.Instrument('missions', 'sgp4', inclination=20, - alt_periapsis=400, **kwargs) + sat = pysat.Instrument( + inst_module=pysatMissions.instruments.missions_sgp4, + **kwargs) # Get last 10 points of day 1 sat.load(2018, 1) From 7af5d655f8a90df5e6232b25844e76ec57f5db3f Mon Sep 17 00:00:00 2001 From: Jeff Klenzing Date: Thu, 23 Dec 2021 15:06:57 -0500 Subject: [PATCH 4/4] STY: simplify --- pysatMissions/tests/test_instruments.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pysatMissions/tests/test_instruments.py b/pysatMissions/tests/test_instruments.py index d5186a1b..3894f242 100644 --- a/pysatMissions/tests/test_instruments.py +++ b/pysatMissions/tests/test_instruments.py @@ -123,9 +123,8 @@ def test_sgp4_data_continuity(self, kwargs): gradient_between_days = day2.iloc[0] - day1.iloc[-1] # Check that the jump between days is within 3 sigma of average gradient - for key in average_gradient.keys(): - del_g = np.abs(average_gradient[key] - gradient_between_days[key]) - assert del_g < 3. * std_gradient[key] + del_g = np.abs(average_gradient - gradient_between_days) + assert np.all(del_g < (3. * std_gradient)) return