Adding PV measurement to estimator

Orekit Python Wrapper usage
1

Hi! sorry to bother you, maybe its a very basic question but Im a new user.

Im trying to add a pvt measurement to an estimator using estimator.addMeasurement(PV(pvt)) and the PV(pvt) part is giving me invalidargserror when applying it to the estimator. The pvt args are the following:
(<class ‘org.orekit.estimation.measurements.PV’>, ‘init’, ([<AbsoluteDate: 2022-04-16T11:52:00.000Z>, JArray[-2768386.5930000003, 5275182.704, 3672200.745], JArray[-759.6318531000001, 4028.2763631999997, -6338.7818408], JArray[1.0, 0.0, 0.0], JArray[1.0, 0.0, 0.0], 1.0],))
could you please indicate me if there is a problem with the way Im providing the measurements or the values? I couldnt find an example for this kind of implementation.
Thanks in advance. Powski

2

Hi @powski

First, welcome to the Orekit forum :slight_smile:

Using estimator.addMeasurement(...) is the good way.
You probably have an error when initializing the PV measurement. I see two possible errors. First, I don’t see the ObservableSatellite in the error message. You maybe forget it. Secondly, I think you use arrays instead of Vector3D to provide the position and velocity data to the PV object.

For the first problem, the ObservableSatellite object is used to represent the satellite for which the orbit is estimated and from which the measurement comes from. If you have only one satellite, you will only have one instance of ObservableSatellite. You can initialize it only once like that:

observableSatellite = ObservableSatellite(0)

Now, I show you an example on how initializing a PV object based on your data.

date = AbsoluteDate(2022, 4, 16, 11, 52, 00.000, TimeScalesFactory.getUTC())
position = Vector3D(-2768386.5930000003, 5275182.704, 3672200.745)
velocity = Vector3D(-759.6318531000001, 4028.2763631999997, -6338.7818408)
pv = PV(date, position, velocity, 1.0, 0.01, 1.0, observableSatellite)

A last remark, I saw you used a sigma of 1.0 meter for position data and 1.0 meter per second for velocity data. I recommend you to reduce the sigma for velocity data compared to the one for position data. Having 2-3 orders of magnitude less for velocity is more realistic. Therefore a value of 1 cm per second (i.e., 0.01) is more relevant.

If you want an example of a full Python-based orbit determination using Orekit, you can have a look on Clément Jonglez’s tutorial on his GitHub repository. He uses PV data stored in a CSV file to perform the orbit determination.

Best regards,
Bryan

1 Like
3

Thank you very much Bryan for the detailed answer, it has all that I needed and more!.
Regarding the observableSatellite param, I included it but later removed it as part of the error testing, thanks also for the recommendation on the sigma for velocity, I will do as indicated.
Last, the suggested tutorial is most welcomed!
Kind regards!

4

dear Bryan, after implementing what you suggested I could call the estimator but its giving me a good number of errors in the java code:

estimatedPropagatorArray = estimator.estimate()

JavaError: <super: <class ‘JavaError’>, >
Java stacktrace:
java.lang.NullPointerException
at org.orekit.time.AbsoluteDate.isEqualTo(AbsoluteDate.java:1210)
at org.orekit.time.AbsoluteDate.isBetweenOrEqualTo(AbsoluteDate.java:1299)
at org.orekit.models.earth.atmosphere.NRLMSISE00.getDensity(NRLMSISE00.java:1156)
at org.orekit.forces.drag.DragForce.acceleration(DragForce.java:80)
at org.orekit.forces.ForceModel.addContribution(ForceModel.java:108)
at org.orekit.propagation.numerical.NumericalPropagator$Main.computeDerivatives(NumericalPropagator.java:899)
at org.orekit.propagation.integration.AbstractIntegratedPropagator$ConvertedMainStateEquations.computeDerivatives(AbstractIntegratedPropagator.java:755)
at org.hipparchus.ode.ExpandableODE.computeDerivatives(ExpandableODE.java:134)
at org.hipparchus.ode.AbstractIntegrator.computeDerivatives(AbstractIntegrator.java:265)
at org.hipparchus.ode.AbstractIntegrator.initIntegration(AbstractIntegrator.java:217)
at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:196)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:474)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:405)
at org.orekit.propagation.PropagatorsParallelizer.propagate(PropagatorsParallelizer.java:140)
at org.orekit.estimation.leastsquares.AbstractBatchLSModel.value(AbstractBatchLSModel.java:311)
at org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresFactory$LocalLeastSquaresProblem.evaluate(LeastSquaresFactory.java:440)
at org.orekit.estimation.leastsquares.BatchLSEstimator$TappedLSProblem.evaluate(BatchLSEstimator.java:615)
at org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer.optimize(GaussNewtonOptimizer.java:163)
at org.orekit.estimation.leastsquares.BatchLSEstimator.estimate(BatchLSEstimator.java:435)

any thoughts on this? I see that many of my implementation is similar to what Clément Jonglez’s tutorial on his GitHub repository included in the repo, only Im using a TLE as initial condition instead of IOD using Gibbs method.
Thanks in advance!

5

Could you show me how do you initialize the NRLMSISE00 model?

Your problem is related to the data used by the model to compute the density. These data are stored in an external file available in the orekit-data. Do you have the orekit-data folder? Inside this folder, do you have a file named SpaceWeather-All-v1.2.txt?

Bryan

6

Yes, thank you Bryan, this is how I initialize the model:

from org.orekit.models.earth.atmosphere import NRLMSISE00
atmosphere = NRLMSISE00(msafe, sun, wgs84Ellipsoid)
from org.orekit.forces.drag import IsotropicDrag
isotropicDrag = IsotropicDrag(area_cd, cd)
from org.orekit.forces.drag import DragForce
dragForce = DragForce(atmosphere, isotropicDrag)
propagatorBuilder.addForceModel(dragForce)

the propagationBuilder is used later for:

from org.orekit.estimation.leastsquares import BatchLSEstimator
estimator = BatchLSEstimator(optimizer, propagatorBuilder)

I have the orekit-data folder and SpaceWeather-All-v1.2.txt file, and I use the “setup_orekit_curdir()” method.

7

Hi @powski,

Looks like Orekit cannot find your solar activity data.

Solar activity is used by the atmosphere models (in your case NRL-MSIS 00) to compute the atmospheric density which is then used to get the drag force.

Could you show use how you initialize your atmosphere model ?

Maxime

8

I have an idea. Looking at your initialization of the NRLMSISE00, it looks like you use the MSAFE data. The data in your orekit-data folder are maybe not up to date.

Could you try to run the update.sh script in your orekit-data folder in order to update the data?

Bryan

9

thank you @MaximeJ and @bcazabonne for your help, I updated the orekit-data with the latest available and the errors remain the same,

this is how I initialize the atmospheric model:

Atmospheric drag

from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation
#from org.orekit.data import DataProvidersManager, ZipJarCrawler, DirectoryCrawler, DataContext
msafe = MarshallSolarActivityFutureEstimation(
‘(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\p{Digit}\p{Digit}\p{Digit}\p{Digit}F10.(?:txt|TXT)’,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
DM = DataContext.getDefault().getDataProvidersManager()
DM.feed(msafe.getSupportedNames(), msafe) # Feeding the F10.7 bulletins to Orekit’s data manager

#from org.orekit.models.earth.atmosphere import NRLMSISE00
#atmosphere = NRLMSISE00(msafe, sun, wgs84Ellipsoid)
from org.orekit.models.earth.atmosphere import DTM2000
atmosphere = DTM2000(msafe, sun, wgs84Ellipsoid)
from org.orekit.forces.drag import IsotropicDrag
isotropicDrag = IsotropicDrag(area_cd, cd)
from org.orekit.forces.drag import DragForce
dragForce = DragForce(atmosphere, isotropicDrag)
propagatorBuilder.addForceModel(dragForce)

Sorry if its not straight forward, is my first time using the tools, its a mix of solutions I found in the examples and Ive tried to use the NRLMSISE00 or the DTM2000 models in different tests and the resulting errors are similar to those indicated yesterday.