From 3d45ea57ab9759c30b96153b94d135cd7c493de5 Mon Sep 17 00:00:00 2001
From: petrush <petrus.hyvonen@sscspace.com>
Date: Wed, 23 Nov 2016 16:38:10 +0100
Subject: [PATCH] Added example of numerical propagator.
---
examples/Example numerical prop.ipynb | 284 ++++++++++++++++++++++++++
1 file changed, 284 insertions(+)
create mode 100644 examples/Example numerical prop.ipynb
diff --git a/examples/Example numerical prop.ipynb b/examples/Example numerical prop.ipynb
new file mode 100644
index 0000000..16d9724
--- /dev/null
+++ b/examples/Example numerical prop.ipynb
@@ -0,0 +1,284 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Numerical Propagation Example"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Populating the interactive namespace from numpy and matplotlib\n"
+ ]
+ }
+ ],
+ "source": [
+ "%pylab inline"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import orekit\n",
+ "orekit.initVM()\n",
+ "\n",
+ "from orekit.pyhelpers import setup_orekit_curdir\n",
+ "setup_orekit_curdir()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "from java.util import Arrays\n",
+ "from orekit import JArray_double"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Do the orekit imports"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "from org.orekit.bodies import OneAxisEllipsoid\n",
+ "from org.orekit.frames import FramesFactory\n",
+ "from org.orekit.data import DataProvidersManager, ZipJarCrawler\n",
+ "from org.orekit.time import TimeScalesFactory, AbsoluteDate\n",
+ "from org.orekit.orbits import KeplerianOrbit\n",
+ "from org.orekit.utils import Constants\n",
+ "from org.orekit.propagation.analytical import EcksteinHechlerPropagator\n",
+ "from org.orekit.propagation.analytical.tle import TLEPropagator\n",
+ "from org.orekit.propagation.conversion import FiniteDifferencePropagatorConverter\n",
+ "from org.orekit.propagation.conversion import TLEPropagatorBuilder\n",
+ "from datetime import datetime\n",
+ "from org.orekit.propagation import SpacecraftState\n",
+ "from org.orekit.orbits import OrbitType, PositionAngle\n",
+ "from org.orekit.propagation.numerical import NumericalPropagator\n",
+ "from org.hipparchus.ode.nonstiff import DormandPrince853Integrator\n",
+ "from org.orekit.forces.gravity.potential import GravityFieldFactory\n",
+ "from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel\n",
+ "from org.orekit.utils import IERSConventions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "import matplotlib.dates as mdates\n",
+ "from math import radians\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "# Some Constants\n",
+ "ae = Constants.WGS84_EARTH_EQUATORIAL_RADIUS\n",
+ "mu = Constants.WGS84_EARTH_MU\n",
+ "utc = TimeScalesFactory.getUTC()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Initial osculating orbit"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "ra = 814 * 1000 # km Apogee!\n",
+ "rp = 786 * 1000 # km Perigee!\n",
+ "i = radians(98.55) # inclination\n",
+ "omega = radians(90.0) # perigee argument\n",
+ "raan = radians(5.1917) # right ascension of ascending node\n",
+ "lv = radians(0.0567634) # True anomaly\n",
+ "epochDate = AbsoluteDate(2012, 01, 26, 16, 00, 00.000, utc)\n",
+ "\n",
+ "a = (rp + ra + 2 * ae) / 2.0 # semi major axis in km\n",
+ "e = 1.0 - (rp + ae) / 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, mu)\n",
+ "satellite_mass = 100.0 # kg\n",
+ "initialState = SpacecraftState(initialOrbit, satellite_mass) "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Set up propagator"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "minStep = 0.001;\n",
+ "maxstep = 1000.0;\n",
+ "initStep = 60.0"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "positionTolerance = 1.0 \n",
+ "orbitType = OrbitType.CARTESIAN\n",
+ "tol = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "integrator = DormandPrince853Integrator(minStep, maxstep, \n",
+ " JArray_double.cast_(tol[0]), # Double array of doubles needs to be casted in Python\n",
+ " JArray_double.cast_(tol[1]))\n",
+ "integrator.setInitialStepSize(initStep)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "propagator_num = NumericalPropagator(integrator)\n",
+ "propagator_num.setOrbitType(orbitType)\n",
+ "propagator_num.setInitialState(initialState)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 12,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed\n",
+ "earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,\n",
+ " Constants.WGS84_EARTH_FLATTENING,\n",
+ " itrf)\n",
+ "gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)\n",
+ "propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Propagate orbit numerically"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "## Create time vector\n",
+ "startDate = AbsoluteDate(2012, 01, 26, 11, 00, 00.000, utc)\n",
+ "\n",
+ "#Overall duration in seconds for extrapolation\n",
+ "duration = 86400 * 3 \n",
+ "step_time = 60 * 3\n",
+ "\n",
+ "# Time array in orekit AbsoluteDate format\n",
+ "t = [startDate.shiftedBy(float(dt)) \\\n",
+ " for dt in np.arange(0, duration, step_time)]\n",
+ "\n",
+ "state = [propagator_num.propagate(tt) for tt in t]\n",
+ "pv_org = [propagator_num.propagate(tt).getPVCoordinates() for tt in t]"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python [default]",
+ "language": "python",
+ "name": "python2"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 2
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython2",
+ "version": "2.7.12"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
--
GitLab