Inconsistency with classical reference

[haji-ali]

I am trying to reproduce the results from a classical reference "Satellite Conjunction Monte Carlo Analysis" by "Salvatore Alfano" in tudat. This reference uses the "two-body analytical propagation" to find the orbit of certain bodies in GEO/MEO and LEO. I wanted to see how different the results from tudat would be, so I modified the example singleSatellitePropagator.cpp to simulate Case 1 in that reference, which is a simple body orbiting the Earth's GEO, by simply setting the following initial data (instead of Asterix's), as taken from the reference (page 16)

Position:  -33552459.274056, -23728303.048015, 0.0
Velocity: -1828.997179397, 2534.1074695609, 0.0

At time 280800, the reference (page 17) has the following values for the mean position and velocity (the co-variance is relatively small, so the answer should not be far off)

Position: 1.53446765e+05,  4.18741559e+07,  0.00000000e+00
Velocity: 3.06687476e+03, -1.13736150e+01,  0.00000000e+00

However, tudat is returning this answer at time 280800 with time-step 10 seconds

Position: 1.22777978e+05,  4.18742582e+07,  0.00000000e+00, 
Velocity: 3.06688226e+03, -9.10041611e+00,  0.00000000e+00

As you can see, the errors in the x-position and y-velocity are relatively large (25% off). I tried other cases as well and saw similar errors in some coordinates.

Do you have any ideas for the reasons of this discrepancy? I did not expect such different results even though reference is using an analytical propagator. I also tried computing the mean (based on the covariance matrix in the reference) and got similar results.

Any help is much appreciated.

[DominicDirkx]

Dear Abdul,

Good to hear that you're using Tudat! Just out of curiosity, could you tell us a bit about how you found, and what you will be using it for? It's always nice to hear from our (new) users.

Your difference may seem a bit large, but remember: linear difference is a better measure for accuracy than a percentage of a single component. The components you are referring to are quite small (123 km for x-position, vs 42000 km for y-position). The difference w.r.t. your reference are about 30 km in x-position, and about 2 m/s in y-velocity. Indeed these values are large, but consider the following possible reasons:

• You are using a different value of the Earth's gravitational parameter in the two cases
• You are using a numerical integrator with a large time step and/or high tolerance

My guess is, this will be the first problem.

The test you are doing has been performed many times (and is part of our unit tests), and with identical physical parameters, and a high-accuracy integrator, the error will be much lower than this,

Best,

Domninic

[haji-ali]

Dear Dominic,

Thanks for you response. I was of course not doubting tudat, just wondering what I was doing wrong. It turns out that the issue was indeed a step-size one. Apologies for the spam.

Can you please elaborate on the "tolerance" parameter that you mentioned? I would like to control the accuracy and I thought the only accuracy parameter is the step size. Is there another parameter I should be controlling?

I found tudat through a post on StackOverflow. I was looking for a software package that I can use as a black-box orbit propagator for my project on uncertainty quantification for collision events.

[DominicDirkx]

Hi Abdul,

No worries, I didn't interpret it as criticism, just wanted to assure you that getting it to come out right is a matter of getting the right settings.

The tolerances are settings that are available in variable step-size integrators, where the step size is autonomously adapted by the integrator, so that the error (approximately) falls in a certain range. This can for instance be useful for elliptical orbits, where the step size at periapsis should be much lower than at apoapsis.

I think that Tudat could be a very useful tool for a project like the one you mention. It may be good to know that we're currently in the process of making a conda installation for Tudat, as well as a Python interface. This will be at least ~3-4 months before we're ready to roll everything out, but may be good to know about.

Our support during the summer period will be very limited, but do feel free to ask some questions here if you cannot find what you're trying to do on our website,

Best,

Dominic