If I give real measurements in the form of topocentric right ascension and declination from a ground station to an Orekit’s Orbit Determination process like the Batch Least Squares, I should NOT correct beforehand for the so-called aberration (a.k.a. light time correction, from the fact that the signal takes a non-zero duration from its reflection on the satellite to its arrival at the sensor, both moving) right?
I am asking as a look to the implementation of AngularRaDec seems to show that this effect is systematically included when computing theoretical measurements, with no way of disabling it.
Yes, you are right.
This is because the effect changes as the orbit is updated. When you want to perform precise orbit determination, you cannot rely on preprocessing removing this effect properly as it does not know the final orbit.
Thanks @Serrof and @luc for discussing this interesting topic. I was actually looking to correct for annual and diurnal aberrations in Orekit (as far as I understand they still need to be corrected for when retrieving topocentric RA-Dec from the astrometric solution of telescope images). I didn’t find how to do it in Orekit, and, while I know how to correct for the annual aberration using other libraries or formulas, I cannot find reliable sources in general on how to correct for the diurnal one. Is there some way to do this in Orekit?
Yes, you have to set up a custom implementation of EstimationModifier for the angular measurement.
We do not provide it yet in Orekit but it would be an interesting contribution.
This effect should remain separated from the raw AngularRaDec measurement and remain implemented as a modifier because angular measurements that are not based on relative direction with respect to background stars must not include this. So applying or not the modifier depends on the method used to perform the measurement.
I am reposting here as I am trying (and failing) to create an observation modifier on these AngularRaDec.
To understand the problem I have created a dummy modifier, that simply recomputes the Ra-Dec, but the result is different to what the estimator has. My question is why?
Here is the sampled code in Python (sorry for the language):
def modify(self, estimated: EstimatedMeasurement) -> None:
# Retrieve data
measurement = AngularRaDec.cast_(estimated.getObservedMeasurement())
ra, dec = estimated.getEstimatedValue()
station = measurement.getStation()
base_frame = station.getBaseFrame()
state = estimated.getStates()[0]
pv_coord = state.getPVCoordinates()
date = state.getDate()
inertial = state.getFrame()
# Get state with light time delay
transform_inertial_to_base = inertial.getStaticTransformTo(base_frame, date)
transform_base_to_inertial = transform_inertial_to_base.getInverse()
position_station = transform_base_to_inertial.transformPosition(Vector3D(0., 0., 0.))
timestamped_state = TimeStampedPVCoordinates(date, pv_coord.getPosition(), pv_coord.getVelocity())
dt = AbstractMeasurement.signalTimeOfFlight(timestamped_state, position_station, date)
position = (state.shiftedBy(-dt)).getPVCoordinates().getPosition()
# Compute elevation in radians
elevation = base_frame.getElevation(position, inertial, date)
# Convert to RA-Dec
azimuth = base_frame.getAzimuth(position, inertial, date)
unit_vec = Vector3D(math.pi / 2. - azimuth, elevation)
transformed = transform_base_to_inertial.transformVector(unit_vec)
right_ascension = transformed.getAlpha()
declination = transformed.getDelta()
delta_ra = right_ascension - ra
delta_dec = declination - dec
estimated.setEstimatedValue([right_ascension, declination])
To be clear, the delta_ra and delta_dec are non-zero, which should not be the case in my understanding.
Here the date will be the one of the state and different from the measurement date, since the light time correction was already applied when constructing the EstimatedMeasurement.
I think the Transform should be computed at measurement date (using measurement.getDate()) instead of using date which is the date of the state (already corrected of light time).
Light time correction was already applied to state so you’re applying it twice.
If you look in Orekit, it’s doing:
// Downlink delay
final Gradient tauD = signalTimeOfFlight(pvaDS, stationDownlink.getPosition(), downlinkDateDS);
// Transit state
final Gradient delta = downlinkDateDS.durationFrom(state.getDate());
final Gradient deltaMTauD = tauD.negate().add(delta);
final SpacecraftState transitState = state.shiftedBy(deltaMTauD.getValue());
It’s using delta: the already-taken-into-account light time correction at previous iteration (of the batch LS).
In your case since your orbit didn’t change, I think you will get deltaMTauD = 0 and the transitState will be the same as State.
So I think that, basically, in your case you can just drop the light time correction part in your code.
so I removed the redundant time delay on the state and used the measurement date for the transformations.
I am hence down to <1e-7 arcsecond difference, yet the OD fails (due to NaN in integration) if I keep the (dummy) modifier, and converges fine without.
Any ideas why?
Ok so the difference modulo a full revolution is small, but the estimated measurement sometimes contains values for the right ascension above 2pi, which the OD cannot handle well.
// Compute right ascension and declination
final Gradient baseRightAscension = staSatReference.getAlpha();
final double twoPiWrap = MathUtils.normalizeAngle(baseRightAscension.getReal(),
getObservedValue()[0]) - baseRightAscension.getReal();
final Gradient rightAscension = baseRightAscension.add(twoPiWrap);
The estimated value baseRightAscension is wrapped in a in a 2\pi interval around the observed value to avoid the case you encountered.
yes and that’s a good thing, but the fact is that I am NOT passing any Right Ascension above 2pi as observed measurements, so am I smelling a bug or is the wrapping a bit unconventional?
Edit: I think I get it, for convergence reasons it chooses the angle closest to the observed one, so if the latter is around 6 for example the estimated would be 6.6 instead of 6.6 - 2pi
Actually a similar trick to allow declination angles in [0, pi/2] U [3pi/2, 2pi] could be good, because if you look at TDM’s definition, it is not imposed to be in [-pi/2, pi/2].