1   /* Copyright 2002-2025 CS GROUP
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.propagation.semianalytical.dsst.utilities.hansen;
18  
19  import java.lang.reflect.Array;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.Field;
23  import org.hipparchus.analysis.differentiation.FieldGradient;
24  import org.hipparchus.analysis.polynomials.PolynomialFunction;
25  import org.hipparchus.util.FastMath;
26  import org.hipparchus.util.MathArrays;
27  import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;
28  
29  /**
30   * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
31   * <p>
32   * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
33   * Collins 4-240 for derivatives. The recursions are transformed into
34   * composition of linear transformations to obtain the associated polynomials
35   * for coefficients and their derivatives - see Petre's paper
36   *
37   * @author Petre Bazavan
38   * @author Lucian Barbulescu
39   * @author Bryan Cazabonne
40   * @param <T> type of the field elements
41   */
42  public class FieldHansenTesseralLinear <T extends CalculusFieldElement<T>> {
43  
44      /** The number of coefficients that will be computed with a set of roots. */
45      private static final int SLICE = 10;
46  
47      /**
48       * The first vector of polynomials associated to Hansen coefficients and
49       * derivatives.
50       */
51      private PolynomialFunction[][] mpvec;
52  
53      /** The second vector of polynomials associated only to derivatives. */
54      private PolynomialFunction[][] mpvecDeriv;
55  
56      /** The Hansen coefficients used as roots. */
57      private final T[][] hansenRoot;
58  
59      /** The derivatives of the Hansen coefficients used as roots. */
60      private final T[][] hansenDerivRoot;
61  
62      /** The minimum value for the order. */
63      private final int Nmin;
64  
65      /** The index of the initial condition, Petre's paper. */
66      private final int N0;
67  
68      /** The number of slices needed to compute the coefficients. */
69      private final int numSlices;
70  
71      /**
72       * The offset used to identify the polynomial that corresponds to a negative.
73       * value of n in the internal array that starts at 0
74       */
75      private final int offset;
76  
77      /** The objects used to calculate initial data by means of Newcomb operators. */
78      private final FieldHansenCoefficientsBySeries<T>[] hansenInit;
79  
80      /**
81       * Constructor.
82       *
83       * @param nMax the maximum (absolute) value of n parameter
84       * @param s s parameter
85       * @param j j parameter
86       * @param n0 the minimum (absolute) value of n
87       * @param maxHansen maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
88       * @param field field used by default
89       */
90      @SuppressWarnings("unchecked")
91      public FieldHansenTesseralLinear(final int nMax, final int s, final int j, final int n0,
92                                       final int maxHansen, final Field<T> field) {
93          //Initialize the fields
94          this.offset = nMax + 1;
95          this.Nmin = -nMax - 1;
96          this.N0 = -n0 - 4;
97  
98          final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
99          //Ensure that only the needed terms are computed
100         this.hansenInit = (FieldHansenCoefficientsBySeries<T>[]) Array.newInstance(FieldHansenCoefficientsBySeries.class, maxRoots);
101         for (int i = 0; i < maxRoots; i++) {
102             this.hansenInit[i] = new FieldHansenCoefficientsBySeries<>(N0 - i + 3, s, j, maxHansen, field);
103         }
104 
105         // The first 4 values are computed with series. No linear combination is needed.
106         final int size = N0 - Nmin;
107         this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
108         hansenRoot = MathArrays.buildArray(field, numSlices, 4);
109         hansenDerivRoot = MathArrays.buildArray(field, numSlices, 4);
110         if (size > 0) {
111             mpvec = new PolynomialFunction[size][];
112             mpvecDeriv = new PolynomialFunction[size][];
113 
114             // Prepare the database of the associated polynomials
115             HansenUtilities.generateTesseralPolynomials(N0, Nmin, offset, SLICE, j, s,
116                                                         mpvec, mpvecDeriv);
117         }
118 
119     }
120 
121     /**
122      * Compute the values for the first four coefficients and their derivatives by means of series.
123      *
124      * @param e2 e²
125      * @param chi &Chi;
126      * @param chi2 &Chi;²
127      */
128     public void computeInitValues(final T e2, final T chi, final T chi2) {
129         // compute the values for n, n+1, n+2 and n+3 by series
130         // See Danielson 2.7.3-(10)
131         //Ensure that only the needed terms are computed
132         final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
133         for (int i = 0; i < maxRoots; i++) {
134             final FieldGradient<T> hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
135             this.hansenRoot[0][i] = hansenKernel.getValue();
136             this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
137         }
138 
139         for (int i = 1; i < numSlices; i++) {
140             for (int k = 0; k < 4; k++) {
141                 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
142                 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];
143 
144                 hansenDerivRoot[i][k] = mv[3].value(chi).multiply(hansenDerivRoot[i - 1][3]).
145                                         add(mv[2].value(chi).multiply(hansenDerivRoot[i - 1][2])).
146                                         add(mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1])).
147                                         add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
148                                         add(sv[3].value(chi).multiply(hansenRoot[i - 1][3])).
149                                         add(sv[2].value(chi).multiply(hansenRoot[i - 1][2])).
150                                         add(sv[1].value(chi).multiply(hansenRoot[i - 1][1])).
151                                         add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]));
152 
153                 hansenRoot[i][k] =  mv[3].value(chi).multiply(hansenRoot[i - 1][3]).
154                                     add(mv[2].value(chi).multiply(hansenRoot[i - 1][2])).
155                                     add(mv[1].value(chi).multiply(hansenRoot[i - 1][1])).
156                                     add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));
157             }
158         }
159     }
160 
161     /**
162      * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
163      *
164      * @param mnm1 -n-1
165      * @param chi χ
166      * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
167      */
168     public T getValue(final int mnm1, final T chi) {
169         //Compute n
170         final int n = -mnm1 - 1;
171 
172         //Compute the potential slice
173         int sliceNo = (n + N0 + 4) / SLICE;
174         if (sliceNo < numSlices) {
175             //Compute the index within the slice
176             final int indexInSlice = (n + N0 + 4) % SLICE;
177 
178             //Check if a root must be returned
179             if (indexInSlice <= 3) {
180                 return hansenRoot[sliceNo][indexInSlice];
181             }
182         } else {
183             //the value was a potential root for a slice, but that slice was not required
184             //Decrease the slice number
185             sliceNo--;
186         }
187 
188         // Computes the coefficient by linear transformation
189         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
190         final PolynomialFunction[] v = mpvec[mnm1 + offset];
191         return v[3].value(chi).multiply(hansenRoot[sliceNo][3]).
192                add(v[2].value(chi).multiply(hansenRoot[sliceNo][2])).
193                add(v[1].value(chi).multiply(hansenRoot[sliceNo][1])).
194                add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));
195 
196     }
197 
198     /**
199      * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
200      *
201      * @param mnm1 -n-1
202      * @param chi χ
203      * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
204      */
205     public T getDerivative(final int mnm1, final T chi) {
206 
207         //Compute n
208         final int n = -mnm1 - 1;
209 
210         //Compute the potential slice
211         int sliceNo = (n + N0 + 4) / SLICE;
212         if (sliceNo < numSlices) {
213             //Compute the index within the slice
214             final int indexInSlice = (n + N0 + 4) % SLICE;
215 
216             //Check if a root must be returned
217             if (indexInSlice <= 3) {
218                 return hansenDerivRoot[sliceNo][indexInSlice];
219             }
220         } else {
221             //the value was a potential root for a slice, but that slice was not required
222             //Decrease the slice number
223             sliceNo--;
224         }
225 
226         // Computes the coefficient by linear transformation
227         // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
228         final PolynomialFunction[] v = mpvec[mnm1 + this.offset];
229         final PolynomialFunction[] vv = mpvecDeriv[mnm1 + this.offset];
230 
231         return v[3].value(chi).multiply(hansenDerivRoot[sliceNo][3]).
232                add(v[2].value(chi).multiply(hansenDerivRoot[sliceNo][2])).
233                add(v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1])).
234                add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0])).
235                add(vv[3].value(chi).multiply(hansenRoot[sliceNo][3])).
236                add(vv[2].value(chi).multiply(hansenRoot[sliceNo][2])).
237                add( vv[1].value(chi).multiply(hansenRoot[sliceNo][1])).
238                add(vv[0].value(chi).multiply(hansenRoot[sliceNo][0]));
239 
240     }
241 
242     /**
243      * Compute a Hansen coefficient with series.
244      * <p>
245      * This class implements the computation of the Hansen kernels
246      * through a power series in e² and that is using
247      * modified Newcomb operators. The reference formulae can be found
248      * in Danielson 2.7.3-10 and 3.3-5
249      * </p>
250      */
251     private static class FieldHansenCoefficientsBySeries <T extends CalculusFieldElement<T>> {
252 
253         /** -n-1 coefficient. */
254         private final int mnm1;
255 
256         /** s coefficient. */
257         private final int s;
258 
259         /** j coefficient. */
260         private final int j;
261 
262         /** Max power in e² for the Newcomb's series expansion. */
263         private final int maxNewcomb;
264 
265         /** Polynomial representing the serie. */
266         private final PolynomialFunction polynomial;
267 
268         /**
269          * Class constructor.
270          *
271          * @param mnm1 -n-1 value
272          * @param s s value
273          * @param j j value
274          * @param maxHansen max power of e² in series expansion
275          * @param field field for the function parameters and value
276          */
277         FieldHansenCoefficientsBySeries(final int mnm1, final int s,
278                                           final int j, final int maxHansen, final Field<T> field) {
279             this.mnm1 = mnm1;
280             this.s = s;
281             this.j = j;
282             this.maxNewcomb = maxHansen;
283             this.polynomial = generatePolynomial();
284         }
285 
286         /** Computes the value of Hansen kernel and its derivative at e².
287          * <p>
288          * The formulae applied are described in Danielson 2.7.3-10 and
289          * 3.3-5
290          * </p>
291          * @param e2 e²
292          * @param chi &Chi;
293          * @param chi2 &Chi;²
294          * @return the value of the Hansen coefficient and its derivative for e²
295          */
296         private FieldGradient<T> getValueGradient(final T e2, final T chi, final T chi2) {
297 
298             final T zero = e2.getField().getZero();
299             //Estimation of the serie expansion at e2
300             final FieldGradient<T> serie = polynomial.value(FieldGradient.variable(1, 0, e2));
301 
302             final T value      =  FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getValue()).divide(chi);
303             final T coef       = zero.subtract(mnm1 + 1.5);
304             final T derivative = coef.multiply(chi2).multiply(value).
305                             add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(0)).divide(chi));
306             return new FieldGradient<>(value, derivative);
307         }
308 
309         /** Generate the serie expansion in e².
310          * <p>
311          * Generate the series expansion in e² used in the formulation
312          * of the Hansen kernel (see Danielson 2.7.3-10):
313          * &Sigma; Y<sup>ns</sup><sub>α+a,α+b</sub>
314          * *e<sup>2α</sup>
315          * </p>
316          * @return polynomial representing the power serie expansion
317          */
318         private PolynomialFunction generatePolynomial() {
319             // Initialization
320             final int aHT = FastMath.max(j - s, 0);
321             final int bHT = FastMath.max(s - j, 0);
322 
323             final double[] coefficients = new double[maxNewcomb + 1];
324 
325             //Loop for getting the Newcomb operators
326             for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
327                 coefficients[alphaHT] =
328                         NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
329             }
330 
331             //Creation of the polynomial
332             return new PolynomialFunction(coefficients);
333         }
334     }
335 
336 }