1 /* Copyright 2002-2022 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.forces;
18
19 import java.lang.reflect.Array;
20 import java.util.ArrayList;
21 import java.util.Arrays;
22 import java.util.Collections;
23 import java.util.HashMap;
24 import java.util.List;
25 import java.util.Map;
26 import java.util.Set;
27 import java.util.TreeMap;
28
29 import org.hipparchus.Field;
30 import org.hipparchus.CalculusFieldElement;
31 import org.hipparchus.analysis.differentiation.FieldGradient;
32 import org.hipparchus.analysis.differentiation.Gradient;
33 import org.hipparchus.util.CombinatoricsUtils;
34 import org.hipparchus.util.FastMath;
35 import org.hipparchus.util.FieldSinCos;
36 import org.hipparchus.util.MathArrays;
37 import org.hipparchus.util.SinCos;
38 import org.orekit.attitudes.AttitudeProvider;
39 import org.orekit.bodies.CelestialBodies;
40 import org.orekit.bodies.CelestialBody;
41 import org.orekit.orbits.FieldOrbit;
42 import org.orekit.orbits.Orbit;
43 import org.orekit.propagation.FieldSpacecraftState;
44 import org.orekit.propagation.PropagationType;
45 import org.orekit.propagation.SpacecraftState;
46 import org.orekit.propagation.events.EventDetector;
47 import org.orekit.propagation.events.FieldEventDetector;
48 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
49 import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
50 import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory;
51 import org.orekit.propagation.semianalytical.dsst.utilities.CoefficientsFactory.NSKey;
52 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
53 import org.orekit.propagation.semianalytical.dsst.utilities.FieldCjSjCoefficient;
54 import org.orekit.propagation.semianalytical.dsst.utilities.FieldShortPeriodicsInterpolatedCoefficient;
55 import org.orekit.propagation.semianalytical.dsst.utilities.JacobiPolynomials;
56 import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
57 import org.orekit.propagation.semianalytical.dsst.utilities.hansen.FieldHansenThirdBodyLinear;
58 import org.orekit.propagation.semianalytical.dsst.utilities.hansen.HansenThirdBodyLinear;
59 import org.orekit.time.AbsoluteDate;
60 import org.orekit.time.FieldAbsoluteDate;
61 import org.orekit.utils.FieldTimeSpanMap;
62 import org.orekit.utils.ParameterDriver;
63 import org.orekit.utils.TimeSpanMap;
64
65 /** Third body attraction perturbation to the
66 * {@link org.orekit.propagation.semianalytical.dsst.DSSTPropagator DSSTPropagator}.
67 *
68 * @author Romain Di Costanzo
69 * @author Pascal Parraud
70 * @author Bryan Cazabonne (field translation)
71 */
72 public class DSSTThirdBody implements DSSTForceModel {
73
74 /** Name of the prefix for short period coefficients keys. */
75 public static final String SHORT_PERIOD_PREFIX = "DSST-3rd-body-";
76
77 /** Name of the single parameter of this model: the attraction coefficient. */
78 public static final String ATTRACTION_COEFFICIENT = " attraction coefficient";
79
80 /** Central attraction scaling factor.
81 * <p>
82 * We use a power of 2 to avoid numeric noise introduction
83 * in the multiplications/divisions sequences.
84 * </p>
85 */
86 private static final double MU_SCALE = FastMath.scalb(1.0, 32);
87
88 /** Retrograde factor I.
89 * <p>
90 * DSST model needs equinoctial orbit as internal representation.
91 * Classical equinoctial elements have discontinuities when inclination
92 * is close to zero. In this representation, I = +1. <br>
93 * To avoid this discontinuity, another representation exists and equinoctial
94 * elements can be expressed in a different way, called "retrograde" orbit.
95 * This implies I = -1. <br>
96 * As Orekit doesn't implement the retrograde orbit, I is always set to +1.
97 * But for the sake of consistency with the theory, the retrograde factor
98 * has been kept in the formulas.
99 * </p>
100 */
101 private static final int I = 1;
102
103 /** Number of points for interpolation. */
104 private static final int INTERPOLATION_POINTS = 3;
105
106 /** Maximum power for eccentricity used in short periodic computation. */
107 private static final int MAX_ECCPOWER_SP = 4;
108
109 /** Max power for summation. */
110 private static final int MAX_POWER = 22;
111
112 /** V<sub>ns</sub> coefficients. */
113 private final TreeMap<NSKey, Double> Vns;
114
115 /** Max frequency of F. */
116 private int maxFreqF;
117
118 /** Max frequency of F for field initialization. */
119 private int maxFieldFreqF;
120
121 /** The 3rd body to consider. */
122 private final CelestialBody body;
123
124 /** Short period terms. */
125 private ThirdBodyShortPeriodicCoefficients shortPeriods;
126
127 /** Short period terms. */
128 private Map<Field<?>, FieldThirdBodyShortPeriodicCoefficients<?>> fieldShortPeriods;
129
130 /** Drivers for third body attraction coefficient and gravitational parameter. */
131 private final List<ParameterDriver> parameterDrivers;
132
133 /** Hansen objects. */
134 private HansenObjects hansen;
135
136 /** Hansen objects for field elements. */
137 private Map<Field<?>, FieldHansenObjects<?>> fieldHansen;
138
139 /** Complete constructor.
140 * @param body the 3rd body to consider
141 * @param mu central attraction coefficient
142 * @see CelestialBodies
143 */
144 public DSSTThirdBody(final CelestialBody body, final double mu) {
145 parameterDrivers = new ArrayList<>(2);
146 parameterDrivers.add(new ParameterDriver(body.getName() + DSSTThirdBody.ATTRACTION_COEFFICIENT,
147 body.getGM(), MU_SCALE,
148 0.0, Double.POSITIVE_INFINITY));
149 parameterDrivers.add(new ParameterDriver(DSSTNewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT,
150 mu, MU_SCALE,
151 0.0, Double.POSITIVE_INFINITY));
152
153 this.body = body;
154 this.Vns = CoefficientsFactory.computeVns(MAX_POWER);
155
156 fieldShortPeriods = new HashMap<>();
157 fieldHansen = new HashMap<>();
158 }
159
160 /** Get third body.
161 * @return third body
162 */
163 public CelestialBody getBody() {
164 return body;
165 }
166
167 /** Computes the highest power of the eccentricity and the highest power
168 * of a/R3 to appear in the truncated analytical power series expansion.
169 * <p>
170 * This method computes the upper value for the 3rd body potential and
171 * determines the maximal powers for the eccentricity and a/R3 producing
172 * potential terms bigger than a defined tolerance.
173 * </p>
174 * @param auxiliaryElements auxiliary elements related to the current orbit
175 * @param type type of the elements used during the propagation
176 * @param parameters values of the force model parameters
177 */
178 @Override
179 public List<ShortPeriodTerms> initializeShortPeriodTerms(final AuxiliaryElements auxiliaryElements,
180 final PropagationType type,
181 final double[] parameters) {
182
183 // Initializes specific parameters.
184 final DSSTThirdBodyContext context = initializeStep(auxiliaryElements, parameters);
185
186 maxFreqF = context.getMaxFreqF();
187
188 hansen = new HansenObjects();
189
190 final int jMax = maxFreqF;
191 shortPeriods = new ThirdBodyShortPeriodicCoefficients(jMax, INTERPOLATION_POINTS,
192 maxFreqF, body.getName(),
193 new TimeSpanMap<Slot>(new Slot(jMax, INTERPOLATION_POINTS)));
194
195 final List<ShortPeriodTerms> list = new ArrayList<ShortPeriodTerms>();
196 list.add(shortPeriods);
197 return list;
198
199 }
200
201 /** {@inheritDoc} */
202 @Override
203 public <T extends CalculusFieldElement<T>> List<FieldShortPeriodTerms<T>> initializeShortPeriodTerms(final FieldAuxiliaryElements<T> auxiliaryElements,
204 final PropagationType type,
205 final T[] parameters) {
206
207 // Field used by default
208 final Field<T> field = auxiliaryElements.getDate().getField();
209
210 // Initializes specific parameters.
211 final FieldDSSTThirdBodyContext<T> context = initializeStep(auxiliaryElements, parameters);
212
213 maxFieldFreqF = context.getMaxFreqF();
214
215 fieldHansen.put(field, new FieldHansenObjects<>(field));
216
217 final int jMax = maxFieldFreqF;
218 final FieldThirdBodyShortPeriodicCoefficients<T> ftbspc =
219 new FieldThirdBodyShortPeriodicCoefficients<>(jMax, INTERPOLATION_POINTS,
220 maxFieldFreqF, body.getName(),
221 new FieldTimeSpanMap<>(new FieldSlot<>(jMax,
222 INTERPOLATION_POINTS),
223 field));
224 fieldShortPeriods.put(field, ftbspc);
225 return Collections.singletonList(ftbspc);
226 }
227
228 /** Performs initialization at each integration step for the current force model.
229 * <p>
230 * This method aims at being called before mean elements rates computation.
231 * </p>
232 * @param auxiliaryElements auxiliary elements related to the current orbit
233 * @param parameters values of the force model parameters
234 * @return new force model context
235 */
236 private DSSTThirdBodyContext initializeStep(final AuxiliaryElements auxiliaryElements, final double[] parameters) {
237 return new DSSTThirdBodyContext(auxiliaryElements, body, parameters);
238 }
239
240 /** Performs initialization at each integration step for the current force model.
241 * <p>
242 * This method aims at being called before mean elements rates computation.
243 * </p>
244 * @param <T> type of the elements
245 * @param auxiliaryElements auxiliary elements related to the current orbit
246 * @param parameters values of the force model parameters
247 * @return new force model context
248 */
249 private <T extends CalculusFieldElement<T>> FieldDSSTThirdBodyContext<T> initializeStep(final FieldAuxiliaryElements<T> auxiliaryElements,
250 final T[] parameters) {
251 return new FieldDSSTThirdBodyContext<>(auxiliaryElements, body, parameters);
252 }
253
254 /** {@inheritDoc} */
255 @Override
256 public double[] getMeanElementRate(final SpacecraftState currentState,
257 final AuxiliaryElements auxiliaryElements, final double[] parameters) {
258
259 // Container for attributes
260 final DSSTThirdBodyContext context = initializeStep(auxiliaryElements, parameters);
261 // Access to potential U derivatives
262 final UAnddU udu = new UAnddU(context, hansen);
263
264 // Compute cross derivatives [Eq. 2.2-(8)]
265 // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
266 final double UAlphaGamma = context.getAlpha() * udu.getdUdGa() - context.getGamma() * udu.getdUdAl();
267 // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
268 final double UBetaGamma = context.getBeta() * udu.getdUdGa() - context.getGamma() * udu.getdUdBe();
269 // Common factor
270 final double pUAGmIqUBGoAB = (auxiliaryElements.getP() * UAlphaGamma - I * auxiliaryElements.getQ() * UBetaGamma) * context.getOoAB();
271
272 // Compute mean elements rates [Eq. 3.1-(1)]
273 final double da = 0.;
274 final double dh = context.getBoA() * udu.getdUdk() + auxiliaryElements.getK() * pUAGmIqUBGoAB;
275 final double dk = -context.getBoA() * udu.getdUdh() - auxiliaryElements.getH() * pUAGmIqUBGoAB;
276 final double dp = context.getMCo2AB() * UBetaGamma;
277 final double dq = context.getMCo2AB() * UAlphaGamma * I;
278 final double dM = context.getM2aoA() * udu.getdUda() + context.getBoABpo() * (auxiliaryElements.getH() * udu.getdUdh() + auxiliaryElements.getK() * udu.getdUdk()) + pUAGmIqUBGoAB;
279
280 return new double[] {da, dk, dh, dq, dp, dM};
281
282 }
283
284 /** {@inheritDoc} */
285 @Override
286 public <T extends CalculusFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> currentState,
287 final FieldAuxiliaryElements<T> auxiliaryElements,
288 final T[] parameters) {
289
290 // Parameters for array building
291 final Field<T> field = currentState.getDate().getField();
292 final T zero = field.getZero();
293
294 // Container for attributes
295 final FieldDSSTThirdBodyContext<T> context = initializeStep(auxiliaryElements, parameters);
296
297 @SuppressWarnings("unchecked")
298 final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
299
300 // Access to potential U derivatives
301 final FieldUAnddU<T> udu = new FieldUAnddU<>(context, fho);
302
303 // Compute cross derivatives [Eq. 2.2-(8)]
304 // U(alpha,gamma) = alpha * dU/dgamma - gamma * dU/dalpha
305 final T UAlphaGamma = udu.getdUdGa().multiply(context.getAlpha()).subtract(udu.getdUdAl().multiply(context.getGamma()));
306 // U(beta,gamma) = beta * dU/dgamma - gamma * dU/dbeta
307 final T UBetaGamma = udu.getdUdGa().multiply(context.getBeta()).subtract(udu.getdUdBe().multiply(context.getGamma()));
308 // Common factor
309 final T pUAGmIqUBGoAB = (UAlphaGamma.multiply(auxiliaryElements.getP()).subtract(UBetaGamma.multiply(auxiliaryElements.getQ()).multiply(I))).multiply(context.getOoAB());
310
311 // Compute mean elements rates [Eq. 3.1-(1)]
312 final T da = zero;
313 final T dh = udu.getdUdk().multiply(context.getBoA()).add(pUAGmIqUBGoAB.multiply(auxiliaryElements.getK()));
314 final T dk = ((udu.getdUdh().multiply(context.getBoA())).negate()).subtract(pUAGmIqUBGoAB.multiply(auxiliaryElements.getH()));
315 final T dp = UBetaGamma.multiply(context.getMCo2AB());
316 final T dq = UAlphaGamma.multiply(I).multiply(context.getMCo2AB());
317 final T dM = pUAGmIqUBGoAB.add(udu.getdUda().multiply(context.getM2aoA())).add((udu.getdUdh().multiply(auxiliaryElements.getH()).add(udu.getdUdk().multiply(auxiliaryElements.getK()))).multiply(context.getBoABpo()));
318
319 final T[] elements = MathArrays.buildArray(field, 6);
320 elements[0] = da;
321 elements[1] = dk;
322 elements[2] = dh;
323 elements[3] = dq;
324 elements[4] = dp;
325 elements[5] = dM;
326
327 return elements;
328
329 }
330
331 /** {@inheritDoc} */
332 @Override
333 public void updateShortPeriodTerms(final double[] parameters, final SpacecraftState... meanStates) {
334
335 final Slot slot = shortPeriods.createSlot(meanStates);
336
337 for (final SpacecraftState meanState : meanStates) {
338
339 // Auxiliary elements related to the current orbit
340 final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(meanState.getOrbit(), I);
341
342 // Container of attributes
343 final DSSTThirdBodyContext context = initializeStep(auxiliaryElements, parameters);
344
345 final GeneratingFunctionCoefficients gfCoefs =
346 new GeneratingFunctionCoefficients(context.getMaxAR3Pow(), MAX_ECCPOWER_SP, context.getMaxAR3Pow() + 1, context, hansen);
347
348 //Compute additional quantities
349 // 2 * a / An
350 final double ax2oAn = -context.getM2aoA() / context.getMeanMotion();
351 // B / An
352 final double BoAn = context.getBoA() / context.getMeanMotion();
353 // 1 / ABn
354 final double ooABn = context.getOoAB() / context.getMeanMotion();
355 // C / 2ABn
356 final double Co2ABn = -context.getMCo2AB() / context.getMeanMotion();
357 // B / (A * (1 + B) * n)
358 final double BoABpon = context.getBoABpo() / context.getMeanMotion();
359 // -3 / n²a² = -3 / nA
360 final double m3onA = -3 / (context.getA() * context.getMeanMotion());
361
362 //Compute the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients.
363 for (int j = 1; j < slot.cij.length; j++) {
364 // First compute the C<sub>i</sub><sup>j</sup> coefficients
365 final double[] currentCij = new double[6];
366
367 // Compute the cross derivatives operator :
368 final double SAlphaGammaCj = context.getAlpha() * gfCoefs.getdSdgammaCj(j) - context.getGamma() * gfCoefs.getdSdalphaCj(j);
369 final double SAlphaBetaCj = context.getAlpha() * gfCoefs.getdSdbetaCj(j) - context.getBeta() * gfCoefs.getdSdalphaCj(j);
370 final double SBetaGammaCj = context.getBeta() * gfCoefs.getdSdgammaCj(j) - context.getGamma() * gfCoefs.getdSdbetaCj(j);
371 final double ShkCj = auxiliaryElements.getH() * gfCoefs.getdSdkCj(j) - auxiliaryElements.getK() * gfCoefs.getdSdhCj(j);
372 final double pSagmIqSbgoABnCj = (auxiliaryElements.getP() * SAlphaGammaCj - I * auxiliaryElements.getQ() * SBetaGammaCj) * ooABn;
373 final double ShkmSabmdSdlCj = ShkCj - SAlphaBetaCj - gfCoefs.getdSdlambdaCj(j);
374
375 currentCij[0] = ax2oAn * gfCoefs.getdSdlambdaCj(j);
376 currentCij[1] = -(BoAn * gfCoefs.getdSdhCj(j) + auxiliaryElements.getH() * pSagmIqSbgoABnCj + auxiliaryElements.getK() * BoABpon * gfCoefs.getdSdlambdaCj(j));
377 currentCij[2] = BoAn * gfCoefs.getdSdkCj(j) + auxiliaryElements.getK() * pSagmIqSbgoABnCj - auxiliaryElements.getH() * BoABpon * gfCoefs.getdSdlambdaCj(j);
378 currentCij[3] = Co2ABn * (auxiliaryElements.getQ() * ShkmSabmdSdlCj - I * SAlphaGammaCj);
379 currentCij[4] = Co2ABn * (auxiliaryElements.getP() * ShkmSabmdSdlCj - SBetaGammaCj);
380 currentCij[5] = -ax2oAn * gfCoefs.getdSdaCj(j) + BoABpon * (auxiliaryElements.getH() * gfCoefs.getdSdhCj(j) + auxiliaryElements.getK() * gfCoefs.getdSdkCj(j)) + pSagmIqSbgoABnCj + m3onA * gfCoefs.getSCj(j);
381
382 // add the computed coefficients to the interpolators
383 slot.cij[j].addGridPoint(meanState.getDate(), currentCij);
384
385 // Compute the S<sub>i</sub><sup>j</sup> coefficients
386 final double[] currentSij = new double[6];
387
388 // Compute the cross derivatives operator :
389 final double SAlphaGammaSj = context.getAlpha() * gfCoefs.getdSdgammaSj(j) - context.getGamma() * gfCoefs.getdSdalphaSj(j);
390 final double SAlphaBetaSj = context.getAlpha() * gfCoefs.getdSdbetaSj(j) - context.getBeta() * gfCoefs.getdSdalphaSj(j);
391 final double SBetaGammaSj = context.getBeta() * gfCoefs.getdSdgammaSj(j) - context.getGamma() * gfCoefs.getdSdbetaSj(j);
392 final double ShkSj = auxiliaryElements.getH() * gfCoefs.getdSdkSj(j) - auxiliaryElements.getK() * gfCoefs.getdSdhSj(j);
393 final double pSagmIqSbgoABnSj = (auxiliaryElements.getP() * SAlphaGammaSj - I * auxiliaryElements.getQ() * SBetaGammaSj) * ooABn;
394 final double ShkmSabmdSdlSj = ShkSj - SAlphaBetaSj - gfCoefs.getdSdlambdaSj(j);
395
396 currentSij[0] = ax2oAn * gfCoefs.getdSdlambdaSj(j);
397 currentSij[1] = -(BoAn * gfCoefs.getdSdhSj(j) + auxiliaryElements.getH() * pSagmIqSbgoABnSj + auxiliaryElements.getK() * BoABpon * gfCoefs.getdSdlambdaSj(j));
398 currentSij[2] = BoAn * gfCoefs.getdSdkSj(j) + auxiliaryElements.getK() * pSagmIqSbgoABnSj - auxiliaryElements.getH() * BoABpon * gfCoefs.getdSdlambdaSj(j);
399 currentSij[3] = Co2ABn * (auxiliaryElements.getQ() * ShkmSabmdSdlSj - I * SAlphaGammaSj);
400 currentSij[4] = Co2ABn * (auxiliaryElements.getP() * ShkmSabmdSdlSj - SBetaGammaSj);
401 currentSij[5] = -ax2oAn * gfCoefs.getdSdaSj(j) + BoABpon * (auxiliaryElements.getH() * gfCoefs.getdSdhSj(j) + auxiliaryElements.getK() * gfCoefs.getdSdkSj(j)) + pSagmIqSbgoABnSj + m3onA * gfCoefs.getSSj(j);
402
403 // add the computed coefficients to the interpolators
404 slot.sij[j].addGridPoint(meanState.getDate(), currentSij);
405
406 if (j == 1) {
407 //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
408 final double[] value = new double[6];
409 for (int i = 0; i < 6; ++i) {
410 value[i] = currentCij[i] * auxiliaryElements.getK() / 2. + currentSij[i] * auxiliaryElements.getH() / 2.;
411 }
412 slot.cij[0].addGridPoint(meanState.getDate(), value);
413 }
414 }
415 }
416 }
417
418 /** {@inheritDoc} */
419 @Override
420 @SuppressWarnings("unchecked")
421 public <T extends CalculusFieldElement<T>> void updateShortPeriodTerms(final T[] parameters,
422 final FieldSpacecraftState<T>... meanStates) {
423
424 // Field used by default
425 final Field<T> field = meanStates[0].getDate().getField();
426
427 final FieldThirdBodyShortPeriodicCoefficients<T> ftbspc = (FieldThirdBodyShortPeriodicCoefficients<T>) fieldShortPeriods.get(field);
428 final FieldSlot<T> slot = ftbspc.createSlot(meanStates);
429 for (final FieldSpacecraftState<T> meanState : meanStates) {
430
431 // Auxiliary elements related to the current orbit
432 final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(meanState.getOrbit(), I);
433
434 // Container of attributes
435 final FieldDSSTThirdBodyContext<T> context = initializeStep(auxiliaryElements, parameters);
436
437 final FieldHansenObjects<T> fho = (FieldHansenObjects<T>) fieldHansen.get(field);
438
439 final FieldGeneratingFunctionCoefficients<T> gfCoefs =
440 new FieldGeneratingFunctionCoefficients<>(context.getMaxAR3Pow(), MAX_ECCPOWER_SP, context.getMaxAR3Pow() + 1, context, fho, field);
441
442 //Compute additional quantities
443 // 2 * a / An
444 final T ax2oAn = context.getM2aoA().negate().divide(context.getMeanMotion());
445 // B / An
446 final T BoAn = context.getBoA().divide(context.getMeanMotion());
447 // 1 / ABn
448 final T ooABn = context.getOoAB().divide(context.getMeanMotion());
449 // C / 2ABn
450 final T Co2ABn = context.getMCo2AB().negate().divide(context.getMeanMotion());
451 // B / (A * (1 + B) * n)
452 final T BoABpon = context.getBoABpo().divide(context.getMeanMotion());
453 // -3 / n²a² = -3 / nA
454 final T m3onA = context.getA().multiply(context.getMeanMotion()).divide(-3.).reciprocal();
455
456 //Compute the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients.
457 for (int j = 1; j < slot.cij.length; j++) {
458 // First compute the C<sub>i</sub><sup>j</sup> coefficients
459 final T[] currentCij = MathArrays.buildArray(field, 6);
460
461 // Compute the cross derivatives operator :
462 final T SAlphaGammaCj = context.getAlpha().multiply(gfCoefs.getdSdgammaCj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdalphaCj(j)));
463 final T SAlphaBetaCj = context.getAlpha().multiply(gfCoefs.getdSdbetaCj(j)).subtract(context.getBeta().multiply(gfCoefs.getdSdalphaCj(j)));
464 final T SBetaGammaCj = context.getBeta().multiply(gfCoefs.getdSdgammaCj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdbetaCj(j)));
465 final T ShkCj = auxiliaryElements.getH().multiply(gfCoefs.getdSdkCj(j)).subtract(auxiliaryElements.getK().multiply(gfCoefs.getdSdhCj(j)));
466 final T pSagmIqSbgoABnCj = ooABn.multiply(auxiliaryElements.getP().multiply(SAlphaGammaCj).subtract(auxiliaryElements.getQ().multiply(SBetaGammaCj).multiply(I)));
467 final T ShkmSabmdSdlCj = ShkCj.subtract(SAlphaBetaCj).subtract(gfCoefs.getdSdlambdaCj(j));
468
469 currentCij[0] = ax2oAn.multiply(gfCoefs.getdSdlambdaCj(j));
470 currentCij[1] = BoAn.multiply(gfCoefs.getdSdhCj(j)).add(auxiliaryElements.getH().multiply(pSagmIqSbgoABnCj)).add(auxiliaryElements.getK().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaCj(j))).negate();
471 currentCij[2] = BoAn.multiply(gfCoefs.getdSdkCj(j)).add(auxiliaryElements.getK().multiply(pSagmIqSbgoABnCj)).subtract(auxiliaryElements.getH().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaCj(j)));
472 currentCij[3] = Co2ABn.multiply(auxiliaryElements.getQ().multiply(ShkmSabmdSdlCj).subtract(SAlphaGammaCj.multiply(I)));
473 currentCij[4] = Co2ABn.multiply(auxiliaryElements.getP().multiply(ShkmSabmdSdlCj).subtract(SBetaGammaCj));
474 currentCij[5] = ax2oAn.negate().multiply(gfCoefs.getdSdaCj(j)).add(BoABpon.multiply(auxiliaryElements.getH().multiply(gfCoefs.getdSdhCj(j)).add(auxiliaryElements.getK().multiply(gfCoefs.getdSdkCj(j))))).add(pSagmIqSbgoABnCj).add(m3onA.multiply(gfCoefs.getSCj(j)));
475
476 // add the computed coefficients to the interpolators
477 slot.cij[j].addGridPoint(meanState.getDate(), currentCij);
478
479 // Compute the S<sub>i</sub><sup>j</sup> coefficients
480 final T[] currentSij = MathArrays.buildArray(field, 6);
481
482 // Compute the cross derivatives operator :
483 final T SAlphaGammaSj = context.getAlpha().multiply(gfCoefs.getdSdgammaSj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdalphaSj(j)));
484 final T SAlphaBetaSj = context.getAlpha().multiply(gfCoefs.getdSdbetaSj(j)).subtract(context.getBeta().multiply(gfCoefs.getdSdalphaSj(j)));
485 final T SBetaGammaSj = context.getBeta().multiply(gfCoefs.getdSdgammaSj(j)).subtract(context.getGamma().multiply(gfCoefs.getdSdbetaSj(j)));
486 final T ShkSj = auxiliaryElements.getH().multiply(gfCoefs.getdSdkSj(j)).subtract(auxiliaryElements.getK().multiply(gfCoefs.getdSdhSj(j)));
487 final T pSagmIqSbgoABnSj = ooABn.multiply(auxiliaryElements.getP().multiply(SAlphaGammaSj).subtract(auxiliaryElements.getQ().multiply(SBetaGammaSj).multiply(I)));
488 final T ShkmSabmdSdlSj = ShkSj.subtract(SAlphaBetaSj).subtract(gfCoefs.getdSdlambdaSj(j));
489
490 currentSij[0] = ax2oAn.multiply(gfCoefs.getdSdlambdaSj(j));
491 currentSij[1] = BoAn.multiply(gfCoefs.getdSdhSj(j)).add(auxiliaryElements.getH().multiply(pSagmIqSbgoABnSj)).add(auxiliaryElements.getK().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaSj(j))).negate();
492 currentSij[2] = BoAn.multiply(gfCoefs.getdSdkSj(j)).add(auxiliaryElements.getK().multiply(pSagmIqSbgoABnSj)).subtract(auxiliaryElements.getH().multiply(BoABpon).multiply(gfCoefs.getdSdlambdaSj(j)));
493 currentSij[3] = Co2ABn.multiply(auxiliaryElements.getQ().multiply(ShkmSabmdSdlSj).subtract(SAlphaGammaSj.multiply(I)));
494 currentSij[4] = Co2ABn.multiply(auxiliaryElements.getP().multiply(ShkmSabmdSdlSj).subtract(SBetaGammaSj));
495 currentSij[5] = ax2oAn.negate().multiply(gfCoefs.getdSdaSj(j)).add(BoABpon.multiply(auxiliaryElements.getH().multiply(gfCoefs.getdSdhSj(j)).add(auxiliaryElements.getK().multiply(gfCoefs.getdSdkSj(j))))).add(pSagmIqSbgoABnSj).add(m3onA.multiply(gfCoefs.getSSj(j)));
496
497 // add the computed coefficients to the interpolators
498 slot.sij[j].addGridPoint(meanState.getDate(), currentSij);
499
500 if (j == 1) {
501 //Compute the C⁰ coefficients using Danielson 2.5.2-15a.
502 final T[] value = MathArrays.buildArray(field, 6);
503 for (int i = 0; i < 6; ++i) {
504 value[i] = currentCij[i].multiply(auxiliaryElements.getK()).divide(2.).add(currentSij[i].multiply(auxiliaryElements.getH()).divide(2.));
505 }
506 slot.cij[0].addGridPoint(meanState.getDate(), value);
507 }
508 }
509 }
510 }
511
512 /** {@inheritDoc} */
513 @Override
514 public EventDetector[] getEventsDetectors() {
515 return null;
516 }
517
518 /** {@inheritDoc} */
519 @Override
520 public <T extends CalculusFieldElement<T>> FieldEventDetector<T>[] getFieldEventsDetectors(final Field<T> field) {
521 return null;
522 }
523
524 /** {@inheritDoc} */
525 @Override
526 public void registerAttitudeProvider(final AttitudeProvider provider) {
527 //nothing is done since this contribution is not sensitive to attitude
528 }
529
530 /** {@inheritDoc} */
531 @Override
532 public List<ParameterDriver> getParametersDrivers() {
533 return Collections.unmodifiableList(parameterDrivers);
534 }
535
536 /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
537 * and their derivatives.
538 * <p>
539 * CS Mathematical Report $3.5.3.2
540 * </p>
541 */
542 private class FourierCjSjCoefficients {
543
544 /** The coefficients G<sub>n, s</sub> and their derivatives. */
545 private final GnsCoefficients gns;
546
547 /** the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k. */
548 private final WnsjEtomjmsCoefficient wnsjEtomjmsCoefficient;
549
550 /** The terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h). */
551 private final CjSjAlphaBetaKH ABDECoefficients;
552
553 /** The Fourier coefficients C<sup>j</sup> and their derivatives.
554 * <p>
555 * Each column of the matrix contains the following values: <br/>
556 * - C<sup>j</sup> <br/>
557 * - dC<sup>j</sup> / da <br/>
558 * - dC<sup>j</sup> / dk <br/>
559 * - dC<sup>j</sup> / dh <br/>
560 * - dC<sup>j</sup> / dα <br/>
561 * - dC<sup>j</sup> / dβ <br/>
562 * - dC<sup>j</sup> / dγ <br/>
563 * </p>
564 */
565 private final double[][] cj;
566
567 /** The S<sup>j</sup> coefficients and their derivatives.
568 * <p>
569 * Each column of the matrix contains the following values: <br/>
570 * - S<sup>j</sup> <br/>
571 * - dS<sup>j</sup> / da <br/>
572 * - dS<sup>j</sup> / dk <br/>
573 * - dS<sup>j</sup> / dh <br/>
574 * - dS<sup>j</sup> / dα <br/>
575 * - dS<sup>j</sup> / dβ <br/>
576 * - dS<sup>j</sup> / dγ <br/>
577 * </p>
578 */
579 private final double[][] sj;
580
581 /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
582 * <p>
583 * See Danielson 4.2-21
584 * </p>
585 */
586 private final double[] cjlambda;
587
588 /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
589 * <p>
590 * See Danielson 4.2-21
591 * </p>
592 */
593 private final double[] sjlambda;
594
595 /** Maximum value for n. */
596 private final int nMax;
597
598 /** Maximum value for s. */
599 private final int sMax;
600
601 /** Maximum value for j. */
602 private final int jMax;
603
604 /**
605 * Private constructor.
606 *
607 * @param nMax maximum value for n index
608 * @param sMax maximum value for s index
609 * @param jMax maximum value for j index
610 * @param context container for attributes
611 */
612 FourierCjSjCoefficients(final int nMax, final int sMax, final int jMax, final DSSTThirdBodyContext context) {
613 //Save parameters
614 this.nMax = nMax;
615 this.sMax = sMax;
616 this.jMax = jMax;
617
618 //Create objects
619 wnsjEtomjmsCoefficient = new WnsjEtomjmsCoefficient(context);
620 ABDECoefficients = new CjSjAlphaBetaKH(context);
621 gns = new GnsCoefficients(nMax, sMax, context);
622
623 //create arays
624 this.cj = new double[7][jMax + 1];
625 this.sj = new double[7][jMax + 1];
626 this.cjlambda = new double[jMax];
627 this.sjlambda = new double[jMax];
628
629 computeCoefficients(context);
630 }
631
632 /**
633 * Compute all coefficients.
634 * @param context container for attributes
635 */
636 private void computeCoefficients(final DSSTThirdBodyContext context) {
637
638 final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
639
640 for (int j = 1; j <= jMax; j++) {
641 // initialise the coefficients
642 for (int i = 0; i <= 6; i++) {
643 cj[i][j] = 0.;
644 sj[i][j] = 0.;
645 }
646 if (j < jMax) {
647 // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
648 cjlambda[j] = 0.;
649 sjlambda[j] = 0.;
650 }
651 for (int s = 0; s <= sMax; s++) {
652
653 // Compute the coefficients A<sub>j, s</sub>, B<sub>j, s</sub>, D<sub>j, s</sub> and E<sub>j, s</sub>
654 ABDECoefficients.computeCoefficients(j, s);
655
656 // compute starting value for n
657 final int minN = FastMath.max(2, FastMath.max(j - 1, s));
658
659 for (int n = minN; n <= nMax; n++) {
660 // check if n-s is even
661 if ((n - s) % 2 == 0) {
662 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n+1, s</sup> and its derivatives
663 final double[] wjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n + 1, context);
664
665 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n+1, s</sup> and its derivatives
666 final double[] wmjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n + 1, context);
667
668 // compute common factors
669 final double coef1 = -(wjnp1semjms[0] * ABDECoefficients.getCoefA() + wmjnp1semjms[0] * ABDECoefficients.getCoefB());
670 final double coef2 = wjnp1semjms[0] * ABDECoefficients.getCoefD() + wmjnp1semjms[0] * ABDECoefficients.getCoefE();
671
672 //Compute C<sup>j</sup>
673 cj[0][j] += gns.getGns(n, s) * coef1;
674 //Compute dC<sup>j</sup> / da
675 cj[1][j] += gns.getdGnsda(n, s) * coef1;
676 //Compute dC<sup>j</sup> / dk
677 cj[2][j] += -gns.getGns(n, s) *
678 (
679 wjnp1semjms[1] * ABDECoefficients.getCoefA() +
680 wjnp1semjms[0] * ABDECoefficients.getdCoefAdk() +
681 wmjnp1semjms[1] * ABDECoefficients.getCoefB() +
682 wmjnp1semjms[0] * ABDECoefficients.getdCoefBdk()
683 );
684 //Compute dC<sup>j</sup> / dh
685 cj[3][j] += -gns.getGns(n, s) *
686 (
687 wjnp1semjms[2] * ABDECoefficients.getCoefA() +
688 wjnp1semjms[0] * ABDECoefficients.getdCoefAdh() +
689 wmjnp1semjms[2] * ABDECoefficients.getCoefB() +
690 wmjnp1semjms[0] * ABDECoefficients.getdCoefBdh()
691 );
692 //Compute dC<sup>j</sup> / dα
693 cj[4][j] += -gns.getGns(n, s) *
694 (
695 wjnp1semjms[0] * ABDECoefficients.getdCoefAdalpha() +
696 wmjnp1semjms[0] * ABDECoefficients.getdCoefBdalpha()
697 );
698 //Compute dC<sup>j</sup> / dβ
699 cj[5][j] += -gns.getGns(n, s) *
700 (
701 wjnp1semjms[0] * ABDECoefficients.getdCoefAdbeta() +
702 wmjnp1semjms[0] * ABDECoefficients.getdCoefBdbeta()
703 );
704 //Compute dC<sup>j</sup> / dγ
705 cj[6][j] += gns.getdGnsdgamma(n, s) * coef1;
706
707 //Compute S<sup>j</sup>
708 sj[0][j] += gns.getGns(n, s) * coef2;
709 //Compute dS<sup>j</sup> / da
710 sj[1][j] += gns.getdGnsda(n, s) * coef2;
711 //Compute dS<sup>j</sup> / dk
712 sj[2][j] += gns.getGns(n, s) *
713 (
714 wjnp1semjms[1] * ABDECoefficients.getCoefD() +
715 wjnp1semjms[0] * ABDECoefficients.getdCoefDdk() +
716 wmjnp1semjms[1] * ABDECoefficients.getCoefE() +
717 wmjnp1semjms[0] * ABDECoefficients.getdCoefEdk()
718 );
719 //Compute dS<sup>j</sup> / dh
720 sj[3][j] += gns.getGns(n, s) *
721 (
722 wjnp1semjms[2] * ABDECoefficients.getCoefD() +
723 wjnp1semjms[0] * ABDECoefficients.getdCoefDdh() +
724 wmjnp1semjms[2] * ABDECoefficients.getCoefE() +
725 wmjnp1semjms[0] * ABDECoefficients.getdCoefEdh()
726 );
727 //Compute dS<sup>j</sup> / dα
728 sj[4][j] += gns.getGns(n, s) *
729 (
730 wjnp1semjms[0] * ABDECoefficients.getdCoefDdalpha() +
731 wmjnp1semjms[0] * ABDECoefficients.getdCoefEdalpha()
732 );
733 //Compute dS<sup>j</sup> / dβ
734 sj[5][j] += gns.getGns(n, s) *
735 (
736 wjnp1semjms[0] * ABDECoefficients.getdCoefDdbeta() +
737 wmjnp1semjms[0] * ABDECoefficients.getdCoefEdbeta()
738 );
739 //Compute dS<sup>j</sup> / dγ
740 sj[6][j] += gns.getdGnsdgamma(n, s) * coef2;
741
742 //Check if n is greater or equal to j and j is at most jMax-1
743 if (n >= j && j < jMax) {
744 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives
745 final double[] wjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n, context);
746
747 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n, s</sup> and its derivatives
748 final double[] wmjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n, context);
749
750 //Compute C<sup>j</sup><sub>,λ</sub>
751 cjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefD() + wmjnsemjms[0] * ABDECoefficients.getCoefE());
752 //Compute S<sup>j</sup><sub>,λ</sub>
753 sjlambda[j] += gns.getGns(n, s) * (wjnsemjms[0] * ABDECoefficients.getCoefA() + wmjnsemjms[0] * ABDECoefficients.getCoefB());
754 }
755 }
756 }
757 }
758 // Divide by j
759 for (int i = 0; i <= 6; i++) {
760 cj[i][j] /= j;
761 sj[i][j] /= j;
762 }
763 }
764 //The C⁰ coefficients are not computed here.
765 //They are evaluated at the final point.
766
767 //C⁰<sub>,λ</sub>
768 cjlambda[0] = auxiliaryElements.getK() * cjlambda[1] / 2. + auxiliaryElements.getH() * sjlambda[1] / 2.;
769 }
770
771 /** Get the Fourier coefficient C<sup>j</sup>.
772 * @param j j index
773 * @return C<sup>j</sup>
774 */
775 public double getCj(final int j) {
776 return cj[0][j];
777 }
778
779 /** Get the derivative dC<sup>j</sup>/da.
780 * @param j j index
781 * @return dC<sup>j</sup>/da
782 */
783 public double getdCjda(final int j) {
784 return cj[1][j];
785 }
786
787 /** Get the derivative dC<sup>j</sup>/dk.
788 * @param j j index
789 * @return dC<sup>j</sup>/dk
790 */
791 public double getdCjdk(final int j) {
792 return cj[2][j];
793 }
794
795 /** Get the derivative dC<sup>j</sup>/dh.
796 * @param j j index
797 * @return dC<sup>j</sup>/dh
798 */
799 public double getdCjdh(final int j) {
800 return cj[3][j];
801 }
802
803 /** Get the derivative dC<sup>j</sup>/dα.
804 * @param j j index
805 * @return dC<sup>j</sup>/dα
806 */
807 public double getdCjdalpha(final int j) {
808 return cj[4][j];
809 }
810
811 /** Get the derivative dC<sup>j</sup>/dβ.
812 * @param j j index
813 * @return dC<sup>j</sup>/dβ
814 */
815 public double getdCjdbeta(final int j) {
816 return cj[5][j];
817 }
818
819 /** Get the derivative dC<sup>j</sup>/dγ.
820 * @param j j index
821 * @return dC<sup>j</sup>/dγ
822 */
823 public double getdCjdgamma(final int j) {
824 return cj[6][j];
825 }
826
827 /** Get the Fourier coefficient S<sup>j</sup>.
828 * @param j j index
829 * @return S<sup>j</sup>
830 */
831 public double getSj(final int j) {
832 return sj[0][j];
833 }
834
835 /** Get the derivative dS<sup>j</sup>/da.
836 * @param j j index
837 * @return dS<sup>j</sup>/da
838 */
839 public double getdSjda(final int j) {
840 return sj[1][j];
841 }
842
843 /** Get the derivative dS<sup>j</sup>/dk.
844 * @param j j index
845 * @return dS<sup>j</sup>/dk
846 */
847 public double getdSjdk(final int j) {
848 return sj[2][j];
849 }
850
851 /** Get the derivative dS<sup>j</sup>/dh.
852 * @param j j index
853 * @return dS<sup>j</sup>/dh
854 */
855 public double getdSjdh(final int j) {
856 return sj[3][j];
857 }
858
859 /** Get the derivative dS<sup>j</sup>/dα.
860 * @param j j index
861 * @return dS<sup>j</sup>/dα
862 */
863 public double getdSjdalpha(final int j) {
864 return sj[4][j];
865 }
866
867 /** Get the derivative dS<sup>j</sup>/dβ.
868 * @param j j index
869 * @return dS<sup>j</sup>/dβ
870 */
871 public double getdSjdbeta(final int j) {
872 return sj[5][j];
873 }
874
875 /** Get the derivative dS<sup>j</sup>/dγ.
876 * @param j j index
877 * @return dS<sup>j</sup>/dγ
878 */
879 public double getdSjdgamma(final int j) {
880 return sj[6][j];
881 }
882
883 /** Get the coefficient C⁰<sub>,λ</sub>.
884 * @return C⁰<sub>,λ</sub>
885 */
886 public double getC0Lambda() {
887 return cjlambda[0];
888 }
889
890 /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
891 * @param j j index
892 * @return C<sup>j</sup><sub>,λ</sub>
893 */
894 public double getCjLambda(final int j) {
895 if (j < 1 || j >= jMax) {
896 return 0.;
897 }
898 return cjlambda[j];
899 }
900
901 /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
902 * @param j j index
903 * @return S<sup>j</sup><sub>,λ</sub>
904 */
905 public double getSjLambda(final int j) {
906 if (j < 1 || j >= jMax) {
907 return 0.;
908 }
909 return sjlambda[j];
910 }
911 }
912
913 /** Computes the C<sup>j</sup> and S<sup>j</sup> coefficients Danielson 4.2-(15,16)
914 * and their derivatives.
915 * <p>
916 * CS Mathematical Report $3.5.3.2
917 * </p>
918 */
919 private class FieldFourierCjSjCoefficients <T extends CalculusFieldElement<T>> {
920
921 /** The coefficients G<sub>n, s</sub> and their derivatives. */
922 private final FieldGnsCoefficients<T> gns;
923
924 /** the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k. */
925 private final FieldWnsjEtomjmsCoefficient<T> wnsjEtomjmsCoefficient;
926
927 /** The terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h). */
928 private final FieldCjSjAlphaBetaKH<T> ABDECoefficients;
929
930 /** The Fourier coefficients C<sup>j</sup> and their derivatives.
931 * <p>
932 * Each column of the matrix contains the following values: <br/>
933 * - C<sup>j</sup> <br/>
934 * - dC<sup>j</sup> / da <br/>
935 * - dC<sup>j</sup> / dk <br/>
936 * - dC<sup>j</sup> / dh <br/>
937 * - dC<sup>j</sup> / dα <br/>
938 * - dC<sup>j</sup> / dβ <br/>
939 * - dC<sup>j</sup> / dγ <br/>
940 * </p>
941 */
942 private final T[][] cj;
943
944 /** The S<sup>j</sup> coefficients and their derivatives.
945 * <p>
946 * Each column of the matrix contains the following values: <br/>
947 * - S<sup>j</sup> <br/>
948 * - dS<sup>j</sup> / da <br/>
949 * - dS<sup>j</sup> / dk <br/>
950 * - dS<sup>j</sup> / dh <br/>
951 * - dS<sup>j</sup> / dα <br/>
952 * - dS<sup>j</sup> / dβ <br/>
953 * - dS<sup>j</sup> / dγ <br/>
954 * </p>
955 */
956 private final T[][] sj;
957
958 /** The Coefficients C<sup>j</sup><sub>,λ</sub>.
959 * <p>
960 * See Danielson 4.2-21
961 * </p>
962 */
963 private final T[] cjlambda;
964
965 /** The Coefficients S<sup>j</sup><sub>,λ</sub>.
966 * <p>
967 * See Danielson 4.2-21
968 * </p>
969 */
970 private final T[] sjlambda;
971
972 /** Zero. */
973 private final T zero;
974
975 /** Maximum value for n. */
976 private final int nMax;
977
978 /** Maximum value for s. */
979 private final int sMax;
980
981 /** Maximum value for j. */
982 private final int jMax;
983
984 /**
985 * Private constructor.
986 *
987 * @param nMax maximum value for n index
988 * @param sMax maximum value for s index
989 * @param jMax maximum value for j index
990 * @param context container for attributes
991 * @param field field used by default
992 */
993 FieldFourierCjSjCoefficients(final int nMax, final int sMax, final int jMax,
994 final FieldDSSTThirdBodyContext<T> context,
995 final Field<T> field) {
996 //Zero
997 this.zero = field.getZero();
998
999 //Save parameters
1000 this.nMax = nMax;
1001 this.sMax = sMax;
1002 this.jMax = jMax;
1003
1004 //Create objects
1005 wnsjEtomjmsCoefficient = new FieldWnsjEtomjmsCoefficient<>(context, field);
1006 ABDECoefficients = new FieldCjSjAlphaBetaKH<>(context, field);
1007 gns = new FieldGnsCoefficients<>(nMax, sMax, context, field);
1008
1009 //create arays
1010 this.cj = MathArrays.buildArray(field, 7, jMax + 1);
1011 this.sj = MathArrays.buildArray(field, 7, jMax + 1);
1012 this.cjlambda = MathArrays.buildArray(field, jMax);
1013 this.sjlambda = MathArrays.buildArray(field, jMax);
1014
1015 computeCoefficients(context, field);
1016 }
1017
1018 /**
1019 * Compute all coefficients.
1020 * @param context container for attributes
1021 * @param field field used by default
1022 */
1023 private void computeCoefficients(final FieldDSSTThirdBodyContext<T> context,
1024 final Field<T> field) {
1025
1026 final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
1027
1028 for (int j = 1; j <= jMax; j++) {
1029 // initialise the coefficients
1030 for (int i = 0; i <= 6; i++) {
1031 cj[i][j] = zero;
1032 sj[i][j] = zero;
1033 }
1034 if (j < jMax) {
1035 // initialise the C<sup>j</sup><sub>,λ</sub> and S<sup>j</sup><sub>,λ</sub> coefficients
1036 cjlambda[j] = zero;
1037 sjlambda[j] = zero;
1038 }
1039 for (int s = 0; s <= sMax; s++) {
1040
1041 // Compute the coefficients A<sub>j, s</sub>, B<sub>j, s</sub>, D<sub>j, s</sub> and E<sub>j, s</sub>
1042 ABDECoefficients.computeCoefficients(j, s);
1043
1044 // compute starting value for n
1045 final int minN = FastMath.max(2, FastMath.max(j - 1, s));
1046
1047 for (int n = minN; n <= nMax; n++) {
1048 // check if n-s is even
1049 if ((n - s) % 2 == 0) {
1050 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n+1, s</sup> and its derivatives
1051 final T[] wjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n + 1, context, field);
1052
1053 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n+1, s</sup> and its derivatives
1054 final T[] wmjnp1semjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n + 1, context, field);
1055
1056 // compute common factors
1057 final T coef1 = (wjnp1semjms[0].multiply(ABDECoefficients.getCoefA()).add(wmjnp1semjms[0].multiply(ABDECoefficients.getCoefB()))).negate();
1058 final T coef2 = wjnp1semjms[0].multiply(ABDECoefficients.getCoefD()).add(wmjnp1semjms[0].multiply(ABDECoefficients.getCoefE()));
1059
1060 //Compute C<sup>j</sup>
1061 cj[0][j] = cj[0][j].add(gns.getGns(n, s).multiply(coef1));
1062 //Compute dC<sup>j</sup> / da
1063 cj[1][j] = cj[1][j].add(gns.getdGnsda(n, s).multiply(coef1));
1064 //Compute dC<sup>j</sup> / dk
1065 cj[2][j] = cj[2][j].add(gns.getGns(n, s).negate().
1066 multiply(
1067 wjnp1semjms[1].multiply(ABDECoefficients.getCoefA()).
1068 add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdk())).
1069 add(wmjnp1semjms[1].multiply(ABDECoefficients.getCoefB())).
1070 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdk()))
1071 ));
1072 //Compute dC<sup>j</sup> / dh
1073 cj[3][j] = cj[3][j].add(gns.getGns(n, s).negate().
1074 multiply(
1075 wjnp1semjms[2].multiply(ABDECoefficients.getCoefA()).
1076 add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdh())).
1077 add(wmjnp1semjms[2].multiply(ABDECoefficients.getCoefB())).
1078 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdh()))
1079 ));
1080 //Compute dC<sup>j</sup> / dα
1081 cj[4][j] = cj[4][j].add(gns.getGns(n, s).negate().
1082 multiply(
1083 wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdalpha()).
1084 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdalpha()))
1085 ));
1086 //Compute dC<sup>j</sup> / dβ
1087 cj[5][j] = cj[5][j].add(gns.getGns(n, s).negate().
1088 multiply(
1089 wjnp1semjms[0].multiply(ABDECoefficients.getdCoefAdbeta()).
1090 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefBdbeta()))
1091 ));
1092 //Compute dC<sup>j</sup> / dγ
1093 cj[6][j] = cj[6][j].add(gns.getdGnsdgamma(n, s).multiply(coef1));
1094
1095 //Compute S<sup>j</sup>
1096 sj[0][j] = sj[0][j].add(gns.getGns(n, s).multiply(coef2));
1097 //Compute dS<sup>j</sup> / da
1098 sj[1][j] = sj[1][j].add(gns.getdGnsda(n, s).multiply(coef2));
1099 //Compute dS<sup>j</sup> / dk
1100 sj[2][j] = sj[2][j].add(gns.getGns(n, s).
1101 multiply(
1102 wjnp1semjms[1].multiply(ABDECoefficients.getCoefD()).
1103 add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdk())).
1104 add(wmjnp1semjms[1].multiply(ABDECoefficients.getCoefE())).
1105 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdk()))
1106 ));
1107 //Compute dS<sup>j</sup> / dh
1108 sj[3][j] = sj[3][j].add(gns.getGns(n, s).
1109 multiply(
1110 wjnp1semjms[2].multiply(ABDECoefficients.getCoefD()).
1111 add(wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdh())).
1112 add(wmjnp1semjms[2].multiply(ABDECoefficients.getCoefE())).
1113 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdh()))
1114 ));
1115 //Compute dS<sup>j</sup> / dα
1116 sj[4][j] = sj[4][j].add(gns.getGns(n, s).
1117 multiply(
1118 wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdalpha()).
1119 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdalpha()))
1120 ));
1121 //Compute dS<sup>j</sup> / dβ
1122 sj[5][j] = sj[5][j].add(gns.getGns(n, s).
1123 multiply(
1124 wjnp1semjms[0].multiply(ABDECoefficients.getdCoefDdbeta()).
1125 add(wmjnp1semjms[0].multiply(ABDECoefficients.getdCoefEdbeta()))
1126 ));
1127 //Compute dS<sup>j</sup> / dγ
1128 sj[6][j] = sj[6][j].add(gns.getdGnsdgamma(n, s).multiply(coef2));
1129
1130 //Check if n is greater or equal to j and j is at most jMax-1
1131 if (n >= j && j < jMax) {
1132 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives
1133 final T[] wjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(j, s, n, context, field);
1134
1135 // compute the coefficient e<sup>-|j-s|</sup>*w<sub>-j</sub><sup>n, s</sup> and its derivatives
1136 final T[] wmjnsemjms = wnsjEtomjmsCoefficient.computeWjnsEmjmsAndDeriv(-j, s, n, context, field);
1137
1138 //Compute C<sup>j</sup><sub>,λ</sub>
1139 cjlambda[j] = cjlambda[j].add(gns.getGns(n, s).multiply(wjnsemjms[0].multiply(ABDECoefficients.getCoefD()).add(wmjnsemjms[0].multiply(ABDECoefficients.getCoefE()))));
1140 //Compute S<sup>j</sup><sub>,λ</sub>
1141 sjlambda[j] = sjlambda[j].add(gns.getGns(n, s).multiply(wjnsemjms[0].multiply(ABDECoefficients.getCoefA()).add(wmjnsemjms[0].multiply(ABDECoefficients.getCoefB()))));
1142 }
1143 }
1144 }
1145 }
1146 // Divide by j
1147 for (int i = 0; i <= 6; i++) {
1148 cj[i][j] = cj[i][j].divide(j);
1149 sj[i][j] = sj[i][j].divide(j);
1150 }
1151 }
1152 //The C⁰ coefficients are not computed here.
1153 //They are evaluated at the final point.
1154
1155 //C⁰<sub>,λ</sub>
1156 cjlambda[0] = auxiliaryElements.getK().multiply(cjlambda[1]).divide(2.).add(auxiliaryElements.getH().multiply(sjlambda[1]).divide(2.));
1157 }
1158
1159 /** Get the Fourier coefficient C<sup>j</sup>.
1160 * @param j j index
1161 * @return C<sup>j</sup>
1162 */
1163 public T getCj(final int j) {
1164 return cj[0][j];
1165 }
1166
1167 /** Get the derivative dC<sup>j</sup>/da.
1168 * @param j j index
1169 * @return dC<sup>j</sup>/da
1170 */
1171 public T getdCjda(final int j) {
1172 return cj[1][j];
1173 }
1174
1175 /** Get the derivative dC<sup>j</sup>/dk.
1176 * @param j j index
1177 * @return dC<sup>j</sup>/dk
1178 */
1179 public T getdCjdk(final int j) {
1180 return cj[2][j];
1181 }
1182
1183 /** Get the derivative dC<sup>j</sup>/dh.
1184 * @param j j index
1185 * @return dC<sup>j</sup>/dh
1186 */
1187 public T getdCjdh(final int j) {
1188 return cj[3][j];
1189 }
1190
1191 /** Get the derivative dC<sup>j</sup>/dα.
1192 * @param j j index
1193 * @return dC<sup>j</sup>/dα
1194 */
1195 public T getdCjdalpha(final int j) {
1196 return cj[4][j];
1197 }
1198
1199 /** Get the derivative dC<sup>j</sup>/dβ.
1200 * @param j j index
1201 * @return dC<sup>j</sup>/dβ
1202 */
1203 public T getdCjdbeta(final int j) {
1204 return cj[5][j];
1205 }
1206
1207 /** Get the derivative dC<sup>j</sup>/dγ.
1208 * @param j j index
1209 * @return dC<sup>j</sup>/dγ
1210 */
1211 public T getdCjdgamma(final int j) {
1212 return cj[6][j];
1213 }
1214
1215 /** Get the Fourier coefficient S<sup>j</sup>.
1216 * @param j j index
1217 * @return S<sup>j</sup>
1218 */
1219 public T getSj(final int j) {
1220 return sj[0][j];
1221 }
1222
1223 /** Get the derivative dS<sup>j</sup>/da.
1224 * @param j j index
1225 * @return dS<sup>j</sup>/da
1226 */
1227 public T getdSjda(final int j) {
1228 return sj[1][j];
1229 }
1230
1231 /** Get the derivative dS<sup>j</sup>/dk.
1232 * @param j j index
1233 * @return dS<sup>j</sup>/dk
1234 */
1235 public T getdSjdk(final int j) {
1236 return sj[2][j];
1237 }
1238
1239 /** Get the derivative dS<sup>j</sup>/dh.
1240 * @param j j index
1241 * @return dS<sup>j</sup>/dh
1242 */
1243 public T getdSjdh(final int j) {
1244 return sj[3][j];
1245 }
1246
1247 /** Get the derivative dS<sup>j</sup>/dα.
1248 * @param j j index
1249 * @return dS<sup>j</sup>/dα
1250 */
1251 public T getdSjdalpha(final int j) {
1252 return sj[4][j];
1253 }
1254
1255 /** Get the derivative dS<sup>j</sup>/dβ.
1256 * @param j j index
1257 * @return dS<sup>j</sup>/dβ
1258 */
1259 public T getdSjdbeta(final int j) {
1260 return sj[5][j];
1261 }
1262
1263 /** Get the derivative dS<sup>j</sup>/dγ.
1264 * @param j j index
1265 * @return dS<sup>j</sup>/dγ
1266 */
1267 public T getdSjdgamma(final int j) {
1268 return sj[6][j];
1269 }
1270
1271 /** Get the coefficient C⁰<sub>,λ</sub>.
1272 * @return C⁰<sub>,λ</sub>
1273 */
1274 public T getC0Lambda() {
1275 return cjlambda[0];
1276 }
1277
1278 /** Get the coefficient C<sup>j</sup><sub>,λ</sub>.
1279 * @param j j index
1280 * @return C<sup>j</sup><sub>,λ</sub>
1281 */
1282 public T getCjLambda(final int j) {
1283 if (j < 1 || j >= jMax) {
1284 return zero;
1285 }
1286 return cjlambda[j];
1287 }
1288
1289 /** Get the coefficient S<sup>j</sup><sub>,λ</sub>.
1290 * @param j j index
1291 * @return S<sup>j</sup><sub>,λ</sub>
1292 */
1293 public T getSjLambda(final int j) {
1294 if (j < 1 || j >= jMax) {
1295 return zero;
1296 }
1297 return sjlambda[j];
1298 }
1299 }
1300
1301 /** This class covers the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k.
1302 *
1303 * <p>
1304 * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
1305 * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
1306 * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
1307 * can be written as: <br/ >
1308 * - for |s| > |j| <br />
1309 * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
1310 * (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
1311 * (-b)<sup>|j-s|</sup> *
1312 * ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
1313 * P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
1314 * <br />
1315 * - for |s| <= |j| <br />
1316 * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
1317 * (-b)<sup>|j-s|</sup> *
1318 * ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
1319 * P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
1320 * </p>
1321 *
1322 * @author Lucian Barbulescu
1323 */
1324 private static class WnsjEtomjmsCoefficient {
1325
1326 /** The value c.
1327 * <p>
1328 * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
1329 * </p>
1330 * */
1331 private final double c;
1332
1333 /** db / dh. */
1334 private final double dbdh;
1335
1336 /** db / dk. */
1337 private final double dbdk;
1338
1339 /** dc / dh = e * db/dh + b * de/dh. */
1340 private final double dcdh;
1341
1342 /** dc / dk = e * db/dk + b * de/dk. */
1343 private final double dcdk;
1344
1345 /** The values (1 - c²)<sup>n</sup>. <br />
1346 * The maximum possible value for the power is N + 1 */
1347 private final double[] omc2tn;
1348
1349 /** The values (1 + c²)<sup>n</sup>. <br />
1350 * The maximum possible value for the power is N + 1 */
1351 private final double[] opc2tn;
1352
1353 /** The values b<sup>|j-s|</sup>. */
1354 private final double[] btjms;
1355
1356 /**
1357 * Standard constructor.
1358 * @param context container for attributes
1359 */
1360 WnsjEtomjmsCoefficient(final DSSTThirdBodyContext context) {
1361
1362 final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
1363
1364 //initialise fields
1365 c = auxiliaryElements.getEcc() * context.getb();
1366 final double c2 = c * c;
1367
1368 //b² * χ
1369 final double b2Chi = context.getb() * context.getb() * context.getX();
1370 //Compute derivatives of b
1371 dbdh = auxiliaryElements.getH() * b2Chi;
1372 dbdk = auxiliaryElements.getK() * b2Chi;
1373
1374 //Compute derivatives of c
1375 if (auxiliaryElements.getEcc() == 0.0) {
1376 // we are at a perfectly circular orbit singularity here
1377 // we arbitrarily consider the perigee is along the X axis,
1378 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
1379 dcdh = auxiliaryElements.getEcc() * dbdh + context.getb();
1380 dcdk = auxiliaryElements.getEcc() * dbdk;
1381 } else {
1382 dcdh = auxiliaryElements.getEcc() * dbdh + context.getb() * auxiliaryElements.getH() / auxiliaryElements.getEcc();
1383 dcdk = auxiliaryElements.getEcc() * dbdk + context.getb() * auxiliaryElements.getK() / auxiliaryElements.getEcc();
1384 }
1385
1386 //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
1387 omc2tn = new double[context.getMaxAR3Pow() + context.getMaxFreqF() + 2];
1388 opc2tn = new double[context.getMaxAR3Pow() + context.getMaxFreqF() + 2];
1389 final double omc2 = 1. - c2;
1390 final double opc2 = 1. + c2;
1391 omc2tn[0] = 1.;
1392 opc2tn[0] = 1.;
1393 for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF() + 1; i++) {
1394 omc2tn[i] = omc2tn[i - 1] * omc2;
1395 opc2tn[i] = opc2tn[i - 1] * opc2;
1396 }
1397
1398 //Compute the powers of b
1399 btjms = new double[context.getMaxAR3Pow() + context.getMaxFreqF() + 1];
1400 btjms[0] = 1.;
1401 for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF(); i++) {
1402 btjms[i] = btjms[i - 1] * context.getb();
1403 }
1404 }
1405
1406 /** Compute the value of the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by h and k. <br />
1407 *
1408 * @param j j index
1409 * @param s s index
1410 * @param n n index
1411 * @param context container for attributes
1412 * @return an array containing the value of the coefficient at index 0, the derivative by k at index 1 and the derivative by h at index 2
1413 */
1414 public double[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n, final DSSTThirdBodyContext context) {
1415 final double[] wjnsemjms = new double[] {0., 0., 0.};
1416
1417 // |j|
1418 final int absJ = FastMath.abs(j);
1419 // |s|
1420 final int absS = FastMath.abs(s);
1421 // |j - s|
1422 final int absJmS = FastMath.abs(j - s);
1423 // |j + s|
1424 final int absJpS = FastMath.abs(j + s);
1425
1426 //The lower index of P. Also the power of (1 - c²)
1427 final int l;
1428 // the factorial ratio coefficient or 1. if |s| <= |j|
1429 final double factCoef;
1430 if (absS > absJ) {
1431 //factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
1432 factCoef = (CombinatoricsUtils.factorialDouble(n + s) / CombinatoricsUtils.factorialDouble(n + j)) * (CombinatoricsUtils.factorialDouble(n - s) / CombinatoricsUtils.factorialDouble(n - j));
1433 l = n - absS;
1434 } else {
1435 factCoef = 1.;
1436 l = n - absJ;
1437 }
1438
1439 // (-1)<sup>|j-s|</sup>
1440 final double sign = absJmS % 2 != 0 ? -1. : 1.;
1441 //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
1442 final double coef1 = omc2tn[l] / opc2tn[n];
1443 //-b<sup>|j-s|</sup>
1444 final double coef2 = sign * btjms[absJmS];
1445 // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
1446 final Gradient jac =
1447 JacobiPolynomials.getValue(l, absJmS, absJpS, Gradient.variable(1, 0, context.getX()));
1448
1449 // the derivative of coef1 by c
1450 final double dcoef1dc = -coef1 * 2. * c * (((double) n) / opc2tn[1] + ((double) l) / omc2tn[1]);
1451 // the derivative of coef1 by h
1452 final double dcoef1dh = dcoef1dc * dcdh;
1453 // the derivative of coef1 by k
1454 final double dcoef1dk = dcoef1dc * dcdk;
1455
1456 // the derivative of coef2 by b
1457 final double dcoef2db = absJmS == 0 ? 0 : sign * (double) absJmS * btjms[absJmS - 1];
1458 // the derivative of coef2 by h
1459 final double dcoef2dh = dcoef2db * dbdh;
1460 // the derivative of coef2 by k
1461 final double dcoef2dk = dcoef2db * dbdk;
1462
1463 // the jacobi polynomial value
1464 final double jacobi = jac.getValue();
1465 // the derivative of the Jacobi polynomial by h
1466 final double djacobidh = jac.getGradient()[0] * context.getHXXX();
1467 // the derivative of the Jacobi polynomial by k
1468 final double djacobidk = jac.getGradient()[0] * context.getKXXX();
1469
1470 //group the above coefficients to limit the mathematical operations
1471 final double term1 = factCoef * coef1 * coef2;
1472 final double term2 = factCoef * coef1 * jacobi;
1473 final double term3 = factCoef * coef2 * jacobi;
1474
1475 //compute e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by k and h
1476 wjnsemjms[0] = term1 * jacobi;
1477 wjnsemjms[1] = dcoef1dk * term3 + dcoef2dk * term2 + djacobidk * term1;
1478 wjnsemjms[2] = dcoef1dh * term3 + dcoef2dh * term2 + djacobidh * term1;
1479
1480 return wjnsemjms;
1481 }
1482 }
1483
1484 /** This class covers the coefficients e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and their derivatives by h and k.
1485 *
1486 * <p>
1487 * Starting from Danielson 4.2-9,10,11 and taking into account that fact that: <br />
1488 * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
1489 * the expression e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup>
1490 * can be written as: <br/ >
1491 * - for |s| > |j| <br />
1492 * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
1493 * (((n + s)!(n - s)!)/((n + j)!(n - j)!)) *
1494 * (-b)<sup>|j-s|</sup> *
1495 * ((1 - c²)<sup>n-|s|</sup>/(1 + c²)<sup>n</sup>) *
1496 * P<sub>n-|s|</sub><sup>|j-s|, |j+s|</sup>(χ) <br />
1497 * <br />
1498 * - for |s| <= |j| <br />
1499 * e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> =
1500 * (-b)<sup>|j-s|</sup> *
1501 * ((1 - c²)<sup>n-|j|</sup>/(1 + c²)<sup>n</sup>) *
1502 * P<sub>n-|j|</sub><sup>|j-s|, |j+s|</sup>(χ)
1503 * </p>
1504 *
1505 * @author Lucian Barbulescu
1506 */
1507 private static class FieldWnsjEtomjmsCoefficient <T extends CalculusFieldElement<T>> {
1508
1509 /** The value c.
1510 * <p>
1511 * c = e / (1 + (1 - e²)<sup>1/2</sup>) = e / (1 + B) = e * b <br/>
1512 * </p>
1513 * */
1514 private final T c;
1515
1516 /** db / dh. */
1517 private final T dbdh;
1518
1519 /** db / dk. */
1520 private final T dbdk;
1521
1522 /** dc / dh = e * db/dh + b * de/dh. */
1523 private final T dcdh;
1524
1525 /** dc / dk = e * db/dk + b * de/dk. */
1526 private final T dcdk;
1527
1528 /** The values (1 - c²)<sup>n</sup>. <br />
1529 * The maximum possible value for the power is N + 1 */
1530 private final T[] omc2tn;
1531
1532 /** The values (1 + c²)<sup>n</sup>. <br />
1533 * The maximum possible value for the power is N + 1 */
1534 private final T[] opc2tn;
1535
1536 /** The values b<sup>|j-s|</sup>. */
1537 private final T[] btjms;
1538
1539 /**
1540 * Standard constructor.
1541 * @param context container for attributes
1542 * @param field field used by default
1543 */
1544 FieldWnsjEtomjmsCoefficient(final FieldDSSTThirdBodyContext<T> context, final Field<T> field) {
1545
1546 final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
1547
1548 //Zero
1549 final T zero = field.getZero();
1550
1551 //initialise fields
1552 c = auxiliaryElements.getEcc().multiply(context.getb());
1553 final T c2 = c.multiply(c);
1554
1555 //b² * χ
1556 final T b2Chi = context.getb().multiply(context.getb()).multiply(context.getX());
1557 //Compute derivatives of b
1558 dbdh = auxiliaryElements.getH().multiply(b2Chi);
1559 dbdk = auxiliaryElements.getK().multiply(b2Chi);
1560
1561 //Compute derivatives of c
1562 if (auxiliaryElements.getEcc().getReal() == 0.0) {
1563 // we are at a perfectly circular orbit singularity here
1564 // we arbitrarily consider the perigee is along the X axis,
1565 // i.e cos(ω + Ω) = h/ecc 1 and sin(ω + Ω) = k/ecc = 0
1566 dcdh = auxiliaryElements.getEcc().multiply(dbdh).add(context.getb());
1567 dcdk = auxiliaryElements.getEcc().multiply(dbdk);
1568 } else {
1569 dcdh = auxiliaryElements.getEcc().multiply(dbdh).add(context.getb().multiply(auxiliaryElements.getH()).divide(auxiliaryElements.getEcc()));
1570 dcdk = auxiliaryElements.getEcc().multiply(dbdk).add(context.getb().multiply(auxiliaryElements.getK()).divide(auxiliaryElements.getEcc()));
1571 }
1572
1573 //Compute the powers (1 - c²)<sup>n</sup> and (1 + c²)<sup>n</sup>
1574 omc2tn = MathArrays.buildArray(field, context.getMaxAR3Pow() + context.getMaxFreqF() + 2);
1575 opc2tn = MathArrays.buildArray(field, context.getMaxAR3Pow() + context.getMaxFreqF() + 2);
1576 final T omc2 = c2.negate().add(1.);
1577 final T opc2 = c2.add(1.);
1578 omc2tn[0] = zero.add(1.);
1579 opc2tn[0] = zero.add(1.);
1580 for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF() + 1; i++) {
1581 omc2tn[i] = omc2tn[i - 1].multiply(omc2);
1582 opc2tn[i] = opc2tn[i - 1].multiply(opc2);
1583 }
1584
1585 //Compute the powers of b
1586 btjms = MathArrays.buildArray(field, context.getMaxAR3Pow() + context.getMaxFreqF() + 1);
1587 btjms[0] = zero.add(1.);
1588 for (int i = 1; i <= context.getMaxAR3Pow() + context.getMaxFreqF(); i++) {
1589 btjms[i] = btjms[i - 1].multiply(context.getb());
1590 }
1591 }
1592
1593 /** Compute the value of the coefficient e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by h and k. <br />
1594 *
1595 * @param j j index
1596 * @param s s index
1597 * @param n n index
1598 * @param context container for attributes
1599 * @param field field used by default
1600 * @return an array containing the value of the coefficient at index 0, the derivative by k at index 1 and the derivative by h at index 2
1601 */
1602 public T[] computeWjnsEmjmsAndDeriv(final int j, final int s, final int n,
1603 final FieldDSSTThirdBodyContext<T> context,
1604 final Field<T> field) {
1605 //Zero
1606 final T zero = field.getZero();
1607
1608 final T[] wjnsemjms = MathArrays.buildArray(field, 3);
1609 Arrays.fill(wjnsemjms, zero);
1610
1611 // |j|
1612 final int absJ = FastMath.abs(j);
1613 // |s|
1614 final int absS = FastMath.abs(s);
1615 // |j - s|
1616 final int absJmS = FastMath.abs(j - s);
1617 // |j + s|
1618 final int absJpS = FastMath.abs(j + s);
1619
1620 //The lower index of P. Also the power of (1 - c²)
1621 final int l;
1622 // the factorial ratio coefficient or 1. if |s| <= |j|
1623 final T factCoef;
1624 if (absS > absJ) {
1625 //factCoef = (fact[n + s] / fact[n + j]) * (fact[n - s] / fact[n - j]);
1626 factCoef = zero.add((CombinatoricsUtils.factorialDouble(n + s) / CombinatoricsUtils.factorialDouble(n + j)) * (CombinatoricsUtils.factorialDouble(n - s) / CombinatoricsUtils.factorialDouble(n - j)));
1627 l = n - absS;
1628 } else {
1629 factCoef = zero.add(1.);
1630 l = n - absJ;
1631 }
1632
1633 // (-1)<sup>|j-s|</sup>
1634 final T sign = absJmS % 2 != 0 ? zero.add(-1.) : zero.add(1.);
1635 //(1 - c²)<sup>n-|s|</sup> / (1 + c²)<sup>n</sup>
1636 final T coef1 = omc2tn[l].divide(opc2tn[n]);
1637 //-b<sup>|j-s|</sup>
1638 final T coef2 = btjms[absJmS].multiply(sign);
1639 // P<sub>l</sub><sup>|j-s|, |j+s|</sup>(χ)
1640 final FieldGradient<T> jac =
1641 JacobiPolynomials.getValue(l, absJmS, absJpS, FieldGradient.variable(1, 0, context.getX()));
1642
1643 // the derivative of coef1 by c
1644 final T dcoef1dc = coef1.negate().multiply(2.).multiply(c).multiply(opc2tn[1].reciprocal().multiply(n).add(omc2tn[1].reciprocal().multiply(l)));
1645 // the derivative of coef1 by h
1646 final T dcoef1dh = dcoef1dc.multiply(dcdh);
1647 // the derivative of coef1 by k
1648 final T dcoef1dk = dcoef1dc.multiply(dcdk);
1649
1650 // the derivative of coef2 by b
1651 final T dcoef2db = absJmS == 0 ? zero : sign.multiply(absJmS).multiply(btjms[absJmS - 1]);
1652 // the derivative of coef2 by h
1653 final T dcoef2dh = dcoef2db.multiply(dbdh);
1654 // the derivative of coef2 by k
1655 final T dcoef2dk = dcoef2db.multiply(dbdk);
1656
1657 // the jacobi polynomial value
1658 final T jacobi = jac.getValue();
1659 // the derivative of the Jacobi polynomial by h
1660 final T djacobidh = jac.getGradient()[0].multiply(context.getHXXX());
1661 // the derivative of the Jacobi polynomial by k
1662 final T djacobidk = jac.getGradient()[0].multiply(context.getKXXX());
1663
1664 //group the above coefficients to limit the mathematical operations
1665 final T term1 = factCoef.multiply(coef1).multiply(coef2);
1666 final T term2 = factCoef.multiply(coef1).multiply(jacobi);
1667 final T term3 = factCoef.multiply(coef2).multiply(jacobi);
1668
1669 //compute e<sup>-|j-s|</sup>*w<sub>j</sub><sup>n, s</sup> and its derivatives by k and h
1670 wjnsemjms[0] = term1.multiply(jacobi);
1671 wjnsemjms[1] = dcoef1dk.multiply(term3).add(dcoef2dk.multiply(term2)).add(djacobidk.multiply(term1));
1672 wjnsemjms[2] = dcoef1dh.multiply(term3).add(dcoef2dh.multiply(term2)).add(djacobidh.multiply(term1));
1673
1674 return wjnsemjms;
1675 }
1676 }
1677
1678 /** The G<sub>n,s</sub> coefficients and their derivatives.
1679 * <p>
1680 * See Danielson, 4.2-17
1681 *
1682 * @author Lucian Barbulescu
1683 */
1684 private class GnsCoefficients {
1685
1686 /** Maximum value for n index. */
1687 private final int nMax;
1688
1689 /** Maximum value for s index. */
1690 private final int sMax;
1691
1692 /** The coefficients G<sub>n,s</sub>. */
1693 private final double gns[][];
1694
1695 /** The derivatives of the coefficients G<sub>n,s</sub> by a. */
1696 private final double dgnsda[][];
1697
1698 /** The derivatives of the coefficients G<sub>n,s</sub> by γ. */
1699 private final double dgnsdgamma[][];
1700
1701 /** Standard constructor.
1702 *
1703 * @param nMax maximum value for n indes
1704 * @param sMax maximum value for s index
1705 * @param context container for attributes
1706 */
1707 GnsCoefficients(final int nMax, final int sMax, final DSSTThirdBodyContext context) {
1708 this.nMax = nMax;
1709 this.sMax = sMax;
1710
1711 final int rows = nMax + 1;
1712 final int columns = sMax + 1;
1713 this.gns = new double[rows][columns];
1714 this.dgnsda = new double[rows][columns];
1715 this.dgnsdgamma = new double[rows][columns];
1716
1717 // Generate the coefficients
1718 generateCoefficients(context);
1719 }
1720 /**
1721 * Compute the coefficient G<sub>n,s</sub> and its derivatives.
1722 * <p>
1723 * Only the derivatives by a and γ are computed as all others are 0
1724 * </p>
1725 * @param context container for attributes
1726 */
1727 private void generateCoefficients(final DSSTThirdBodyContext context) {
1728
1729 final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
1730
1731 for (int s = 0; s <= sMax; s++) {
1732 // The n index is always at least the maximum between 2 and s
1733 final int minN = FastMath.max(2, s);
1734 for (int n = minN; n <= nMax; n++) {
1735 // compute the coefficients only if (n - s) % 2 == 0
1736 if ( (n - s) % 2 == 0 ) {
1737 // Kronecker symbol (2 - delta(0,s))
1738 final double delta0s = (s == 0) ? 1. : 2.;
1739 final double vns = Vns.get(new NSKey(n, s));
1740 final double coef0 = delta0s * context.getAoR3Pow()[n] * vns * context.getMuoR3();
1741 final double coef1 = coef0 * context.getQns()[n][s];
1742 // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
1743 // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
1744 final double dqns = (n == s) ? 0. : context.getQns()[n][s + 1];
1745
1746 //Compute the coefficient and its derivatives.
1747 this.gns[n][s] = coef1;
1748 this.dgnsda[n][s] = coef1 * n / auxiliaryElements.getSma();
1749 this.dgnsdgamma[n][s] = coef0 * dqns;
1750 } else {
1751 // the coefficient and its derivatives is 0
1752 this.gns[n][s] = 0.;
1753 this.dgnsda[n][s] = 0.;
1754 this.dgnsdgamma[n][s] = 0.;
1755 }
1756 }
1757 }
1758 }
1759
1760 /** Get the coefficient G<sub>n,s</sub>.
1761 *
1762 * @param n n index
1763 * @param s s index
1764 * @return the coefficient G<sub>n,s</sub>
1765 */
1766 public double getGns(final int n, final int s) {
1767 return this.gns[n][s];
1768 }
1769
1770 /** Get the derivative dG<sub>n,s</sub> / da.
1771 *
1772 * @param n n index
1773 * @param s s index
1774 * @return the derivative dG<sub>n,s</sub> / da
1775 */
1776 public double getdGnsda(final int n, final int s) {
1777 return this.dgnsda[n][s];
1778 }
1779
1780 /** Get the derivative dG<sub>n,s</sub> / dγ.
1781 *
1782 * @param n n index
1783 * @param s s index
1784 * @return the derivative dG<sub>n,s</sub> / dγ
1785 */
1786 public double getdGnsdgamma(final int n, final int s) {
1787 return this.dgnsdgamma[n][s];
1788 }
1789 }
1790
1791 /** The G<sub>n,s</sub> coefficients and their derivatives.
1792 * <p>
1793 * See Danielson, 4.2-17
1794 *
1795 * @author Lucian Barbulescu
1796 */
1797 private class FieldGnsCoefficients <T extends CalculusFieldElement<T>> {
1798
1799 /** Maximum value for n index. */
1800 private final int nMax;
1801
1802 /** Maximum value for s index. */
1803 private final int sMax;
1804
1805 /** The coefficients G<sub>n,s</sub>. */
1806 private final T gns[][];
1807
1808 /** The derivatives of the coefficients G<sub>n,s</sub> by a. */
1809 private final T dgnsda[][];
1810
1811 /** The derivatives of the coefficients G<sub>n,s</sub> by γ. */
1812 private final T dgnsdgamma[][];
1813
1814 /** Standard constructor.
1815 *
1816 * @param nMax maximum value for n indes
1817 * @param sMax maximum value for s index
1818 * @param context container for attributes
1819 * @param field field used by default
1820 */
1821 FieldGnsCoefficients(final int nMax, final int sMax,
1822 final FieldDSSTThirdBodyContext<T> context,
1823 final Field<T> field) {
1824 this.nMax = nMax;
1825 this.sMax = sMax;
1826
1827 final int rows = nMax + 1;
1828 final int columns = sMax + 1;
1829 this.gns = MathArrays.buildArray(field, rows, columns);
1830 this.dgnsda = MathArrays.buildArray(field, rows, columns);
1831 this.dgnsdgamma = MathArrays.buildArray(field, rows, columns);
1832
1833 // Generate the coefficients
1834 generateCoefficients(context, field);
1835 }
1836 /**
1837 * Compute the coefficient G<sub>n,s</sub> and its derivatives.
1838 * <p>
1839 * Only the derivatives by a and γ are computed as all others are 0
1840 * </p>
1841 * @param context container for attributes
1842 * @param field field used by default
1843 */
1844 private void generateCoefficients(final FieldDSSTThirdBodyContext<T> context,
1845 final Field<T> field) {
1846
1847 //Zero
1848 final T zero = field.getZero();
1849
1850 final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
1851
1852 for (int s = 0; s <= sMax; s++) {
1853 // The n index is always at least the maximum between 2 and s
1854 final int minN = FastMath.max(2, s);
1855 for (int n = minN; n <= nMax; n++) {
1856 // compute the coefficients only if (n - s) % 2 == 0
1857 if ( (n - s) % 2 == 0 ) {
1858 // Kronecker symbol (2 - delta(0,s))
1859 final T delta0s = (s == 0) ? zero.add(1.) : zero.add(2.);
1860 final double vns = Vns.get(new NSKey(n, s));
1861 final T coef0 = context.getAoR3Pow()[n].multiply(vns).multiply(context.getMuoR3()).multiply(delta0s);
1862 final T coef1 = coef0.multiply(context.getQns()[n][s]);
1863 // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
1864 // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
1865 final T dqns = (n == s) ? zero : context.getQns()[n][s + 1];
1866
1867 //Compute the coefficient and its derivatives.
1868 this.gns[n][s] = coef1;
1869 this.dgnsda[n][s] = coef1.multiply(n).divide(auxiliaryElements.getSma());
1870 this.dgnsdgamma[n][s] = coef0.multiply(dqns);
1871 } else {
1872 // the coefficient and its derivatives is 0
1873 this.gns[n][s] = zero;
1874 this.dgnsda[n][s] = zero;
1875 this.dgnsdgamma[n][s] = zero;
1876 }
1877 }
1878 }
1879 }
1880
1881 /** Get the coefficient G<sub>n,s</sub>.
1882 *
1883 * @param n n index
1884 * @param s s index
1885 * @return the coefficient G<sub>n,s</sub>
1886 */
1887 public T getGns(final int n, final int s) {
1888 return this.gns[n][s];
1889 }
1890
1891 /** Get the derivative dG<sub>n,s</sub> / da.
1892 *
1893 * @param n n index
1894 * @param s s index
1895 * @return the derivative dG<sub>n,s</sub> / da
1896 */
1897 public T getdGnsda(final int n, final int s) {
1898 return this.dgnsda[n][s];
1899 }
1900
1901 /** Get the derivative dG<sub>n,s</sub> / dγ.
1902 *
1903 * @param n n index
1904 * @param s s index
1905 * @return the derivative dG<sub>n,s</sub> / dγ
1906 */
1907 public T getdGnsdgamma(final int n, final int s) {
1908 return this.dgnsdgamma[n][s];
1909 }
1910 }
1911
1912 /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
1913 *
1914 * <p>
1915 * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
1916 * - sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) <br />
1917 * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
1918 * - C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) <br />
1919 * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
1920 * For the ease of usage the above terms are renamed A<sub>js</sub>, B<sub>js</sub>, D<sub>js</sub> and E<sub>js</sub> respectively <br />
1921 * See the CS Mathematical report $3.5.3.2 for more details
1922 * </p>
1923 * @author Lucian Barbulescu
1924 */
1925 private static class CjSjAlphaBetaKH {
1926
1927 /** The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) series. */
1928 private final CjSjCoefficient cjsjkh;
1929
1930 /** The C<sub>j</sub>(α, β) and the S<sub>j</sub>(α, β) series. */
1931 private final CjSjCoefficient cjsjalbe;
1932
1933 /** The coeficient sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h)
1934 * and its derivative by k, h, α and β. */
1935 private final double coefAandDeriv[];
1936
1937 /** The coeficient C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h)
1938 * and its derivative by k, h, α and β. */
1939 private final double coefBandDeriv[];
1940
1941 /** The coeficient C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h)
1942 * and its derivative by k, h, α and β. */
1943 private final double coefDandDeriv[];
1944
1945 /** The coeficient C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h)
1946 * and its derivative by k, h, α and β. */
1947 private final double coefEandDeriv[];
1948
1949 /**
1950 * Standard constructor.
1951 * @param context container for attributes
1952 */
1953 CjSjAlphaBetaKH(final DSSTThirdBodyContext context) {
1954
1955 final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
1956
1957 cjsjkh = new CjSjCoefficient(auxiliaryElements.getK(), auxiliaryElements.getH());
1958 cjsjalbe = new CjSjCoefficient(context.getAlpha(), context.getBeta());
1959
1960 coefAandDeriv = new double[5];
1961 coefBandDeriv = new double[5];
1962 coefDandDeriv = new double[5];
1963 coefEandDeriv = new double[5];
1964 }
1965
1966 /** Compute the coefficients and their derivatives for a given (j,s) pair.
1967 * @param j j index
1968 * @param s s index
1969 */
1970 public void computeCoefficients(final int j, final int s) {
1971 // sign of j-s
1972 final int sign = j < s ? -1 : 1;
1973
1974 //|j-s|
1975 final int absJmS = FastMath.abs(j - s);
1976
1977 //j+s
1978 final int jps = j + s;
1979
1980 //Compute the coefficient A and its derivatives
1981 coefAandDeriv[0] = sign * cjsjalbe.getCj(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getCj(absJmS);
1982 coefAandDeriv[1] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDk(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDk(absJmS);
1983 coefAandDeriv[2] = sign * cjsjalbe.getCj(s) * cjsjkh.getDsjDh(absJmS) + cjsjalbe.getSj(s) * cjsjkh.getDcjDh(absJmS);
1984 coefAandDeriv[3] = sign * cjsjalbe.getDcjDk(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDk(s) * cjsjkh.getCj(absJmS);
1985 coefAandDeriv[4] = sign * cjsjalbe.getDcjDh(s) * cjsjkh.getSj(absJmS) + cjsjalbe.getDsjDh(s) * cjsjkh.getCj(absJmS);
1986
1987 //Compute the coefficient B and its derivatives
1988 coefBandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getSj(jps) - cjsjalbe.getSj(s) * cjsjkh.getCj(jps);
1989 coefBandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDsjDk(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDk(jps);
1990 coefBandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDsjDh(jps) - cjsjalbe.getSj(s) * cjsjkh.getDcjDh(jps);
1991 coefBandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDk(s) * cjsjkh.getCj(jps);
1992 coefBandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getSj(jps) - cjsjalbe.getDsjDh(s) * cjsjkh.getCj(jps);
1993
1994 //Compute the coefficient D and its derivatives
1995 coefDandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getSj(absJmS);
1996 coefDandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDk(absJmS);
1997 coefDandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(absJmS) - sign * cjsjalbe.getSj(s) * cjsjkh.getDsjDh(absJmS);
1998 coefDandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDk(s) * cjsjkh.getSj(absJmS);
1999 coefDandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(absJmS) - sign * cjsjalbe.getDsjDh(s) * cjsjkh.getSj(absJmS);
2000
2001 //Compute the coefficient E and its derivatives
2002 coefEandDeriv[0] = cjsjalbe.getCj(s) * cjsjkh.getCj(jps) + cjsjalbe.getSj(s) * cjsjkh.getSj(jps);
2003 coefEandDeriv[1] = cjsjalbe.getCj(s) * cjsjkh.getDcjDk(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDk(jps);
2004 coefEandDeriv[2] = cjsjalbe.getCj(s) * cjsjkh.getDcjDh(jps) + cjsjalbe.getSj(s) * cjsjkh.getDsjDh(jps);
2005 coefEandDeriv[3] = cjsjalbe.getDcjDk(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDk(s) * cjsjkh.getSj(jps);
2006 coefEandDeriv[4] = cjsjalbe.getDcjDh(s) * cjsjkh.getCj(jps) + cjsjalbe.getDsjDh(s) * cjsjkh.getSj(jps);
2007 }
2008
2009 /** Get the value of coefficient A<sub>j,s</sub>.
2010 *
2011 * @return the coefficient A<sub>j,s</sub>
2012 */
2013 public double getCoefA() {
2014 return coefAandDeriv[0];
2015 }
2016
2017 /** Get the value of coefficient dA<sub>j,s</sub>/dk.
2018 *
2019 * @return the coefficient dA<sub>j,s</sub>/dk
2020 */
2021 public double getdCoefAdk() {
2022 return coefAandDeriv[1];
2023 }
2024
2025 /** Get the value of coefficient dA<sub>j,s</sub>/dh.
2026 *
2027 * @return the coefficient dA<sub>j,s</sub>/dh
2028 */
2029 public double getdCoefAdh() {
2030 return coefAandDeriv[2];
2031 }
2032
2033 /** Get the value of coefficient dA<sub>j,s</sub>/dα.
2034 *
2035 * @return the coefficient dA<sub>j,s</sub>/dα
2036 */
2037 public double getdCoefAdalpha() {
2038 return coefAandDeriv[3];
2039 }
2040
2041 /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
2042 *
2043 * @return the coefficient dA<sub>j,s</sub>/dβ
2044 */
2045 public double getdCoefAdbeta() {
2046 return coefAandDeriv[4];
2047 }
2048
2049 /** Get the value of coefficient B<sub>j,s</sub>.
2050 *
2051 * @return the coefficient B<sub>j,s</sub>
2052 */
2053 public double getCoefB() {
2054 return coefBandDeriv[0];
2055 }
2056
2057 /** Get the value of coefficient dB<sub>j,s</sub>/dk.
2058 *
2059 * @return the coefficient dB<sub>j,s</sub>/dk
2060 */
2061 public double getdCoefBdk() {
2062 return coefBandDeriv[1];
2063 }
2064
2065 /** Get the value of coefficient dB<sub>j,s</sub>/dh.
2066 *
2067 * @return the coefficient dB<sub>j,s</sub>/dh
2068 */
2069 public double getdCoefBdh() {
2070 return coefBandDeriv[2];
2071 }
2072
2073 /** Get the value of coefficient dB<sub>j,s</sub>/dα.
2074 *
2075 * @return the coefficient dB<sub>j,s</sub>/dα
2076 */
2077 public double getdCoefBdalpha() {
2078 return coefBandDeriv[3];
2079 }
2080
2081 /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
2082 *
2083 * @return the coefficient dB<sub>j,s</sub>/dβ
2084 */
2085 public double getdCoefBdbeta() {
2086 return coefBandDeriv[4];
2087 }
2088
2089 /** Get the value of coefficient D<sub>j,s</sub>.
2090 *
2091 * @return the coefficient D<sub>j,s</sub>
2092 */
2093 public double getCoefD() {
2094 return coefDandDeriv[0];
2095 }
2096
2097 /** Get the value of coefficient dD<sub>j,s</sub>/dk.
2098 *
2099 * @return the coefficient dD<sub>j,s</sub>/dk
2100 */
2101 public double getdCoefDdk() {
2102 return coefDandDeriv[1];
2103 }
2104
2105 /** Get the value of coefficient dD<sub>j,s</sub>/dh.
2106 *
2107 * @return the coefficient dD<sub>j,s</sub>/dh
2108 */
2109 public double getdCoefDdh() {
2110 return coefDandDeriv[2];
2111 }
2112
2113 /** Get the value of coefficient dD<sub>j,s</sub>/dα.
2114 *
2115 * @return the coefficient dD<sub>j,s</sub>/dα
2116 */
2117 public double getdCoefDdalpha() {
2118 return coefDandDeriv[3];
2119 }
2120
2121 /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
2122 *
2123 * @return the coefficient dD<sub>j,s</sub>/dβ
2124 */
2125 public double getdCoefDdbeta() {
2126 return coefDandDeriv[4];
2127 }
2128
2129 /** Get the value of coefficient E<sub>j,s</sub>.
2130 *
2131 * @return the coefficient E<sub>j,s</sub>
2132 */
2133 public double getCoefE() {
2134 return coefEandDeriv[0];
2135 }
2136
2137 /** Get the value of coefficient dE<sub>j,s</sub>/dk.
2138 *
2139 * @return the coefficient dE<sub>j,s</sub>/dk
2140 */
2141 public double getdCoefEdk() {
2142 return coefEandDeriv[1];
2143 }
2144
2145 /** Get the value of coefficient dE<sub>j,s</sub>/dh.
2146 *
2147 * @return the coefficient dE<sub>j,s</sub>/dh
2148 */
2149 public double getdCoefEdh() {
2150 return coefEandDeriv[2];
2151 }
2152
2153 /** Get the value of coefficient dE<sub>j,s</sub>/dα.
2154 *
2155 * @return the coefficient dE<sub>j,s</sub>/dα
2156 */
2157 public double getdCoefEdalpha() {
2158 return coefEandDeriv[3];
2159 }
2160
2161 /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
2162 *
2163 * @return the coefficient dE<sub>j,s</sub>/dβ
2164 */
2165 public double getdCoefEdbeta() {
2166 return coefEandDeriv[4];
2167 }
2168 }
2169
2170 /** This class computes the terms containing the coefficients C<sub>j</sub> and S<sub>j</sub> of (α, β) or (k, h).
2171 *
2172 * <p>
2173 * The following terms and their derivatives by k, h, alpha and beta are considered: <br/ >
2174 * - sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) <br />
2175 * - C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) <br />
2176 * - C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) <br />
2177 * - C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) <br />
2178 * For the ease of usage the above terms are renamed A<sub>js</sub>, B<sub>js</sub>, D<sub>js</sub> and E<sub>js</sub> respectively <br />
2179 * See the CS Mathematical report $3.5.3.2 for more details
2180 * </p>
2181 * @author Lucian Barbulescu
2182 */
2183 private static class FieldCjSjAlphaBetaKH <T extends CalculusFieldElement<T>> {
2184
2185 /** The C<sub>j</sub>(k, h) and the S<sub>j</sub>(k, h) series. */
2186 private final FieldCjSjCoefficient<T> cjsjkh;
2187
2188 /** The C<sub>j</sub>(α, β) and the S<sub>j</sub>(α, β) series. */
2189 private final FieldCjSjCoefficient<T> cjsjalbe;
2190
2191 /** The coeficient sign(j-s) * C<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h) + S<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h)
2192 * and its derivative by k, h, α and β. */
2193 private final T coefAandDeriv[];
2194
2195 /** The coeficient C<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h) - S<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h)
2196 * and its derivative by k, h, α and β. */
2197 private final T coefBandDeriv[];
2198
2199 /** The coeficient C<sub>s</sub>(α, β) * C<sub>|j-s|</sub>(k, h) - sign(j-s) * S<sub>s</sub>(α, β) * S<sub>|j-s|</sub>(k, h)
2200 * and its derivative by k, h, α and β. */
2201 private final T coefDandDeriv[];
2202
2203 /** The coeficient C<sub>s</sub>(α, β) * C<sub>j+s</sub>(k, h) + S<sub>s</sub>(α, β) * S<sub>j+s</sub>(k, h)
2204 * and its derivative by k, h, α and β. */
2205 private final T coefEandDeriv[];
2206
2207 /**
2208 * Standard constructor.
2209 * @param context container for attributes
2210 * @param field field used by default
2211 */
2212 FieldCjSjAlphaBetaKH(final FieldDSSTThirdBodyContext<T> context, final Field<T> field) {
2213
2214 final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
2215
2216 cjsjkh = new FieldCjSjCoefficient<>(auxiliaryElements.getK(), auxiliaryElements.getH(), field);
2217 cjsjalbe = new FieldCjSjCoefficient<>(context.getAlpha(), context.getBeta(), field);
2218
2219 coefAandDeriv = MathArrays.buildArray(field, 5);
2220 coefBandDeriv = MathArrays.buildArray(field, 5);
2221 coefDandDeriv = MathArrays.buildArray(field, 5);
2222 coefEandDeriv = MathArrays.buildArray(field, 5);
2223 }
2224
2225 /** Compute the coefficients and their derivatives for a given (j,s) pair.
2226 * @param j j index
2227 * @param s s index
2228 */
2229 public void computeCoefficients(final int j, final int s) {
2230 // sign of j-s
2231 final int sign = j < s ? -1 : 1;
2232
2233 //|j-s|
2234 final int absJmS = FastMath.abs(j - s);
2235
2236 //j+s
2237 final int jps = j + s;
2238
2239 //Compute the coefficient A and its derivatives
2240 coefAandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getCj(absJmS)));
2241 coefAandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDk(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDk(absJmS)));
2242 coefAandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDh(absJmS)).multiply(sign).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDh(absJmS)));
2243 coefAandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getCj(absJmS)));
2244 coefAandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign).add(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getCj(absJmS)));
2245
2246 //Compute the coefficient B and its derivatives
2247 coefBandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getCj(jps)));
2248 coefBandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDk(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDk(jps)));
2249 coefBandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDsjDh(jps)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDcjDh(jps)));
2250 coefBandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getCj(jps)));
2251 coefBandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getSj(jps)).subtract(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getCj(jps)));
2252
2253 //Compute the coefficient D and its derivatives
2254 coefDandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
2255 coefDandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDk(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDk(absJmS)).multiply(sign));
2256 coefDandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDh(absJmS)).subtract(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDh(absJmS)).multiply(sign));
2257 coefDandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
2258 coefDandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getCj(absJmS)).subtract(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getSj(absJmS)).multiply(sign));
2259
2260 //Compute the coefficient E and its derivatives
2261 coefEandDeriv[0] = cjsjalbe.getCj(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getSj(jps)));
2262 coefEandDeriv[1] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDk(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDk(jps)));
2263 coefEandDeriv[2] = cjsjalbe.getCj(s).multiply(cjsjkh.getDcjDh(jps)).add(cjsjalbe.getSj(s).multiply(cjsjkh.getDsjDh(jps)));
2264 coefEandDeriv[3] = cjsjalbe.getDcjDk(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getDsjDk(s).multiply(cjsjkh.getSj(jps)));
2265 coefEandDeriv[4] = cjsjalbe.getDcjDh(s).multiply(cjsjkh.getCj(jps)).add(cjsjalbe.getDsjDh(s).multiply(cjsjkh.getSj(jps)));
2266 }
2267
2268 /** Get the value of coefficient A<sub>j,s</sub>.
2269 *
2270 * @return the coefficient A<sub>j,s</sub>
2271 */
2272 public T getCoefA() {
2273 return coefAandDeriv[0];
2274 }
2275
2276 /** Get the value of coefficient dA<sub>j,s</sub>/dk.
2277 *
2278 * @return the coefficient dA<sub>j,s</sub>/dk
2279 */
2280 public T getdCoefAdk() {
2281 return coefAandDeriv[1];
2282 }
2283
2284 /** Get the value of coefficient dA<sub>j,s</sub>/dh.
2285 *
2286 * @return the coefficient dA<sub>j,s</sub>/dh
2287 */
2288 public T getdCoefAdh() {
2289 return coefAandDeriv[2];
2290 }
2291
2292 /** Get the value of coefficient dA<sub>j,s</sub>/dα.
2293 *
2294 * @return the coefficient dA<sub>j,s</sub>/dα
2295 */
2296 public T getdCoefAdalpha() {
2297 return coefAandDeriv[3];
2298 }
2299
2300 /** Get the value of coefficient dA<sub>j,s</sub>/dβ.
2301 *
2302 * @return the coefficient dA<sub>j,s</sub>/dβ
2303 */
2304 public T getdCoefAdbeta() {
2305 return coefAandDeriv[4];
2306 }
2307
2308 /** Get the value of coefficient B<sub>j,s</sub>.
2309 *
2310 * @return the coefficient B<sub>j,s</sub>
2311 */
2312 public T getCoefB() {
2313 return coefBandDeriv[0];
2314 }
2315
2316 /** Get the value of coefficient dB<sub>j,s</sub>/dk.
2317 *
2318 * @return the coefficient dB<sub>j,s</sub>/dk
2319 */
2320 public T getdCoefBdk() {
2321 return coefBandDeriv[1];
2322 }
2323
2324 /** Get the value of coefficient dB<sub>j,s</sub>/dh.
2325 *
2326 * @return the coefficient dB<sub>j,s</sub>/dh
2327 */
2328 public T getdCoefBdh() {
2329 return coefBandDeriv[2];
2330 }
2331
2332 /** Get the value of coefficient dB<sub>j,s</sub>/dα.
2333 *
2334 * @return the coefficient dB<sub>j,s</sub>/dα
2335 */
2336 public T getdCoefBdalpha() {
2337 return coefBandDeriv[3];
2338 }
2339
2340 /** Get the value of coefficient dB<sub>j,s</sub>/dβ.
2341 *
2342 * @return the coefficient dB<sub>j,s</sub>/dβ
2343 */
2344 public T getdCoefBdbeta() {
2345 return coefBandDeriv[4];
2346 }
2347
2348 /** Get the value of coefficient D<sub>j,s</sub>.
2349 *
2350 * @return the coefficient D<sub>j,s</sub>
2351 */
2352 public T getCoefD() {
2353 return coefDandDeriv[0];
2354 }
2355
2356 /** Get the value of coefficient dD<sub>j,s</sub>/dk.
2357 *
2358 * @return the coefficient dD<sub>j,s</sub>/dk
2359 */
2360 public T getdCoefDdk() {
2361 return coefDandDeriv[1];
2362 }
2363
2364 /** Get the value of coefficient dD<sub>j,s</sub>/dh.
2365 *
2366 * @return the coefficient dD<sub>j,s</sub>/dh
2367 */
2368 public T getdCoefDdh() {
2369 return coefDandDeriv[2];
2370 }
2371
2372 /** Get the value of coefficient dD<sub>j,s</sub>/dα.
2373 *
2374 * @return the coefficient dD<sub>j,s</sub>/dα
2375 */
2376 public T getdCoefDdalpha() {
2377 return coefDandDeriv[3];
2378 }
2379
2380 /** Get the value of coefficient dD<sub>j,s</sub>/dβ.
2381 *
2382 * @return the coefficient dD<sub>j,s</sub>/dβ
2383 */
2384 public T getdCoefDdbeta() {
2385 return coefDandDeriv[4];
2386 }
2387
2388 /** Get the value of coefficient E<sub>j,s</sub>.
2389 *
2390 * @return the coefficient E<sub>j,s</sub>
2391 */
2392 public T getCoefE() {
2393 return coefEandDeriv[0];
2394 }
2395
2396 /** Get the value of coefficient dE<sub>j,s</sub>/dk.
2397 *
2398 * @return the coefficient dE<sub>j,s</sub>/dk
2399 */
2400 public T getdCoefEdk() {
2401 return coefEandDeriv[1];
2402 }
2403
2404 /** Get the value of coefficient dE<sub>j,s</sub>/dh.
2405 *
2406 * @return the coefficient dE<sub>j,s</sub>/dh
2407 */
2408 public T getdCoefEdh() {
2409 return coefEandDeriv[2];
2410 }
2411
2412 /** Get the value of coefficient dE<sub>j,s</sub>/dα.
2413 *
2414 * @return the coefficient dE<sub>j,s</sub>/dα
2415 */
2416 public T getdCoefEdalpha() {
2417 return coefEandDeriv[3];
2418 }
2419
2420 /** Get the value of coefficient dE<sub>j,s</sub>/dβ.
2421 *
2422 * @return the coefficient dE<sub>j,s</sub>/dβ
2423 */
2424 public T getdCoefEdbeta() {
2425 return coefEandDeriv[4];
2426 }
2427 }
2428
2429 /** This class computes the coefficients for the generating function S and its derivatives.
2430 * <p>
2431 * The form of the generating functions is: <br>
2432 * S = C⁰ + Σ<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
2433 * The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
2434 * presented in Danielson 4.2-14,15 except for the case j=1 where
2435 * C¹ = C¹<sub>Fourier</sub> - hU and
2436 * S¹ = S¹<sub>Fourier</sub> + kU <br>
2437 * Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
2438 * are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
2439 * </p>
2440 * @author Lucian Barbulescu
2441 */
2442 private class GeneratingFunctionCoefficients {
2443
2444 /** The Fourier coefficients as presented in Danielson 4.2-14,15. */
2445 private final FourierCjSjCoefficients cjsjFourier;
2446
2447 /** Maximum value of j index. */
2448 private final int jMax;
2449
2450 /** The coefficients C<sup>j</sup> of the function S and its derivatives.
2451 * <p>
2452 * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
2453 * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
2454 * - S <br/>
2455 * - dS / da <br/>
2456 * - dS / dk <br/>
2457 * - dS / dh <br/>
2458 * - dS / dα <br/>
2459 * - dS / dβ <br/>
2460 * - dS / dγ <br/>
2461 * - dS / dλ
2462 * </p>
2463 */
2464 private final double[][] cjCoefs;
2465
2466 /** The coefficients S<sup>j</sup> of the function S and its derivatives.
2467 * <p>
2468 * The index j belongs to the interval [0,jMax].<br>
2469 * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
2470 * - S <br/>
2471 * - dS / da <br/>
2472 * - dS / dk <br/>
2473 * - dS / dh <br/>
2474 * - dS / dα <br/>
2475 * - dS / dβ <br/>
2476 * - dS / dγ <br/>
2477 * - dS / dλ
2478 * </p>
2479 */
2480 private final double[][] sjCoefs;
2481
2482 /**
2483 * Standard constructor.
2484 *
2485 * @param nMax maximum value of n index
2486 * @param sMax maximum value of s index
2487 * @param jMax maximum value of j index
2488 * @param context container for attributes
2489 * @param hansen hansen objects
2490 */
2491 GeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax,
2492 final DSSTThirdBodyContext context, final HansenObjects hansen) {
2493 this.jMax = jMax;
2494 this.cjsjFourier = new FourierCjSjCoefficients(nMax, sMax, jMax, context);
2495 this.cjCoefs = new double[8][jMax + 1];
2496 this.sjCoefs = new double[8][jMax + 1];
2497
2498 computeGeneratingFunctionCoefficients(context, hansen);
2499 }
2500
2501 /**
2502 * Compute the coefficients for the generating function S and its derivatives.
2503 * @param context container for attributes
2504 * @param hansenObjects hansen objects
2505 */
2506 private void computeGeneratingFunctionCoefficients(final DSSTThirdBodyContext context, final HansenObjects hansenObjects) {
2507
2508 final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
2509
2510 // Access to potential U derivatives
2511 final UAnddU udu = new UAnddU(context, hansenObjects);
2512
2513 //Compute the C<sup>j</sup> coefficients
2514 for (int j = 1; j <= jMax; j++) {
2515 //Compute the C<sup>j</sup> coefficients
2516 cjCoefs[0][j] = cjsjFourier.getCj(j);
2517 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
2518 cjCoefs[2][j] = cjsjFourier.getdCjdk(j) - (cjsjFourier.getSjLambda(j - 1) - cjsjFourier.getSjLambda(j + 1)) / 2;
2519 cjCoefs[3][j] = cjsjFourier.getdCjdh(j) - (cjsjFourier.getCjLambda(j - 1) + cjsjFourier.getCjLambda(j + 1)) / 2;
2520 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
2521 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
2522 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
2523 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);
2524
2525 //Compute the S<sup>j</sup> coefficients
2526 sjCoefs[0][j] = cjsjFourier.getSj(j);
2527 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
2528 sjCoefs[2][j] = cjsjFourier.getdSjdk(j) + (cjsjFourier.getCjLambda(j - 1) - cjsjFourier.getCjLambda(j + 1)) / 2;
2529 sjCoefs[3][j] = cjsjFourier.getdSjdh(j) - (cjsjFourier.getSjLambda(j - 1) + cjsjFourier.getSjLambda(j + 1)) / 2;
2530 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
2531 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
2532 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
2533 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);
2534
2535 //In the special case j == 1 there are some additional terms to be added
2536 if (j == 1) {
2537 //Additional terms for C<sup>j</sup> coefficients
2538 cjCoefs[0][j] += -auxiliaryElements.getH() * udu.getU();
2539 cjCoefs[1][j] += -auxiliaryElements.getH() * udu.getdUda();
2540 cjCoefs[2][j] += -auxiliaryElements.getH() * udu.getdUdk();
2541 cjCoefs[3][j] += -(auxiliaryElements.getH() * udu.getdUdh() + udu.getU() + cjsjFourier.getC0Lambda());
2542 cjCoefs[4][j] += -auxiliaryElements.getH() * udu.getdUdAl();
2543 cjCoefs[5][j] += -auxiliaryElements.getH() * udu.getdUdBe();
2544 cjCoefs[6][j] += -auxiliaryElements.getH() * udu.getdUdGa();
2545
2546 //Additional terms for S<sup>j</sup> coefficients
2547 sjCoefs[0][j] += auxiliaryElements.getK() * udu.getU();
2548 sjCoefs[1][j] += auxiliaryElements.getK() * udu.getdUda();
2549 sjCoefs[2][j] += auxiliaryElements.getK() * udu.getdUdk() + udu.getU() + cjsjFourier.getC0Lambda();
2550 sjCoefs[3][j] += auxiliaryElements.getK() * udu.getdUdh();
2551 sjCoefs[4][j] += auxiliaryElements.getK() * udu.getdUdAl();
2552 sjCoefs[5][j] += auxiliaryElements.getK() * udu.getdUdBe();
2553 sjCoefs[6][j] += auxiliaryElements.getK() * udu.getdUdGa();
2554 }
2555 }
2556 }
2557
2558 /** Get the coefficient C<sup>j</sup> for the function S.
2559 * <br>
2560 * Possible values for j are within the interval [0,jMax].
2561 * The value 0 is used to obtain the free coefficient C⁰
2562 * @param j j index
2563 * @return C<sup>j</sup> for the function S
2564 */
2565 public double getSCj(final int j) {
2566 return cjCoefs[0][j];
2567 }
2568
2569 /** Get the coefficient S<sup>j</sup> for the function S.
2570 * <br>
2571 * Possible values for j are within the interval [1,jMax].
2572 * @param j j index
2573 * @return S<sup>j</sup> for the function S
2574 */
2575 public double getSSj(final int j) {
2576 return sjCoefs[0][j];
2577 }
2578
2579 /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
2580 * <br>
2581 * Possible values for j are within the interval [0,jMax].
2582 * The value 0 is used to obtain the free coefficient C⁰
2583 * @param j j index
2584 * @return C<sup>j</sup> for the function dS/da
2585 */
2586 public double getdSdaCj(final int j) {
2587 return cjCoefs[1][j];
2588 }
2589
2590 /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
2591 * <br>
2592 * Possible values for j are within the interval [1,jMax].
2593 * @param j j index
2594 * @return S<sup>j</sup> for the derivative dS/da
2595 */
2596 public double getdSdaSj(final int j) {
2597 return sjCoefs[1][j];
2598 }
2599
2600 /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
2601 * <br>
2602 * Possible values for j are within the interval [0,jMax].
2603 * The value 0 is used to obtain the free coefficient C⁰
2604 * @param j j index
2605 * @return C<sup>j</sup> for the function dS/dk
2606 */
2607 public double getdSdkCj(final int j) {
2608 return cjCoefs[2][j];
2609 }
2610
2611 /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
2612 * <br>
2613 * Possible values for j are within the interval [1,jMax].
2614 * @param j j index
2615 * @return S<sup>j</sup> for the derivative dS/dk
2616 */
2617 public double getdSdkSj(final int j) {
2618 return sjCoefs[2][j];
2619 }
2620
2621 /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
2622 * <br>
2623 * Possible values for j are within the interval [0,jMax].
2624 * The value 0 is used to obtain the free coefficient C⁰
2625 * @param j j index
2626 * @return C<sup>j</sup> for the function dS/dh
2627 */
2628 public double getdSdhCj(final int j) {
2629 return cjCoefs[3][j];
2630 }
2631
2632 /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
2633 * <br>
2634 * Possible values for j are within the interval [1,jMax].
2635 * @param j j index
2636 * @return S<sup>j</sup> for the derivative dS/dh
2637 */
2638 public double getdSdhSj(final int j) {
2639 return sjCoefs[3][j];
2640 }
2641
2642 /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
2643 * <br>
2644 * Possible values for j are within the interval [0,jMax].
2645 * The value 0 is used to obtain the free coefficient C⁰
2646 * @param j j index
2647 * @return C<sup>j</sup> for the function dS/dα
2648 */
2649 public double getdSdalphaCj(final int j) {
2650 return cjCoefs[4][j];
2651 }
2652
2653 /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
2654 * <br>
2655 * Possible values for j are within the interval [1,jMax].
2656 * @param j j index
2657 * @return S<sup>j</sup> for the derivative dS/dα
2658 */
2659 public double getdSdalphaSj(final int j) {
2660 return sjCoefs[4][j];
2661 }
2662
2663 /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
2664 * <br>
2665 * Possible values for j are within the interval [0,jMax].
2666 * The value 0 is used to obtain the free coefficient C⁰
2667 * @param j j index
2668 * @return C<sup>j</sup> for the function dS/dβ
2669 */
2670 public double getdSdbetaCj(final int j) {
2671 return cjCoefs[5][j];
2672 }
2673
2674 /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
2675 * <br>
2676 * Possible values for j are within the interval [1,jMax].
2677 * @param j j index
2678 * @return S<sup>j</sup> for the derivative dS/dβ
2679 */
2680 public double getdSdbetaSj(final int j) {
2681 return sjCoefs[5][j];
2682 }
2683
2684 /** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
2685 * <br>
2686 * Possible values for j are within the interval [0,jMax].
2687 * The value 0 is used to obtain the free coefficient C⁰
2688 * @param j j index
2689 * @return C<sup>j</sup> for the function dS/dγ
2690 */
2691 public double getdSdgammaCj(final int j) {
2692 return cjCoefs[6][j];
2693 }
2694
2695 /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
2696 * <br>
2697 * Possible values for j are within the interval [1,jMax].
2698 * @param j j index
2699 * @return S<sup>j</sup> for the derivative dS/dγ
2700 */
2701 public double getdSdgammaSj(final int j) {
2702 return sjCoefs[6][j];
2703 }
2704
2705 /** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
2706 * <br>
2707 * Possible values for j are within the interval [0,jMax].
2708 * The value 0 is used to obtain the free coefficient C⁰
2709 * @param j j index
2710 * @return C<sup>j</sup> for the function dS/dλ
2711 */
2712 public double getdSdlambdaCj(final int j) {
2713 return cjCoefs[7][j];
2714 }
2715
2716 /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
2717 * <br>
2718 * Possible values for j are within the interval [1,jMax].
2719 * @param j j index
2720 * @return S<sup>j</sup> for the derivative dS/dλ
2721 */
2722 public double getdSdlambdaSj(final int j) {
2723 return sjCoefs[7][j];
2724 }
2725 }
2726
2727 /** This class computes the coefficients for the generating function S and its derivatives.
2728 * <p>
2729 * The form of the generating functions is: <br>
2730 * S = C⁰ + Σ<sub>j=1</sub><sup>N+1</sup>(C<sup>j</sup> * cos(jF) + S<sup>j</sup> * sin(jF)) <br>
2731 * The coefficients C⁰, C<sup>j</sup>, S<sup>j</sup> are the Fourrier coefficients
2732 * presented in Danielson 4.2-14,15 except for the case j=1 where
2733 * C¹ = C¹<sub>Fourier</sub> - hU and
2734 * S¹ = S¹<sub>Fourier</sub> + kU <br>
2735 * Also the coefficients of the derivatives of S by a, k, h, α, β, γ and λ
2736 * are computed end expressed in a similar manner. The formulas used are 4.2-19, 20, 23, 24
2737 * </p>
2738 * @author Lucian Barbulescu
2739 */
2740 private class FieldGeneratingFunctionCoefficients <T extends CalculusFieldElement<T>> {
2741
2742 /** The Fourier coefficients as presented in Danielson 4.2-14,15. */
2743 private final FieldFourierCjSjCoefficients<T> cjsjFourier;
2744
2745 /** Maximum value of j index. */
2746 private final int jMax;
2747
2748 /** The coefficients C<sup>j</sup> of the function S and its derivatives.
2749 * <p>
2750 * The index j belongs to the interval [0,jMax]. The coefficient C⁰ is the free coefficient.<br>
2751 * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
2752 * - S <br/>
2753 * - dS / da <br/>
2754 * - dS / dk <br/>
2755 * - dS / dh <br/>
2756 * - dS / dα <br/>
2757 * - dS / dβ <br/>
2758 * - dS / dγ <br/>
2759 * - dS / dλ
2760 * </p>
2761 */
2762 private final T[][] cjCoefs;
2763
2764 /** The coefficients S<sup>j</sup> of the function S and its derivatives.
2765 * <p>
2766 * The index j belongs to the interval [0,jMax].<br>
2767 * Each column of the matrix contains the coefficient corresponding to the following functions: <br/>
2768 * - S <br/>
2769 * - dS / da <br/>
2770 * - dS / dk <br/>
2771 * - dS / dh <br/>
2772 * - dS / dα <br/>
2773 * - dS / dβ <br/>
2774 * - dS / dγ <br/>
2775 * - dS / dλ
2776 * </p>
2777 */
2778 private final T[][] sjCoefs;
2779
2780 /**
2781 * Standard constructor.
2782 *
2783 * @param nMax maximum value of n index
2784 * @param sMax maximum value of s index
2785 * @param jMax maximum value of j index
2786 * @param context container for attributes
2787 * @param hansen hansen objects
2788 * @param field field used by default
2789 */
2790 FieldGeneratingFunctionCoefficients(final int nMax, final int sMax, final int jMax,
2791 final FieldDSSTThirdBodyContext<T> context,
2792 final FieldHansenObjects<T> hansen,
2793 final Field<T> field) {
2794 this.jMax = jMax;
2795 this.cjsjFourier = new FieldFourierCjSjCoefficients<>(nMax, sMax, jMax, context, field);
2796 this.cjCoefs = MathArrays.buildArray(field, 8, jMax + 1);
2797 this.sjCoefs = MathArrays.buildArray(field, 8, jMax + 1);
2798
2799 computeGeneratingFunctionCoefficients(context, hansen);
2800 }
2801
2802 /**
2803 * Compute the coefficients for the generating function S and its derivatives.
2804 * @param context container for attributes
2805 * @param hansenObjects hansen objects
2806 */
2807 private void computeGeneratingFunctionCoefficients(final FieldDSSTThirdBodyContext<T> context,
2808 final FieldHansenObjects<T> hansenObjects) {
2809
2810 final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
2811
2812 // Access to potential U derivatives
2813 final FieldUAnddU<T> udu = new FieldUAnddU<>(context, hansenObjects);
2814
2815 //Compute the C<sup>j</sup> coefficients
2816 for (int j = 1; j <= jMax; j++) {
2817 //Compute the C<sup>j</sup> coefficients
2818 cjCoefs[0][j] = cjsjFourier.getCj(j);
2819 cjCoefs[1][j] = cjsjFourier.getdCjda(j);
2820 cjCoefs[2][j] = cjsjFourier.getdCjdk(j).subtract((cjsjFourier.getSjLambda(j - 1).subtract(cjsjFourier.getSjLambda(j + 1))).divide(2.));
2821 cjCoefs[3][j] = cjsjFourier.getdCjdh(j).subtract((cjsjFourier.getCjLambda(j - 1).add(cjsjFourier.getCjLambda(j + 1))).divide(2.));
2822 cjCoefs[4][j] = cjsjFourier.getdCjdalpha(j);
2823 cjCoefs[5][j] = cjsjFourier.getdCjdbeta(j);
2824 cjCoefs[6][j] = cjsjFourier.getdCjdgamma(j);
2825 cjCoefs[7][j] = cjsjFourier.getCjLambda(j);
2826
2827 //Compute the S<sup>j</sup> coefficients
2828 sjCoefs[0][j] = cjsjFourier.getSj(j);
2829 sjCoefs[1][j] = cjsjFourier.getdSjda(j);
2830 sjCoefs[2][j] = cjsjFourier.getdSjdk(j).add((cjsjFourier.getCjLambda(j - 1).subtract(cjsjFourier.getCjLambda(j + 1))).divide(2.));
2831 sjCoefs[3][j] = cjsjFourier.getdSjdh(j).subtract((cjsjFourier.getSjLambda(j - 1).add(cjsjFourier.getSjLambda(j + 1))).divide(2.));
2832 sjCoefs[4][j] = cjsjFourier.getdSjdalpha(j);
2833 sjCoefs[5][j] = cjsjFourier.getdSjdbeta(j);
2834 sjCoefs[6][j] = cjsjFourier.getdSjdgamma(j);
2835 sjCoefs[7][j] = cjsjFourier.getSjLambda(j);
2836
2837 //In the special case j == 1 there are some additional terms to be added
2838 if (j == 1) {
2839 //Additional terms for C<sup>j</sup> coefficients
2840 cjCoefs[0][j] = cjCoefs[0][j].add(auxiliaryElements.getH().negate().multiply(udu.getU()));
2841 cjCoefs[1][j] = cjCoefs[1][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUda()));
2842 cjCoefs[2][j] = cjCoefs[2][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdk()));
2843 cjCoefs[3][j] = cjCoefs[3][j].add(auxiliaryElements.getH().multiply(udu.getdUdh()).add(udu.getU()).add(cjsjFourier.getC0Lambda()).negate());
2844 cjCoefs[4][j] = cjCoefs[4][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdAl()));
2845 cjCoefs[5][j] = cjCoefs[5][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdBe()));
2846 cjCoefs[6][j] = cjCoefs[6][j].add(auxiliaryElements.getH().negate().multiply(udu.getdUdGa()));
2847
2848 //Additional terms for S<sup>j</sup> coefficients
2849 sjCoefs[0][j] = sjCoefs[0][j].add(auxiliaryElements.getK().multiply(udu.getU()));
2850 sjCoefs[1][j] = sjCoefs[1][j].add(auxiliaryElements.getK().multiply(udu.getdUda()));
2851 sjCoefs[2][j] = sjCoefs[2][j].add(auxiliaryElements.getK().multiply(udu.getdUdk()).add(udu.getU()).add(cjsjFourier.getC0Lambda()));
2852 sjCoefs[3][j] = sjCoefs[3][j].add(auxiliaryElements.getK().multiply(udu.getdUdh()));
2853 sjCoefs[4][j] = sjCoefs[4][j].add(auxiliaryElements.getK().multiply(udu.getdUdAl()));
2854 sjCoefs[5][j] = sjCoefs[5][j].add(auxiliaryElements.getK().multiply(udu.getdUdBe()));
2855 sjCoefs[6][j] = sjCoefs[6][j].add(auxiliaryElements.getK().multiply(udu.getdUdGa()));
2856 }
2857 }
2858 }
2859
2860 /** Get the coefficient C<sup>j</sup> for the function S.
2861 * <br>
2862 * Possible values for j are within the interval [0,jMax].
2863 * The value 0 is used to obtain the free coefficient C⁰
2864 * @param j j index
2865 * @return C<sup>j</sup> for the function S
2866 */
2867 public T getSCj(final int j) {
2868 return cjCoefs[0][j];
2869 }
2870
2871 /** Get the coefficient S<sup>j</sup> for the function S.
2872 * <br>
2873 * Possible values for j are within the interval [1,jMax].
2874 * @param j j index
2875 * @return S<sup>j</sup> for the function S
2876 */
2877 public T getSSj(final int j) {
2878 return sjCoefs[0][j];
2879 }
2880
2881 /** Get the coefficient C<sup>j</sup> for the derivative dS/da.
2882 * <br>
2883 * Possible values for j are within the interval [0,jMax].
2884 * The value 0 is used to obtain the free coefficient C⁰
2885 * @param j j index
2886 * @return C<sup>j</sup> for the function dS/da
2887 */
2888 public T getdSdaCj(final int j) {
2889 return cjCoefs[1][j];
2890 }
2891
2892 /** Get the coefficient S<sup>j</sup> for the derivative dS/da.
2893 * <br>
2894 * Possible values for j are within the interval [1,jMax].
2895 * @param j j index
2896 * @return S<sup>j</sup> for the derivative dS/da
2897 */
2898 public T getdSdaSj(final int j) {
2899 return sjCoefs[1][j];
2900 }
2901
2902 /** Get the coefficient C<sup>j</sup> for the derivative dS/dk
2903 * <br>
2904 * Possible values for j are within the interval [0,jMax].
2905 * The value 0 is used to obtain the free coefficient C⁰
2906 * @param j j index
2907 * @return C<sup>j</sup> for the function dS/dk
2908 */
2909 public T getdSdkCj(final int j) {
2910 return cjCoefs[2][j];
2911 }
2912
2913 /** Get the coefficient S<sup>j</sup> for the derivative dS/dk.
2914 * <br>
2915 * Possible values for j are within the interval [1,jMax].
2916 * @param j j index
2917 * @return S<sup>j</sup> for the derivative dS/dk
2918 */
2919 public T getdSdkSj(final int j) {
2920 return sjCoefs[2][j];
2921 }
2922
2923 /** Get the coefficient C<sup>j</sup> for the derivative dS/dh
2924 * <br>
2925 * Possible values for j are within the interval [0,jMax].
2926 * The value 0 is used to obtain the free coefficient C⁰
2927 * @param j j index
2928 * @return C<sup>j</sup> for the function dS/dh
2929 */
2930 public T getdSdhCj(final int j) {
2931 return cjCoefs[3][j];
2932 }
2933
2934 /** Get the coefficient S<sup>j</sup> for the derivative dS/dh.
2935 * <br>
2936 * Possible values for j are within the interval [1,jMax].
2937 * @param j j index
2938 * @return S<sup>j</sup> for the derivative dS/dh
2939 */
2940 public T getdSdhSj(final int j) {
2941 return sjCoefs[3][j];
2942 }
2943
2944 /** Get the coefficient C<sup>j</sup> for the derivative dS/dα
2945 * <br>
2946 * Possible values for j are within the interval [0,jMax].
2947 * The value 0 is used to obtain the free coefficient C⁰
2948 * @param j j index
2949 * @return C<sup>j</sup> for the function dS/dα
2950 */
2951 public T getdSdalphaCj(final int j) {
2952 return cjCoefs[4][j];
2953 }
2954
2955 /** Get the coefficient S<sup>j</sup> for the derivative dS/dα.
2956 * <br>
2957 * Possible values for j are within the interval [1,jMax].
2958 * @param j j index
2959 * @return S<sup>j</sup> for the derivative dS/dα
2960 */
2961 public T getdSdalphaSj(final int j) {
2962 return sjCoefs[4][j];
2963 }
2964
2965 /** Get the coefficient C<sup>j</sup> for the derivative dS/dβ
2966 * <br>
2967 * Possible values for j are within the interval [0,jMax].
2968 * The value 0 is used to obtain the free coefficient C⁰
2969 * @param j j index
2970 * @return C<sup>j</sup> for the function dS/dβ
2971 */
2972 public T getdSdbetaCj(final int j) {
2973 return cjCoefs[5][j];
2974 }
2975
2976 /** Get the coefficient S<sup>j</sup> for the derivative dS/dβ.
2977 * <br>
2978 * Possible values for j are within the interval [1,jMax].
2979 * @param j j index
2980 * @return S<sup>j</sup> for the derivative dS/dβ
2981 */
2982 public T getdSdbetaSj(final int j) {
2983 return sjCoefs[5][j];
2984 }
2985
2986 /** Get the coefficient C<sup>j</sup> for the derivative dS/dγ
2987 * <br>
2988 * Possible values for j are within the interval [0,jMax].
2989 * The value 0 is used to obtain the free coefficient C⁰
2990 * @param j j index
2991 * @return C<sup>j</sup> for the function dS/dγ
2992 */
2993 public T getdSdgammaCj(final int j) {
2994 return cjCoefs[6][j];
2995 }
2996
2997 /** Get the coefficient S<sup>j</sup> for the derivative dS/dγ.
2998 * <br>
2999 * Possible values for j are within the interval [1,jMax].
3000 * @param j j index
3001 * @return S<sup>j</sup> for the derivative dS/dγ
3002 */
3003 public T getdSdgammaSj(final int j) {
3004 return sjCoefs[6][j];
3005 }
3006
3007 /** Get the coefficient C<sup>j</sup> for the derivative dS/dλ
3008 * <br>
3009 * Possible values for j are within the interval [0,jMax].
3010 * The value 0 is used to obtain the free coefficient C⁰
3011 * @param j j index
3012 * @return C<sup>j</sup> for the function dS/dλ
3013 */
3014 public T getdSdlambdaCj(final int j) {
3015 return cjCoefs[7][j];
3016 }
3017
3018 /** Get the coefficient S<sup>j</sup> for the derivative dS/dλ.
3019 * <br>
3020 * Possible values for j are within the interval [1,jMax].
3021 * @param j j index
3022 * @return S<sup>j</sup> for the derivative dS/dλ
3023 */
3024 public T getdSdlambdaSj(final int j) {
3025 return sjCoefs[7][j];
3026 }
3027 }
3028
3029 /**
3030 * The coefficients used to compute the short periodic contribution for the Third body perturbation.
3031 * <p>
3032 * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
3033 * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
3034 * are computed by replacing the corresponding values in formula 2.5.5-10.
3035 * </p>
3036 * @author Lucian Barbulescu
3037 */
3038 private static class ThirdBodyShortPeriodicCoefficients implements ShortPeriodTerms {
3039
3040 /** Maximal value for j. */
3041 private final int jMax;
3042
3043 /** Number of points used in the interpolation process. */
3044 private final int interpolationPoints;
3045
3046 /** Max frequency of F. */
3047 private final int maxFreqF;
3048
3049 /** Coefficients prefix. */
3050 private final String prefix;
3051
3052 /** All coefficients slots. */
3053 private final transient TimeSpanMap<Slot> slots;
3054
3055 /**
3056 * Standard constructor.
3057 * @param interpolationPoints number of points used in the interpolation process
3058 * @param jMax maximal value for j
3059 * @param maxFreqF Max frequency of F
3060 * @param bodyName third body name
3061 * @param slots all coefficients slots
3062 */
3063 ThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
3064 final int maxFreqF, final String bodyName,
3065 final TimeSpanMap<Slot> slots) {
3066 this.jMax = jMax;
3067 this.interpolationPoints = interpolationPoints;
3068 this.maxFreqF = maxFreqF;
3069 this.prefix = DSSTThirdBody.SHORT_PERIOD_PREFIX + bodyName + "-";
3070 this.slots = slots;
3071 }
3072
3073 /** Get the slot valid for some date.
3074 * @param meanStates mean states defining the slot
3075 * @return slot valid at the specified date
3076 */
3077 public Slot createSlot(final SpacecraftState... meanStates) {
3078 final Slot slot = new Slot(jMax, interpolationPoints);
3079 final AbsoluteDate first = meanStates[0].getDate();
3080 final AbsoluteDate last = meanStates[meanStates.length - 1].getDate();
3081 final int compare = first.compareTo(last);
3082 if (compare < 0) {
3083 slots.addValidAfter(slot, first, false);
3084 } else if (compare > 0) {
3085 slots.addValidBefore(slot, first, false);
3086 } else {
3087 // single date, valid for all time
3088 slots.addValidAfter(slot, AbsoluteDate.PAST_INFINITY, false);
3089 }
3090 return slot;
3091 }
3092
3093 /** {@inheritDoc} */
3094 @Override
3095 public double[] value(final Orbit meanOrbit) {
3096
3097 // select the coefficients slot
3098 final Slot slot = slots.get(meanOrbit.getDate());
3099
3100 // the current eccentric longitude
3101 final double F = meanOrbit.getLE();
3102
3103 //initialize the short periodic contribution with the corresponding C⁰ coeficient
3104 final double[] shortPeriodic = slot.cij[0].value(meanOrbit.getDate());
3105
3106 // Add the cos and sin dependent terms
3107 for (int j = 1; j <= maxFreqF; j++) {
3108 //compute cos and sin
3109 final SinCos scjF = FastMath.sinCos(j * F);
3110
3111 final double[] c = slot.cij[j].value(meanOrbit.getDate());
3112 final double[] s = slot.sij[j].value(meanOrbit.getDate());
3113 for (int i = 0; i < 6; i++) {
3114 shortPeriodic[i] += c[i] * scjF.cos() + s[i] * scjF.sin();
3115 }
3116 }
3117
3118 return shortPeriodic;
3119
3120 }
3121
3122 /** {@inheritDoc} */
3123 @Override
3124 public String getCoefficientsKeyPrefix() {
3125 return prefix;
3126 }
3127
3128 /** {@inheritDoc}
3129 * <p>
3130 * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
3131 * maxFreqF sj coefficients where maxFreqF depends on the orbit.
3132 * The j index is the integer multiplier for the eccentric longitude argument
3133 * in the cj and sj coefficients.
3134 * </p>
3135 */
3136 @Override
3137 public Map<String, double[]> getCoefficients(final AbsoluteDate date, final Set<String> selected) {
3138
3139 // select the coefficients slot
3140 final Slot slot = slots.get(date);
3141
3142 final Map<String, double[]> coefficients = new HashMap<String, double[]>(2 * maxFreqF + 1);
3143 storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
3144 for (int j = 1; j <= maxFreqF; j++) {
3145 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
3146 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
3147 }
3148 return coefficients;
3149
3150 }
3151
3152 /** Put a coefficient in a map if selected.
3153 * @param map map to populate
3154 * @param selected set of coefficients that should be put in the map
3155 * (empty set means all coefficients are selected)
3156 * @param value coefficient value
3157 * @param id coefficient identifier
3158 * @param indices list of coefficient indices
3159 */
3160 private void storeIfSelected(final Map<String, double[]> map, final Set<String> selected,
3161 final double[] value, final String id, final int... indices) {
3162 final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
3163 keyBuilder.append(id);
3164 for (int index : indices) {
3165 keyBuilder.append('[').append(index).append(']');
3166 }
3167 final String key = keyBuilder.toString();
3168 if (selected.isEmpty() || selected.contains(key)) {
3169 map.put(key, value);
3170 }
3171 }
3172
3173 }
3174
3175 /**
3176 * The coefficients used to compute the short periodic contribution for the Third body perturbation.
3177 * <p>
3178 * The short periodic contribution for the Third Body is expressed in Danielson 4.2-25.<br>
3179 * The coefficients C<sub>i</sub>⁰, C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup>
3180 * are computed by replacing the corresponding values in formula 2.5.5-10.
3181 * </p>
3182 * @author Lucian Barbulescu
3183 */
3184 private static class FieldThirdBodyShortPeriodicCoefficients <T extends CalculusFieldElement<T>> implements FieldShortPeriodTerms<T> {
3185
3186 /** Maximal value for j. */
3187 private final int jMax;
3188
3189 /** Number of points used in the interpolation process. */
3190 private final int interpolationPoints;
3191
3192 /** Max frequency of F. */
3193 private final int maxFreqF;
3194
3195 /** Coefficients prefix. */
3196 private final String prefix;
3197
3198 /** All coefficients slots. */
3199 private final transient FieldTimeSpanMap<FieldSlot<T>, T> slots;
3200
3201 /**
3202 * Standard constructor.
3203 * @param interpolationPoints number of points used in the interpolation process
3204 * @param jMax maximal value for j
3205 * @param maxFreqF Max frequency of F
3206 * @param bodyName third body name
3207 * @param slots all coefficients slots
3208 */
3209 FieldThirdBodyShortPeriodicCoefficients(final int jMax, final int interpolationPoints,
3210 final int maxFreqF, final String bodyName,
3211 final FieldTimeSpanMap<FieldSlot<T>, T> slots) {
3212 this.jMax = jMax;
3213 this.interpolationPoints = interpolationPoints;
3214 this.maxFreqF = maxFreqF;
3215 this.prefix = DSSTThirdBody.SHORT_PERIOD_PREFIX + bodyName + "-";
3216 this.slots = slots;
3217 }
3218
3219 /** Get the slot valid for some date.
3220 * @param meanStates mean states defining the slot
3221 * @return slot valid at the specified date
3222 */
3223 @SuppressWarnings("unchecked")
3224 public FieldSlot<T> createSlot(final FieldSpacecraftState<T>... meanStates) {
3225 final FieldSlot<T> slot = new FieldSlot<>(jMax, interpolationPoints);
3226 final FieldAbsoluteDate<T> first = meanStates[0].getDate();
3227 final FieldAbsoluteDate<T> last = meanStates[meanStates.length - 1].getDate();
3228 if (first.compareTo(last) <= 0) {
3229 slots.addValidAfter(slot, first);
3230 } else {
3231 slots.addValidBefore(slot, first);
3232 }
3233 return slot;
3234 }
3235
3236 /** {@inheritDoc} */
3237 @Override
3238 public T[] value(final FieldOrbit<T> meanOrbit) {
3239
3240 // select the coefficients slot
3241 final FieldSlot<T> slot = slots.get(meanOrbit.getDate());
3242
3243 // the current eccentric longitude
3244 final T F = meanOrbit.getLE();
3245
3246 //initialize the short periodic contribution with the corresponding C⁰ coeficient
3247 final T[] shortPeriodic = (T[]) slot.cij[0].value(meanOrbit.getDate());
3248
3249 // Add the cos and sin dependent terms
3250 for (int j = 1; j <= maxFreqF; j++) {
3251 //compute cos and sin
3252 final FieldSinCos<T> scjF = FastMath.sinCos(F.multiply(j));
3253
3254 final T[] c = (T[]) slot.cij[j].value(meanOrbit.getDate());
3255 final T[] s = (T[]) slot.sij[j].value(meanOrbit.getDate());
3256 for (int i = 0; i < 6; i++) {
3257 shortPeriodic[i] = shortPeriodic[i].add(c[i].multiply(scjF.cos()).add(s[i].multiply(scjF.sin())));
3258 }
3259 }
3260
3261 return shortPeriodic;
3262
3263 }
3264
3265 /** {@inheritDoc} */
3266 @Override
3267 public String getCoefficientsKeyPrefix() {
3268 return prefix;
3269 }
3270
3271 /** {@inheritDoc}
3272 * <p>
3273 * For third body attraction forces,there are maxFreqF + 1 cj coefficients,
3274 * maxFreqF sj coefficients where maxFreqF depends on the orbit.
3275 * The j index is the integer multiplier for the eccentric longitude argument
3276 * in the cj and sj coefficients.
3277 * </p>
3278 */
3279 @Override
3280 public Map<String, T[]> getCoefficients(final FieldAbsoluteDate<T> date, final Set<String> selected) {
3281
3282 // select the coefficients slot
3283 final FieldSlot<T> slot = slots.get(date);
3284
3285 final Map<String, T[]> coefficients = new HashMap<String, T[]>(2 * maxFreqF + 1);
3286 storeIfSelected(coefficients, selected, slot.cij[0].value(date), "c", 0);
3287 for (int j = 1; j <= maxFreqF; j++) {
3288 storeIfSelected(coefficients, selected, slot.cij[j].value(date), "c", j);
3289 storeIfSelected(coefficients, selected, slot.sij[j].value(date), "s", j);
3290 }
3291 return coefficients;
3292
3293 }
3294
3295 /** Put a coefficient in a map if selected.
3296 * @param map map to populate
3297 * @param selected set of coefficients that should be put in the map
3298 * (empty set means all coefficients are selected)
3299 * @param value coefficient value
3300 * @param id coefficient identifier
3301 * @param indices list of coefficient indices
3302 */
3303 private void storeIfSelected(final Map<String, T[]> map, final Set<String> selected,
3304 final T[] value, final String id, final int... indices) {
3305 final StringBuilder keyBuilder = new StringBuilder(getCoefficientsKeyPrefix());
3306 keyBuilder.append(id);
3307 for (int index : indices) {
3308 keyBuilder.append('[').append(index).append(']');
3309 }
3310 final String key = keyBuilder.toString();
3311 if (selected.isEmpty() || selected.contains(key)) {
3312 map.put(key, value);
3313 }
3314 }
3315
3316 }
3317
3318 /** Coefficients valid for one time slot. */
3319 private static class Slot {
3320
3321 /** The coefficients C<sub>i</sub><sup>j</sup>.
3322 * <p>
3323 * The index order is cij[j][i] <br/>
3324 * i corresponds to the equinoctial element, as follows: <br/>
3325 * - i=0 for a <br/>
3326 * - i=1 for k <br/>
3327 * - i=2 for h <br/>
3328 * - i=3 for q <br/>
3329 * - i=4 for p <br/>
3330 * - i=5 for λ <br/>
3331 * </p>
3332 */
3333 private final ShortPeriodicsInterpolatedCoefficient[] cij;
3334
3335 /** The coefficients S<sub>i</sub><sup>j</sup>.
3336 * <p>
3337 * The index order is sij[j][i] <br/>
3338 * i corresponds to the equinoctial element, as follows: <br/>
3339 * - i=0 for a <br/>
3340 * - i=1 for k <br/>
3341 * - i=2 for h <br/>
3342 * - i=3 for q <br/>
3343 * - i=4 for p <br/>
3344 * - i=5 for λ <br/>
3345 * </p>
3346 */
3347 private final ShortPeriodicsInterpolatedCoefficient[] sij;
3348
3349 /** Simple constructor.
3350 * @param jMax maximum value for j index
3351 * @param interpolationPoints number of points used in the interpolation process
3352 */
3353 Slot(final int jMax, final int interpolationPoints) {
3354 // allocate the coefficients arrays
3355 cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
3356 sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1];
3357 for (int j = 0; j <= jMax; j++) {
3358 cij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
3359 sij[j] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
3360 }
3361
3362
3363 }
3364 }
3365
3366 /** Coefficients valid for one time slot. */
3367 private static class FieldSlot <T extends CalculusFieldElement<T>> {
3368
3369 /** The coefficients C<sub>i</sub><sup>j</sup>.
3370 * <p>
3371 * The index order is cij[j][i] <br/>
3372 * i corresponds to the equinoctial element, as follows: <br/>
3373 * - i=0 for a <br/>
3374 * - i=1 for k <br/>
3375 * - i=2 for h <br/>
3376 * - i=3 for q <br/>
3377 * - i=4 for p <br/>
3378 * - i=5 for λ <br/>
3379 * </p>
3380 */
3381 private final FieldShortPeriodicsInterpolatedCoefficient<T>[] cij;
3382
3383 /** The coefficients S<sub>i</sub><sup>j</sup>.
3384 * <p>
3385 * The index order is sij[j][i] <br/>
3386 * i corresponds to the equinoctial element, as follows: <br/>
3387 * - i=0 for a <br/>
3388 * - i=1 for k <br/>
3389 * - i=2 for h <br/>
3390 * - i=3 for q <br/>
3391 * - i=4 for p <br/>
3392 * - i=5 for λ <br/>
3393 * </p>
3394 */
3395 private final FieldShortPeriodicsInterpolatedCoefficient<T>[] sij;
3396
3397 /** Simple constructor.
3398 * @param jMax maximum value for j index
3399 * @param interpolationPoints number of points used in the interpolation process
3400 */
3401 @SuppressWarnings("unchecked")
3402 FieldSlot(final int jMax, final int interpolationPoints) {
3403 // allocate the coefficients arrays
3404 cij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
3405 sij = (FieldShortPeriodicsInterpolatedCoefficient<T>[]) Array.newInstance(FieldShortPeriodicsInterpolatedCoefficient.class, jMax + 1);
3406 for (int j = 0; j <= jMax; j++) {
3407 cij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
3408 sij[j] = new FieldShortPeriodicsInterpolatedCoefficient<>(interpolationPoints);
3409 }
3410
3411
3412 }
3413 }
3414
3415 /** Compute potential and potential derivatives with respect to orbital parameters. */
3416 private class UAnddU {
3417
3418 /** The current value of the U function. <br/>
3419 * Needed for the short periodic contribution */
3420 private double U;
3421
3422 /** dU / da. */
3423 private double dUda;
3424
3425 /** dU / dk. */
3426 private double dUdk;
3427
3428 /** dU / dh. */
3429 private double dUdh;
3430
3431 /** dU / dAlpha. */
3432 private double dUdAl;
3433
3434 /** dU / dBeta. */
3435 private double dUdBe;
3436
3437 /** dU / dGamma. */
3438 private double dUdGa;
3439
3440 /** Simple constuctor.
3441 * @param context container for attributes
3442 * @param hansen hansen objects
3443 */
3444 UAnddU(final DSSTThirdBodyContext context,
3445 final HansenObjects hansen) {
3446 // Auxiliary elements related to the current orbit
3447 final AuxiliaryElements auxiliaryElements = context.getAuxiliaryElements();
3448
3449 // Gs and Hs coefficients
3450 final double[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta(), context.getMaxEccPow());
3451
3452 // Initialise U.
3453 U = 0.;
3454
3455 // Potential derivatives
3456 dUda = 0.;
3457 dUdk = 0.;
3458 dUdh = 0.;
3459 dUdAl = 0.;
3460 dUdBe = 0.;
3461 dUdGa = 0.;
3462
3463 for (int s = 0; s <= context.getMaxEccPow(); s++) {
3464
3465 // initialise the Hansen roots
3466 hansen.computeHansenObjectsInitValues(context, auxiliaryElements.getB(), s);
3467
3468 // Get the current Gs coefficient
3469 final double gs = GsHs[0][s];
3470
3471 // Compute Gs partial derivatives from 3.1-(9)
3472 double dGsdh = 0.;
3473 double dGsdk = 0.;
3474 double dGsdAl = 0.;
3475 double dGsdBe = 0.;
3476 if (s > 0) {
3477 // First get the G(s-1) and the H(s-1) coefficients
3478 final double sxGsm1 = s * GsHs[0][s - 1];
3479 final double sxHsm1 = s * GsHs[1][s - 1];
3480 // Then compute derivatives
3481 dGsdh = context.getBeta() * sxGsm1 - context.getAlpha() * sxHsm1;
3482 dGsdk = context.getAlpha() * sxGsm1 + context.getBeta() * sxHsm1;
3483 dGsdAl = auxiliaryElements.getK() * sxGsm1 - auxiliaryElements.getH() * sxHsm1;
3484 dGsdBe = auxiliaryElements.getH() * sxGsm1 + auxiliaryElements.getK() * sxHsm1;
3485 }
3486
3487 // Kronecker symbol (2 - delta(0,s))
3488 final double delta0s = (s == 0) ? 1. : 2.;
3489
3490 for (int n = FastMath.max(2, s); n <= context.getMaxAR3Pow(); n++) {
3491 // (n - s) must be even
3492 if ((n - s) % 2 == 0) {
3493 // Extract data from previous computation :
3494 final double kns = hansen.getHansenObjects()[s].getValue(n, auxiliaryElements.getB());
3495 final double dkns = hansen.getHansenObjects()[s].getDerivative(n, auxiliaryElements.getB());
3496
3497 final double vns = Vns.get(new NSKey(n, s));
3498 final double coef0 = delta0s * context.getAoR3Pow()[n] * vns;
3499 final double coef1 = coef0 * context.getQns()[n][s];
3500 final double coef2 = coef1 * kns;
3501 // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
3502 // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
3503 final double dqns = (n == s) ? 0. : context.getQns()[n][s + 1];
3504
3505 //Compute U:
3506 U += coef2 * gs;
3507
3508 // Compute dU / da :
3509 dUda += coef2 * n * gs;
3510 // Compute dU / dh
3511 dUdh += coef1 * (kns * dGsdh + context.getHXXX() * gs * dkns);
3512 // Compute dU / dk
3513 dUdk += coef1 * (kns * dGsdk + context.getKXXX() * gs * dkns);
3514 // Compute dU / dAlpha
3515 dUdAl += coef2 * dGsdAl;
3516 // Compute dU / dBeta
3517 dUdBe += coef2 * dGsdBe;
3518 // Compute dU / dGamma
3519 dUdGa += coef0 * kns * dqns * gs;
3520 }
3521 }
3522 }
3523
3524 // multiply by mu3 / R3
3525 this.U = U * context.getMuoR3();
3526
3527 this.dUda = dUda * context.getMuoR3() / auxiliaryElements.getSma();
3528 this.dUdk = dUdk * context.getMuoR3();
3529 this.dUdh = dUdh * context.getMuoR3();
3530 this.dUdAl = dUdAl * context.getMuoR3();
3531 this.dUdBe = dUdBe * context.getMuoR3();
3532 this.dUdGa = dUdGa * context.getMuoR3();
3533
3534 }
3535
3536 /** Return value of U.
3537 * @return U
3538 */
3539 public double getU() {
3540 return U;
3541 }
3542
3543 /** Return value of dU / da.
3544 * @return dUda
3545 */
3546 public double getdUda() {
3547 return dUda;
3548 }
3549
3550 /** Return value of dU / dk.
3551 * @return dUdk
3552 */
3553 public double getdUdk() {
3554 return dUdk;
3555 }
3556
3557 /** Return value of dU / dh.
3558 * @return dUdh
3559 */
3560 public double getdUdh() {
3561 return dUdh;
3562 }
3563
3564 /** Return value of dU / dAlpha.
3565 * @return dUdAl
3566 */
3567 public double getdUdAl() {
3568 return dUdAl;
3569 }
3570
3571 /** Return value of dU / dBeta.
3572 * @return dUdBe
3573 */
3574 public double getdUdBe() {
3575 return dUdBe;
3576 }
3577
3578 /** Return value of dU / dGamma.
3579 * @return dUdGa
3580 */
3581 public double getdUdGa() {
3582 return dUdGa;
3583 }
3584
3585 }
3586
3587 /** Compute potential and potential derivatives with respect to orbital parameters. */
3588 private class FieldUAnddU <T extends CalculusFieldElement<T>> {
3589
3590 /** The current value of the U function. <br/>
3591 * Needed for the short periodic contribution */
3592 private T U;
3593
3594 /** dU / da. */
3595 private T dUda;
3596
3597 /** dU / dk. */
3598 private T dUdk;
3599
3600 /** dU / dh. */
3601 private T dUdh;
3602
3603 /** dU / dAlpha. */
3604 private T dUdAl;
3605
3606 /** dU / dBeta. */
3607 private T dUdBe;
3608
3609 /** dU / dGamma. */
3610 private T dUdGa;
3611
3612 /** Simple constuctor.
3613 * @param context container for attributes
3614 * @param hansen hansen objects
3615 */
3616 FieldUAnddU(final FieldDSSTThirdBodyContext<T> context,
3617 final FieldHansenObjects<T> hansen) {
3618
3619 // Auxiliary elements related to the current orbit
3620 final FieldAuxiliaryElements<T> auxiliaryElements = context.getFieldAuxiliaryElements();
3621
3622 // Field for array building
3623 final Field<T> field = auxiliaryElements.getDate().getField();
3624 // Zero for initialization
3625 final T zero = field.getZero();
3626
3627 // Gs and Hs coefficients
3628 final T[][] GsHs = CoefficientsFactory.computeGsHs(auxiliaryElements.getK(), auxiliaryElements.getH(), context.getAlpha(), context.getBeta(), context.getMaxEccPow(), field);
3629
3630 // Initialise U.
3631 U = zero;
3632
3633 // Potential derivatives
3634 dUda = zero;
3635 dUdk = zero;
3636 dUdh = zero;
3637 dUdAl = zero;
3638 dUdBe = zero;
3639 dUdGa = zero;
3640
3641 for (int s = 0; s <= context.getMaxEccPow(); s++) {
3642 // initialise the Hansen roots
3643 hansen.computeHansenObjectsInitValues(context, auxiliaryElements.getB(), s);
3644
3645 // Get the current Gs coefficient
3646 final T gs = GsHs[0][s];
3647
3648 // Compute Gs partial derivatives from 3.1-(9)
3649 T dGsdh = zero;
3650 T dGsdk = zero;
3651 T dGsdAl = zero;
3652 T dGsdBe = zero;
3653 if (s > 0) {
3654 // First get the G(s-1) and the H(s-1) coefficients
3655 final T sxGsm1 = GsHs[0][s - 1].multiply(s);
3656 final T sxHsm1 = GsHs[1][s - 1].multiply(s);
3657 // Then compute derivatives
3658 dGsdh = sxGsm1.multiply(context.getBeta()).subtract(sxHsm1.multiply(context.getAlpha()));
3659 dGsdk = sxGsm1.multiply(context.getAlpha()).add(sxHsm1.multiply(context.getBeta()));
3660 dGsdAl = sxGsm1.multiply(auxiliaryElements.getK()).subtract(sxHsm1.multiply(auxiliaryElements.getH()));
3661 dGsdBe = sxGsm1.multiply(auxiliaryElements.getH()).add(sxHsm1.multiply(auxiliaryElements.getK()));
3662 }
3663
3664 // Kronecker symbol (2 - delta(0,s))
3665 final T delta0s = zero.add((s == 0) ? 1. : 2.);
3666
3667 for (int n = FastMath.max(2, s); n <= context.getMaxAR3Pow(); n++) {
3668 // (n - s) must be even
3669 if ((n - s) % 2 == 0) {
3670 // Extract data from previous computation :
3671 final T kns = (T) hansen.getHansenObjects()[s].getValue(n, auxiliaryElements.getB());
3672 final T dkns = (T) hansen.getHansenObjects()[s].getDerivative(n, auxiliaryElements.getB());
3673
3674 final double vns = Vns.get(new NSKey(n, s));
3675 final T coef0 = delta0s.multiply(vns).multiply(context.getAoR3Pow()[n]);
3676 final T coef1 = coef0.multiply(context.getQns()[n][s]);
3677 final T coef2 = coef1.multiply(kns);
3678 // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
3679 // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
3680 final T dqns = (n == s) ? zero : context.getQns()[n][s + 1];
3681
3682 //Compute U:
3683 U = U.add(coef2.multiply(gs));
3684
3685 // Compute dU / da :
3686 dUda = dUda.add(coef2.multiply(n).multiply(gs));
3687 // Compute dU / dh
3688 dUdh = dUdh.add(coef1.multiply(dGsdh.multiply(kns).add(context.getHXXX().multiply(gs).multiply(dkns))));
3689 // Compute dU / dk
3690 dUdk = dUdk.add(coef1.multiply(dGsdk.multiply(kns).add(context.getKXXX().multiply(gs).multiply(dkns))));
3691 // Compute dU / dAlpha
3692 dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
3693 // Compute dU / dBeta
3694 dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
3695 // Compute dU / dGamma
3696 dUdGa = dUdGa.add(coef0.multiply(kns).multiply(dqns).multiply(gs));
3697 }
3698 }
3699 }
3700
3701 // multiply by mu3 / R3
3702 this.U = U.multiply(context.getMuoR3());
3703
3704 this.dUda = dUda.multiply(context.getMuoR3().divide(auxiliaryElements.getSma()));
3705 this.dUdk = dUdk.multiply(context.getMuoR3());
3706 this.dUdh = dUdh.multiply(context.getMuoR3());
3707 this.dUdAl = dUdAl.multiply(context.getMuoR3());
3708 this.dUdBe = dUdBe.multiply(context.getMuoR3());
3709 this.dUdGa = dUdGa.multiply(context.getMuoR3());
3710
3711 }
3712
3713 /** Return value of U.
3714 * @return U
3715 */
3716 public T getU() {
3717 return U;
3718 }
3719
3720 /** Return value of dU / da.
3721 * @return dUda
3722 */
3723 public T getdUda() {
3724 return dUda;
3725 }
3726
3727 /** Return value of dU / dk.
3728 * @return dUdk
3729 */
3730 public T getdUdk() {
3731 return dUdk;
3732 }
3733
3734 /** Return value of dU / dh.
3735 * @return dUdh
3736 */
3737 public T getdUdh() {
3738 return dUdh;
3739 }
3740
3741 /** Return value of dU / dAlpha.
3742 * @return dUdAl
3743 */
3744 public T getdUdAl() {
3745 return dUdAl;
3746 }
3747
3748 /** Return value of dU / dBeta.
3749 * @return dUdBe
3750 */
3751 public T getdUdBe() {
3752 return dUdBe;
3753 }
3754
3755 /** Return value of dU / dGamma.
3756 * @return dUdGa
3757 */
3758 public T getdUdGa() {
3759 return dUdGa;
3760 }
3761
3762 }
3763
3764 /** Computes init values of the Hansen Objects. */
3765 private static class HansenObjects {
3766
3767 /** Max power for summation. */
3768 private static final int MAX_POWER = 22;
3769
3770 /** An array that contains the objects needed to build the Hansen coefficients. <br/>
3771 * The index is s */
3772 private final HansenThirdBodyLinear[] hansenObjects;
3773
3774 /** Simple constructor. */
3775 HansenObjects() {
3776 this.hansenObjects = new HansenThirdBodyLinear[MAX_POWER + 1];
3777 for (int s = 0; s <= MAX_POWER; s++) {
3778 this.hansenObjects[s] = new HansenThirdBodyLinear(MAX_POWER, s);
3779 }
3780 }
3781
3782 /** Compute init values for hansen objects.
3783 * @param context container for attributes
3784 * @param B = sqrt(1 - e²).
3785 * @param element element of the array to compute the init values
3786 */
3787 public void computeHansenObjectsInitValues(final DSSTThirdBodyContext context, final double B, final int element) {
3788 hansenObjects[element].computeInitValues(B, context.getBB(), context.getBBB());
3789 }
3790
3791 /** Get the Hansen Objects.
3792 * @return hansenObjects
3793 */
3794 public HansenThirdBodyLinear[] getHansenObjects() {
3795 return hansenObjects;
3796 }
3797
3798 }
3799
3800 /** Computes init values of the Hansen Objects. */
3801 private static class FieldHansenObjects<T extends CalculusFieldElement<T>> {
3802
3803 /** Max power for summation. */
3804 private static final int MAX_POWER = 22;
3805
3806 /** An array that contains the objects needed to build the Hansen coefficients. <br/>
3807 * The index is s */
3808 private final FieldHansenThirdBodyLinear<T>[] hansenObjects;
3809
3810 /** Simple constructor.
3811 * @param field field used by default
3812 */
3813 @SuppressWarnings("unchecked")
3814 FieldHansenObjects(final Field<T> field) {
3815 this.hansenObjects = (FieldHansenThirdBodyLinear<T>[]) Array.newInstance(FieldHansenThirdBodyLinear.class, MAX_POWER + 1);
3816 for (int s = 0; s <= MAX_POWER; s++) {
3817 this.hansenObjects[s] = new FieldHansenThirdBodyLinear<>(MAX_POWER, s, field);
3818 }
3819 }
3820
3821 /** Initialise the Hansen roots for third body problem.
3822 * @param context container for attributes
3823 * @param B = sqrt(1 - e²).
3824 * @param element element of the array to compute the init values
3825 */
3826 public void computeHansenObjectsInitValues(final FieldDSSTThirdBodyContext<T> context,
3827 final T B, final int element) {
3828 hansenObjects[element].computeInitValues(B, context.getBB(), context.getBBB());
3829 }
3830
3831 /** Get the Hansen Objects.
3832 * @return hansenObjects
3833 */
3834 public FieldHansenThirdBodyLinear<T>[] getHansenObjects() {
3835 return hansenObjects;
3836 }
3837
3838 }
3839
3840 }