CartesianAdjointJ2Term.java

  1. /* Copyright 2022-2025 Romain Serra
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.control.indirect.adjoint;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.analysis.differentiation.FieldGradient;
  22. import org.hipparchus.analysis.differentiation.FieldGradientField;
  23. import org.hipparchus.analysis.differentiation.Gradient;
  24. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  25. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  26. import org.hipparchus.util.MathArrays;
  27. import org.orekit.forces.gravity.J2OnlyPerturbation;
  28. import org.orekit.frames.FieldStaticTransform;
  29. import org.orekit.frames.Frame;
  30. import org.orekit.frames.StaticTransform;
  31. import org.orekit.time.AbsoluteDate;
  32. import org.orekit.time.FieldAbsoluteDate;

  34. /**
  35.  * Class defining a (constant) J2 contributions in the adjoint equations for Cartesian coordinates.
  36.  * If present, then the propagator should also include a constant J2 term (oblateness) of the central body.
  37.  * @author Romain Serra
  38.  * @see CartesianAdjointEquationTerm
  39.  * @see org.orekit.forces.gravity.J2OnlyPerturbation
  40.  * @since 12.2
  41.  */
  42. public class CartesianAdjointJ2Term extends AbstractCartesianAdjointGravitationalTerm {

  44.     /** J2 coefficient. */
  45.     private final double j2;

  47.     /** Equatorial radius of central body. */
  48.     private final double rEq;

  50.     /** Frame where J2 applies. */
  51.     private final Frame j2Frame;

  53.     /**
  54.      * Constructor.
  55.      * @param mu central body gravitational parameter.
  56.      * @param rEq equatorial radius
  57.      * @param j2 J2 coefficient
  58.      * @param j2Frame J2 frame
  59.      */
  60.     public CartesianAdjointJ2Term(final double mu, final double rEq, final double j2,
  61.                                   final Frame j2Frame) {
  62.         super(mu);
  63.         this.j2 = j2;
  64.         this.rEq = rEq;
  65.         this.j2Frame = j2Frame;
  66.     }

  68.     /**
  69.      * Getter for central body equatorial radius.
  70.      * @return equatorial radius
  71.      */
  72.     public double getrEq() {
  73.         return rEq;
  74.     }

  76.     /**
  77.      * Getter for J2.
  78.      * @return J2 coefficient
  79.      */
  80.     public double getJ2() {
  81.         return j2;
  82.     }

  84.     /** {@inheritDoc} */
  85.     @Override
  86.     public double[] getPositionAdjointContribution(final AbsoluteDate date, final double[] stateVariables,
  87.                                                    final double[] adjointVariables, final Frame frame) {
  88.         final double[] contribution = new double[3];
  89.         final int numberOfGradientVariables = 3;
  90.         final FieldVector3D<Gradient> position = new FieldVector3D<>(Gradient.variable(numberOfGradientVariables, 0, stateVariables[0]),
  91.             Gradient.variable(numberOfGradientVariables, 1, stateVariables[1]),
  92.             Gradient.variable(numberOfGradientVariables, 2, stateVariables[2]));
  93.         final StaticTransform transform = frame.getStaticTransformTo(j2Frame, date);
  94.         final FieldVector3D<Gradient> positionInJ2Frame = transform.transformPosition(position);
  95.         final Gradient fieldJ2 = Gradient.constant(numberOfGradientVariables, j2);
  96.         final FieldVector3D<Gradient> accelerationInJ2Frame = J2OnlyPerturbation.computeAccelerationInJ2Frame(positionInJ2Frame,
  97.                 getMu(), rEq, fieldJ2);
  98.         final FieldVector3D<Gradient> acceleration = transform.getStaticInverse().transformVector(accelerationInJ2Frame);
  99.         final double pvx = adjointVariables[3];
  100.         final double pvy = adjointVariables[4];
  101.         final double pvz = adjointVariables[5];
  102.         final double[] gradientAccelerationX = acceleration.getX().getGradient();
  103.         final double[] gradientAccelerationY = acceleration.getY().getGradient();
  104.         final double[] gradientAccelerationZ = acceleration.getZ().getGradient();
  105.         contribution[0] = -(gradientAccelerationX[0] * pvx + gradientAccelerationY[0] * pvy + gradientAccelerationZ[0] * pvz);
  106.         contribution[1] = -(gradientAccelerationX[1] * pvx + gradientAccelerationY[1] * pvy + gradientAccelerationZ[1] * pvz);
  107.         contribution[2] = -(gradientAccelerationX[2] * pvx + gradientAccelerationY[2] * pvy + gradientAccelerationZ[2] * pvz);
  108.         return contribution;
  109.     }

  111.     /** {@inheritDoc} */
  112.     @Override
  113.     public <T extends CalculusFieldElement<T>> T[] getPositionAdjointFieldContribution(final FieldAbsoluteDate<T> date,
  114.                                                                                        final T[] stateVariables,
  115.                                                                                        final T[] adjointVariables,
  116.                                                                                        final Frame frame) {
  117.         final Field<T> field = adjointVariables[0].getField();
  118.         final T[] contribution = MathArrays.buildArray(field, 3);
  119.         final int numberOfGradientVariables = 3;
  120.         final FieldVector3D<FieldGradient<T>> position = new FieldVector3D<>(FieldGradient.variable(numberOfGradientVariables, 0, stateVariables[0]),
  121.             FieldGradient.variable(numberOfGradientVariables, 1, stateVariables[1]),
  122.             FieldGradient.variable(numberOfGradientVariables, 2, stateVariables[2]));
  123.         final T shift = date.durationFrom(date.toAbsoluteDate());
  124.         final FieldGradientField<T> gradientField = FieldGradientField.getField(field, 3);
  125.         final FieldAbsoluteDate<FieldGradient<T>> gradientDate = new FieldAbsoluteDate<>(gradientField, date.toAbsoluteDate())
  126.             .shiftedBy(FieldGradient.constant(numberOfGradientVariables, shift));
  127.         final FieldStaticTransform<FieldGradient<T>> transform = frame.getStaticTransformTo(j2Frame, gradientDate);
  128.         final FieldVector3D<FieldGradient<T>> positionInJ2Frame = transform.transformPosition(position);
  129.         final FieldGradient<T> fieldJ2 = FieldGradient.constant(numberOfGradientVariables, field.getZero().newInstance(j2));
  130.         final FieldVector3D<FieldGradient<T>> accelerationInJ2Frame = J2OnlyPerturbation.computeAccelerationInJ2Frame(positionInJ2Frame,
  131.                 getMu(), rEq, fieldJ2);
  132.         final FieldVector3D<FieldGradient<T>> acceleration = transform.getStaticInverse().transformVector(accelerationInJ2Frame);
  133.         final T pvx = adjointVariables[3];
  134.         final T pvy = adjointVariables[4];
  135.         final T pvz = adjointVariables[5];
  136.         final T[] gradientAccelerationX = acceleration.getX().getGradient();
  137.         final T[] gradientAccelerationY = acceleration.getY().getGradient();
  138.         final T[] gradientAccelerationZ = acceleration.getZ().getGradient();
  139.         contribution[0] = gradientAccelerationX[0].multiply(pvx).add(gradientAccelerationY[0].multiply(pvy)).add(gradientAccelerationZ[0].multiply(pvz));
  140.         contribution[1] = gradientAccelerationX[1].multiply(pvx).add(gradientAccelerationY[1].multiply(pvy)).add(gradientAccelerationZ[1].multiply(pvz));
  141.         contribution[2] = gradientAccelerationX[2].multiply(pvx).add(gradientAccelerationY[2].multiply(pvy)).add(gradientAccelerationZ[2].multiply(pvz));
  142.         contribution[0] = contribution[0].negate();
  143.         contribution[1] = contribution[1].negate();
  144.         contribution[2] = contribution[2].negate();
  145.         return contribution;
  146.     }

  148.     /** {@inheritDoc} */
  149.     @Override
  150.     public Vector3D getAcceleration(final AbsoluteDate date, final double[] stateVariables,
  151.                                     final Frame frame) {
  152.         final StaticTransform transform = frame.getStaticTransformTo(j2Frame, date);
  153.         final Vector3D positionInJ2Frame = transform.transformPosition(new Vector3D(stateVariables[0], stateVariables[1], stateVariables[2]));
  154.         final Vector3D accelerationInJ2Frame = J2OnlyPerturbation.computeAccelerationInJ2Frame(positionInJ2Frame,
  155.                 getMu(), rEq, getJ2());
  156.         return transform.getStaticInverse().transformVector(accelerationInJ2Frame);
  157.     }

  159.     /** {@inheritDoc} */
  160.     @Override
  161.     public <T extends CalculusFieldElement<T>> FieldVector3D<T> getFieldAcceleration(final FieldAbsoluteDate<T> date,
  162.                                                                                      final T[] stateVariables,
  163.                                                                                      final Frame frame) {
  164.         final FieldStaticTransform<T> transform = frame.getStaticTransformTo(j2Frame, date);
  165.         final FieldVector3D<T> positionInJ2Frame = transform.transformPosition(new FieldVector3D<>(stateVariables[0], stateVariables[1], stateVariables[2]));
  166.         final FieldVector3D<T> accelerationInJ2Frame = J2OnlyPerturbation.computeAccelerationInJ2Frame(positionInJ2Frame,
  167.                 getMu(), rEq, date.getField().getZero().newInstance(getJ2()));
  168.         return transform.getStaticInverse().transformVector(accelerationInJ2Frame);
  169.     }
  170. }
Created with JaCoCo 0.8.12.202403310830