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.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-242 or echivalently, Danielson 2.7.3-(6) for Hansen Coefficients and
26 * Collins 4-245 or Danielson 3.1-(7) 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 HansenZonalLinear {
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 s coefficient. */
54 private final int s;
55
56 /**
57 * The offset used to identify the polynomial that corresponds to a negative
58 * value of n in the internal array that starts at 0.
59 */
60 private final int offset;
61
62 /** The number of slices needed to compute the coefficients. */
63 private final int numSlices;
64
65 /** 2<sup>s</sup>. */
66 private final double twots;
67
68 /** 2*s+1. */
69 private final int twosp1;
70
71 /** 2*s. */
72 private final int twos;
73
74 /** (2*s+1) / 2<sup>s</sup>. */
75 private final double twosp1otwots;
76
77 /**
78 * Constructor.
79 *
80 * @param nMax the maximum (absolute) value of n coefficient
81 * @param s s coefficient
82 */
83 public HansenZonalLinear(final int nMax, final int s) {
84
85 //Initialize fields
86 final int Nmin = -nMax - 1;
87 final int N0 = -(s + 2);
88 this.offset = nMax + 1;
89 this.s = s;
90 this.twots = FastMath.pow(2., s);
91 this.twos = 2 * s;
92 this.twosp1 = this.twos + 1;
93 this.twosp1otwots = (double) this.twosp1 / this.twots;
94
95 // prepare structures for stored data
96 final int size = nMax - s - 1;
97 mpvec = new PolynomialFunction[size][];
98 mpvecDeriv = new PolynomialFunction[size][];
99
100 this.numSlices = FastMath.max((int) FastMath.ceil(((double) size) / SLICE), 1);
101 hansenRoot = new double[numSlices][2];
102 hansenDerivRoot = new double[numSlices][2];
103
104 // Prepare the data base of associated polynomials
105 HansenUtilities.generateZonalPolynomials(N0, Nmin, offset, SLICE, s,
106 mpvec, mpvecDeriv);
107
108 }
109
110 /**
111 * Compute the roots for the Hansen coefficients and their derivatives.
112 *
113 * @param chi 1 / sqrt(1 - e²)
114 */
115 public void computeInitValues(final double chi) {
116 // compute the values for n=s and n=s+1
117 // See Danielson 2.7.3-(6a,b)
118 hansenRoot[0][0] = 0;
119 hansenRoot[0][1] = FastMath.pow(chi, this.twosp1) / this.twots;
120 hansenDerivRoot[0][0] = 0;
121 hansenDerivRoot[0][1] = this.twosp1otwots * FastMath.pow(chi, this.twos);
122
123 final int st = -s - 1;
124 for (int i = 1; i < numSlices; i++) {
125 for (int j = 0; j < 2; j++) {
126 // Get the required polynomials
127 final PolynomialFunction[] mv = mpvec[st - (i * SLICE) - j + offset];
128 final PolynomialFunction[] sv = mpvecDeriv[st - (i * SLICE) - j + offset];
129
130 //Compute the root derivatives
131 hansenDerivRoot[i][j] = mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
132 mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
133 (sv[1].value(chi) * hansenRoot[i - 1][1] +
134 sv[0].value(chi) * hansenRoot[i - 1][0]
135 ) / chi;
136 hansenRoot[i][j] = mv[1].value(chi) * hansenRoot[i - 1][1] +
137 mv[0].value(chi) * hansenRoot[i - 1][0];
138
139 }
140
141 }
142 }
143
144 /**
145 * Get the K₀<sup>-n-1,s</sup> coefficient value.
146 *
147 * <p> The s value is given in the class constructor
148 *
149 * @param mnm1 (-n-1) coefficient
150 * @param chi The value of χ
151 * @return K₀<sup>-n-1,s</sup>
152 */
153 public double getValue(final int mnm1, final double chi) {
154
155 //Compute n
156 final int n = -mnm1 - 1;
157
158 //Compute the potential slice
159 int sliceNo = (n - s) / SLICE;
160 if (sliceNo < numSlices) {
161 //Compute the index within the slice
162 final int indexInSlice = (n - s) % SLICE;
163
164 //Check if a root must be returned
165 if (indexInSlice <= 1) {
166 return hansenRoot[sliceNo][indexInSlice];
167 }
168 } else {
169 //the value was a potential root for a slice, but that slice was not required
170 //Decrease the slice number
171 sliceNo--;
172 }
173
174 // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
175 final PolynomialFunction[] v = mpvec[mnm1 + offset];
176 double ret = v[1].value(chi) * hansenRoot[sliceNo][1];
177 if (hansenRoot[sliceNo][0] != 0) {
178 ret += v[0].value(chi) * hansenRoot[sliceNo][0];
179 }
180 return ret;
181 }
182
183 /**
184 * Get the dK₀<sup>-n-1,s</sup> / dΧ coefficient derivative.
185 *
186 * <p> The s value is given in the class constructor.
187 *
188 * @param mnm1 (-n-1) coefficient
189 * @param chi The value of χ
190 * @return dK₀<sup>-n-1,s</sup> / dΧ
191 */
192 public double getDerivative(final int mnm1, final double chi) {
193
194 //Compute n
195 final int n = -mnm1 - 1;
196
197 //Compute the potential slice
198 int sliceNo = (n - s) / SLICE;
199 if (sliceNo < numSlices) {
200 //Compute the index within the slice
201 final int indexInSlice = (n - s) % SLICE;
202
203 //Check if a root must be returned
204 if (indexInSlice <= 1) {
205 return hansenDerivRoot[sliceNo][indexInSlice];
206 }
207 } else {
208 //the value was a potential root for a slice, but that slice was not required
209 //Decrease the slice number
210 sliceNo--;
211 }
212
213 // Danielson 3.1-(7c) and Petre's paper
214 final PolynomialFunction[] v = mpvec[mnm1 + offset];
215 double ret = v[1].value(chi) * hansenDerivRoot[sliceNo][1];
216 if (hansenDerivRoot[sliceNo][0] != 0) {
217 ret += v[0].value(chi) * hansenDerivRoot[sliceNo][0];
218 }
219
220 // Danielson 2.7.3-(6b)
221 final PolynomialFunction[] v1 = mpvecDeriv[mnm1 + offset];
222 double hret = v1[1].value(chi) * hansenRoot[sliceNo][1];
223 if (hansenRoot[sliceNo][0] != 0) {
224 hret += v1[0].value(chi) * hansenRoot[sliceNo][0];
225 }
226 ret += hret / chi;
227
228 return ret;
229
230 }
231
232 }