Hi all again.
After several days of testing, I come back to you with what I think is the same problem, but this time I’ve investigated it a bit more. I think it’s worth it to re-formulate the problem since last time I hadn’t checked against other derivative computation methods. I have also written a small piece of code that clearly shows the difference between the two and is easily modifiable to see the problem appear and disappear as the configuration is changed (propagation time, initial altitude, etc).
So again, I’m propagating from a low altitude (initial semi-major axis 6550 km and nearly zero eccentricity) for 6000 seconds with a FieldNumericalPropagator. As in the previous case, the integrator I’m using is a DormandPrince853 (field version), but the problem also happens with other variable stepsize integrators available in Hipparchus. Other orbit parameters are basically left at random. As for propagation parameters, tolerances and stepsizes limits are as described in my first post, but changing them to more “reasonable” values doesn’t seem to fix the problem (this is easy to change in the example provided down below). Force models applied are the following:
- Newtonian grav. attraction
- HolmesFeatherstone Attraction Model (18,18).
- Sun and Moon as third bodies.
- Isotropic Drag configured with a HarrisPriester atmosphere (S = 10m2, Cd = 2.0).
- SolarRadiationPressure configured with IsotropicRadiationCNES95Convention.
To monitor the value of the derivatives, the ones from the first component of the PV array (x position) are extracted from the DerivativeStructure object (in the Field version). For comparison, a PartialDerivativeEquations object is used in the Non-Field propagation and the first row of the Jacobian obtained after propagation is checked against the extracted derivatives in the Field version.
The basic problem is that after propagating for a certain period of time with the Field integrator with an adaptive stepsize, the derivatives of the PV variables with respect to the initial PV diverge. Trying to isolate the problem has been unsuccessful for me, since it seems to depend on a lot of different parameters. I’ll try to list them as thoroughly as I can:
- First of all, the problem dissapears when the Drag and SolarRadiationPressure are not added. This is already a big hint that the problem has to do with at least one of these ForceModels.
- Most of the times, leaving the SolarRadiationPressure on and taking Drag off the propagator fixes the problem (but it’s obsviously not well suited to re-entry cases).
- Derivatives diverge when the Drag is active, under the cited conditions, but if SolarRadiationPressure is also active the resulting derivatives are even bigger (by some orders of magnitude). This is why I don’t completely rule out SRP as a source of problems.
- As said in my first post, either making the tolerance big enough in the variable stepsize integrator or integrating with the fixed stepsize one (with no regards to tolerances and a fixed step of around 1 second) makes the derivatives “not diverge”, but the results are also not the same as the ones obtained though the PDE object (which seem to be the best guess of the three).
- Altitude is a key factor here. Contrary to what I thought when I started this post, it’s not Altitude = 0 that creates the problem. This happens even before the HarrisPriester model runs out of altitude (100 km ). At that altitude, all other models should be comfortable enough to not arrive at any singularities. However, at a higher altitude (initial semi-major axis = 6850 km) the problem goes away.
- I have also tested this with a SimpleExponentialAtmosphere, and the problem is still there, it only happens at a different altitude.
This quickly-made piece of code should illustrate the issue well enough. What it does is create both propagators (field and non-field) basically side-by-side, with the same input and observing the “same” output. Orbital parameters are set at the beggining and the propagation time (which can be lowered to see the effect on the derivatives) is the same for both cases.
FieldBugTest.java (9.7 KB)
(It will look for the orekit-data folder in the Home directory, but appart from that it should work just fine.)
The output of this code are two printed vectors, showing the values of the previously mentioned derivatives side-by-side. If run as it is, the code will produce vastly different results from both methods. If the initial semi-major axis is increased 100 km, the problem virtually dissapears and both methods produce very similar results.
I’m still not sure this is not something I’ve done rather than a bug, but in that case I guess you can tell me what’s wrong with my whole setup.
Thank you very much in advance. Best regards,