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