Using Fields for time-derivatives in measurement class

I’m writing a measurement class for Ra/Dec and Ra/Dec-rates. I’m comfortable calculating the time-derivatives by hand, but thought it might be a good to try out the auto-diff methods that are available.

But I don’t know enough about these techniques and I’m stuck. Comparing to the AngularRaDec class, instead of

   // Station-satellite vector expressed in inertial frame
   final FieldVector3D<Gradient> staSatInertial = transitStateDS.getPosition().subtract(stationDownlink.getPosition());

I can compute

    // Station-satellite vector expressed in inertial frame
   final FieldVector3D<FieldUnivariateDerivative1<Gradient>> staSatInertial =
           transitStateDS.toUnivariateDerivative1Vector()
                   .subtract(stationDownlink.toUnivariateDerivative1Vector());

Which means that

    final Gradient baseRightAscension = staSatReference.getAlpha().getValue();
    final Gradient rightAscensionRate = staSatReference.getAlpha().getFirstDerivative();

However, in between the calculation of staSatInertial and the angles, there is a transformation from the frame of the spacecraft state into the frame of the measurements (which produces staSatReference). It currently looks like this

    // Field transform from inertial to reference frame at station's reception date
    final FieldTransform<Gradient> inertialToReferenceDownlink =
                    state.getFrame().getTransformTo(referenceFrame, downlinkDateDS);

but I would like this

   final FieldTransform<FieldUnivariateDerivative1<Gradient>> inertialToReferenceDownlink =
           state.getFrame().getTransformTo(referenceFrame, downlinkDateDS);

The issue is that the downlinkDateDS is an FieldAbsoluteDate<Gradient>, not FieldAbsoluteDate<FieldUnivariateDerivative1<Gradient>>. I don’t know how to include the time derivative structure into the date as well as the spatial derivatives (or even if that makes sense)?

This “works”, but is it correct?

   // Field transform from inertial to reference frame at station's reception date
   final Gradient delta1 = downlinkDateDS.durationFrom(getDate());
   final FieldAbsoluteDate<FieldUnivariateDerivative1<Gradient>> downlinkDateDS1 =
           new FieldAbsoluteDate<>(getDate(), new FieldUnivariateDerivative1<>(delta1, GradientField.getField(nbParams).getZero()));
   final FieldTransform<FieldUnivariateDerivative1<Gradient>> inertialToReferenceDownlink =
           state.getFrame().getTransformTo(referenceFrame, downlinkDateDS1);

Hi @markrutten ,

That’s an interesting implementation of the time-derivatives of Ra-Dec values.

In order to have the delta1 value in terms of FieldUnivariateDerivative1<Gradient> you can directly do the following instruction:


final FieldUnivariateDerivative1<Gradient> delta1 = new FieldUnivariateDerivative1<>(delta, GradientField.getField(nbParams).getOne());

It means that dt/dt = 1.0

This is something close to what we do to initialize the UnivariateDerivative2 instances in analytical orbit propagation models like GNSS, Eckstein-Hechler and Brouwer-Lyddane.

After that, you can continue the implementation of the method. It will compile. I’m curious to know if at the end you have good values for the time-derivatives of right-ascension and declination angles. Maybe you could do a cross validation between the auto-diff and the “by hand” methods.

Will it be possible to share the results?

Best regards,
Bryan

Thanks @bcazabonne!

I have some promising results. It’s a test on a single GEO object, so not thorough, but the results are consistent.

I compared the rates calculated with auto-diff against rates calculated “by hand” within the measurement class

    // Compute RA and Dec rates
    final Gradient xySqSum = (staSatReference.getX().multiply(staSatReference.getX()))
            .add(staSatReference.getY().multiply(staSatReference.getY()));
    final Gradient rightAscensionRate = ((staSatReference.getX().multiply(staSatReferenceVel.getY()))
            .subtract(staSatReference.getY().multiply(staSatReferenceVel.getX()))).divide(xySqSum);

    final Gradient declinationRate = (staSatReferenceVel.getZ().multiply(xySqSum)
            .subtract(staSatReference.getX().multiply(staSatReference.getZ()).multiply(staSatReferenceVel.getX()))
            .subtract(staSatReference.getY().multiply(staSatReference.getZ()).multiply(staSatReferenceVel.getY())))
            .divide(staSatReference.getNormSq()).divide(xySqSum.sqrt());

against values calculated outside the measurement classes

                // Angles accounting for light-time
                TimeStampedPVCoordinates objectPV = state.getPVCoordinates(inertialFrame);
                TimeStampedPVCoordinates sitePV = groundStation.getBaseFrame()
                        .getPVCoordinates(state.getDate(), inertialFrame);
                double tof = AbstractMeasurement.signalTimeOfFlight(objectPV, sitePV.getPosition(), state.getDate());
                TimeStampedPVCoordinates obsPV = objectPV.shiftedBy(-tof);

                Vector3D relPosition = obsPV.getPosition().subtract(sitePV.getPosition());
                Vector3D relVelocity = obsPV.getVelocity().subtract(sitePV.getVelocity());

                // Right ascension
                double rightAscension = relPosition.getAlpha();

                double xySqSum = relPosition.getX() * relPosition.getX() +
                        relPosition.getY() * relPosition.getY();
                double rightAscensionRate = ((relPosition.getX() * relVelocity.getY())
                        - (relPosition.getY() * relVelocity.getX())) / xySqSum;

                // Declination
                double declination = relPosition.getDelta();

                double declinationRate = (relVelocity.getZ() * xySqSum -
                        relPosition.getX() * relPosition.getZ() * relVelocity.getX() -
                        relPosition.getY() * relPosition.getZ() * relVelocity.getY()) /
                        relPosition.getNormSq() / Math.sqrt(xySqSum);

An example of the comparison is

[0.25777076215066574, 0.048354465779878196, 7.290022488455892E-5, 7.936399906260245E-7]
[0.25777076215066574, 0.048354465779878196, 7.290023186254899E-5, 7.936405219582532E-7]
RA: 0.2577707621506837 :: 6.955832266661889E-14 :: 6.955832266661889E-14
Dec: 0.048354465779515896 :: 7.492595803953079E-12 :: 7.492595803953079E-12
RA-rate: 7.290023200947462E-5 :: 9.773516766510257E-8 :: 2.0154343346443987E-9
Dec-rate: 7.936429678314021E-7 :: 3.7513298382235697E-6 :: 3.0818400538887233E-6

The first line is the theoretical measurement calculated using the derivative structures, the second line is with rates calculated explicitly within the measurement class. The last 4 lines compares those values against the values I computed without any derivative structures (as above). The first column are the calculated values, the second and third are relative errors against the measurement classes.

So RA and Dec are computed within numerical precision. The rates are within 6 or so digits of precision of each other.

The spacecraft state I was passing to theoreticalEvaluation was produced from a TLE, so in TEME coordinates. If I explicitly convert that state to one in the same coordinates as the measurements (I am using GCRF), then the results from all three methods match to within numerical precision:

[0.2577707621506837, 0.04835446577951594, 7.29002320094746E-5, 7.93642967831405E-7]
[0.2577707621506837, 0.04835446577951594, 7.29002320094746E-5, 7.936429678314052E-7]
RA: 0.2577707621506837 :: 0.0 :: 0.0
Dec: 0.048354465779515896 :: 8.610034823521971E-16 :: 8.610034823521971E-16
RA-rate: 7.290023200947462E-5 :: 1.8590513064906882E-16 :: 1.8590513064906882E-16
Dec-rate: 7.936429678314021E-7 :: 3.602043126261972E-15 :: 3.86886113561471E-15

I’m happy to see that the auto-diff method produces results close to the “by hand” method. Congratulations for the implementation!

I have a quick question. I’m not familiar with time-derivatives angular measurements. Is it something you use often (or other people you know)? If so, is there an interest in adding an AngularRaDecRate class in Orekit?

Bryan

Good question! I haven’t used angle rates before and I don’t think they would be widely used. We’re trying it out as a way of summarising a sequence of very high time-rate angular measurements, in order to reduce comms bandwidth. If our experiment is successful I’ll ask if it’s OK to contribute the angle rates class.

1 Like