Skip to content
Snippets Groups Projects
Commit 77b0c239 authored by Petrus Hyvönen's avatar Petrus Hyvönen
Browse files

Added example on general propagation

parent a2fc4e62
No related branches found
No related tags found
No related merge requests found
%% Cell type:markdown id: tags:
# Propagation
%% Cell type:markdown id: tags:
## Authors
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
## Learning Goals
* *What is an orbit propagator*: What is the background, what types are there, and why
* *How do I propagate my satellite*: How is it implemented in Orekit
## Keywords
orekit, propagation
%% Cell type:code id: tags:
``` python
%matplotlib inline
from math import radians, degrees
```
%% Cell type:markdown id: tags:
Initialize orkit and bring up the python-java interface
%% Cell type:code id: tags:
``` python
import orekit
vm = orekit.initVM()
```
%% Cell type:markdown id: tags:
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 id: tags:
``` python
from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()
```
%% Cell type:markdown id: tags:
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 id: tags:
# Propagation
%% Cell type:markdown id: tags:
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 id: tags:
The method of propagating a satellite orbit can be divided into three categories:
- 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
- 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
- Semianalytical: Semianalytical combines features of numerical and analytical method to get a good mix of accuracy and efficency.
%% Cell type:markdown id: tags:
# Analytical Propagators
%% Cell type:markdown id: tags:
In orekit there are a number of analytical propagators.
%% Cell type:markdown id: tags:
## Keplerian Propagator
This is a simple propagator that models a Keplerian orbit around a planet, based on the mass of the central body, µ= GM.
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 id: tags:
``` python
from org.orekit.orbits import KeplerianOrbit, PositionAngle
from org.orekit.propagation.analytical import KeplerianPropagator
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants
from org.orekit.frames import FramesFactory
```
%% Cell type:code id: tags:
``` python
utc = TimeScalesFactory.getUTC()
```
%% Cell type:code id: tags:
``` python
ra = 400 * 1000 # Apogee
rp = 500 * 1000 # Perigee
i = radians(87.0) # inclination
omega = radians(20.0) # perigee argument
raan = radians(10.0) # right ascension of ascending node
lv = radians(0.0) # True anomaly
epochDate = AbsoluteDate(2020, 1, 1, 0, 0, 00.000, utc)
initialDate = epochDate
a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a
## Inertial frame where the satellite is defined
inertialFrame = FramesFactory.getEME2000()
## Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
PositionAngle.TRUE,
inertialFrame, epochDate, Constants.WGS84_EARTH_MU)
initialOrbit
```
%% Output
<KeplerianOrbit: Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 0.0;}>
%% Cell type:code id: tags:
``` python
propagator = KeplerianPropagator(initialOrbit)
```
%% Cell type:markdown id: tags:
We can show the initial state that the propagator will start from:
%% Cell type:code id: tags:
``` python
propagator.getInitialState()
```
%% Output
<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={}}>
%% Cell type:markdown id: tags:
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 id: tags:
``` python
propagator.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))
```
%% Output
<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={}}>
%% Cell type:markdown id: tags:
Note that only one variable changed, which?
%% Cell type:markdown id: tags:
## Eckstein-Hechler Propagator
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.
The EH propagator is only applicable for near circular orbits, typically used for LEO satellites.
The [orekit documentation for the EH propagator]() gives more details.
%% Cell type:code id: tags:
``` python
from org.orekit.propagation.analytical import EcksteinHechlerPropagator
from org.orekit.orbits import OrbitType
```
%% Cell type:code id: tags:
``` python
propagator_eh = EcksteinHechlerPropagator(initialOrbit,
Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20,
Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40,
Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60)
```
%% Cell type:code id: tags:
``` python
propagator_eh.getInitialState()
```
%% Output
<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={}}>
%% Cell type:code id: tags:
``` python
end_state = propagator_eh.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))
end_state
```
%% Output
<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={}}>
%% Cell type:code id: tags:
``` python
OrbitType.KEPLERIAN.convertType(end_state.getOrbit())
```
%% Output
<Orbit: Keplerian parameters: {a: 6816647.287563898; e: 0.008498145471038656; i: 86.9974042555557; pa: -172.94021521488222; raan: 9.17350883469523; v: 114.45253754444634;}>
%% Cell type:markdown id: tags:
## SGP4 / SDP4 Propagator
%% Cell type:markdown id: tags:
This analytical propagator is dedicated to propagation of Two-Line Elements (TLE).
%% Cell type:markdown id: tags:
See separate example
%% Cell type:markdown id: tags:
# Numerical Propagators
%% Cell type:markdown id: tags:
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 id: tags:
## Simple Propagation of Equation of Motion
%% Cell type:markdown id: tags:
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.
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.
The following code snippet shows a typical setting for Low Earth Orbit propagation:
%% Cell type:markdown id: tags:
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 id: tags:
``` python
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.propagation import SpacecraftState
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.utils import IERSConventions
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from orekit import JArray_double
```
%% Cell type:code id: tags:
``` python
minStep = 0.001
maxstep = 1000.0
initStep = 60.0
```
%% Cell type:markdown id: tags:
The spatial tolerance can be specified (meters)
%% Cell type:code id: tags:
``` python
positionTolerance = 1.0
```
%% Cell type:code id: tags:
``` python
tolerances = NumericalPropagator.tolerances(positionTolerance,
initialOrbit,
initialOrbit.getType())
```
%% Cell type:markdown id: tags:
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 id: tags:
``` python
integrator = DormandPrince853Integrator(minStep, maxstep,
JArray_double.cast_(tolerances[0]), # Double array of doubles needs to be casted in Python
JArray_double.cast_(tolerances[1]))
integrator.setInitialStepSize(initStep)
```
%% Cell type:code id: tags:
``` python
satellite_mass = 100.0 # The models need a spacecraft mass, unit kg.
initialState = SpacecraftState(initialOrbit, satellite_mass)
```
%% Cell type:code id: tags:
``` python
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(OrbitType.CARTESIAN)
propagator_num.setInitialState(initialState)
```
%% Cell type:markdown id: tags:
For the propagator to make sense it needs some forces acting on the satellite. Here we are adding a gravity field model.
For a more detailed propagation, other force models can be added.
%% Cell type:code id: tags:
``` python
gravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)
propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider))
```
%% Cell type:code id: tags:
``` python
OrbitType.KEPLERIAN.convertType(propagator_num.getInitialState().getOrbit())
```
%% Output
<Orbit: Keplerian parameters: {a: 6828137.0; e: -0.007322641593160872; i: 86.99999999999999; pa: 20.0; raan: 10.0; v: 0.0;}>
%% Cell type:code id: tags:
``` python
end_state = propagator_num.propagate(initialDate, initialDate.shiftedBy(3600.0 * 48))
end_state
```
%% Output
<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={}}>
%% Cell type:code id: tags:
``` python
OrbitType.KEPLERIAN.convertType(end_state.getOrbit()) # Note that this is the Osculating orbit!
```
%% Output
<Orbit: Keplerian parameters: {a: 6816573.271804827; e: 0.008516228214247105; i: 86.99611217388401; pa: -173.06144237396225; raan: 9.169665558677691; v: 114.5987064928954;}>
%% Cell type:code id: tags:
``` python
```
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment