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