Released v2.0.1

.gitignore

@@ -1,9 +1,9 @@
-target/
 __pycache__/
-output/
 build/
 dist/
 orbdetpy.egg-info/
+orbdetpy/orekit-data/
+orbdetpy/target/
 *.json
 *.oem
 *.tdm

README.md

@@ -1,14 +1,10 @@
-Introduction
-------------
-
-This is orbdetpy, a library of Python and Java routines for orbit
-determination. It is a thin Python wrapper for our estimation tools
-and Orekit, which are both written in Java. 
+This is orbdetpy, a Python and Java library for orbit determination.
+It is a thin Python wrapper for our Java estimation tools and Orekit.
 
 Features
 --------
 
-The force model for orbit propagation currently includes:
+The force model for orbit propagation can be configured to include:
 
 1. EGM96 gravity field up to degree and order 360.
 2. Earth solid tides due to the influence of the Sun and Moon.
@@ -26,56 +22,59 @@ acclerations that result from unmodeled dynamics, maneuvers etc.
 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 point to your JDK installation. The `java`
-   executable must also be in your system path.
+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
+   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
-   orbdetpy and other package dependencies. If you wish to use the `develop`
-   or other experimental branches from GitHub, `git clone` them and run
-   `pip install -e .` from the top level `orbdetpy` folder.
-
-3. Source code, example programs and data files can be downloaded from 
-   <https://github.com/ut-astria/orbdetpy>.
+   orbdetpy and other package dependencies.
 
-4. Update the astrodynamics data in `orbdetpy/data` periodically by running
-   the following. You might need to do so as root on Unix-like systems.
+3. Update the astrodynamics data under `orbdetpy/orekit-data` periodically by
+   running the following. You might need to be the `root` user on some systems.
 
    `python -c "from orbdetpy.astro_data import update_data; update_data();"`
 
-5. Apache Maven 3+ is needed if you hack the Java code and need to
-   rebuild the JAR files. Switch to the `orbdetpy/` folder and run the
-   following, where `os_cpu_type` is `linux-x86_64`, `linux-x86_32`,
-   `windows-x86_64`, `windows-x86_32`, `osx-x86_64`, or `osx-x86_32`,
-   depending on your CPU architecture and OS.
-
-   `mvn -e -Dos.detected.classifier=os_cpu_type package`
-
-   If you are on Intel/AMD 64-bit Linux the command-line simplifies to:
-
-   `mvn -e package`
-
 Examples
 --------
 
-1. `predict_passes.py` : Predict satellite passes for ground stations or
+1. `fit_radec.py` : Run OD with real angles measurements.
+
+2. `predict_passes.py` : Predict satellite passes for ground stations or
    geographic regions using TLEs. Current elements may be obtained from
    sites such as <http://www.celestrak.com>.
 
-2. `propagate_tle.py` : Propagate TLEs given by command-line arguments.
+3. `propagate_tle.py` : Propagate TLEs given by command-line arguments.
 
-3. `test_conversion.py` : Test reference frame and other conversion functions.
+4. `test_conversion.py` : Test reference frame and other conversion functions.
 
-4. `test_estimation.py` : Demonstrates measurement simulation and orbit
+5. `test_estimation.py` : Demonstrates measurement simulation and orbit
    determination functions.
 
-5. `test_interpolation.py` : Interpolate state vectors.
+6. `test_interpolation.py` : Interpolate state vectors.
+
+Development
+-----------
+
+1. Developers will need Apache Maven 3+ to build the Java library. Build
+   using the following from the `orbdetpy/` sub-folder, where `os_cpu_type` is
+   `linux-x86_64`, `linux-x86_32`, `windows-x86_64`, `windows-x86_32`,
+   `osx-x86_64`, or `osx-x86_32` depending on your CPU and OS:
+
+   `mvn -e -Dos.detected.classifier=os_cpu_type package`
+
+   The command-line is simpler on Intel/AMD 64-bit Linux:
+
+   `mvn -e package`
+
+2. Download and extract <https://github.com/ut-astria/orbdetpy/releases/download/2.0.1/orekit-data.tar.gz>
+   under the `orbdetpy/` sub-folder.
 
 Known Issues
 ------------
 
-1. In Microsoft Windows, you might receive warnings from the "Windows
-   Defender Firewall" when you import `orbdetpy`. You can safely allow
-   `orbdetpy` access to the network because this only involves a single
-   TCP port on `localhost`.
+1. You might receive warnings from the Windows Defender Firewall on Microsoft
+   Windows. Grant `orbdetpy` network access permissions.
+
+2. If you use the `multiprocessing` Python package, imports and calls into
+   `orbdetpy` must not span `multiprocessing` function calls. That is, `orbdetpy`
+   can be used in the parent process or the spawned child processes, but not both.

examples/fit_radec.py

@@ -0,0 +1,90 @@
+# fit_radec.py - Run OD on real angles data.
+# Copyright (C) 2020 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
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# 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 configure, add_station, build_measurement, Filter, MeasurementType, Constant
+from orbdetpy.conversion import get_J2000_epoch_offset, get_UTC_string
+from orbdetpy.estimation import determine_orbit
+from orbdetpy.plotting.estimation import plot
+
+# Real measurements: UTC, RA (rad), dec (rad)
+real_obs = [
+    ["2020-07-24T03:19:36.035Z", 1.073579, 1.381180],
+    ["2020-07-24T03:19:40.551Z", 0.950022, 1.386485],
+    ["2020-07-24T03:19:45.226Z", 0.813143, 1.388843],
+    ["2020-07-24T03:20:08.080Z", 0.177733, 1.348062],
+    ["2020-07-24T03:20:12.593Z", 0.083064, 1.330833],
+    ["2020-07-24T03:20:17.135Z", 6.282924, 1.311199],
+    ["2020-07-24T03:20:35.638Z", 6.043317, 1.213738],
+    ["2020-07-24T03:20:40.159Z", 6.002756, 1.186840],
+    ["2020-07-24T03:20:44.824Z", 5.966159, 1.158158],
+    ["2020-07-24T03:20:49.538Z", 5.933827, 1.128431],
+    ["2020-07-24T03:20:54.236Z", 5.905503, 1.098115],
+    ["2020-07-24T03:20:58.770Z", 5.881280, 1.068338],
+    ["2020-07-24T03:21:07.908Z", 5.840156, 1.007179],
+    ["2020-07-24T03:21:12.645Z", 5.822224, 0.975141],
+    ["2020-07-24T03:21:17.240Z", 5.806584, 0.943872],
+    ["2020-07-24T03:21:26.480Z", 5.779668, 0.880952],
+    ["2020-07-24T03:21:31.131Z", 5.768043, 0.849405],
+    ["2020-07-24T03:21:35.713Z", 5.757657, 0.818491],
+    ["2020-07-24T03:21:40.249Z", 5.748296, 0.788056],
+    ["2020-07-24T03:21:44.927Z", 5.739520, 0.757038],
+    ["2020-07-24T03:21:54.265Z", 5.724259, 0.696125],
+    ["2020-07-24T03:21:59.018Z", 5.717508, 0.665763],
+    ["2020-07-24T03:22:03.679Z", 5.711479, 0.636510],
+    ["2020-07-24T03:22:08.283Z", 5.706017, 0.608080],
+    ["2020-07-24T03:22:12.904Z", 5.701007, 0.580095],
+    ["2020-07-24T03:22:17.432Z", 5.696508, 0.553190],
+    ["2020-07-24T03:22:22.073Z", 5.692283, 0.526236],
+    ["2020-07-24T03:22:26.745Z", 5.688394, 0.499645],
+    ["2020-07-24T03:22:31.417Z", 5.684841, 0.473698],
+    ["2020-07-24T03:22:35.997Z", 5.681670, 0.448829],
+    ["2020-07-24T03:22:40.696Z", 5.678703, 0.423978],
+    ["2020-07-24T03:22:45.390Z", 5.675990, 0.399725],
+    ["2020-07-24T03:22:49.930Z", 5.673622, 0.376921]
+]
+
+# Configure orbit determination
+config = configure(prop_start=get_J2000_epoch_offset("2020-07-24T03:15:01.526Z"),
+                   prop_initial_state=[-3150507.866, -3090624.805, 5330071.203, 3600.098, -6482.160, -1626.428],
+                   prop_end=get_J2000_epoch_offset(real_obs[-1][0]))
+
+# Define ground station(s)
+add_station(config, "UTA-ASTRIA-NMSkies", 0.5743, -1.8418, 2225.0)
+
+# Define measurement types; RA/dec used here
+config.measurements[MeasurementType.RIGHT_ASCENSION].error[:] = [2.0*Constant.ARC_SECOND_TO_RAD]
+config.measurements[MeasurementType.DECLINATION].error[:] = [2.0*Constant.ARC_SECOND_TO_RAD]
+
+# Use Filter.EXTENDED_KALMAN for the EKF
+config.estm_filter = Filter.UNSCENTED_KALMAN
+
+# Build a list of Measurement objects, one for each RA/dec pair
+meas_obj = []
+for o in real_obs:
+    meas_obj.append(build_measurement(get_J2000_epoch_offset(o[0]), "UTA-ASTRIA-NMSkies", o[1:]))
+
+# Run OD. Fitting single object and hence the [0]
+fit = determine_orbit([config], [meas_obj])[0]
+# Check for estimation errors
+if (isinstance(fit, str)):
+    print(fit)
+    exit(1)
+
+# Plot OD results
+plot(config, meas_obj, fit, interactive=True, estim_param=False)
+for f in fit:
+    # print(f) to dump pre-fits/post-fits/covariances
+    print(get_UTC_string(f.time), f.estimated_state)

examples/propagate_tle.py

@@ -17,6 +17,7 @@
 import sys
 import argparse
 from datetime import datetime, timedelta
+from orbdetpy.conversion import get_UTC_string
 
 day0 = datetime.today()
 day1 = day0 + timedelta(days=1)
@@ -53,6 +54,6 @@
         print("\nObject {}:".format(elements[i][2:7]))
         i += 2
         for m in o.array:
-            print(m.time, m.true_state)
+            print(get_UTC_string(m.time), m.true_state)
 except Exception as exc:
     print(exc)

examples/test_conversion.py

@@ -14,7 +14,7 @@
 # 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
+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)
 
@@ -36,21 +36,21 @@
 
 # Position => Lat/Lon/Alt
 lla = pos_to_lla(Frame.GCRF, time, gcrf_pv[:3])
-print("pos_to_lla: Lat={:.2f}deg, Lon={:.2f}deg, Alt={:.2f}m".format(
+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)
-print("pv_to_elem: a={:.2f}m, e={:.4f}, i={:.2f}deg, O={:.2f}deg, w={:.2f}deg, M={:.2f}deg, theta={:.2f}deg".format(
+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[:-1])))
+    *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)
-print("radec_to_azel: Azimuth={:.3f}deg, Elevation={:.3f}deg".format(az/Constant.DEGREE_TO_RAD, el/Constant.DEGREE_TO_RAD))
+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)
-print("azel_to_radec: RA={:.3f}deg, Dec={:.3f}deg\n".format(r/Constant.DEGREE_TO_RAD, d/Constant.DEGREE_TO_RAD))
+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) 2019-2020 University of Texas
+# Copyright (C) 2020 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,21 +16,19 @@
 
 from numpy import diag
 from numpy.random import multivariate_normal
-from orbdetpy import (configure, add_facet, add_maneuver, add_station, AttitudeType, Filter,
-                      ManeuverTrigger, ManeuverType, MeasurementType, Constant)
+from orbdetpy import (configure, add_facet, add_maneuver, add_station, AttitudeType,
+                      Filter, ManeuverTrigger, ManeuverType, MeasurementType, Constant)
 from orbdetpy.ccsds import export_OEM, export_TDM
 from orbdetpy.conversion import get_J2000_epoch_offset
 from orbdetpy.estimation import determine_orbit
 from orbdetpy.propagation import propagate_orbits
-from orbdetpy.plotting.estimation import plot
+from orbdetpy.plotting.estimation import plot as estimation_plot
+from orbdetpy.plotting.simulation import plot as simulation_plot
 
 # Set up configuration for measurement generation
 cfg = configure(prop_start=get_J2000_epoch_offset("2019-05-01T00:00:00Z"),
-                prop_initial_state=[-23183898.259, 35170229.755, 43425.075,
-                                    -2566.938, -1692.19, 138.948],
-                prop_end=get_J2000_epoch_offset("2019-05-01T01:00:00Z"),
-                prop_step=300.0,
-                sim_measurements=True)
+                prop_initial_state=[-23183898.259, 35170229.755, 43425.075, -2566.938, -1692.19, 138.948],
+                prop_end=get_J2000_epoch_offset("2019-05-01T01:00:00Z"), prop_step=300.0, sim_measurements=True)
 
 # Uncomment to define box-wing model
 #add_facet(cfg, Constant.PLUS_I, 3.0)
@@ -57,33 +55,50 @@
             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
 cfg.measurements[MeasurementType.AZIMUTH].error[:] = [Constant.ARC_SECOND_TO_RAD]
 cfg.measurements[MeasurementType.ELEVATION].error[:] = [Constant.ARC_SECOND_TO_RAD]
 
+# You can provide RANGE+RANGE_RATE or RANGE alone or RANGE_RATE alone
 #cfg.measurements[MeasurementType.RANGE].error[:] = [10.0]
 #cfg.measurements[MeasurementType.RANGE_RATE].error[:] = [0.1]
 
+# RIGHT_ASCENSION and DECLINATION must be given
 #cfg.measurements[MeasurementType.RIGHT_ASCENSION].error[:] = [Constant.ARC_SECOND_TO_RAD]
 #cfg.measurements[MeasurementType.DECLINATION].error[:] = [Constant.ARC_SECOND_TO_RAD]
 
+# Inertial position as "measurements"
+#cfg.measurements[MeasurementType.POSITION].error[:] = [100.0, 200.0, 300.0]
+
+# Inertial position/velocity as "measurements"
+#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
 
+# Uncomment to plot orbital elements etc. from the simulation
+#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"))
 
 # Perturb truth initial state before running OD
 cfg.prop_initial_state[:] = multivariate_normal(cfg.prop_initial_state, diag([1E6, 1E6, 1E6, 1E2, 1E2, 1E2]))
-
 # Specify non-zero step to get propagated states/covariances between measurement updates
 cfg.prop_step = 0.0
-cfg.estm_filter = Filter.EXTENDED_KALMAN
+cfg.estm_filter = Filter.UNSCENTED_KALMAN
 
-# Run OD on simulated measurements and plot residuals
+# Run OD on simulated measurements
 fit = determine_orbit([cfg], [meas])[0]
-plot(cfg, meas, fit, interactive=True, estim_param=False)
+# Check for estimation errors
+if (isinstance(fit, str)):
+    print(fit)
+    exit(1)
+
+# 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", "COVID-19"))
+#    fp.write(export_OEM(cfg, fit, "2020-001A", "COVID-19"))

examples/test_interpolation.py

@@ -15,7 +15,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 from orbdetpy import configure, Frame
-from orbdetpy.conversion import get_J2000_epoch_offset
+from orbdetpy.conversion import get_J2000_epoch_offset, get_UTC_string
 from orbdetpy.propagation import propagate_orbits
 from orbdetpy.utilities import interpolate_ephemeris
 
@@ -32,4 +32,4 @@
 # Interpolate over the same period at 1 minute intervals
 for i in interpolate_ephemeris(Frame.GCRF, times, states, len(times),
                                Frame.GCRF, cfg.prop_start, cfg.prop_end, 60.0):
-    print(i.time, i.true_state)
+    print(get_UTC_string(i.time), i.true_state)