diff --git a/examples/Propagation.ipynb b/examples/Propagation.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..efba6f8207e3adf33019530d7457938e5a0f2ad3 --- /dev/null +++ b/examples/Propagation.ipynb @@ -0,0 +1,622 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Propagation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Authors\n", + "\n", + "Lots of parts are directly from the orekit documentation on [propagation](https://www.orekit.org/site-orekit-latest/architecture/propagation.html), with some updates, simplifications and Pythonification by Petrus Hyvönen, SSC\n", + "\n", + "## Learning Goals\n", + "* *What is an orbit propagator*: What is the background, what types are there, and why\n", + "* *How do I propagate my satellite*: How is it implemented in Orekit\n", + "\n", + "## Keywords\n", + "orekit, propagation" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "from math import radians, degrees" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Initialize orkit and bring up the python-java interface" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import orekit\n", + "vm = orekit.initVM()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now set up the pointer to the orekit-data.zip file, using one of the helper files. The file should be in current directory if not specified otherwise." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime\n", + "setup_orekit_curdir()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we are set up to import and use objects from the orekit library. Packages can be imported as they were native Python packages" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Propagation " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Propagation is the prediction of the evolution of a system from an initial state. In Orekit, this initial state is represented by a SpacecraftState, which is a simple container for all needed information : orbit, mass, kinematics, attitude, date, frame etc." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The method of propagating a satellite orbit can be divided into three categories:\n", + "\n", + "- Analytical Propagators: These are based on mathematical analytical models, which commonly does not need so much computing power and are genereally fast but not neccessary precise in complex environments\n", + "- Numerical Propagators: These propagators are based on a numerical models where forces are integrated over time by a large number of calculations. Can handle complex models of different forces acting on a spacecraft\n", + "- Semianalytical: Semianalytical combines features of numerical and analytical method to get a good mix of accuracy and efficency." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Analytical Propagators " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In orekit there are a number of analytical propagators.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Keplerian Propagator\n", + "\n", + "This is a simple propagator that models a Keplerian orbit around a planet, based on the mass of the central body, µ= GM.\n", + "\n", + "The [Keplerian Orbit](https://www.orekit.org/site-orekit-latest/apidocs/org/orekit/propagation/analytical/KeplerianPropagator.html) at the orekit documentation API shows the usage. A basic example is:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "from org.orekit.orbits import KeplerianOrbit, PositionAngle\n", + "from org.orekit.propagation.analytical import KeplerianPropagator\n", + "from org.orekit.time import AbsoluteDate, TimeScalesFactory\n", + "from org.orekit.utils import Constants\n", + "from org.orekit.frames import FramesFactory" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "utc = TimeScalesFactory.getUTC()" + ] + }, + { + "cell_type": "code", + "execution_count": 159, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<KeplerianOrbit: Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 0.0;}>" + ] + }, + "execution_count": 159, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ra = 400 * 1000 # Apogee\n", + "rp = 500 * 1000 # Perigee\n", + "i = radians(87.0) # inclination\n", + "omega = radians(20.0) # perigee argument\n", + "raan = radians(10.0) # right ascension of ascending node\n", + "lv = radians(0.0) # True anomaly\n", + "\n", + "epochDate = AbsoluteDate(2020, 1, 1, 0, 0, 00.000, utc)\n", + "initialDate = epochDate\n", + "\n", + "a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0 \n", + "e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a\n", + "\n", + "## Inertial frame where the satellite is defined\n", + "inertialFrame = FramesFactory.getEME2000()\n", + "\n", + "## Orbit construction as Keplerian\n", + "initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,\n", + " PositionAngle.TRUE,\n", + " inertialFrame, epochDate, Constants.WGS84_EARTH_MU)\n", + "initialOrbit" + ] + }, + { + "cell_type": "code", + "execution_count": 160, + "metadata": {}, + "outputs": [], + "source": [ + "propagator = KeplerianPropagator(initialOrbit)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can show the initial state that the propagator will start from:" + ] + }, + { + "cell_type": "code", + "execution_count": 164, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<SpacecraftState: SpacecraftState{orbit=Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 0.0;}, attitude=org.orekit.attitudes.Attitude@db57326, mass=1000.0, additional={}}>" + ] + }, + "execution_count": 164, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "propagator.getInitialState()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A basic way to execute the propagator is through the propagate(start, end) method. In this example we propagate the orbit for 48 hours from initialDate." + ] + }, + { + "cell_type": "code", + "execution_count": 165, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<SpacecraftState: SpacecraftState{orbit=Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 11079.353102947607;}, attitude=org.orekit.attitudes.Attitude@18a136ac, mass=1000.0, additional={}}>" + ] + }, + "execution_count": 165, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "propagator.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that only one variable changed, which?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Eckstein-Hechler Propagator\n", + "\n", + "The Eckstein-Hechler propagator is an analytical propagator that can use a significant more elaborated model of the gravity field, including the J2 to J6 potential zonal coefficients. It uses mean orbital parameters to compute the new position.\n", + "\n", + "The EH propagator is only applicable for near circular orbits, typically used for LEO satellites.\n", + "\n", + "The [orekit documentation for the EH propagator]() gives more details.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 166, + "metadata": {}, + "outputs": [], + "source": [ + "from org.orekit.propagation.analytical import EcksteinHechlerPropagator\n", + "from org.orekit.orbits import OrbitType" + ] + }, + { + "cell_type": "code", + "execution_count": 167, + "metadata": {}, + "outputs": [], + "source": [ + "propagator_eh = EcksteinHechlerPropagator(initialOrbit, \n", + " Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,\n", + " Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20,\n", + " Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40,\n", + " Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60)" + ] + }, + { + "cell_type": "code", + "execution_count": 168, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<SpacecraftState: SpacecraftState{orbit=Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 0.0;}, attitude=org.orekit.attitudes.Attitude@176b3f44, mass=1000.0, additional={}}>" + ] + }, + "execution_count": 168, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "propagator_eh.getInitialState()" + ] + }, + { + "cell_type": "code", + "execution_count": 170, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<SpacecraftState: SpacecraftState{orbit=Cartesian parameters: {P(3578228.3421075065, 268431.7711005597, -5823469.25383565), V(6411.063408954351, 1243.983278874011, 3927.1303963350906)}, attitude=org.orekit.attitudes.Attitude@5d12a356, mass=1000.0, additional={}}>" + ] + }, + "execution_count": 170, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "end_state = propagator_eh.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))\n", + "end_state" + ] + }, + { + "cell_type": "code", + "execution_count": 171, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<Orbit: Keplerian parameters: {a: 6816647.287563898; e: 0.008498145471038656; i: 86.9974042555557; pa: -172.94021521488222; raan: 9.17350883469523; v: 114.45253754444634;}>" + ] + }, + "execution_count": 171, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "OrbitType.KEPLERIAN.convertType(end_state.getOrbit())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## SGP4 / SDP4 Propagator" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This analytical propagator is dedicated to propagation of Two-Line Elements (TLE)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "See separate example" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Numerical Propagators" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Numerical propagation is one of the most important parts of the Orekit project. Based on Hipparchus ordinary differential equations integrators, the NumericalPropagator class realizes the interface between space mechanics and mathematical resolutions. Despite its utilization seems daunting on first sight, it is in fact quite straigthforward to use." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simple Propagation of Equation of Motion " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The mathematical problem to integrate is a dimension-seven time-derivative equations system. The six first elements of the state vector are the orbital parameters, which may be any orbit type (KeplerianOrbit, CircularOrbit, EquinoctialOrbit or CartesianOrbit) in meters and radians, and the last element is the mass in kilograms. It is possible to have more elements in the state vector if AdditionalEquations have been added (typically PartialDerivativesEquations which is an implementation of AdditionalEquations devoted to integration of Jacobian matrices). The time derivatives are computed automatically by the Orekit using the Gauss equations for the first parameters corresponding to the selected orbit type and the flow rate for mass evolution during maneuvers. The user only needs to register the various force models needed for the simulation. Various force models are already available in the library and specialized ones can be added by users easily for specific needs.\n", + "\n", + "The integrators (first order integrators) provided by Hipparchus need the state vector at t0, the state vector first time derivative at t0, and then calculates the next step state vector, and asks for the next first time derivative, etc. until it reaches the final asked date. These underlying numerical integrators can also be configured. Typical tuning parameters for adaptive stepsize integrators are the min, max and perhaps start step size as well as the absolute and/or relative errors thresholds. \n", + "\n", + "The following code snippet shows a typical setting for Low Earth Orbit propagation:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The numerical propagation is based on an integrator with variable step size. These are specified, as other time units in Orekit, in seconds." + ] + }, + { + "cell_type": "code", + "execution_count": 172, + "metadata": {}, + "outputs": [], + "source": [ + "from org.orekit.propagation.numerical import NumericalPropagator\n", + "from org.hipparchus.ode.nonstiff import DormandPrince853Integrator\n", + "from org.orekit.propagation import SpacecraftState\n", + "from org.orekit.bodies import OneAxisEllipsoid\n", + "from org.orekit.utils import IERSConventions\n", + "from org.orekit.forces.gravity.potential import GravityFieldFactory\n", + "from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel\n", + "\n", + "from orekit import JArray_double" + ] + }, + { + "cell_type": "code", + "execution_count": 185, + "metadata": {}, + "outputs": [], + "source": [ + "minStep = 0.001\n", + "maxstep = 1000.0\n", + "initStep = 60.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The spatial tolerance can be specified (meters)" + ] + }, + { + "cell_type": "code", + "execution_count": 186, + "metadata": {}, + "outputs": [], + "source": [ + "positionTolerance = 1.0 " + ] + }, + { + "cell_type": "code", + "execution_count": 187, + "metadata": {}, + "outputs": [], + "source": [ + "tolerances = NumericalPropagator.tolerances(positionTolerance, \n", + " initialOrbit, \n", + " initialOrbit.getType())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The actual integrator, in this case DormandPrince853, is part of the Hipparchos library. Note that the tolerances needs casting in Python to an array of doubles (floats)." + ] + }, + { + "cell_type": "code", + "execution_count": 188, + "metadata": {}, + "outputs": [], + "source": [ + "integrator = DormandPrince853Integrator(minStep, maxstep, \n", + " JArray_double.cast_(tolerances[0]), # Double array of doubles needs to be casted in Python\n", + " JArray_double.cast_(tolerances[1]))\n", + "integrator.setInitialStepSize(initStep)" + ] + }, + { + "cell_type": "code", + "execution_count": 189, + "metadata": {}, + "outputs": [], + "source": [ + "satellite_mass = 100.0 # The models need a spacecraft mass, unit kg.\n", + "initialState = SpacecraftState(initialOrbit, satellite_mass) " + ] + }, + { + "cell_type": "code", + "execution_count": 202, + "metadata": {}, + "outputs": [], + "source": [ + "propagator_num = NumericalPropagator(integrator)\n", + "propagator_num.setOrbitType(OrbitType.CARTESIAN)\n", + "propagator_num.setInitialState(initialState)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the propagator to make sense it needs some forces acting on the satellite. Here we are adding a gravity field model.\n", + "\n", + "For a more detailed propagation, other force models can be added." + ] + }, + { + "cell_type": "code", + "execution_count": 203, + "metadata": {}, + "outputs": [], + "source": [ + "gravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)\n", + "propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider))" + ] + }, + { + "cell_type": "code", + "execution_count": 204, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<Orbit: Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 0.0;}>" + ] + }, + "execution_count": 204, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "OrbitType.KEPLERIAN.convertType(propagator_num.getInitialState().getOrbit())" + ] + }, + { + "cell_type": "code", + "execution_count": 215, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<SpacecraftState: SpacecraftState{orbit=Cartesian parameters: {P(3580817.970383157, 268551.4783581181, -5822001.791514523), V(6409.3095167843, 1243.4899776182563, 3929.8276034988135)}, attitude=org.orekit.attitudes.Attitude@79c7532f, mass=100.0, additional={}}>" + ] + }, + "execution_count": 215, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "end_state = propagator_num.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))\n", + "end_state" + ] + }, + { + "cell_type": "code", + "execution_count": 217, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "<Orbit: Keplerian parameters: {a: 6816573.271804827; e: 0.008516228214247105; i: 86.99611217388401; pa: -173.06144237396225; raan: 9.169665558677691; v: 114.5987064928954;}>" + ] + }, + "execution_count": 217, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "OrbitType.KEPLERIAN.convertType(end_state.getOrbit()) # Note that this is the Osculating orbit!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}