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