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