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;
18  
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;
33  
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 {
43  
44      /** J2 coefficient. */
45      private final double j2;
46  
47      /** Equatorial radius of central body. */
48      private final double rEq;
49  
50      /** Frame where J2 applies. */
51      private final Frame j2Frame;
52  
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      }
67  
68      /**
69       * Getter for central body equatorial radius.
70       * @return equatorial radius
71       */
72      public double getrEq() {
73          return rEq;
74      }
75  
76      /**
77       * Getter for J2.
78       * @return J2 coefficient
79       */
80      public double getJ2() {
81          return j2;
82      }
83  
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     }
110 
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     }
147 
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     }
158 
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 }