Re: [Orekit Users] Numerical Propagator using the Python Orekit wrapper

Thanks so much for the hint Evan! I was able to access the JArray casting using the following:

orekit.JArray('double').cast_(tolerances[0])

As per the tutorial, I was able to successfully add the gravity field to the propagator. But, I am now stuck on getting the atmosphere models to work.

I want to simulate the orbit of a CubeSat under the influence of Drag (Spherical model) and

a 10 x 10 field. The spacecraft mass and model are as follow:

initialDate = AbsoluteDate(2016, 01, 01, 23, 30, 00.000,utc);
mu = 3.986004415e+14;
# Initial orbit
a = 6878.0e3; # semi major axis in meters
e = 0.000002; # eccentricity
i = 0.122173; # inclination
omega = 3.14159; # perigee argument
raan = 4.55531; # right ascention of ascending node
lM = 0.0; # mean anomaly
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu)
# Initial state definition
initialState = SpacecraftState(initialOrbit, 5.0) # 5 kg
satmodel = drag.IsotropicDrag(0.039, 2.2) # Cross sectional area and the drag coefficient

The above describes a satellite in an almost circular orbit at 500 km altitude. The mass is 5 kg and the inverse Ballistic coefficient is 58.3 kg/m^2 (m = 5 kg, CD = 2.2, surface area = 0.039 m^2)

I am trying to set up the DTM 2000 atmosphere model as:

# Initialize the drag
sun = CelestialBodyFactory.getSun() # sun position for drag
earthBody = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, True))
## DTM2000
avg = drag.MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
sflux = drag.MarshallSolarActivityFutureEstimation("Jun2016F10.txt", avg)
atmosphere = drag.DTM2000(sflux, sun, earthBody)
dragForce = drag.DragForce(atmosphere, satmodel)

I am wondering if this is the correct way to use the DTM model. I am doubting my implementation because using the Average strength level results in an altitude of 487.5 km after 30 days of propagation. However, using either a weak or a strong strength level crashes the propagation saying that the altitude went below 120 km. I could maybe understand the High level doing that, but I am surprised by the weak level's behavior.

Also, I can't figure out how to get the JB2006 model to work. Especially, loading the JB2006 input parameters.

I would really appreciate it if someone can show a simple example (Python or Java) where drag (with the DTM and the JB2006 atmospheres) is used on a spherical spacecraft.

Best Regards,

Vivek

On Tue, Jul 26, 2016 at 5:18 AM, Evan Ward <evan.ward@nrl.navy.mil> wrote:

Hi Vivek,

It appears the python wrapper believes tolerances[0] and tolerances[1] are of type Object instead of double[]. The python wrapper has trouble with the Java type system and overloaded methods in some cases. Have a look at [1]. It suggests using JArray('double').cast_(tolerances[0]) to tell python that it is actually a double[].

Best Regards,

Evan

[1] https://lucene.apache.org/pylucene/jcc/features.html#handling-arrays

On 07/25/2016 06:12 PM, Vivek Vittaldev wrote:

Hi everyone,

I am trying to use the Python Orekit wrapper to numerically integrate orbits and I'm running into a few problems.

I am following the tutorial here:

https://www.orekit.org/static/tutorial/propagation.html

I've gotten the Slave mode and the Ephemeris mode part of the tutorial to work. But I get an error when I try to set up the integrator, i.e., in the following part of the tutorial:

// Adaptive step integrator
// with a minimum step of 0.001 and a maximum step of 1000
double minStep = 0.001;
double maxstep = 1000.0;
double positionTolerance = 10.0;
OrbitType propagationType = OrbitType.KEPLERIAN;
double[][] tolerances =
NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
AdaptiveStepsizeIntegrator integrator =
new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);

The error I get is:

InvalidArgsError: (<type 'DormandPrince853Integrator'>, '__init__', (0.001, 1000.0, <Object: [D@e50a6f6>, <Object: [D@14ec4505>))

It looks like it is something to do with the handling of tolerances. I have attached my IPython notebook in case anyone can help me. The code I have for running the example is:

# Import all the required libraries
#Initialize orekit and JVM
import orekit
orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
#Apache math
from org.orekit.bodies import CelestialBodyFactory, CelestialBody
from org.orekit.errors import OrekitException
from org.orekit.frames import LOFType, Frame, FramesFactory, Transform
from org.orekit.orbits import KeplerianOrbit, Orbit, CartesianOrbit, OrbitType, PositionAngle
from org.orekit.propagation import Propagator, SpacecraftState
from org.orekit.propagation.analytical import tle, KeplerianPropagator
from org.orekit.propagation.numerical import NumericalPropagator, Jacobianizer, TimeDerivativesEquations
from org.orekit.time import TimeScalesFactory, AbsoluteDate, DateComponents, TimeComponents
from org.orekit.utils import IERSConventions, Constants, PVCoordinates, PVCoordinatesProvider
from org.orekit.python import PythonEventHandler, PythonOrekitFixedStepHandler
from org.orekit.utils import IERSConventions
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.geometry.euclidean.threed import RotationOrder
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator, ClassicalRungeKuttaIntegrator

# Inertial frame

inertialFrame = FramesFactory.getEME2000();

# Initial date

utc = TimeScalesFactory.getUTC();

initialDate = AbsoluteDate(2004, 01, 01, 23, 30, 00.000,utc);

# Central attraction coefficient

mu = 3.986004415e+14;

# Initial orbit

a = 24396159.0; # semi major axis in meters

e = 0.72831215; # eccentricity

i = 0.122173; # inclination

omega = 3.14159; # perigee argument

raan = 4.55531; # right ascention of ascending node

lM = 0.0; # mean anomaly

initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu)

# Initial state definition

initialState = SpacecraftState(initialOrbit);

# Adaptive step integrator

# with a minimum step of 0.001 and a maximum step of 1000

minStep = 0.001;

maxstep = 1000.0;

positionTolerance = 10.0;

propagationType = OrbitType.KEPLERIAN;

tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);

integrator = DormandPrince853Integrator(minStep,

maxstep,

tolerances[0],

tolerances[1]) # Error here

propagator = NumericalPropagator(integrator)

propagator.setOrbitType(propagationType)

I wonder if I am missing any imports. I would greatly appreciate any help!

Best Regards,

Vivek