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.geometry.euclidean.threed.FieldVector3D;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.util.FastMath;
23  import org.hipparchus.util.MathArrays;
24  
25  /**
26   * Abstract class for common computations regarding adjoint dynamics and Newtonian gravity for Cartesian coordinates.
27   *
28   * @author Romain Serra
29   * @see CartesianAdjointEquationTerm
30   * @since 12.2
31   */
32  public abstract class AbstractCartesianAdjointNewtonianTerm extends AbstractCartesianAdjointGravitationalTerm {
33  
34      /** Minus three. */
35      private static final double MINUS_THREE = -3;
36  
37      /**
38       * Constructor.
39       * @param mu body gravitational parameter
40       */
41      protected AbstractCartesianAdjointNewtonianTerm(final double mu) {
42          super(mu);
43      }
44  
45      /**
46       * Computes the generic term of a Newtonian attraction (central or not).
47       * @param relativePosition relative position w.r.t. the body
48       * @param adjointVariables adjoint variables
49       * @return contribution to velocity adjoint derivatives
50       */
51      protected double[] getNewtonianVelocityAdjointContribution(final double[] relativePosition,
52                                                                 final double[] adjointVariables) {
53          final double[] contribution = new double[3];
54          final double x = relativePosition[0];
55          final double y = relativePosition[1];
56          final double z = relativePosition[2];
57          final double x2 = x * x;
58          final double y2 = y * y;
59          final double z2 = z * z;
60          final double r2 = x2 + y2 + z2;
61          final double r = FastMath.sqrt(r2);
62          final double factor = getMu() / (r2 * r2 * r);
63          final double xy = x * y;
64          final double xz = x * z;
65          final double yz = y * z;
66          final double pvx = adjointVariables[3];
67          final double pvy = adjointVariables[4];
68          final double pvz = adjointVariables[5];
69          contribution[0] = ((x2 * MINUS_THREE + r2) * pvx + xy * MINUS_THREE * pvy + xz * MINUS_THREE * pvz) * factor;
70          contribution[1] = ((y2 * MINUS_THREE + r2) * pvy + xy * MINUS_THREE * pvx + yz * MINUS_THREE * pvz) * factor;
71          contribution[2] = ((z2 * MINUS_THREE + r2) * pvz + yz * MINUS_THREE * pvy + xz * MINUS_THREE * pvx) * factor;
72          return contribution;
73      }
74  
75      /**
76       * Computes the generic term of a Newtonian attraction (central or not).
77       * @param relativePosition relative position w.r.t. the body
78       * @param adjointVariables adjoint variables
79       * @param <T> field type
80       * @return contribution to velocity adjoint derivatives
81       */
82      protected <T extends CalculusFieldElement<T>> T[] getFieldNewtonianVelocityAdjointContribution(final T[] relativePosition,
83                                                                                                      final T[] adjointVariables) {
84          final T[] contribution = MathArrays.buildArray(adjointVariables[0].getField(), 3);
85          final T x = relativePosition[0];
86          final T y = relativePosition[1];
87          final T z = relativePosition[2];
88          final T x2 = x.multiply(x);
89          final T y2 = y.multiply(y);
90          final T z2 = z.multiply(z);
91          final T r2 = x2.add(y2).add(z2);
92          final T r = r2.sqrt();
93          final T factor = (r2.multiply(r2).multiply(r)).reciprocal().multiply(getMu());
94          final T xy = x.multiply(y);
95          final T xz = x.multiply(z);
96          final T yz = y.multiply(z);
97          final T pvx = adjointVariables[3];
98          final T pvy = adjointVariables[4];
99          final T pvz = adjointVariables[5];
100         contribution[0] = ((x2.multiply(MINUS_THREE).add(r2)).multiply(pvx).
101                 add((xy.multiply(MINUS_THREE)).multiply(pvy)).
102                 add((xz.multiply(MINUS_THREE)).multiply(pvz))).multiply(factor);
103         contribution[1] = ((xy.multiply(MINUS_THREE)).multiply(pvx).
104                 add((y2.multiply(MINUS_THREE).add(r2)).multiply(pvy)).
105                 add((yz.multiply(MINUS_THREE)).multiply(pvz))).multiply(factor);
106         contribution[2] = ((xz.multiply(MINUS_THREE)).multiply(pvx).
107                 add((yz.multiply(MINUS_THREE)).multiply(pvy)).
108                 add((z2.multiply(MINUS_THREE).add(r2)).multiply(pvz))).multiply(factor);
109         return contribution;
110     }
111 
112     /**
113      * Compute the Newtonian acceleration.
114      * @param position position vector as array
115      * @return Newtonian acceleration vector
116      */
117     protected Vector3D getNewtonianAcceleration(final double[] position) {
118         final Vector3D positionVector = new Vector3D(position[0], position[1], position[2]);
119         final double squaredRadius = positionVector.getNormSq();
120         final double factor = -getMu() / (squaredRadius * FastMath.sqrt(squaredRadius));
121         return positionVector.scalarMultiply(factor);
122     }
123 
124     /**
125      * Compute the Newtonian acceleration.
126      * @param position position vector as array
127      * @param <T> field type
128      * @return Newtonian acceleration vector
129      */
130     protected <T extends CalculusFieldElement<T>> FieldVector3D<T> getFieldNewtonianAcceleration(final T[] position) {
131         final FieldVector3D<T> positionVector = new FieldVector3D<>(position[0], position[1], position[2]);
132         final T squaredRadius = positionVector.getNormSq();
133         final T factor = (squaredRadius.multiply(FastMath.sqrt(squaredRadius))).reciprocal().multiply(-getMu());
134         return positionVector.scalarMultiply(factor);
135     }
136 }