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