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