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 org.hipparchus.Field;
20 import org.hipparchus.CalculusFieldElement;
21 import org.hipparchus.analysis.polynomials.PolynomialFunction;
22 import org.hipparchus.util.FastMath;
23 import org.hipparchus.util.MathArrays;
24
25 /**
26 * Hansen coefficients K(t,n,s) for t=0 and n < 0.
27 * <p>
28 *Implements Collins 4-242 or echivalently, Danielson 2.7.3-(6) for Hansen Coefficients and
29 * Collins 4-245 or Danielson 3.1-(7) for derivatives. The recursions are transformed into
30 * composition of linear transformations to obtain the associated polynomials
31 * for coefficients and their derivatives - see Petre's paper
32 *
33 * @author Petre Bazavan
34 * @author Lucian Barbulescu
35 * @author Bryan Cazabonne
36 * @param <T> type of the field elements
37 */
38 public class FieldHansenZonalLinear <T extends CalculusFieldElement<T>> {
39
40 /** The number of coefficients that will be computed with a set of roots. */
41 private static final int SLICE = 10;
42
43 /**
44 * The first vector of polynomials associated to Hansen coefficients and
45 * derivatives.
46 */
47 private final PolynomialFunction[][] mpvec;
48
49 /** The second vector of polynomials associated only to derivatives. */
50 private final PolynomialFunction[][] mpvecDeriv;
51
52 /** The Hansen coefficients used as roots. */
53 private final T[][] hansenRoot;
54
55 /** The derivatives of the Hansen coefficients used as roots. */
56 private final T[][] hansenDerivRoot;
57
58 /** The s coefficient. */
59 private final int s;
60
61 /**
62 * The offset used to identify the polynomial that corresponds to a negative
63 * value of n in the internal array that starts at 0.
64 */
65 private final int offset;
66
67 /** The number of slices needed to compute the coefficients. */
68 private final int numSlices;
69
70 /** 2<sup>s</sup>. */
71 private final double twots;
72
73 /** 2*s+1. */
74 private final int twosp1;
75
76 /** 2*s. */
77 private final int twos;
78
79 /** (2*s+1) / 2<sup>s</sup>. */
80 private final double twosp1otwots;
81
82 /**
83 * Constructor.
84 *
85 * @param nMax the maximum (absolute) value of n coefficient
86 * @param s s coefficient
87 * @param field field used by default
88 */
89 public FieldHansenZonalLinear(final int nMax, final int s, final Field<T> field) {
90 //Initialize fields
91 final int Nmin = -nMax - 1;
92 final int N0 = -(s + 2);
93 this.offset = nMax + 1;
94 this.s = s;
95 this.twots = FastMath.pow(2., s);
96 this.twos = 2 * s;
97 this.twosp1 = this.twos + 1;
98 this.twosp1otwots = (double) this.twosp1 / this.twots;
99
100 // prepare structures for stored data
101 final int size = nMax - s - 1;
102 mpvec = new PolynomialFunction[size][];
103 mpvecDeriv = new PolynomialFunction[size][];
104
105 this.numSlices = FastMath.max((int) FastMath.ceil(((double) size) / SLICE), 1);
106 hansenRoot = MathArrays.buildArray(field, numSlices, 2);
107 hansenDerivRoot = MathArrays.buildArray(field, numSlices, 2);
108
109 // Prepare the data base of associated polynomials
110 HansenUtilities.generateZonalPolynomials(N0, Nmin, offset, SLICE, s,
111 mpvec, mpvecDeriv);
112
113 }
114
115 /**
116 * Compute the roots for the Hansen coefficients and their derivatives.
117 *
118 * @param chi 1 / sqrt(1 - e²)
119 */
120 public void computeInitValues(final T chi) {
121 final Field<T> field = chi.getField();
122 final T zero = field.getZero();
123 // compute the values for n=s and n=s+1
124 // See Danielson 2.7.3-(6a,b)
125 hansenRoot[0][0] = zero;
126 hansenRoot[0][1] = FastMath.pow(chi, this.twosp1).divide(this.twots);
127 hansenDerivRoot[0][0] = zero;
128 hansenDerivRoot[0][1] = FastMath.pow(chi, this.twos).multiply(this.twosp1otwots);
129
130 final int st = -s - 1;
131 for (int i = 1; i < numSlices; i++) {
132 for (int j = 0; j < 2; j++) {
133 // Get the required polynomials
134 final PolynomialFunction[] mv = mpvec[st - (i * SLICE) - j + offset];
135 final PolynomialFunction[] sv = mpvecDeriv[st - (i * SLICE) - j + offset];
136
137 //Compute the root derivatives
138 hansenDerivRoot[i][j] = mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1]).
139 add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
140 add((sv[1].value(chi).multiply(hansenRoot[i - 1][1]).
141 add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]))
142 ).divide(chi));
143 hansenRoot[i][j] = mv[1].value(chi).multiply(hansenRoot[i - 1][1]).
144 add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));
145
146 }
147
148 }
149 }
150
151 /**
152 * Get the K₀<sup>-n-1,s</sup> coefficient value.
153 *
154 * <p> The s value is given in the class constructor
155 *
156 * @param mnm1 (-n-1) coefficient
157 * @param chi The value of χ
158 * @return K₀<sup>-n-1,s</sup>
159 */
160 public T getValue(final int mnm1, final T chi) {
161 //Compute n
162 final int n = -mnm1 - 1;
163
164 //Compute the potential slice
165 int sliceNo = (n - s) / SLICE;
166 if (sliceNo < numSlices) {
167 //Compute the index within the slice
168 final int indexInSlice = (n - s) % SLICE;
169
170 //Check if a root must be returned
171 if (indexInSlice <= 1) {
172 return hansenRoot[sliceNo][indexInSlice];
173 }
174 } else {
175 //the value was a potential root for a slice, but that slice was not required
176 //Decrease the slice number
177 sliceNo--;
178 }
179
180 // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
181 final PolynomialFunction[] v = mpvec[mnm1 + offset];
182 T ret = v[1].value(chi).multiply(hansenRoot[sliceNo][1]);
183 if (hansenRoot[sliceNo][0].getReal() != 0) {
184 ret = ret.add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));
185 }
186 return ret;
187 }
188
189 /**
190 * Get the dK₀<sup>-n-1,s</sup> / dΧ coefficient derivative.
191 *
192 * <p> The s value is given in the class constructor.
193 *
194 * @param mnm1 (-n-1) coefficient
195 * @param chi The value of χ
196 * @return dK₀<sup>-n-1,s</sup> / dΧ
197 */
198 public T getDerivative(final int mnm1, final T chi) {
199 //Compute n
200 final int n = -mnm1 - 1;
201
202 //Compute the potential slice
203 int sliceNo = (n - s) / SLICE;
204 if (sliceNo < numSlices) {
205 //Compute the index within the slice
206 final int indexInSlice = (n - s) % SLICE;
207
208 //Check if a root must be returned
209 if (indexInSlice <= 1) {
210 return hansenDerivRoot[sliceNo][indexInSlice];
211 }
212 } else {
213 //the value was a potential root for a slice, but that slice was not required
214 //Decrease the slice number
215 sliceNo--;
216 }
217
218 // Danielson 3.1-(7c) and Petre's paper
219 final PolynomialFunction[] v = mpvec[mnm1 + offset];
220 T ret = v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1]);
221 if (hansenDerivRoot[sliceNo][0].getReal() != 0) {
222 ret = ret.add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0]));
223 }
224
225 // Danielson 2.7.3-(6b)
226 final PolynomialFunction[] v1 = mpvecDeriv[mnm1 + offset];
227 T hret = v1[1].value(chi).multiply(hansenRoot[sliceNo][1]);
228 if (hansenRoot[sliceNo][0].getReal() != 0) {
229 hret = hret.add(v1[0].value(chi).multiply(hansenRoot[sliceNo][0]));
230 }
231 ret = ret.add(hret.divide(chi));
232
233 return ret;
234
235 }
236
237 }