Use of DerivativeStructure to obtain a rotation

Orekit usage
1

Hello Orekit community,

I have a question regarding how to use FieldVector3D to obtain a Transform. For giving a bit more of context I am trying to create a new eventDetector that will trigger whenever the sun is in the field of view of the satellite’s sensor.

Prior to start working with derivatives, the method worked as expected. However, I did not get the readings related to the rotation rate and acceleration. That is why I am trying to implement the derivatives to the model.

Now, I am not able to apply properly the FieldRotation method, so the resulting transform does not have the expected values of rotation. Furthermore, the rotation rate and acceleration obtained are still zero.

Is there an error in the way I am defining the respective FieldVector3D for the rotation? The goal is to have the solar arrays always pointing in the sun direction.

Here is an example of the test I am running to compare the obtained rotation:

public class FieldRotationTest{

    @Test
    public void testFieldRotation() {

        AbsoluteDate date = new AbsoluteDate();

        Frame eme2000 = FramesFactory.getEME2000();

        KeplerianOrbit kep = new KeplerianOrbit(7000.e3, 1e-3, 1.60, 1, 2, 0, PositionAngle.MEAN,
                eme2000, date, Constants.EGM96_EARTH_MU);

        // relative sun pvCoord
        // in inertial frame
        Vector3D satSunDirection = computeSunDirection(kep, date, eme2000);

        FieldVector3D<DerivativeStructure> pointingSat = kep.getPVCoordinates(date, eme2000).toDerivativeStructureVector(2).normalize();
        DerivativeStructure zeroWithDerivatives = pointingSat.getX().getFactory().variable(0, 0.0);
        FieldAbsoluteDate fDate = new FieldAbsoluteDate(date, zeroWithDerivatives);

        FieldVector3D<DerivativeStructure> solarPanelNormalSatFieldVector = new FieldVector3D<>(fDate.getField(), Vector3D.PLUS_J).normalize();
        Vector3D solarPanelNormalSatVector = solarPanelNormalSatFieldVector.toVector3D();

        Rotation expectedRotation = buildRotation(satSunDirection, solarPanelNormalSatVector);
        Transform derivativeFieldTransform = buildRotationWithFieldDerivatives(satSunDirection, fDate, solarPanelNormalSatFieldVector);
        boolean equal = compareRotations(expectedRotation, derivativeFieldTransform.getRotation());

        System.out.println("Both cases represent the same rotation? " + equal);
        System.out.println("Rotation rate: " + ArrayUtils.toString(derivativeFieldTransform.getRotationRate()));
        System.out.println("Acceleration rate: " + ArrayUtils.toString(derivativeFieldTransform.getRotationAcceleration()));
    }

    public Vector3D computeSunDirection(PVCoordinatesProvider pvProv, AbsoluteDate date, Frame frame) {
        PVCoordinates sunPVCoord = new PVCoordinates(pvProv.getPVCoordinates(date, frame),
                CelestialBodyFactory.getSun().getPVCoordinates(date, frame));
        return sunPVCoord.getPosition();
    }

    public Rotation buildRotation(Vector3D satSunDirection, Vector3D solarPanelNormalSat) {

        Rotation inertial2BodyRotation = new Rotation(satSunDirection, solarPanelNormalSat);

        return inertial2BodyRotation;

    }

    public Transform buildRotationWithFieldDerivatives(Vector3D satSunDirection, FieldAbsoluteDate fDate, FieldVector3D<DerivativeStructure> solarPanelNormalSatFieldVector) {

        FieldVector3D<DerivativeStructure> pointedBody = new FieldVector3D<>(fDate.getField(), satSunDirection).normalize();
        FieldVector3D<DerivativeStructure> solarPanelSat = new FieldVector3D<>(fDate.getField(), Vector3D.PLUS_K).normalize();

        FieldRotation<DerivativeStructure> celToSatRotationDS =
                new FieldRotation<DerivativeStructure>(pointedBody, solarPanelSat, solarPanelNormalSatFieldVector, solarPanelSat);

        Transform transform = new Transform(fDate.toAbsoluteDate(), new AngularCoordinates(celToSatRotationDS));

        return transform;

    }

    public boolean compareRotations(Rotation r1, Rotation r2) {

        // Compare the rotation angles
        double angleTolerance = 1e-12; // Tolerance for comparing angles
        double angleDifference = r1.getAngle() - r2.getAngle();
        if (Math.abs(angleDifference) > angleTolerance) {
            return false;
        }

        // Compare the rotation axes
        Vector3D axis1 = r1.getAxis(RotationConvention.VECTOR_OPERATOR);
        Vector3D axis2 = r2.getAxis(RotationConvention.VECTOR_OPERATOR);
        if (!axis1.equals(axis2)) {
            return false;
        }

        return true;
    }


}

Thank you very much for the help.
Adrian

2

Hi @adri

Welcome to the Orekit forum!

First of all, I think the FieldPropagation.java tutorial could help you a lot.

In addition, please fin below some comments:

  • In the code below, I recommend you to start your DSFactory as soon as possible.
        AbsoluteDate date = new AbsoluteDate();

        Frame eme2000 = FramesFactory.getEME2000();

        KeplerianOrbit kep = new KeplerianOrbit(7000.e3, 1e-3, 1.60, 1, 2, 0, PositionAngle.MEAN,
                eme2000, date, Constants.EGM96_EARTH_MU);

For instance:

        final DSFactory factory = new DSFactory(params, order);
        final DerivativeStructure a    = factory.variable(0, 7000.e3);
        final DerivativeStructure e    = factory.variable(1, 1e-3);
        final DerivativeStructure i    = factory.variable(2, 1.60);
        final DerivativeStructure pa   = factory.variable(3, 1);
        final DerivativeStructure raan = factory.variable(4,  2);
        final DerivativeStructure ni   = factory.variable(5, 0);
        final DerivativeStructure  mu = factory.constant(Constants.EGM96_EARTH_MU);

        FieldAbsoluteDate<DerivativeStructure> date = new FieldAbsoluteDate(a.getField());
        FieldKeplerianOrbit<DerivativeStructure> kep = new FieldKeplerianOrbit(a, e, i, pa, raan, ni, PositionAngle.MEAN, eme2000, date, mu);

Using field a soon a possible is very important to be sure to don’t lost any derivative.

  • For this

DerivativeStructure zeroWithDerivatives = pointingSat.getX().getFactory().variable(0, 0.0);

I recommend you to do

DerivativeStructure zeroWithDerivatives = a.getZero();

  • The method below shall be implemented using Field:

      // relative sun pvCoord
  // in inertial frame
  Vector3D satSunDirection = computeSunDirection(kep, date, eme2000);

  • Based on the first comment, update

 FieldVector3D<DerivativeStructure> pointingSat = kep.getPVCoordinates(date, eme2000).toDerivativeStructureVector(2).normalize();

With

FieldVector3D<DerivativeStructure> pointingSat = kep.getPVCoordinates(date, eme2000).normalize();

  • You can remove this:

FieldAbsoluteDate fDate = new FieldAbsoluteDate(date, zeroWithDerivatives);

  • Try to not use the constructor FieldVector3D(field, Vector3D) but directly FieldVector3D parameters. In other words, satSunDirection shall be directly a FieldVector3D.

To summarize, Field parameters shall be initialized as soon as possible and used in the code.

I hope these comments will help you!

Best regards,
Bryan

3

Thank you very much for your answer @bcazabonne , I have adapted the definition of all the variables according to the derivativeStructure, but I still have a doubt. My goal is to get a rotation with this definition using derivativeStructures, but also provide me with the rotationRate and rotationAcceleration a part of the quaternions.

I used to achieve this with the following line of code:

FieldVector3D<DerivativeStructure> pointingSat = kep.getPVCoordinates(date, eme2000).toDerivativeStructureVector(2).normalize();

Now how can I get it with this new definition, because obviously when using DSFactory I get a lot more derivatives with respect the different keplerian elements and the dimensions end up being inconsistent.

Transform derivativeFieldTransform = buildRotationWithFieldDerivatives(satSunDirection, date, solarPanelNormalSatFieldVector);

My goal was to only have 2 derivatives with respect to time so that I could extract the previously mentioned rotation data.

Thanks in advance,

Adrian

4

In order to have time derivatives, you could look at UnivariateDerivatives2 class instead of DerivativeStructure.

Bryan