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