Skip to content
GitLab
Commit c01a858b authored by Petrus Hyvönen's avatar Petrus Hyvönen
Browse files

Minor updates to examples

parent 3936c4b6
Expand all Hide whitespace changes
Inline Side-by-side
examples/Example EarthObservation - Attitude Sequence.ipynb
View file @ c01a858b
%% Cell type:markdown id: tags:
# Attitude Sequence
THIS EXAMPLE IS CURRENTLY NOT WORKING - NEEDS UPDATE FOR OREKIT 8.0
%% Cell type:markdown id: tags:
This example is based on the attitude tutorial in the orekit library and webpage https://www.orekit.org/static/tutorial/attitude.html by Luc Maisonobe written in Java. This version of the example is based on ipython notebook, using a python-wrapped version of the orekit library, as a demonstration of the python capabilities and as a tutorial. Conversion made by Petrus Hyvönen, SSC, 2014.
%% Cell type:markdown id: tags:
This example uses AttitudesSequence that enables easy switching between attitude laws on event occurrences when propagating some SpacecraftState.
%% Cell type:code id: tags:
``` python
%pylab inline
figsize(12,12)
```
%% Output
Populating the interactive namespace from numpy and matplotlib
%% Cell type:code id: tags:
``` python
#initialize orekit and JVM
import orekit
orekit.initVM()
```
%% Output
<jcc.JCCEnv at 0x278fb40>
<jcc.JCCEnv at 0x108337228>
%% Cell type:code id: tags:
``` python
# setup the orekit data loading, the file orekit-data.zip shall be in same directory as notebook.
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
```
%% Cell type:code id: tags:
``` python
from org.hipparchus.geometry.euclidean.threed import RotationOrder;
from org.hipparchus.geometry.euclidean.threed import Vector3D;
```
%% Cell type:code id: tags:
``` python
from org.orekit.attitudes import AttitudeProvider;
from org.orekit.attitudes import AttitudesSequence;
from org.orekit.attitudes import LofOffset;
from org.orekit.bodies import CelestialBodyFactory;
from org.orekit.errors import OrekitException;
from org.orekit.errors import PropagationException;
from org.orekit.frames import FramesFactory;
from org.orekit.frames import LOFType;
from org.orekit.orbits import KeplerianOrbit;
from org.orekit.orbits import Orbit;
from org.orekit.propagation import Propagator;
from org.orekit.propagation import SpacecraftState;
from org.orekit.propagation.analytical import EcksteinHechlerPropagator;
from org.orekit.propagation.events import EclipseDetector;
from org.orekit.propagation.events import EventDetector;
from org.orekit.propagation.events.handlers import EventHandler;
from org.orekit.propagation.sampling import OrekitFixedStepHandler;
from org.orekit.time import AbsoluteDate;
from org.orekit.time import TimeScalesFactory;
from org.orekit.utils import Constants;
from org.orekit.utils import PVCoordinates;
from org.orekit.utils import PVCoordinatesProvider;
```
%% Cell type:code id: tags:
``` python
from org.orekit.python import PythonEventHandler, PythonOrekitFixedStepHandler
```
%% Cell type:markdown id: tags:
Let's set up an initial state as:
- a date in UTC time scale
- an orbit defined by the position and the velocity of the spacecraft in the EME2000 inertial frame and an associated central attraction coefficient chosen among many physical constants available in Orekit.
The initial orbit is here defined as a KeplerianOrbit.
%% Cell type:code id: tags:
``` python
initialDate = AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC())
initialDate = AbsoluteDate(2004, 1, 1, 23, 30, 00.000, TimeScalesFactory.getUTC())
position = Vector3D(-6142438.668, 3492467.560, -25767.25680);
velocity = Vector3D(505.8479685, 942.7809215, 7435.922231);
initialOrbit = KeplerianOrbit(PVCoordinates(position, velocity),
FramesFactory.getEME2000(), initialDate,
Constants.EIGEN5C_EARTH_MU);
```
%% Cell type:code id: tags:
``` python
print initialOrbit
print(initialOrbit)
```
%% Output
keplerian parameters: {a: 7069220.386682823; e: 4.777356060557311E-4; i: 98.18525099174988; pa: 13.741061002484528; raan: 150.34825333049; v: -13.952151446378437;}
%% Cell type:markdown id: tags:
More details on the orbit representation can be found in the orbits section of the library architecture documentation. Note that in the python wrapping of orekit the format of numbers is important, such that floats needs to be specified with a decimal point such as in the AbsoluteDate above where seconds are a float.
Let's define a couple of AttitudeLaw for our satellite, one for daytime and one for nighttime, built upon Local Orbital Frame (Lof) Offset laws. VVLH is a coodrinate system based on Vehicle Velocity, Local Horizontal frame (Z axis aligned with nadir vector, X-Y is in horizontal plane, where the X-axis is along the velocity vector). The offset is defined in Cardan angles (RotationOrder.XYZ).
%% Cell type:code id: tags:
``` python
dayObservationLaw = LofOffset(initialOrbit.getFrame(),
LOFType.VVLH,
RotationOrder.XYZ,
math.radians(20), math.radians(40), 0.0)
```
%% Cell type:markdown id: tags:
In night we want the pointing to be in LOF aligned.
%% Cell type:code id: tags:
``` python
nightRestingLaw = LofOffset(initialOrbit.getFrame(), LOFType.VVLH)
```
%% Cell type:markdown id: tags:
Let's also define some EventDetectors for the switching. For this tutorial's requirements, two specialized event detectors are built
upon an EclipseDetector:
one, dayNightEvent, to detect the day to night transition,
another, nightDayEvent, to detect the night to day transition.
This is done by overriding the eventOccurred method of the standard EclipseDetector. To override methods in a java class in Python,
a special wrapped version of the class is used, called PythonElevationDetector.
%% Cell type:code id: tags:
``` python
class myNightEclipseDetector(PythonEventHandler):
def eventOccurred(self, s, detector, increasing):
if not increasing:
print s.getDate()," : event occurred, entering eclipse => switching to night law"
print(s.getDate()," : event occurred, entering eclipse => switching to night law")
return EventHandler.Action.CONTINUE
def resetState(self, detector, oldState):
return oldState;
```
%% Cell type:code id: tags:
``` python
class myDayEclipseDetector(PythonEventHandler):
def eventOccurred(self, s, detector, increasing):
if increasing:
print s.getDate()," : event occurred, exiting eclipse => switching to day law"
print(s.getDate()," : event occurred, exiting eclipse => switching to day law")
return EventHandler.Action.CONTINUE
def resetState(self, detector, oldState):
return oldState
```
%% Cell type:markdown id: tags:
For the eclipse calculation we want to use the Sun and the Earth.
%% Cell type:code id: tags:
``` python
sun = CelestialBodyFactory.getSun()
earth = CelestialBodyFactory.getEarth()
```
%% Cell type:markdown id: tags:
Create the actual detectors with the the handlers specified above.
%% Cell type:code id: tags:
``` python
dayNightEvent = EclipseDetector(sun, 696000000., earth, Constants.WGS84_EARTH_EQUATORIAL_RADIUS)
dayNightEvent = dayNightEvent.withHandler(myNightEclipseDetector().of_(EclipseDetector))
```
%% Cell type:code id: tags:
``` python
nightDayEvent = EclipseDetector(sun, 696000000., earth, Constants.WGS84_EARTH_EQUATORIAL_RADIUS)
nightDayEvent = nightDayEvent.withHandler(myDayEclipseDetector().of_(EclipseDetector))
```
%% Cell type:markdown id: tags:
More details on event detectors can be found in the propagation section of the library architecture documentation.
An AttitudesSequence is then defined, for the sake of this tutorial, by adding two switching conditions acting as a simple loop:
- the first one enables the transition from dayObservationLaw to nightRestingLaw when a decreasing dayNightEvent occurs,
- the second one enables the transition from nightRestingLaw to dayObservationLaw when an increasing nightDayEvent occurs.
%% Cell type:code id: tags:
``` python
attitudesSequence = AttitudesSequence()
attitudesSequence.addSwitchingCondition(dayObservationLaw, dayNightEvent, False, True, nightRestingLaw)
attitudesSequence.addSwitchingCondition(nightRestingLaw, nightDayEvent, True, False, dayObservationLaw)
attitudesSequence.addSwitchingCondition(dayObservationLaw, nightRestingLaw, dayNightEvent, False, True, 10.0)
attitudesSequence.addSwitchingCondition(nightRestingLaw, dayObservationLaw, nightDayEvent, True, False, 10.0)
```
%% Output
---------------------------------------------------------------------------
InvalidArgsError Traceback (most recent call last)
<ipython-input-21-320c3d02c46d> in <module>()
<ipython-input-21-075de152e259> in <module>()
1 attitudesSequence = AttitudesSequence()
----> 2 attitudesSequence.addSwitchingCondition(dayObservationLaw, dayNightEvent, False, True, nightRestingLaw)
3 attitudesSequence.addSwitchingCondition(nightRestingLaw, nightDayEvent, True, False, dayObservationLaw)
InvalidArgsError: (<type 'AttitudesSequence'>, 'addSwitchingCondition', (<LofOffset: org.orekit.attitudes.LofOffset@17c386de>, <EclipseDetector: org.orekit.propagation.events.EclipseDetector@45dd4eda>, False, True, <LofOffset: org.orekit.attitudes.LofOffset@5af97850>))
----> 2 attitudesSequence.addSwitchingCondition(dayObservationLaw, nightRestingLaw, dayNightEvent, False, True, 10.0)
3 #attitudesSequence.addSwitchingCondition(nightRestingLaw, dayObservationLaw, nightDayEvent, True, False, 10.0)
InvalidArgsError: (<type 'AttitudesSequence'>, 'addSwitchingCondition', (<LofOffset: org.orekit.attitudes.LofOffset@17c386de>, <LofOffset: org.orekit.attitudes.LofOffset@5af97850>, <EclipseDetector: org.orekit.propagation.events.EclipseDetector@45dd4eda>, False, True, 10.0))
%% Cell type:markdown id: tags:
An AttitudesSequence needs at least one switching condition to be meaningful, but there is no upper limit.
An active AttitudeLaw may have several switch events and next law settings, leading to different activation patterns depending on which event is triggered first.
So, don't forget to set the current active law according to the current state:
%% Cell type:code id: tags:
``` python
if (dayNightEvent.g(SpacecraftState(initialOrbit)) >= 0):
# initial position is in daytime
attitudesSequence.resetActiveProvider(dayObservationLaw);
else:
# initial position is in nighttime
attitudesSequence.resetActiveProvider(nightRestingLaw);
```
%% Cell type:markdown id: tags:
Now, let's choose some propagator to compute the spacecraft motion. We will use an EcksteinHechlerPropagator based on the analytical Eckstein-Hechler model. The propagator is built upon the initialOrbit, the attitudeSequence and physical constants for the gravity field potential.
%% Cell type:code id: tags:
``` python
propagator = EcksteinHechlerPropagator(initialOrbit, attitudesSequence,
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:markdown id: tags:
The attitudeSequence must register all the switching events before propagation.
%% Cell type:code id: tags:
``` python
# Register the switching events to the propagator
attitudesSequence.registerSwitchEvents(propagator)
```
%% Cell type:markdown id: tags:
The propagator operating mode is set to master mode with fixed step. The implementation of the interface OrekitFixedStepHandler aims to define the handleStep method called within the loop. For the purpose of this tutorial, the handleStep method will print at the current date two angles, the first one indicates if the spacecraft is eclipsed while the second informs about the current attitude law.
%% Cell type:code id: tags:
``` python
class mystephandler(PythonOrekitFixedStepHandler):
eclipseAngles = []
pointingOffsets = []
dates = []
def init(self,s0, t):
pass
def handleStep(self,currentState, isLast):
# the Earth position in spacecraft frame should be along spacecraft Z axis
# during nigthtime and away from it during daytime due to roll and pitch offsets
earth = currentState.toTransform().transformPosition(Vector3D.ZERO)
pointingOffset = Vector3D.angle(earth, Vector3D.PLUS_K)
# the g function is the eclipse indicator, its an angle between Sun and Earth limb,
# positive when Sun is outside of Earth limb, negative when Sun is hidden by Earth limb
eclipseAngle = dayNightEvent.g(currentState)
print "%s %6.3f %6.1f" % (currentState.getDate(), eclipseAngle, math.degrees(pointingOffset))
self.eclipseAngles.append(eclipseAngle)
self.pointingOffsets.append(math.degrees(pointingOffset))
self.dates.append(currentState.getDate())
```
%% Output
File "<ipython-input-24-7076b8b616e7>", line 19
print "%s %6.3f %6.1f" % (currentState.getDate(), eclipseAngle, math.degrees(pointingOffset))
^
SyntaxError: invalid syntax
%% Cell type:code id: tags:
``` python
handler = mystephandler()
propagator.setMasterMode(180.0, handler)
```
%% Cell type:markdown id: tags:
More details on propagation modes can be found in the propagation section of the library architecture documentation.
Finally, the propagator is just asked to propagate for a given duration.
%% Cell type:code id: tags:
``` python
# Propagate from the initial date for the fixed duration
finalState = propagator.propagate(initialDate.shiftedBy(12600.))
```
%% Output
2004-01-01T23:30:00.000 -0.203 0.0
2004-01-01T23:33:00.000 -0.311 0.0
2004-01-01T23:36:00.000 -0.392 0.0
2004-01-01T23:39:00.000 -0.437 0.0
2004-01-01T23:42:00.000 -0.439 0.0
2004-01-01T23:45:00.000 -0.397 0.0
2004-01-01T23:48:00.000 -0.318 0.0
2004-01-01T23:51:00.000 -0.212 0.0
2004-01-01T23:55:57.968 : event occurred, exiting eclipse => switching to day law
2004-01-01T23:54:00.000 -0.088 0.0
2004-01-01T23:57:00.000 0.048 44.0
2004-01-02T00:00:00.000 0.191 44.0
2004-01-02T00:03:00.000 0.338 44.0
2004-01-02T00:06:00.000 0.487 44.0
2004-01-02T00:09:00.000 0.636 44.0
2004-01-02T00:12:00.000 0.782 44.0
2004-01-02T00:15:00.000 0.922 44.0
2004-01-02T00:18:00.000 1.052 44.0
2004-01-02T00:21:00.000 1.168 44.0
2004-01-02T00:24:00.000 1.261 44.0
2004-01-02T00:27:00.000 1.322 44.0
2004-01-02T00:30:00.000 1.342 44.0
2004-01-02T00:33:00.000 1.317 44.0
2004-01-02T00:36:00.000 1.252 44.0
2004-01-02T00:39:00.000 1.156 44.0
2004-01-02T00:42:00.000 1.039 44.0
2004-01-02T00:45:00.000 0.908 44.0
2004-01-02T00:48:00.000 0.767 44.0
2004-01-02T00:51:00.000 0.621 44.0
2004-01-02T00:54:00.000 0.473 44.0
2004-01-02T00:57:00.000 0.325 44.0
2004-01-02T01:00:00.000 0.178 44.0
2004-01-02T01:03:45.919 : event occurred, entering eclipse => switching to night law
2004-01-02T01:03:00.000 0.035 44.0
2004-01-02T01:06:00.000 -0.100 0.0
2004-01-02T01:09:00.000 -0.222 0.0
2004-01-02T01:12:00.000 -0.326 0.0
2004-01-02T01:15:00.000 -0.402 0.0
2004-01-02T01:18:00.000 -0.441 0.0
2004-01-02T01:21:00.000 -0.435 0.0
2004-01-02T01:24:00.000 -0.386 0.0
2004-01-02T01:27:00.000 -0.302 0.0
2004-01-02T01:30:00.000 -0.193 0.0
2004-01-02T01:34:28.690 : event occurred, exiting eclipse => switching to day law
2004-01-02T01:33:00.000 -0.067 0.0
2004-01-02T01:36:00.000 0.071 44.0
2004-01-02T01:39:00.000 0.215 44.0
2004-01-02T01:42:00.000 0.363 44.0
2004-01-02T01:45:00.000 0.512 44.0
2004-01-02T01:48:00.000 0.660 44.0
2004-01-02T01:51:00.000 0.805 44.0
2004-01-02T01:54:00.000 0.944 44.0
2004-01-02T01:57:00.000 1.072 44.0
2004-01-02T02:00:00.000 1.185 44.0
2004-01-02T02:03:00.000 1.273 44.0
2004-01-02T02:06:00.000 1.328 44.0
2004-01-02T02:09:00.000 1.341 44.0
2004-01-02T02:12:00.000 1.309 44.0
2004-01-02T02:15:00.000 1.239 44.0
2004-01-02T02:18:00.000 1.139 44.0
2004-01-02T02:21:00.000 1.018 44.0
2004-01-02T02:24:00.000 0.885 44.0
2004-01-02T02:27:00.000 0.744 44.0
2004-01-02T02:30:00.000 0.597 44.0
2004-01-02T02:33:00.000 0.449 44.0
2004-01-02T02:36:00.000 0.300 44.0
2004-01-02T02:39:00.000 0.154 44.0
2004-01-02T02:42:16.591 : event occurred, entering eclipse => switching to night law
2004-01-02T02:42:00.000 0.013 44.0
2004-01-02T02:45:00.000 -0.121 0.0
2004-01-02T02:48:00.000 -0.241 0.0
2004-01-02T02:51:00.000 -0.341 0.0
2004-01-02T02:54:00.000 -0.411 0.0
2004-01-02T02:57:00.000 -0.443 0.0
2004-01-02T03:00:00.000 -0.430 0.0
%% Cell type:code id: tags:
``` python
print("Propagation ended at " + finalState.getDate().toString())
```
%% Output
Propagation ended at 2004-01-02T03:00:00.000
%% Cell type:markdown id: tags:
To plot nicely in Python, we should convert the orekit AbsoluteDate to python DateTime objects.
%% Cell type:code id: tags:
``` python
from datetime import datetime
def absolutedate_to_datetime(orekit_date):
utc = TimeScalesFactory.getUTC()
or_comp = orekit_date.getComponents(utc)
or_date = or_comp.getDate()
or_time = or_comp.getTime()
seconds = or_time.getSecond()
return datetime(or_date.getYear(),
or_date.getMonth(),
or_date.getDay(),
or_time.getHour(),
or_time.getMinute(),
int(math.floor(seconds)),
int(1000.0 * (seconds - math.floor(seconds))))
```
%% Cell type:code id: tags:
``` python
pydates = [absolutedate_to_datetime(t) for t in handler.dates]
```
%% Cell type:code id: tags:
``` python
plt.xticks(rotation=90)
ax = plt.gca()
ax2 = ax.twinx()
ax.plot(pydates,handler.eclipseAngles, color='blue')
ax.set_ylabel('Eclipse Angle', color='blue')
ax2.plot(pydates, handler.pointingOffsets, color='green')
ax2.set_ylabel('Pointing offset', color='green');
```
%% Output
%% Cell type:markdown id: tags:
The java code for similar example can be found in the source tree of the library, in file:
src/tutorials/fr/cs/examples/attitude/EarthObservation.java
......
examples/Example numerical prop.ipynb
View file @ c01a858b
%% Cell type:markdown id: tags:
# Numerical Propagation Example
%% Cell type:code id: tags:
``` python
%pylab inline
```
%% Output
Populating the interactive namespace from numpy and matplotlib