There seems to be weird behaviors of the following method when elevation angle is very low (near zero) : org.orekit.models.earth.troposphere.SaastamoinenModel#pathDelay(double, double, double[], org.orekit.time.AbsoluteDate)
The Saastamoinen model is not defined for altitudes below 0.0. for continuity reasons, we use the value for h = 0 when altitude is negative.
However, in my experience, there can still be non-physical path delays computed when elevation=0. I suspect the tangent computation final double tan = FastMath.tan(z);
I have wrapped the class in a custom BoundedSaastamoinenModel that forces elevation to be higher than 1e-3 (I have taken a comfortable margin since it has no impact for my use case) and all my problems are gone.
Have I misused the SaastamoinenModel in some way ? Or is there a numerical issue here that could be fixed ?
I understand your problem and I don’t think it comes from your utilization of the model. So, we have to fix it!
A way to fix the problem is, as you proposed, to impose a minimum elevation value in the model. I think that the minimum value shall be an attribute of the class in order to be defined by a user. However, I am not sure if the attribute has to be in the model (SaastamoinenModel) or in the the tropospheric delay modifiers (RangeTroposphericDelayModifier, AngularTroposphericDelayModifier, etc).
In reality, I think it shall be in the measurement modifiers
Currently the measurement modifiers only compute the tropospheric delay if the elevation angle is greater than 0.0
// only consider measures above the horizon
if (elevation > 0) {
// altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// delay in meters
final double delay = tropoModel.pathDelay(elevation, height, tropoModel.getParameters(), state.getDate());
return delay;
}
return 0;
A good fix is to replace the 0.0 value by an attribute that can be defined by the user.
// only consider measures above the threshold
if (elevation > threshold) {
// altitude AMSL in meters
final double height = getStationHeightAMSL(station);
// delay in meters
final double delay = tropoModel.pathDelay(elevation, height, tropoModel.getParameters(), state.getDate());
return delay;
}
return 0;
Where threshold is the user-defined attribute (1e-3 in your example). This architecture imposes to have 2 class constructors: one with the default margin (i.e. 0.0) and one where the user can defined the threshold value.
This modification has to be performed for all the tropospheric delay modifiers and for the field and non-field methods in the modifiers.
What do you think ?
Can you open an issue on the issue tracker ?
Thank you for your answer. I like the idea of having the threshold as a user-defined attribute. However, I am not really comfortable having it that far from the numerical quirk. Since SaastamoinenModel is a public class, any user of Orekit can use it directly, without relying on one of the existing measurements modifiers. I feel like this class should know its own limitations, and protect the rest of the code against those limitations.
I’ve done some research on the Saastamoinen modified model, but could not find a recommended elevation range. So unless I have overlooked some paper, I guess we will have to decide our own threshold (and leaving it customizable for the user is probably a good idea).
I have also done some more tests on the pathDelay() method. After a debug session I can confirm that the divergence comes only from the computation of tan(z), with a zenith angle of \pi/2 when elevation is zero. All other terms seem fine. After a quick dichotomy it seems that the path delay is non-physical up to about 0.04 rad. For lower elevations, the return value is negative (which seems wrong to me) with a limit to -\infty when elevation goes towards zero.
A generic fix could be to have pathDelay() throw an exception when the path delay is negative. This way, we would not have to bother with a threshold. However, since the path delay is supposed to increase at low elevations, this would mean that some values that are barely above zero would be able to slip through despite being non-physical. Another solution could be to set up an arbitrary low elevation value (such as 0.05rad) and leave setters in case the user wants to change it.
With the first proposal (throw an exception), some difficulties can occur during orbit determination processes. Indeed, the elevation angle is computed inside the modifier. So, the user cannot control the computation of this elevation angle. In that respect, if an exception is thrown, the user can have some difficulties to fix the problem. The best way is to let the user to control the divergence problem of the model. That’s why I like the second proposal.
As you proposed, we can define a new attribute in the model: threshold. This attribute is not declare as final in order to be modifiable. It can be initialized to 0.05 as you proposed. Finally, a new method can be implemented inside the model in order to let the user to set a new value for the elevation threshold.
I have checked that the SaastamoinenModel class remains accessible for the user even though it is usually called by modifiers. It seems to be the case: all examples that I found accept the model as a constructor parameter. The value usually is SaastamoinenModel.getStandardModel() which the user could change to SaastamoinenModel.getStandardModel().setElevationThreshold(xxx).
I have opened a merge request. @bcazabonne if you have the time, could you please review my changes ? There is no rush, this is a much lower priority than the release of 10.1
