Custom Atmosphere Implementation

mia · November 14, 2024, 4:32pm

Hi! I’m using the orekit python wrapper to create a custom atmosphere field. I’ve followed this post to create a CustomAtmosphere2 class which implements the PythonAtmosphere class, and I filled in getDensity and getVelocity. To verify my implementation, I tried to replicate the JB2008 model. To do that, I called the JB2008 getDensity() method in my custom atmosphere’s getDensity() method. I have found that the propagated orbit using my custom atmosphere and the JB2008 model produce much different orbits. I suspect that there might be a frame transformation issue. What am I missing here?

Apologies for the large amount of code. Please direct your attention to the CustomAtmosphere2 class. Thanks!!

# Initialize orekit and JVM
import orekit
vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()

from org.orekit.orbits import KeplerianOrbit, EquinoctialOrbit, PositionAngle, OrbitType, CartesianOrbit
from org.orekit.frames import FramesFactory, LOFType, Frame
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, PVCoordinates
from org.orekit.time import AbsoluteDate
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, CelestialBodyFactory, CelestialBody
from org.orekit.utils import IERSConventions
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel, ThirdBodyAttraction, OceanTides, SolidTides
from orekit import JArray_double, JArray
from java.util import ArrayList

from org.orekit.forces.radiation import SolarRadiationPressure, IsotropicRadiationSingleCoefficient, RadiationSensitive
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData, JB2008SpaceEnvironmentData
from org.orekit.forces.drag import IsotropicDrag, DragForce
from org.orekit.models.earth.atmosphere import DTM2000, HarrisPriester, JB2008, NRLMSISE00, Atmosphere, SimpleExponentialAtmosphere
# from org.orekit.data import DataProvidersManager

# from org.orekit.attitudes import LofOffset, InertialProvider, AttitudeProvider
from org.orekit.propagation.events import DateDetector
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.propagation.semianalytical.dsst import DSSTPropagator
from org.orekit.propagation import PropagationType
from org.orekit.propagation.semianalytical.dsst.forces import DSSTThirdBody
from org.orekit.propagation.semianalytical.dsst.forces import DSSTSolarRadiationPressure
from org.orekit.models.earth import ReferenceEllipsoid
from org.orekit.propagation.semianalytical.dsst.forces import DSSTZonal   
from org.orekit.propagation.semianalytical.dsst.forces import DSSTTesseral
from org.orekit.propagation.semianalytical.dsst.forces import DSSTAtmosphericDrag
from org.orekit.propagation.conversion import NumericalPropagatorBuilder, DormandPrince853IntegratorBuilder, DormandPrince54IntegratorBuilder, OsculatingToMeanElementsConverter, ClassicalRungeKuttaIntegratorBuilder
from org.hipparchus.linear import RealMatrix, DiagonalMatrix
from org.orekit.forces import ForceModel

from math import radians, degrees, pi
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
import pandas as pd
import pickle
import time
import matplotlib.ticker as ticker
from scipy.io import loadmat

plt.rcParams.update({'font.size': 12})

from org.orekit.models.earth.atmosphere import PythonAtmosphere
from org.orekit.models.earth.atmosphere import JB2008
from org.orekit.time import AbsoluteDate
from org.orekit.frames import Frame
from org.orekit.utils import PVCoordinates
from org.hipparchus.geometry.euclidean.threed import Vector3D

class CustomAtmosphere2(PythonAtmosphere):
    def __init__(self, cswl, sun, earth):
        super().__init__()
        self.atm = JB2008(cswl, sun, earth)
        self.earth = earth
        self.earth_rotation_rate = 7.2921159e-5  # Earth's rotation rate in rad/s

    def getDensity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        transform = frame.getTransformTo(self.earth.getBodyFrame(), date)  # Transformation to Earth frame
        pv_coordinates = transform.transformPVCoordinates(PVCoordinates(position, Vector3D.ZERO))
        position_in_eci = pv_coordinates.getPosition()
        return self.atm.getDensity(date, position_in_eci, self.earth.getBodyFrame())
    
    def getVelocity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        try:
            if frame == self.earth.getBodyFrame():
                # If in Earth-fixed frame, return zero velocity
                velocity = Vector3D.ZERO
            else:
                # Transform position and velocity to the Earth-centered inertial frame
                transform = frame.getTransformTo(self.earth.getBodyFrame(), date)
                pv_coordinates = transform.transformPVCoordinates(PVCoordinates(position, Vector3D.ZERO))
                position_in_eci = pv_coordinates.getPosition()
                
                # Cross product omega x r to compute velocity
                velocity = Vector3D.crossProduct(Vector3D(0.0, 0.0, self.earth_rotation_rate), position_in_eci)
                    
            return velocity
        except Exception as e:
            print(f"Error in getVelocity: {e}")
            raise RuntimeError(f"Error in getVelocity: {e}")
    

def prop_orbit(mass,BC,srpArea,pos,vel,tss,sel_d, inertialFrame,mu,earth,gravity_force,drag_force,sun_force,threebody_force,tides_force):

    ts = tss

    # Orbital perturbations
    # Earth non spehrical gravity field
    degree = 70 # 70
    torder = 70 
    # Drag Force
    dragCoeff = 2.2
    crossSection = float(( mass*BC ) / dragCoeff) # m^2
    # Solar radiation pressure
    cr = 1.0
    
    
    # =============================================================================
    #% Orbit setup
    # =============================================================================

    # Reference Orbit (km to m)
    POS = Vector3D([float(1e3*pos[0]), float(1e3*pos[1]), float(1e3*pos[2])])
    VEL = Vector3D([float(1e3*vel[0]), float(1e3*vel[1]), float(1e3*vel[2])])

    initialOrbit = CartesianOrbit(PVCoordinates(POS, VEL),
                    inertialFrame, # The frame in which the parameters are defined (must be a pseudo-inertial frame)
                    ts[0],   # Sets the date of the orbital parameters
                    mu)   # Sets the central attraction coefficient (m³/s²)

    # =============================================================================
    #% Propagator setup
    # =============================================================================

    orbitType = initialOrbit.getType()
    initialState = SpacecraftState(initialOrbit, mass)
    minStep = 1.0e-10
    maxstep = 1.0e4
    initStep = 0.1
    positionTolerance = 1e-4
    tol = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)
    integrator = DormandPrince853Integrator(minStep, maxstep, 
    JArray_double.cast_(tol[0]),  # Double array of doubles needs to be casted in Python
    JArray_double.cast_(tol[1]))
    integrator.setInitialStepSize(initStep)
    propagator_num = NumericalPropagator(integrator)
    propagator_num.setOrbitType(orbitType)
    propagator_num.setInitialState(initialState)

    # freezingDate = AbsoluteDate.J2000_EPOCH;
    gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(degree, torder)
    gravityForce = HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider)
    if gravity_force == 1:
        propagator_num.addForceModel(gravityForce)

    sun = CelestialBodyFactory.getSun()
    moon = CelestialBodyFactory.getMoon()

    # Solar radiation pressure
    spacecraft = IsotropicRadiationSingleCoefficient(srpArea, cr)
#     print(isinstance(spacecraft, RadiationSensitive))
    srpProvider = SolarRadiationPressure(sun, earth.getEquatorialRadius(), spacecraft)
    if sun_force == 1:
        propagator_num.addForceModel(srpProvider)
        
    #Atmosphere model
    satmodel = IsotropicDrag(crossSection, dragCoeff) # Cross sectional area and the drag coefficient

    # Drag Force
    
    if sel_d == 1:
        cswl = CssiSpaceWeatherData("SpaceWeather-All-v1.2_ORIG.txt") # cswl.getDailyFlux(date)
        atmosphere = DTM2000(cswl, sun, earth)  
    elif sel_d == 2:
        cswl = JB2008SpaceEnvironmentData("SOLFSMY.TXT","DTCFILE.TXT") 
        atmosphere = JB2008(cswl, sun, earth)
    elif sel_d == 3:
        cswl = JB2008SpaceEnvironmentData("SOLFSMY.TXT","DTCFILE.TXT") 
        atmosphere = CustomAtmosphere(cswl, sun, earth, initialState, inertialFrame)
    
    dragForce = DragForce(atmosphere, satmodel)
    if drag_force == 1:
        propagator_num.addForceModel(dragForce)   

    # Solid Tides
    solidTidesBodies = ArrayList().of_(CelestialBody)
    solidTidesBodies.add(sun)
    solidTidesBodies.add(moon)
    solidTidesBodies = solidTidesBodies.toArray()

    solidTides = SolidTides(earth.getBodyFrame(), 
                            gravityProvider.getAe(), gravityProvider.getMu(),
                            gravityProvider.getTideSystem(), 
                            IERSConventions.IERS_2010,
                            TimeScalesFactory.getUT1(IERSConventions.IERS_2010, True), 
                            solidTidesBodies)
    if tides_force == 1:
        propagator_num.addForceModel(solidTides)

    ThirdBodyAttraction
    if threebody_force == 1:
        propagator_num.addForceModel(ThirdBodyAttraction(sun))
        propagator_num.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getMoon())) 

    # =============================================================================
    #% Orbit propagation
    # =============================================================================
#     print("Forces:", propagator_num.getAllForceModels())
    
    states = [propagator_num.propagate(tt)  for tt in ts]
    
    return states

# =============================================================================
#% Constants/Basics
# =============================================================================
r_Earth = Constants.IERS2010_EARTH_EQUATORIAL_RADIUS #m
itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed
inertialFrame = FramesFactory.getEME2000()
earth = OneAxisEllipsoid(r_Earth,
                         Constants.IERS2010_EARTH_FLATTENING,
                         itrf)
mu = Constants.IERS2010_EARTH_MU #m^3/s^2

utc = TimeScalesFactory.getUTC()

# Selection

sel_orig = 1
sel_custom = 1
#--------------------------------------------------------------------------------------
# Reference Orbit

rp0 = r_Earth + 320 * 1e3 # perigee radius (m)
ra0 = r_Earth + 1505 * 1e3 # apogee radius (m)
secs = 60.0
mins = 10
days = 10
step = secs*mins # propagation step time [s]
dur = days*86400.0 # propagation duration time [s]
satellite_mass = 260.0
crossSection = 3.2*1.6 # m^2
srpArea = 30.0 # m^2
#--------------------------------------------------------------------------------------

deg = pi/180
a0 = (rp0 + ra0)/2
e0 = (ra0 - rp0)/(ra0 + rp0)
w0 = 30 * deg
i0 = 45 * deg
ra0 = 0 * deg
M0 = 0 * deg #true anomaly

degree = 70 
torder = 70 
dragCoeff = 2.2
cr = 1.0

epochDate = AbsoluteDate(2022, 1, 1, 0, 0, 00.000, utc)
initialDate = epochDate
start_time = initialDate
tspan1 = [start_time.shiftedBy(float(dt)) for dt in np.linspace(0, dur, int(dur/step))]
tspan2 = [tspan1[0],tspan1[-1]]

minStep = 1e-6
maxstep = 100.0
initStep = 1.0
positionTolerance = 1e-4
    
cswl = JB2008SpaceEnvironmentData("SOLFSMY.TXT","DTCFILE.TXT") 
sun = CelestialBodyFactory.getSun()


satmodel = IsotropicDrag(crossSection, dragCoeff) # Cross sectional area and the drag coefficient

# Propagation setup for original JB2008
if sel_orig == 1:
    initialOrbit = KeplerianOrbit(a0, e0, i0, w0, ra0, M0, PositionAngle.TRUE, inertialFrame, initialDate, mu)
    orbitType = initialOrbit.getType()
    initialState = SpacecraftState(initialOrbit, satellite_mass)
    tol = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)

    integrator = DormandPrince853Integrator(minStep, maxstep, JArray_double.cast_(tol[0]), JArray_double.cast_(tol[1]))
    integrator.setInitialStepSize(initStep)

    propagator_num = NumericalPropagator(integrator)
    propagator_num.setOrbitType(orbitType)
    propagator_num.setInitialState(initialState)

    gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(degree, torder)
    gravityForce = HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider)
    propagator_num.addForceModel(gravityForce)

    print("correct-JB2008")
    atmosphere = JB2008(cswl, sun, earth)
    dragForce = DragForce(atmosphere, satmodel)
    propagator_num.addForceModel(dragForce)

    states = [initialState]
    tic = time.time()
    for i1 in range(len(tspan1) - 1):
        states.append(propagator_num.propagate(tspan1[i1], tspan1[i1 + 1]))
    toc = time.time()

    posvel = [state.getPVCoordinates() for state in states]
    poss = [state.getPosition() for state in posvel]
    vels = [state.getVelocity() for state in posvel]
    px = [pos.getX() * 1e-3 for pos in poss]
    py = [pos.getY() * 1e-3 for pos in poss]
    pz = [pos.getZ() * 1e-3 for pos in poss]
    vx = [vel.getX() * 1e-3 for vel in vels]
    vy = [vel.getY() * 1e-3 for vel in vels]
    vz = [vel.getZ() * 1e-3 for vel in vels]
    stat_list = [dur, toc - tic, px[-1], py[-1], pz[-1], vx[-1], vy[-1], vz[-1], step]
    print("Time interval [s]:", stat_list[0])
    print("Time step [s]:", stat_list[8])
    print("CPU time [s]:", stat_list[1])
    print("Final Pos [km]:", np.linalg.norm([px[-1], py[-1], pz[-1]]))
    print("Final Vel [km]:", np.linalg.norm([vx[-1], vy[-1], vz[-1]]))
    print("\n")
    poss_orig = poss
    vels_orig = vels

# Propagation setup for custom JB2008
if sel_custom == 1:
    initialOrbit = KeplerianOrbit(a0, e0, i0, w0, ra0, M0, PositionAngle.TRUE, inertialFrame, initialDate, mu)
    orbitType = initialOrbit.getType()
    initialState = SpacecraftState(initialOrbit, satellite_mass)
    tol = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)

    integrator = DormandPrince853Integrator(minStep, maxstep, JArray_double.cast_(tol[0]), JArray_double.cast_(tol[1]))
    integrator.setInitialStepSize(initStep)

    propagator_num = NumericalPropagator(integrator)
    propagator_num.setOrbitType(orbitType)
    propagator_num.setInitialState(initialState)

    gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(degree, torder)
    gravityForce = HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider)
    propagator_num.addForceModel(gravityForce)

    print("custom-JB2008")
    atmosphere = CustomAtmosphere2(cswl, sun, earth)
    dragForce = DragForce(atmosphere, satmodel)
    propagator_num.addForceModel(dragForce)

    states = [initialState]
    tic = time.time()
    for i1 in range(len(tspan1) - 1):
        states.append(propagator_num.propagate(tspan1[i1], tspan1[i1 + 1]))
    toc = time.time()

    posvel = [state.getPVCoordinates() for state in states]
    poss = [state.getPosition() for state in posvel]
    vels = [state.getVelocity() for state in posvel]
    px = [pos.getX() * 1e-3 for pos in poss]
    py = [pos.getY() * 1e-3 for pos in poss]
    pz = [pos.getZ() * 1e-3 for pos in poss]
    vx = [vel.getX() * 1e-3 for vel in vels]
    vy = [vel.getY() * 1e-3 for vel in vels]
    vz = [vel.getZ() * 1e-3 for vel in vels]
    stat_list = [dur, toc - tic, px[-1], py[-1], pz[-1], vx[-1], vy[-1], vz[-1], step]
    print("Time interval [s]:", stat_list[0])
    print("Time step [s]:", stat_list[8])
    print("CPU time [s]:", stat_list[1])
    print("Final Pos [km]:", np.linalg.norm([px[-1], py[-1], pz[-1]]))
    print("Final Vel [km]:", np.linalg.norm([vx[-1], vy[-1], vz[-1]]))
    print("\n")
    poss_custom = poss
    vels_custom = vels
#--------------------------------------------------------------------------------------
colors = [[0, 0, 0], [0, 0.4470, 0.7410], [0.8500, 0.3250, 0.0980], 
          [0.9290, 0.6940, 0.1250], [0.4940, 0.1840, 0.5560], [0.4660, 0.6740, 0.1880]]
sel_LineWidth = 3
sel_LineWidth2 = 1.5
sel_MarkerWidth = 10
sel_LineWidthAxis = 1
sel_FontSize = 14  

step_p = 1
tspan3 = np.linspace(0, dur, int(dur/step))/60

px_o = np.array([pos.getX()*1e-3 for pos in poss_orig])
py_o = np.array([pos.getY()*1e-3 for pos in poss_orig])
pz_o = np.array([pos.getZ()*1e-3 for pos in poss_orig])
vx_o = np.array([vel.getX()*1e-3 for vel in vels_orig])
vy_o = np.array([vel.getY()*1e-3 for vel in vels_orig])
vz_o = np.array([vel.getZ()*1e-3 for vel in vels_orig])

px_c = np.array([pos.getX()*1e-3 for pos in poss_custom])
py_c = np.array([pos.getY()*1e-3 for pos in poss_custom])
pz_c = np.array([pos.getZ()*1e-3 for pos in poss_custom])
vx_c = np.array([vel.getX()*1e-3 for vel in vels_custom])
vy_c = np.array([vel.getY()*1e-3 for vel in vels_custom])
vz_c = np.array([vel.getZ()*1e-3 for vel in vels_custom])

err_x = abs(px_o - px_c)
err_y = abs(py_o - py_c)
err_z = abs(pz_o - pz_c)
err_vx = abs(vx_o - vx_c)
err_vy = abs(vy_o - vy_c)
err_vz = abs(vz_o - vz_c)

err_pos = []
err_vel = []
for i1 in range(len(tspan3)):
    err_pos.append( abs( np.linalg.norm( np.array([px_o[i1],py_o[i1],pz_o[i1]]) - np.array([px_c[i1],py_c[i1],pz_c[i1]]) ) ) )
    err_vel.append( abs( np.linalg.norm( np.array([vx_o[i1],vy_o[i1],vz_o[i1]]) - np.array([vx_c[i1],vy_c[i1],vz_c[i1]]) ) ) )
#     err_pos.append( abs( np.linalg.norm([px_o[i1],py_o[i1],pz_o[i1]]) - np.linalg.norm([px_c[i1],py_c[i1],pz_c[i1]]) ) )
#     err_vel.append( abs( np.linalg.norm([vx_o[i1],vy_o[i1],vz_o[i1]]) - np.linalg.norm([vx_c[i1],vy_c[i1],vz_c[i1]]) ) )
err_pos = np.array(err_pos)
err_vel = np.array(err_vel)

plt.figure(facecolor='white') # ,figsize=(10, 6)
plt.grid(True)
plt.plot(tspan3[::step_p], err_pos[::step_p], linewidth=sel_LineWidth, color=colors[1], linestyle='-')
plt.xlabel('Time (min)', fontsize=sel_FontSize)
plt.ylabel('Error - Position norm (km)', fontsize=sel_FontSize)
plt.xticks(fontsize=sel_FontSize)
plt.yticks(fontsize=sel_FontSize)
plt.tight_layout()
plt.gca().tick_params(width=sel_LineWidthAxis)
plt.show()

plt.figure(facecolor='white')
plt.grid(True)
plt.plot(tspan3[::step_p], err_vel[::step_p], linewidth=sel_LineWidth, color=colors[1], linestyle='-')
plt.xlabel('Time (min)', fontsize=sel_FontSize)
plt.ylabel('Error - Velocity norm (km/s)', fontsize=sel_FontSize)
plt.xticks(fontsize=sel_FontSize)
plt.yticks(fontsize=sel_FontSize)
plt.tight_layout()
plt.gca().tick_params(width=sel_LineWidthAxis)
plt.show()

plt.figure(facecolor='white')
plt.grid(True)
plt.plot(tspan3[::step_p], err_x[::step_p], linewidth=sel_LineWidth, label="x", color=colors[1], linestyle='-')
plt.plot(tspan3[::step_p], err_y[::step_p], linewidth=sel_LineWidth, label="y", color=colors[2], linestyle='-')
plt.plot(tspan3[::step_p], err_z[::step_p], linewidth=sel_LineWidth, label="z", color=colors[3], linestyle='-')
plt.xlabel('Time (min)', fontsize=sel_FontSize)
plt.ylabel('Error - Position (km)', fontsize=sel_FontSize)
plt.xticks(fontsize=sel_FontSize)
plt.yticks(fontsize=sel_FontSize)
plt.tight_layout()
plt.gca().tick_params(width=sel_LineWidthAxis)
plt.legend(loc='best')
plt.show()

plt.figure(facecolor='white')
plt.grid(True)
plt.plot(tspan3[::step_p], err_vx[::step_p], linewidth=sel_LineWidth, label="vx", color=colors[1], linestyle='-')
plt.plot(tspan3[::step_p], err_vy[::step_p], linewidth=sel_LineWidth, label="vy", color=colors[2], linestyle='-')
plt.plot(tspan3[::step_p], err_vz[::step_p], linewidth=sel_LineWidth, label="vz", color=colors[3], linestyle='-')
plt.xlabel('Time (min)', fontsize=sel_FontSize)
plt.ylabel('Error - Velocity (km/s)', fontsize=sel_FontSize)
plt.xticks(fontsize=sel_FontSize)
plt.yticks(fontsize=sel_FontSize)
plt.tight_layout()
plt.gca().tick_params(width=sel_LineWidthAxis)
plt.legend(loc='best')
plt.show()

baubin · November 14, 2024, 6:37pm

Apologies if I’m misunderstanding, but are you saying that when you call JB2008 from your custom class, you get a different result than when you call JB2008 directly?

mia · November 14, 2024, 7:18pm

Yes. I am giving the same inputs, but somehow my custom class (which is just calling JB2008) is not resulting in the same orbit as JB2008.

baubin · November 14, 2024, 7:22pm

I would check a couple of things straight off:

Do you get the same density from your custom class vs. when you use the orekit model directly?
Are you entering the same position value in ECI as well as the same datetime value into both models?

mia · November 15, 2024, 5:57pm

Hi, I think we figured it out. The big change is that we removed the frame transformation in getDensity and implemented getVelocity based on the atmosphere interface Here’s my custom atmosphere class:

class CustomAtmosphere2(PythonAtmosphere):
    def __init__(self, cswl, sun, earth):
        super().__init__()
        self.atm = JB2008(cswl, sun, earth)
        self.earth = earth
    def getDensity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        return self.atm.getDensity(date, position, frame)
    def getVelocity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        bodyToFrame = frame.getTransformTo(self.earth.getBodyFrame(), date)
        bodyToFrame_inv = bodyToFrame.getInverse()
        pv_coordinates = bodyToFrame.transformPVCoordinates(PVCoordinates(position, Vector3D.ZERO))
        posInBody = pv_coordinates.getPosition()
        pvBody = PVCoordinates(posInBody, Vector3D.ZERO)
        pvFrame = bodyToFrame_inv.transformPVCoordinates(pvBody)
        return pvFrame.getVelocity()

There is still a very slight error compared to the JB2008 class:

correct-JB2008
Time interval [s]: 1728000.0
Time step [s]: 600.0
CPU time [s]: 55.52110171318054
Final Pos [km]: 6890.477134085638
Final Vel [km]: 7.597406306978231

custom-JB2008
Time interval [s]: 1728000.0
Time step [s]: 600.0
CPU time [s]: 67.58265495300293
Final Pos [km]: 6890.477133735982
Final Vel [km]: 7.597406307300406

The position and velocity errors are quite small, but what’s confusing to me is the difference in CPU time. Does anyone have thoughts here?

Serrof · November 15, 2024, 6:20pm

Hi there,

I bet the performance penalty comes from the frame transform. Use getKinematicTransform and the transformPVOnly method. If you don’t, the code will compute everything it needs to convert acceleration vectors as well, which is useless in your case.

Cheers,
Romain.

GiovanniLavezzi · November 16, 2024, 11:41am

Hi, I’m working with Mia on the custom atmosphere implementation in Orekit. At the moment we are using Orekit v.11.3.2, whereas (if understood correctly) getKinematicTransform is available starting from v.12. That’s why I was not able to replicate the same implementation in atmosphere interface. However the original JB2008 implementation is getting a lower cpu time.
My questions are:

Is there a different way to get the same cpu time as the original JB2008 in v.11.3.2? And in general, is that the correct/best implementation?
How much would be the improvement in cpu time if using v.12?

Thanks,
Giovanni

Serrof · November 16, 2024, 2:44pm

Hi,

If you’re using the JCC Python wrapper, you’ll get an overhead from custom implementations.
About Orekit Java, there has been performance improvements across the board on propagation since 11.0, but the effects really depend on the applications (perturbations, attitude, etc)

Cheers,
Romain.

GiovanniLavezzi · November 18, 2024, 8:49am

Hi, yes, I’m using Orekit v.11.3.2 in python. The goal is to have a high-fidelity propagator (inclusive of Earth’s non-spherical gravity field, atmospheric drag, solar radiation pressure, lunisolar third body accelerations, and Earth solid tides) to test custom atmosphere models.
Based on your answer, I didn’t get it if there is a more efficient/standard way to code the custom class, specifically the getVelocity part, since I based my implementation on the link above which is Orekit v.12:

class CustomAtmosphere2(PythonAtmosphere):
    def __init__(self, cswl, sun, earth):
        super().__init__()
        self.atm = JB2008(cswl, sun, earth)
        self.earth = earth
    def getDensity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        return self.atm.getDensity(date, position, frame)
    def getVelocity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        bodyToFrame = frame.getTransformTo(self.earth.getBodyFrame(), date)
        bodyToFrame_inv = bodyToFrame.getInverse()
        pv_coordinates = bodyToFrame.transformPVCoordinates(PVCoordinates(position, Vector3D.ZERO))
        posInBody = pv_coordinates.getPosition()
        pvBody = PVCoordinates(posInBody, Vector3D.ZERO)
        pvFrame = bodyToFrame_inv.transformPVCoordinates(pvBody)
        return pvFrame.getVelocity()

Thank you,
Giovanni

Serrof · November 18, 2024, 4:17pm

What I meant about the JCC-based wrapper is that when you implement methods in python, the back and forth calls with java slow down the run. To avoid that you can try the jpype technology. @yzokras and @petrus.hyvonen might be able to give you some pointers.

Cheers,
Romain

yzokras · November 25, 2024, 11:21am

Unfortunately the new jpype-based Python wrapper is not faster than the JCC one, the performance will still suffer badly from back-and-forth calls between Python and Java.

But fortunately, the Jpype-based Python wrapper makes it much easier to wrap your custom Java code into Python. So I highly recommend you write your custom atmosphere model in Java.

You can use the Ephemerista project as an example which is using Orekit Jpype under the hood for propagation:

Example of custom java code
a Maven pom.xml file for building your custom Java code into a .jar file
Use mvn package to build your JAR file
Then you need to add the path of your JAR file to the java classpath

PS: the Orekit jpype wrapper is only available for Orekit 12, but if you are stuck with a large legacy codebase based on Orekit 11, we could consider releasing a version based on Orekit 11.3 on PyPi.

mia · December 11, 2024, 4:38pm

Thank you!

We did end up just using the python wrapper. For anyone curious, here is the implementation that worked for us in version 12.1:

class CustomAtmosphere(PythonAtmosphere):
    
    def __init__(self, cswl, sun, earth):
        super().__init__()
        self.atm = JB2008(cswl, sun, earth)
        self.earth = earth
    def getDensity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        return self.atm.getDensity(date, position, frame)
    def getVelocity(self, date: AbsoluteDate, position: Vector3D, frame: Frame):
        bodyToFrame = self.earth.getBodyFrame().getKinematicTransformTo(frame, date)
        posInBody = bodyToFrame.getStaticInverse().transformPosition(position)
        pv_body = PVCoordinates(posInBody, Vector3D.ZERO)
        pvFrame = bodyToFrame.transformOnlyPV(pv_body)
        return pvFrame.getVelocity()

Vincent · December 11, 2024, 4:50pm

Thank you for providing a working code sample !

Cheers,
Vincent