Released v2.0.5

README.md

@@ -22,11 +22,11 @@ filter to estimate unmodeled accelerations.
 Installation
 ------------
 
-1. Install the Java Development Kit 8 (1.8) from <http://openjdk.java.net/install/index.html>.
-   Set the `JAVA_HOME` environment variable to the JDK installation
+1. Install Java SE 11 (11.0.10) from <https://www.oracle.com/javadownload>.
+   Set the `JAVA_HOME` environment variable to the Java installation
    folder. The `java` executable must be added to the system path.
 
-2. Install Python 3.7.0 or higher and run `pip install orbdetpy` to install
+2. Install Python 3.8.0 or higher and run `pip install orbdetpy` to install
    orbdetpy and other package dependencies.
 
 3. Update the astrodynamics data under `orbdetpy/orekit-data` periodically by
@@ -37,7 +37,7 @@ Installation
 Development
 -----------
 
-1. Download and extract <https://github.com/ut-astria/orbdetpy/releases/download/2.0.4/orekit-data.tar.gz>
+1. Download and extract <https://github.com/ut-astria/orbdetpy/releases/download/2.0.5/orekit-data.tar.gz>
    under the `orbdetpy/` sub-folder.
 
 2. Developers will need Apache Maven 3+ to build the Java library. Build

docs/orbdetpy/index.html → docs/index.html

@@ -43,12 +43,12 @@ <h2 id="features">Features</h2>
 <h2 id="installation">Installation</h2>
 <ol>
 <li>
-<p>Install the Java Development Kit 8 (1.8) from <a href="http://openjdk.java.net/install/index.html">http://openjdk.java.net/install/index.html</a>.
-Set the <code>JAVA_HOME</code> environment variable to the JDK installation
+<p>Install Java SE 11 (11.0.10) from <a href="https://www.oracle.com/javadownload">https://www.oracle.com/javadownload</a>.
+Set the <code>JAVA_HOME</code> environment variable to the Java installation
 folder. The <code>java</code> executable must be added to the system path.</p>
 </li>
 <li>
-<p>Install Python 3.7.0 or higher and run <code>pip install <a title="orbdetpy" href="#orbdetpy">orbdetpy</a></code> to install
+<p>Install Python 3.8.0 or higher and run <code>pip install <a title="orbdetpy" href="#orbdetpy">orbdetpy</a></code> to install
 orbdetpy and other package dependencies.</p>
 </li>
 <li>
@@ -60,7 +60,7 @@ <h2 id="installation">Installation</h2>
 <h2 id="development">Development</h2>
 <ol>
 <li>
-<p>Download and extract <a href="https://github.com/ut-astria/orbdetpy/releases/download/2.0.4/orekit-data.tar.gz">https://github.com/ut-astria/orbdetpy/releases/download/2.0.4/orekit-data.tar.gz</a>
+<p>Download and extract <a href="https://github.com/ut-astria/orbdetpy/releases/download/2.0.5/orekit-data.tar.gz">https://github.com/ut-astria/orbdetpy/releases/download/2.0.5/orekit-data.tar.gz</a>
 under the <code>orbdetpy/</code> sub-folder.</p>
 </li>
 <li>

examples/fit_radec.py

@@ -1,5 +1,5 @@
 # fit_radec.py - Run OD on real angles data.
-# 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
@@ -97,5 +97,6 @@
 for f in fit:
     # print(f) to dump pre-fits/post-fits/covariances
     print(get_UTC_string(f.time), f.station, f.estimated_state[:6])
+
 # Plot OD results
-plot(config, meas_obj, fit, interactive=True, estim_param=True)
+plot(config, meas_obj, fit, interactive=True, estim_param=False)

examples/test_conversion.py

@@ -1,5 +1,5 @@
 # test_conversion.py - Test conversion functions.
-# 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
@@ -14,43 +14,49 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
-from orbdetpy import Constant, Frame, PositionAngle
-from orbdetpy.conversion import (azel_to_radec, elem_to_pv, get_J2000_epoch_offset, get_UTC_string,
-                                 pos_to_lla, pv_to_elem, radec_to_azel, transform_frame)
+from datetime import datetime
+import orbdetpy.conversion as conv
+from orbdetpy import Constant, Epoch, Frame, PositionAngle
 
 # TT <=> UTC
-# NOTE: J2000.0 epoch is "2000-01-01T12:00:00 TT" and NOT UTC!
-print("UTC at J2000.0 epoch = {}".format(get_UTC_string(0.0)))
-print("Time offset of UTC 2000-01-01T12:00:00Z from J2000.0 epoch = {}s\n".format(
-    get_J2000_epoch_offset("2000-01-01T12:00:00Z")))
-
-time = get_J2000_epoch_offset("1998-08-10T23:10:00Z")
+now = datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%S.%fZ")
+print(f"J2000_EPOCH = {conv.get_UTC_string(0.0)}")
+print(f"""J2000_EPOCH => 2000-01-01T12:00:00Z = {conv.get_J2000_epoch_offset("2000-01-01T12:00:00Z")}s""")
+print(f"""J2000_EPOCH => {now} (now) = {conv.get_J2000_epoch_offset(now)}s""")
+
+# Time difference between pairs of epochs
+for key, val in vars(Epoch).items():
+    if (key.endswith("_EPOCH") and key != "J2000_EPOCH"):
+        print(f"J2000_EPOCH => {key} = {conv.get_epoch_difference(Epoch.J2000_EPOCH, val)}s")
+print()
+
+time = conv.get_J2000_epoch_offset("1998-08-10T23:10:00Z")
 # Inertial GCRF state vector
 gcrf_pv = [-6045E3, -3490E3, 2500E3, -3.457E3, 6.618E3, 2.533E3]
 
 # Frame transformations
-itrf_pv = transform_frame(Frame.GCRF, time, gcrf_pv, Frame.ITRF_CIO_CONV_2010_ACCURATE_EOP)
+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(
-    *transform_frame(Frame.ITRF_CIO_CONV_2010_ACCURATE_EOP, time, itrf_pv, Frame.GCRF)))
+    *conv.transform_frame(Frame.ITRF_CIO_CONV_2010_ACCURATE_EOP, time, itrf_pv, Frame.GCRF)))
 
 # Position => Lat/Lon/Alt
-lla = pos_to_lla(Frame.GCRF, time, gcrf_pv[:3])
+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]))
 
 # State vector <=> orbital elements
-elem = pv_to_elem(Frame.GCRF, time, gcrf_pv)
+elem = conv.pv_to_elem(Frame.GCRF, time, gcrf_pv)
 print("pv_to_elem: a={:.2f}m, e={:.4f}, i={:.2f}d, O={:.2f}d, w={:.2f}d, M={:.2f}d, theta={:.2f}d, E={:.2f}d".format(
     elem[0], elem[1], *[e/Constant.DEGREE_TO_RAD for e in elem[2:]]))
 print("elem_to_pv: {:.2f}m, {:.2f}m, {:.2f}m, {:.2f}m/s, {:.2f}m/s, {:.2f}m/s\n".format(
-    *elem_to_pv(Frame.GCRF, 0.0, *elem[:-2], PositionAngle.MEAN)))
+    *conv.elem_to_pv(Frame.GCRF, 0.0, *elem[:-2], PositionAngle.MEAN)))
 
 # RA/Dec <=> Az/El
 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 = radec_to_azel(Frame.GCRF, time, ra, dec, lat, lon, alt)
+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))
 
-r, d = azel_to_radec(time, az, el, lat, lon, alt, Frame.GCRF)
+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))

examples/test_estimation.py

@@ -1,5 +1,5 @@
 # test_estimation.py - Run measurement simulation and OD tests.
-# 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
@@ -16,7 +16,7 @@
 
 from numpy import diag
 from numpy.random import multivariate_normal
-from orbdetpy import (configure, add_facet, add_maneuver, add_station, AttitudeType,
+from orbdetpy import (configure, add_facet, add_maneuver, add_station, build_measurement, AttitudeType,
                       Filter, ManeuverTrigger, ManeuverType, MeasurementType, Constant)
 from orbdetpy.ccsds import export_OEM, export_TDM
 from orbdetpy.conversion import get_J2000_epoch_offset, get_UTC_string
@@ -79,10 +79,19 @@
 #cfg.measurements[MeasurementType.POSITION_VELOCITY].error[:] = [100.0, 200.0, 300.0, 3.0, 2.0, 1.0]
 
 # Propagate orbits and generate measurements
-meas = propagate_orbits([cfg])[0].array
+meas = list(m for m in propagate_orbits([cfg])[0].array)
+
+# Zero valued measurements are ignored but they force the filters to output
+# propagated states/covariances at arbitrary intermediate times
+t0 = get_J2000_epoch_offset("2019-05-01T00:26:02.853Z")
+for i in range(5):
+    meas.append(build_measurement(t0 + i*60.0, "Maui", [0.0, 0.0]))
+
+# Measurements must be sorted in ascending chronological order
+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:
@@ -106,7 +115,7 @@
     print(get_UTC_string(f.time), f.station, f.estimated_state[:6])
 
 # Plot OD results
-estimation_plot(cfg, meas, fit, interactive=True, estim_param=True)
+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:

orbdetpy/__init__.py

@@ -1,5 +1,5 @@
 # __init__.py - orbdetpy package initialization.
-# Copyright (C) 2018-2020 University of Texas
+# Copyright (C) 2018-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 @@
 .. include:: ../README.md
 """
 
-from os import path
+from os import path, stat
 from math import pi, sqrt
 from typing import List, Optional, Tuple
 from .version import __version__
@@ -360,7 +360,9 @@ def ltr_to_matrix(lower_triangle: List[float])->List[List[float]]:
     __pdoc__ = {m: False for m in ("Facet", "Maneuver", "Measurement", "Parameter", "Settings",
                                    "Station", "TDMFileFormat", "rpc", "version")}
     _root_dir = path.dirname(path.abspath(path.realpath(__file__)))
-    _data_dir = path.join(_root_dir, "orekit-data")
     _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/ccsds.py

@@ -39,6 +39,7 @@ def export_OEM(cfg: Settings, obs, obj_id: str, obj_name: str, time_list: List[s
     Ephemerides in OEM format.
     """
 
+    utc_list = get_UTC_string([o.time for o in obs])
     oem_header = f"""CCSDS_OEM_VERS = 2.0
 CREATION_DATE = {datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%S")}
 ORIGINATOR = UT-Austin
@@ -49,8 +50,8 @@ def export_OEM(cfg: Settings, obs, obj_id: str, obj_name: str, time_list: List[s
 CENTER_NAME = EARTH
 REF_FRAME = {cfg.prop_inertial_frame}
 TIME_SYSTEM = UTC
-START_TIME = {time_list[0] if (time_list) else get_UTC_string(obs[0].time)[:-1]}
-STOP_TIME = {time_list[-1] if (time_list) else get_UTC_string(obs[-1].time)[:-1]}
+START_TIME = {time_list[0] if (time_list) else utc_list[0][:-1]}
+STOP_TIME = {time_list[-1] if (time_list) else utc_list[-1][:-1]}
 META_STOP
 
 """
@@ -59,11 +60,11 @@ def export_OEM(cfg: Settings, obs, obj_id: str, obj_name: str, time_list: List[s
     is_estm = (hasattr(obs[0], "estimated_state") and hasattr(obs[0], "estimated_covariance")
                and hasattr(obs[0], "propagated_covariance"))
     eph_key = "estimated_state" if (is_estm) else "true_state"
-    for o in obs:
+    for utc, o in zip(utc_list, obs):
+        utc = utc[:-1]
         if (o.time in added):
             continue
         added.add(o.time)
-        utc = get_UTC_string(o.time)[:-1]
         if (time_list is not None and utc not in time_list):
             continue
         X = [x/1000.0 for x in getattr(o, eph_key)[:6]]
@@ -103,6 +104,7 @@ def export_TDM(cfg: Settings, obs, obj_id: str, station_list: List[str]=None)->s
     """
 
     miter = cfg.measurements.keys()
+    utc_list = get_UTC_string([o.time for o in obs])
     if (MeasurementType.RIGHT_ASCENSION in miter and MeasurementType.DECLINATION in miter):
         obstype = f"ANGLE_TYPE = RADEC\nREFERENCE_FRAME = {cfg.prop_inertial_frame}"
         obspath = "1,2"
@@ -137,18 +139,20 @@ def export_TDM(cfg: Settings, obs, obj_id: str, station_list: List[str]=None)->s
 META_STOP
 """)
         blocks.append("DATA_START")
-        for o in obs:
+
+        for utc, o in zip(utc_list, obs):
+            utc = utc[:-1]
             if (o.station != sname):
                 continue
-            utc = get_UTC_string(o.time)[:-1]
             if (MeasurementType.RANGE in miter or MeasurementType.RANGE_RATE in miter):
                 if (MeasurementType.RANGE in miter):
                     blocks.append(f"RANGE = {utc} {o.values[0]/1000.0}")
                 if (MeasurementType.RANGE_RATE in miter):
                     blocks.append(f"DOPPLER_INSTANTANEOUS = {utc} {o.values[-1]/1000.0}")
             if ((MeasurementType.AZIMUTH in miter and MeasurementType.ELEVATION in miter) or
                 (MeasurementType.RIGHT_ASCENSION in miter and MeasurementType.DECLINATION in miter)):
-                blocks.append(f"""ANGLE_1 = {utc} {o.values[0]/Constant.DEGREE_TO_RAD}\nANGLE_2 = {utc} {o.values[1]/Constant.DEGREE_TO_RAD}""")
+                blocks.append((f"""ANGLE_1 = {utc} {o.values[0]/Constant.DEGREE_TO_RAD}\n"""
+                               f"""ANGLE_2 = {utc} {o.values[1]/Constant.DEGREE_TO_RAD}"""))
         blocks.append(f"DATA_STOP\n")
 
     return(tdm_header + "\n".join(blocks))

orbdetpy/plotting/estimation.py

@@ -1,5 +1,5 @@
 # estimation.py - Plot orbit determination results.
-# Copyright (C) 2018-2020 University of Texas
+# Copyright (C) 2018-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
@@ -43,7 +43,7 @@ def plot(cfg: Settings, measurements, orbit_fit, interactive: bool=False,
     key = list(cfg.measurements.keys())
     key.sort()
     meas_names = {0: "Azimuth", 1: "Elevation", 2: "Range", 3: "Range Rate", 4: "Right Ascension", 5: "Declination"}
-    inp = [m for m in measurements if (len(m.station) > 0 or len(m.values) >= 3)]
+    inp = [m for m in measurements if (any(x != 0.0 for x in m.values) and (len(m.station) > 0 or len(m.values) >= 3))]
     out = [f for f in orbit_fit if (len(f.pre_fit)*len(f.post_fit) > 0)]
     dmcidx, dmcrun = 6, (cfg.estm_DMC_corr_time > 0.0 and cfg.estm_DMC_sigma_pert > 0.0)
 
@@ -70,22 +70,23 @@ def plot(cfg: Settings, measurements, orbit_fit, interactive: bool=False,
             for m in key:
                 parnames.append(f"{sk}-{meas_names[m]}")
 
-    tstamp, prefit, posfit, inocov, params, estmacc, estmcov, colors = [], [], [], [], [], [], [], []
+    diag_pos = [0, 2, 5, 9, 14, 20]
+    has_truth = len(inp) > 0 and hasattr(inp[0], "true_state") and len(inp[0].true_state) >= 6
+    tstamp,prefit,posfit,inocov,params,estmacc,estmcov,colors,state_err,state_cov = [],[],[],[],[],[],[],[],[],[]
     for i, o in zip(inp, out):
         tstamp.append(i.time)
         colors.append(cmap[o.station if len(o.station) else None])
         prefit.append([ix - ox for ix, ox in zip(i.values, o.pre_fit)])
         posfit.append([ix - ox for ix, ox in zip(i.values, o.post_fit)])
-
-        p = []
-        for m in range(len(o.innovation_covariance)):
-            if (m in [0, 2, 5, 9, 14, 20]):
-                p.append(3.0*numpy.sqrt(o.innovation_covariance[m]))
-        if (len(p) > 0):
-            inocov.append(p)
+        if (len(o.innovation_covariance) > 0):
+            inocov.append([3.0*numpy.sqrt(o.innovation_covariance[m]) for m in diag_pos if (m < len(o.innovation_covariance))])
         else:
             inocov.append([0.0]*len(prefit[0]))
 
+        if (has_truth):
+            state_err.append([ix - ox for ix, ox in zip(i.true_state, o.estimated_state)])
+            state_cov.append([3.0*numpy.sqrt(o.estimated_covariance[m]) for m in diag_pos if (m < len(o.estimated_covariance))])
+
         if (estim_param and len(o.estimated_state) > 6):
             if (dmcrun):
                 params.append(o.estimated_state[6:dmcidx] + o.estimated_state[dmcidx+3:])
@@ -148,13 +149,13 @@ def plot(cfg: Settings, measurements, orbit_fit, interactive: bool=False,
 
     plt.figure(1)
     plt.suptitle("Post-Fit Residuals")
-    for i in range(pre.shape[-1]):
+    for i in range(pos.shape[-1]):
         if (MeasurementType.POSITION in key):
             plt.subplot(3, 1, order[i])
         elif (MeasurementType.POSITION_VELOCITY in key):
             plt.subplot(3, 2, order[i])
         else:
-            plt.subplot(pre.shape[-1], 1, i + 1)
+            plt.subplot(pos.shape[-1], 1, i + 1)
         plt.scatter(tim, pos[:,i], color=colors, marker="o", s=7)
         plt.legend(handles=patches, loc="best")
         plt.plot(tim, -cov[:,i], "-r")
@@ -201,6 +202,28 @@ def plot(cfg: Settings, measurements, orbit_fit, interactive: bool=False,
                 outfiles.append(output_file_path + "_estacc.png")
                 plt.savefig(outfiles[-1], format="png")
 
+    if (has_truth):
+        state_err = numpy.array(state_err)
+        state_cov = numpy.array(state_cov)
+        order = [1, 3, 5, 2, 4, 6]
+        units = ["m", "m", "m", "m/s", "m/s", "m/s"]
+        ylabs = [r"$\Delta x$", r"$\Delta y$", r"$\Delta z$", r"$\Delta v_x$", r"$\Delta v_y$", r"$\Delta v_z$"]
+        plt.figure(4)
+        plt.suptitle("Position and Velocity Errors")
+        for i in range(6):
+            plt.subplot(3, 2, order[i])
+            plt.scatter(tim, state_err[:,i], color=colors, marker="o", s=7)
+            plt.plot(tim, -state_cov[:,i], "-r")
+            plt.plot(tim,  state_cov[:,i], "-r")
+            plt.xlabel("Time [hr]")
+            plt.ylabel(f"{ylabs[i]} [{units[i]}]")
+            plt.grid(True)
+
+        plt.tight_layout(rect=[0, 0.03, 1, 0.95])
+        if (output_file_path):
+            outfiles.append(output_file_path + "_posvel_error.png")
+            plt.savefig(outfiles[-1], format="png")
+
     if (interactive):
         plt.show()
     plt.close("all")

orbdetpy/plotting/simulation.py

@@ -1,5 +1,5 @@
 # simulation.py - Plot simulation results.
-# Copyright (C) 2018-2020 University of Texas
+# Copyright (C) 2018-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
@@ -40,6 +40,8 @@ def plot(sim_data, interactive: bool=False, output_file_path: str=None):
             rv = [x/1000.0 for x in o.true_state[:6]]
         else:
             rv = [x/1000.0 for x in o.estimated_state[:6]]
+        if (not rv):
+            continue
         rn = norm(rv[:3])
         vn = norm(rv[3:])
         h = cross(rv[:3], rv[3:])