[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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



I was able to get the DTM atmosphere to work. The reason for not working earlier is that the position tolerance and the maximum step size were too high. The numerical integrations are working after reducing both of them.

However, I am still stuck on getting the JB2006 atmosphere to work. I can't seem to figure out what how I can read/create the JB2006InputParameters.  I used MarshallSolarActivityFutureEstimation for DTM2000InputParameters, which I hope is correct. How can I load the JB2006InputParameters? I would appreciate any help in this matter.

Best Regards,

Vivek

On Wed, Jul 27, 2016 at 2:50 PM, Vivek Vittaldev <v.vittaldev@utexas.edu> wrote:
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:

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