SRP effect in propagation

[deschku]

Hello,

I’ve just started using the Orekit Python Wrapper to build a propagator for satellite orbits. I’m using STK to compare and match the specifications of the forces in the model, but solar radiation pressure I haven’t been able to tackle. It adds periodic differences with a max of 7500 meters after 24 hours while STK adds about 10-12 meters. Making Cr=0 or reducing the area_srp does nothing to change the effect. Making the mass of the satellite 1kg instead of 1000kg does make a difference. I have concluded that it has something to do with how i set up my force model, the code for which is below. I use RK4 in both propagators.
I’m new to this, thank you in advance.

from org.orekit.time import TimeScalesFactory, AbsoluteDate
from org.orekit.frames import FramesFactory
from org.orekit.bodies import CelestialBodyFactory
from org.orekit.utils import Constants as orekit_constants, IERSConventions

utc = TimeScalesFactory.getUTC() 
epoch = AbsoluteDate(2025, 3, 4, 9, 0, 0.0, utc)  # start time of measurements (taken from file name)
ecef = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
gcrf = FramesFactory.getGCRF()

mu = orekit_constants.EGM96_EARTH_MU

from org.orekit.propagation.conversion import ClassicalRungeKuttaIntegratorBuilder
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation import SpacecraftState
from org.orekit.orbits import OrbitType

sigma_position = 10
# sigma_velocity = 0.1 # m/s

# estimator parameters
estimator_position_scale = 1.0 # m

# orbit propagator parameters
prop_min_step = 0.1 # s
prop_max_step = 60.0 # s
prop_position_error = 10.0 # m

# integrator/propagator builder
integrator_builder = ClassicalRungeKuttaIntegratorBuilder(prop_max_step)
integrator = integrator_builder.buildIntegrator(ic, OrbitType.CARTESIAN)

propagator = NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN) 

satellite_mass = 1000.0 #can be set explicitly
initial_state = SpacecraftState(ic, satellite_mass)
propagator.setInitialState(initial_state)

# setup force model

from org.orekit.forces.gravity import ThirdBodyAttraction, HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.radiation import SolarRadiationPressure, IsotropicRadiationSingleCoefficient
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.utils import Constants, IERSConventions

# Earth geopotential model
earth_frame = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth_shape = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                               Constants.WGS84_EARTH_FLATTENING,
                               earth_frame)

# Gravity field
gravity_provider = GravityFieldFactory.getNormalizedProvider(21, 21)
propagator.addForceModel(
    HolmesFeatherstoneAttractionModel(earth_frame, gravity_provider)
)

# Third bodies
sun = CelestialBodyFactory.getSun()
moon = CelestialBodyFactory.getMoon()
propagator.addForceModel(ThirdBodyAttraction(sun))
propagator.addForceModel(ThirdBodyAttraction(moon))
                         
# Drag Force
from org.orekit.models.earth.atmosphere import NRLMSISE00, PythonNRLMSISE00InputParameters
from org.orekit.forces.drag import DragForce, IsotropicDrag

class ConstantSolarActivity(PythonNRLMSISE00InputParameters):
    def __init__(self, f10, f10_avg, ap):
        super(ConstantSolarActivity, self).__init__()
        self.f10 = float(f10)
        self.f10_avg = float(f10_avg)
        self.ap = float(ap)
        
    def getMinDate(self):
        return AbsoluteDate.PAST_INFINITY
        
    def getMaxDate(self):
        return AbsoluteDate.FUTURE_INFINITY
        
    def getDailyFlux(self, date):
        return self.f10
        
    def getAverageFlux(self, date):
        return self.f10_avg
        
    def getAp(self, date):

        return [self.ap] * 7

f10 = 150.0
f10_avg = 150.0
ap = 15.0

solar_params = ConstantSolarActivity(f10, f10_avg, ap)
atm = NRLMSISE00(solar_params, sun, earth_shape)

cross_section = 8.5194   # m^2 (given that orekit sets a default mass of 1000kg, this corresponds to a 0.0085194 m^2/kg area/mass ratio))
drag_coef = 2.2

iso_drag = IsotropicDrag(cross_section, drag_coef)
drag_force = DragForce(atm, iso_drag)
propagator.addForceModel(drag_force)

# Solar Radiation Pressure
cr = 1.2
area_srp = 8.5194 # m^2 (corresponds to 0.0085194 m^2/kg area/mass ratio with default mass)
srp_model = IsotropicRadiationSingleCoefficient(area_srp, cr)
srp = SolarRadiationPressure(sun, earth_shape, srp_model) # (uses the Earth as occulting body and a dual cone shadow model)

propagator.addForceModel(srp)

[MaximeJ]

Hi @deschku,

Welcome to the forum !!

@deschku

but solar radiation pressure I haven’t been able to tackle. It adds periodic differences with a max of 7500 meters after 24 hours while STK adds about 10-12 meters. Making Cr=0 or reducing the area_srp does nothing to change the effect. Making the mass of the satellite 1kg instead of 1000kg does make a difference.

What exactly are you trying to achieve ? Matching Orekit’s results with STK’s ?
SRP acceleration depends on the ratio Cr * Surface / Mass (see here) so it’s weird that setting Cr=0 or reducing srp_area doesn’t change anything.
Could you explain a bit more what you don’t understand in the usage of SRP ?

Cheers,
Maxim

[deschku]

Hello,

@MaximeJ

What exactly are you trying to achieve ? Matching Orekit’s results with STK’s ?

Yes, I’m trying to achieve the highest match-up between STK’s and my propagator’s results possible. I understand it has been done before but when I check similar code there doesn’t seem to be a difference when it comes to the SRP set-up.

@MaximeJ

Could you explain a bit more what you don’t understand in the usage of SRP ?

While I understand how the SRP force works, I want to see whether the way I set up the force model introduces these extra many meters to the propagation. For example, in the plots below I performed two propagations, one with a simple 21x21 harmonics model, and the other where I turned on SRP with Cr=0. Simply adding the force creates a huge difference. I thought maybe it was a frame mismatch but switching from ITRF to GCRF changes nothing. Since I’m new to Orekit maybe there is something i set up wrong and causes this difference.

[MaximeJ]

Thank you for the details.

@deschku

For example, in the plots below I performed two propagations, one with a simple 21x21 harmonics model, and the other where I turned on SRP with Cr=0

ok, that should not happen. Setting Cr=0 should cancel the SRP. There is something weird here. I didn’t see anything obviously wrong in your code so I’m a bit at a loss…
I should try on my side and see if I get any difference.

[deschku]

Hello Maxim,

as it turns out a more robust integrator did the trick. I used DorPri with low position error instead of RK4, subtracted positions with SRP off from SRP on and got the desired 10-12 meter difference. I’m guessing that such a small difference when positions are in the millions of meters is a bit of a hard task for a simple Runge-Kutta integrator?
The reason I thought it didn’t have to do with the integrator is because I used the STK propagation for comparison and while it calculated the SRP correctly it completely threw off the harmonics model (it uses RK4).

Thank you for your help and interest and sorry for the mix-up!

[MaximeJ]

Hello @deschku,

@deschku

Thank you for your help and interest and sorry for the mix-up!

You’re welcome and no worries! I’m actually happy things are working on your side.