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