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