Merge branch 'develop'

README.md

@@ -37,7 +37,7 @@ Installation
 Development
 -----------
 
-1. Download and extract <https://github.com/ut-astria/orbdetpy/releases/download/2.0.6/orekit-data.tar.gz>
+1. Download and extract <https://github.com/ut-astria/orbdetpy/releases/download/2.0.7/orekit-data.tar.gz>
    under the `orbdetpy/` sub-folder.
 
 2. Developers will need Apache Maven 3+ to build the Java library. Build
@@ -51,6 +51,8 @@ Development
 
    `mvn -e package`
 
+3. Run `pip install -e ./` from the top-level folder containing `setup.py`.
+
 Examples
 --------
 

examples/fit_radec.py

@@ -69,8 +69,8 @@
 initial_state = iod_laplace(Frame.EME2000, lat, lon, alt, times, ra, dec)
 
 # Configure orbit determination
-config = configure(prop_start=initial_epoch, prop_initial_state=initial_state,
-                   prop_end=get_J2000_epoch_offset(real_obs[-1][0]))
+config = configure(prop_start=initial_epoch, prop_initial_state=initial_state, prop_end=get_J2000_epoch_offset(real_obs[-1][0]),
+                   rso_mass=250.0, rso_area=10.0, estm_DMC_corr_time=15.0, estm_DMC_sigma_pert=3.0)
 
 # Define ground station(s)
 add_station(config, "UTA-ASTRIA-NMSkies", lat, lon, alt)
@@ -81,6 +81,8 @@
 
 # Use Filter.EXTENDED_KALMAN for the EKF
 config.estm_filter = Filter.UNSCENTED_KALMAN
+# Set the initial covariance; diagonal entries for the estimated state vector and parameters 
+config.estm_covariance[:] = [25E6, 25E6, 25E6, 1E2, 1E2, 1E2, 1.00, 0.25, 1E-6, 1E-6, 1E-6]
 
 # Build a list of Measurement objects, one for each RA/dec pair
 meas_obj = []

examples/predict_passes.py

@@ -1,5 +1,5 @@
 # predict_passes.py - Predict satellite passes for ground stations.
-# Copyright (C) 2020 University of Texas
+# Copyright (C) 2020-2021 University of Texas
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -18,7 +18,7 @@
 from os import path
 import datetime as dt
 import tkinter as tk
-from orbdetpy import add_station, configure, Constant, Frame, MeasurementType
+from orbdetpy import add_station, configure, Constant, DragModel, Frame, MeasurementType, OutputFlag
 from orbdetpy.conversion import get_J2000_epoch_offset, get_UTC_string, pos_to_lla
 from orbdetpy.propagation import propagate_orbits
 
@@ -113,8 +113,7 @@ def load_settings(self):
 
     def save_settings(self):
         cfg = {}
-        for f in ["latitude", "longitude", "altitude", "fov_azimuth", "fov_elevation",
-                  "fov_aperture", "tle"]:
+        for f in ["latitude", "longitude", "altitude", "fov_azimuth", "fov_elevation", "fov_aperture", "tle"]:
             cfg[f] = getattr(self, f).get("0.0", "end-1c") if (f == "tle") else getattr(self, f).get()
 
         self.master.destroy()
@@ -125,53 +124,48 @@ def predict(self):
         self.predict["state"] = "disabled"
         start = get_J2000_epoch_offset(self.start_time.get())
         end = get_J2000_epoch_offset(self.end_time.get())
-        data, tle = {}, [l for l in self.tle.get("0.0", "end-1c").splitlines()
-                         if l.startswith("1") or l.startswith("2")]
+        data, tle = {}, [l for l in self.tle.get("0.0", "end-1c").splitlines() if l.startswith("1") or l.startswith("2")]
 
-        for f in ["latitude", "longitude", "altitude", "fov_azimuth", "fov_elevation",
-                  "fov_aperture", "step_size"]:
+        for f in ["latitude", "longitude", "altitude", "fov_azimuth", "fov_elevation", "fov_aperture", "step_size"]:
             data[f] = [float(t.strip()) for t in getattr(self, f).get().split(",")]
             if (f not in ["altitude", "step_size"]):
                 data[f] = [d*Constant.DEGREE_TO_RAD for d in data[f]]
 
         cfg_list = []
         sim_meas = len(data["latitude"]) <= 2 or len(data["latitude"]) != len(data["longitude"])
         for i in range(0, len(tle), 2):
-            cfg_list.append(configure(prop_start=start, prop_initial_TLE=tle[i:i+2],
-                                      prop_end=end, prop_step=data["step_size"][0],
-                                      sim_measurements=sim_meas))
-            if (not sim_meas):
-                cfg_list[-1].geo_zone_lat_lon[:] = [l for ll in zip(data["latitude"], data["longitude"])
-                                                    for l in ll]
-                continue
-
-            add_station(cfg_list[-1], "Sensor", data["latitude"][0], data["longitude"][0],
-                        data["altitude"][0], data["fov_azimuth"][0], data["fov_elevation"][0],
-                        data["fov_aperture"][0])
-            cfg_list[-1].measurements[MeasurementType.AZIMUTH].error[:] = [0.0]
-            cfg_list[-1].measurements[MeasurementType.ELEVATION].error[:] = [0.0]
+            cfg_list.append(configure(prop_start=start, prop_initial_TLE=tle[i:i+2], prop_end=end, prop_step=data["step_size"][0],
+                                      sim_measurements=sim_meas, gravity_degree=-1, gravity_order=-1, ocean_tides_degree=-1,
+                                      ocean_tides_order=-1, third_body_sun=False, third_body_moon=False, solid_tides_sun=False,
+                                      solid_tides_moon=False, drag_model=DragModel.UNDEFINED, rp_sun=False))
+            cfg_list[-1].output_flags |= OutputFlag.OUTPUT_ECLIPSE
+            if (sim_meas):
+                add_station(cfg_list[-1], "Sensor", data["latitude"][0], data["longitude"][0], data["altitude"][0],
+                            data["fov_azimuth"][0], data["fov_elevation"][0], data["fov_aperture"][0])
+                cfg_list[-1].measurements[MeasurementType.AZIMUTH].error[:] = [0.0]
+                cfg_list[-1].measurements[MeasurementType.ELEVATION].error[:] = [0.0]
+            else:
+                cfg_list[-1].geo_zone_lat_lon[:] = [l for ll in zip(data["latitude"], data["longitude"]) for l in ll]
 
         if (len(cfg_list)):
             i = 0
+            lookup = {0.0: "Sunlit", 0.5: "Penumbra", 1.0: "Umbra"}
             self.output.delete("0.0", tk.END)
             if (sim_meas):
                 self.output_label["text"] = "UTC, Azimuth [deg], Elevation [deg]"
             else:
                 self.output_label["text"] = "UTC, Latitude [deg], Longitude [deg], Altitude [m]"
 
             for o in propagate_orbits(cfg_list):
-                self.output.insert(tk.END, "\nObject {}:\n".format(tle[i][2:7]))
+                self.output.insert(tk.END, f"\nObject {tle[i][2:7]}:\n")
                 i += 2
                 for m in o.array:
                     if (sim_meas):
-                        self.output.insert(tk.END, "{}: {:.5f}, {:.5f}\n".format(
-                            get_UTC_string(m.time), (m.values[0]/Constant.DEGREE_TO_RAD + 360)%360,
-                            m.values[1]/Constant.DEGREE_TO_RAD))
+                        self.output.insert(tk.END, "{}: {:.5f}, {:.5f} ({})\n".format(get_UTC_string(m.time),                                                                                      (m.values[0]/Constant.DEGREE_TO_RAD + 360)%360, m.values[1]/Constant.DEGREE_TO_RAD, lookup[m.true_state[-1]]))
                     else:
                         lla = pos_to_lla(Frame.GCRF, m.time, m.true_state)
-                        self.output.insert(tk.END, "{}: {:.5f}, {:.5f}, {:.2f}\n".format(
-                            get_UTC_string(m.time), lla[0]/Constant.DEGREE_TO_RAD, lla[1]/Constant.DEGREE_TO_RAD, lla[2]))
-
+                        self.output.insert(tk.END, "{}: {:.5f}, {:.5f}, {:.2f} ({})\n".format(get_UTC_string(m.time),
+                            lla[0]/Constant.DEGREE_TO_RAD, lla[1]/Constant.DEGREE_TO_RAD, lla[2], lookup[m.true_state[-1]]))
         self.predict["state"] = "normal"
 
 Application(master=tk.Tk()).mainloop()

examples/propagate_tle.py

@@ -1,5 +1,5 @@
 # propagate_tle.py - Propagate many TLEs in parallel.
-# Copyright (C) 2019-2020 University of Texas
+# Copyright (C) 2019-2021 University of Texas
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -17,7 +17,8 @@
 import sys
 import argparse
 from datetime import datetime, timedelta
-from orbdetpy import configure
+from orbdetpy import configure, DragModel
+from orbdetpy.ccsds import export_OEM
 from orbdetpy.conversion import get_J2000_epoch_offset, get_UTC_string
 from orbdetpy.propagation import propagate_orbits
 
@@ -32,25 +33,29 @@
 parser.add_argument("--end-time", help="propagation end time",
                     default=datetime(day1.year, day1.month, day1.day).strftime(timefmt))
 parser.add_argument("--step-size", help="step size [sec]", type=float, default=900.0)
+parser.add_argument("--export-oem", help="Export output to CCSDS OEM files", action="store_true", default=False)
 if (len(sys.argv) == 1):
     parser.print_help()
     exit(1)
 
 args = parser.parse_args()
 start = get_J2000_epoch_offset(args.start_time)
 end = get_J2000_epoch_offset(args.end_time)
-config, elements = [], [l for l in getattr(args, "tle-file").read().splitlines()
-                        if l.startswith("1") or l.startswith("2")]
+config, elements = [], [l for l in getattr(args, "tle-file").read().splitlines() if l.startswith("1") or l.startswith("2")]
 for i in range(0, len(elements), 2):
-    config.append(configure(prop_start=start, prop_initial_TLE=elements[i:i+2],
-                            prop_end=end, prop_step=args.step_size))
+    config.append(configure(prop_start=start, prop_initial_TLE=elements[i:i+2], prop_end=end, prop_step=args.step_size, gravity_degree=-1,
+                            gravity_order=-1, ocean_tides_degree=-1, ocean_tides_order=-1, third_body_sun=False, third_body_moon=False,
+                            solid_tides_sun=False, solid_tides_moon=False, drag_model=DragModel.UNDEFINED, rp_sun=False))
 
 try:
-    i = 0
-    for o in propagate_orbits(config):
-        print("\nObject {}:".format(elements[i][2:7]))
-        i += 2
-        for m in o.array:
+    for idx, obj in enumerate(propagate_orbits(config)):
+        obj_id = elements[idx*2][2:7]
+        print(f"\nObject {obj_id}:")
+        for m in obj.array:
             print(get_UTC_string(m.time), m.true_state)
+
+        if (args.export_oem):
+            with open(obj_id + ".oem", "w") as fp:
+                fp.write(export_OEM(config[idx], obj.array, obj_id, obj_id))
 except Exception as exc:
     print(exc)

examples/test_conversion.py

@@ -36,14 +36,16 @@
 
 # Frame transformations
 itrf_pv = conv.transform_frame(Frame.GCRF, time, gcrf_pv, Frame.ITRF_CIO_CONV_2010_ACCURATE_EOP)
-print("GCRF->ITRF: {:.2f}m, {:.2f}m, {:.2f}m, {:.2f}m/s, {:.2f}m/s, {:.2f}m/s".format(*itrf_pv))
-print("ITRF->GCRF: {:.2f}m, {:.2f}m, {:.2f}m, {:.2f}m/s, {:.2f}m/s, {:.2f}m/s\n".format(
+print("GCRF=>ITRF: {:.2f}m, {:.2f}m, {:.2f}m, {:.2f}m/s, {:.2f}m/s, {:.2f}m/s".format(*itrf_pv))
+print("ITRF=>GCRF: {:.2f}m, {:.2f}m, {:.2f}m, {:.2f}m/s, {:.2f}m/s, {:.2f}m/s\n".format(
     *conv.transform_frame(Frame.ITRF_CIO_CONV_2010_ACCURATE_EOP, time, itrf_pv, Frame.GCRF)))
 
 # Position => Lat/Lon/Alt
 lla = conv.pos_to_lla(Frame.GCRF, time, gcrf_pv[:3])
-print("pos_to_lla: Lat={:.2f}d, Lon={:.2f}d, Alt={:.2f}m\n".format(
-    lla[0]/Constant.DEGREE_TO_RAD, lla[1]/Constant.DEGREE_TO_RAD, lla[2]))
+print(f"pos_to_lla: Lat={lla[0]/Constant.DEGREE_TO_RAD:.2f}d, Lon={lla[1]/Constant.DEGREE_TO_RAD:.2f}d, Alt={lla[2]:.2f}m")
+
+# Lat/Lon/Alt => Position
+print("lla_to_pos: {:.2f}m, {:.2f}m, {:.2f}m\n".format(*conv.lla_to_pos(time, lla)))
 
 # State vector <=> orbital elements
 elem = conv.pv_to_elem(Frame.GCRF, time, gcrf_pv)
@@ -56,7 +58,7 @@
 lat, lon, alt = 52.5*Constant.DEGREE_TO_RAD, -1.917*Constant.DEGREE_TO_RAD, 0.0
 ra, dec = 250.425*Constant.DEGREE_TO_RAD, 36.467*Constant.DEGREE_TO_RAD
 az, el = conv.radec_to_azel(Frame.GCRF, time, ra, dec, lat, lon, alt)
-print("radec_to_azel: Azimuth={:.3f}d, Elevation={:.3f}d".format(az/Constant.DEGREE_TO_RAD, el/Constant.DEGREE_TO_RAD))
+print(f"radec_to_azel: Azimuth={az/Constant.DEGREE_TO_RAD:.3f}d, Elevation={el/Constant.DEGREE_TO_RAD:.3f}d")
 
 r, d = conv.azel_to_radec(time, az, el, lat, lon, alt, Frame.GCRF)
-print("azel_to_radec: RA={:.3f}d, Dec={:.3f}d\n".format(r/Constant.DEGREE_TO_RAD, d/Constant.DEGREE_TO_RAD))
+print(f"azel_to_radec: RA={r/Constant.DEGREE_TO_RAD:.3f}d, Dec={d/Constant.DEGREE_TO_RAD:.3f}d\n")

examples/test_estimation.py

@@ -51,13 +51,11 @@
 
 # ManeuverType.*_CHANGE maneuvers take a single "delta" argument representing change in the
 # corresponding element; values are in meters for distances and radians for angles
-#add_maneuver(cfg, get_J2000_epoch_offset("2019-05-01T00:10:00Z"), ManeuverTrigger.DATE_TIME,
-#             None, ManeuverType.PERIGEE_CHANGE, [50000.0])
+#add_maneuver(cfg, get_J2000_epoch_offset("2019-05-01T00:10:00Z"), ManeuverTrigger.DATE_TIME, None, ManeuverType.PERIGEE_CHANGE, [5E4])
 
 # Define ground stations
 add_station(cfg, "Maui", 0.3614, -2.7272, 3059.0)
-add_station(cfg, "Millstone", 0.7438, -1.2652, 100.0, fov_azimuth=-2.75,
-            fov_elevation=0.62, fov_aperture=15*Constant.DEGREE_TO_RAD)
+add_station(cfg, "Millstone", 0.7438, -1.2652, 100.0, fov_azimuth=-2.75, fov_elevation=0.62, fov_aperture=15*Constant.DEGREE_TO_RAD)
 
 # Uncomment mutually exclusive sections below to choose different measurements 
 # AZIMUTH and ELEVATION must be given
@@ -91,16 +89,18 @@
 meas.sort(key=lambda m: m.time)
 
 # Plot orbital elements from the simulation
-#simulation_plot(meas, interactive=True)
+simulation_plot(meas, interactive=True)
 
-# Uncomment to export simulated measurements to a CCSDS TDM file
-#with open("orbdetpy_sim.tdm", "w") as fp:
-#    fp.write(export_TDM(cfg, meas, "COVID-19"))
+# Export simulated measurements to a CCSDS TDM file
+with open("orbdetpy_sim.tdm", "w") as fp:
+    fp.write(export_TDM(cfg, meas, "DEBRIS"))
 
 # Perturb truth initial state before running OD
 cfg.prop_initial_state[:] = multivariate_normal(cfg.prop_initial_state, diag([1E6, 1E6, 1E6, 1E2, 1E2, 1E2]))
 # Use Filter.EXTENDED_KALMAN for the EKF
 cfg.estm_filter = Filter.UNSCENTED_KALMAN
+# Set the initial covariance; diagonal entries for the estimated state vector and parameters 
+cfg.estm_covariance[:] = [1E6, 1E6, 1E6, 1E2, 1E2, 1E2, 1.00, 0.25, 1E-6, 1E-6, 1E-6]
 
 # Run OD on simulated measurements
 fit = determine_orbit([cfg], [meas])[0]
@@ -116,6 +116,6 @@
 # Plot OD results
 estimation_plot(cfg, meas, fit, interactive=True, estim_param=False)
 
-# Uncomment to export ephemeris to a CCSDS OEM file
-#with open("orbdetpy_fit.oem", "w") as fp:
-#    fp.write(export_OEM(cfg, fit, "2020-001A", "COVID-19"))
+# Export ephemeris to a CCSDS OEM file
+with open("orbdetpy_fit.oem", "w") as fp:
+    fp.write(export_OEM(cfg, fit, "2021-001A", "DEBRIS"))

orbdetpy/__init__.py

@@ -21,9 +21,8 @@
 from os import path, stat
 from math import pi, sqrt
 from typing import List, Optional, Tuple
-from .version import __version__
+from orbdetpy.version import __version__
 from orbdetpy.rpc.messages_pb2 import Facet, Maneuver, Measurement, Parameter, Settings, Station
-from orbdetpy.rpc.server import RemoteServer
 
 class AttitudeType():
     """Orekit attitude providers.
@@ -74,6 +73,7 @@ class Filter():
 
     EXTENDED_KALMAN = 0
     UNSCENTED_KALMAN = 1
+    BATCH_LEAST_SQUARES = 2
 
 class Frame():
     """Orekit reference frames.
@@ -139,6 +139,10 @@ class ManeuverTrigger():
     DATE_TIME = 1
     LONGITUDE_CROSSING = 2
     APSIDE_CROSSING = 3
+    LONGITUDE_EXTREMUM = 4
+    LATITUDE_CROSSING = 5
+    LATITUDE_EXTREMUM = 6
+    NODE_CROSSING = 7
 
 class ManeuverType():
     """Maneuver types.
@@ -194,6 +198,7 @@ class OutputFlag():
     OUTPUT_INNO_COV = 4
     OUTPUT_RESIDUALS = 8
     OUTPUT_DENSITY = 16
+    OUTPUT_ECLIPSE = 32
 
 class Constant():
     """Miscellaneous constants.
@@ -363,6 +368,3 @@ def ltr_to_matrix(lower_triangle: List[float])->List[List[float]]:
     _libs_dir = path.join(_root_dir, "target")
     _jar_file = path.join(_libs_dir, "orbdetpy-server-{}.jar".format(__version__))
     _data_dir = path.join(_root_dir, "orekit-data")
-    stat(_jar_file)
-    stat(_data_dir)
-    RemoteServer.connect(_data_dir, _jar_file)

orbdetpy/astro_data.py

@@ -1,5 +1,5 @@
 # astro_data.py - Update Orekit astrodynamics data files.
-# Copyright (C) 2019-2020 University of Texas
+# Copyright (C) 2019-2021 University of Texas
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -14,9 +14,10 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
+import tarfile
 import requests
-from os import path
-from orbdetpy import _data_dir
+from os import path, remove
+from orbdetpy import _root_dir, _data_dir
 
 def format_weather(lines: str)->str:
     """Re-format space weather data into a more efficient form.
@@ -30,14 +31,29 @@ def format_weather(lines: str)->str:
     for line in lines.splitlines():
         if (line == "END DAILY_PREDICTED"):
             break
-        if (len(line) > 0 and line[0].isnumeric()):
-            output.append(",".join([line[i:j] for i, j in zip(c1, c2)]))
+        if (line and line[0].isnumeric()):
+            output.append(",".join((line[i:j].strip() for i, j in zip(c1, c2))))
     return("\n".join(output))
 
 def update_data()->None:
     """Download and re-format astrodynamics data from multiple sources.
     """
 
+    if (not path.isdir(_data_dir)):
+        uri = "https://github.com/ut-astria/orbdetpy/releases/download/2.0.6/orekit-data.tar.gz"
+        print(f"Downloading {uri}")
+        resp = requests.get(uri, timeout=10.0, stream=True)
+        if (resp.status_code == requests.codes.ok):
+            tgz = path.join(_root_dir, "orekit-data.tar.gz")
+            with open(tgz, "wb") as fp:
+                fp.write(resp.raw.read())
+            tar = tarfile.open(tgz, "r:gz")
+            tar.extractall()
+            tar.close()
+            remove(tgz)
+        else:
+            print(f"HTTP error: {resp.status_code}")
+
     updates = [["https://datacenter.iers.org/data/latestVersion/7_FINALS.ALL_IAU1980_V2013_017.txt",
                 path.join(_data_dir, "Earth-Orientation-Parameters", "IAU-1980", "finals.all"), None],
                ["https://datacenter.iers.org/data/latestVersion/9_FINALS.ALL_IAU2000_V2013_019.txt",
@@ -51,8 +67,11 @@ def update_data()->None:
             resp = requests.get(u[0], timeout=10.0)
             if (resp.status_code == requests.codes.ok):
                 with open(u[1], "w") as fp:
-                    fp.write(u[2](resp.text) if (u[2] is not None) else resp.text)
+                    fp.write(u[2](resp.text) if (u[2]) else resp.text)
             else:
                 print(f"HTTP error: {resp.status_code}")
         except Exception as exc:
             print(exc)
+
+if (__name__ == '__main__'):
+    update_data()