1 /* Copyright 2002-2015 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.semianalytical.dsst.forces;
18
19 import org.apache.commons.math3.analysis.UnivariateVectorFunction;
20 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
21 import org.apache.commons.math3.util.FastMath;
22 import org.orekit.attitudes.Attitude;
23 import org.orekit.attitudes.AttitudeProvider;
24 import org.orekit.errors.OrekitException;
25 import org.orekit.errors.OrekitExceptionWrapper;
26 import org.orekit.forces.ForceModel;
27 import org.orekit.frames.Frame;
28 import org.orekit.orbits.EquinoctialOrbit;
29 import org.orekit.orbits.Orbit;
30 import org.orekit.orbits.OrbitType;
31 import org.orekit.orbits.PositionAngle;
32 import org.orekit.propagation.SpacecraftState;
33 import org.orekit.propagation.numerical.TimeDerivativesEquations;
34 import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
35 import org.orekit.propagation.semianalytical.dsst.utilities.CjSjCoefficient;
36 import org.orekit.propagation.semianalytical.dsst.utilities.ShortPeriodicsInterpolatedCoefficient;
37 import org.orekit.time.AbsoluteDate;
38
39 /** Common handling of {@link DSSTForceModel} methods for Gaussian contributions to DSST propagation.
40 * <p>
41 * This abstract class allows to provide easily a subset of {@link DSSTForceModel} methods
42 * for specific Gaussian contributions.
43 * </p><p>
44 * This class implements the notion of numerical averaging of the DSST theory.
45 * Numerical averaging is mainly used for non-conservative disturbing forces such as
46 * atmospheric drag and solar radiation pressure.
47 * </p><p>
48 * Gaussian contributions can be expressed as: da<sub>i</sub>/dt = δa<sub>i</sub>/δv . q<br>
49 * where:
50 * <ul>
51 * <li>a<sub>i</sub> are the six equinoctial elements</li>
52 * <li>v is the velocity vector</li>
53 * <li>q is the perturbing acceleration due to the considered force</li>
54 * </ul>
55 * The averaging process and other considerations lead to integrate this contribution
56 * over the true longitude L possibly taking into account some limits.
57 * </p><p>
58 * To create a numerically averaged contribution, one needs only to provide a
59 * {@link ForceModel} and to implement in the derived class the method:
60 * {@link #getLLimits(SpacecraftState)}.
61 * </p>
62 * @author Pascal Parraud
63 */
64 public abstract class AbstractGaussianContribution implements DSSTForceModel {
65
66 /** Available orders for Gauss quadrature. */
67 private static final int[] GAUSS_ORDER = {12, 16, 20, 24, 32, 40, 48};
68
69 /** Max rank in Gauss quadrature orders array. */
70 private static final int MAX_ORDER_RANK = GAUSS_ORDER.length - 1;
71
72 /** Number of points for interpolation. */
73 private static final int INTERPOLATION_POINTS = 3;
74
75 /** Maximum value for j index. */
76 private static final int JMAX = 12;
77 // CHECKSTYLE: stop VisibilityModifierCheck
78
79 /** Retrograde factor. */
80 protected int I;
81
82 /** a. */
83 protected double a;
84 /** e<sub>x</sub>. */
85 protected double k;
86 /** e<sub>y</sub>. */
87 protected double h;
88 /** h<sub>x</sub>. */
89 protected double q;
90 /** h<sub>y</sub>. */
91 protected double p;
92
93 /** Eccentricity. */
94 protected double ecc;
95
96 /** Kepler mean motion: n = sqrt(μ / a³). */
97 protected double n;
98
99 /** Mean longitude. */
100 protected double lm;
101
102 /** Equinoctial frame f vector. */
103 protected Vector3D f;
104 /** Equinoctial frame g vector. */
105 protected Vector3D g;
106 /** Equinoctial frame w vector. */
107 protected Vector3D w;
108
109 /** A = sqrt(μ * a). */
110 protected double A;
111 /** B = sqrt(1 - h² - k²). */
112 protected double B;
113 /** C = 1 + p² + q². */
114 protected double C;
115
116 /** 2 / (n² * a) . */
117 protected double ton2a;
118 /** 1 / A .*/
119 protected double ooA;
120 /** 1 / (A * B) .*/
121 protected double ooAB;
122 /** C / (2 * A * B) .*/
123 protected double co2AB;
124 /** 1 / (1 + B) .*/
125 protected double ooBpo;
126 /** 1 / μ .*/
127 protected double ooMu;
128 /** μ .*/
129 protected double mu;
130
131 // CHECKSTYLE: resume VisibilityModifierCheck
132
133 /** Contribution to be numerically averaged. */
134 private final ForceModel contribution;
135
136 /** Gauss integrator. */
137 private final double threshold;
138
139 /** Gauss integrator. */
140 private GaussQuadrature integrator;
141
142 /** Flag for Gauss order computation. */
143 private boolean isDirty;
144
145 /** Attitude provider. */
146 private AttitudeProvider attitudeProvider;
147
148 /** The C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients used to compute
149 * the short-periodic gaussian contribution. */
150 private GaussianShortPeriodicCoefficients gaussianSPCoefs;
151
152 /** The frame used to describe the orbits. */
153 private Frame frame;
154
155 /** Build a new instance.
156 *
157 * @param threshold tolerance for the choice of the Gauss quadrature order
158 * @param contribution the {@link ForceModel} to be numerically averaged
159 */
160 protected AbstractGaussianContribution(final double threshold,
161 final ForceModel contribution) {
162 this.contribution = contribution;
163 this.threshold = threshold;
164 this.integrator = new GaussQuadrature(GAUSS_ORDER[MAX_ORDER_RANK]);
165 this.isDirty = true;
166 }
167
168 /** {@inheritDoc} */
169 @Override
170 public void initialize(final AuxiliaryElements aux, final boolean meanOnly)
171 throws OrekitException {
172
173 // save the frame
174 this.frame = aux.getFrame();
175
176 if (!meanOnly) {
177 // initialize the short periodic coefficient generator, if needed.
178 this.gaussianSPCoefs = new GaussianShortPeriodicCoefficients(JMAX, INTERPOLATION_POINTS);
179 }
180 }
181
182 /** {@inheritDoc} */
183 @Override
184 public void initializeStep(final AuxiliaryElements aux)
185 throws OrekitException {
186
187 // Equinoctial elements
188 a = aux.getSma();
189 k = aux.getK();
190 h = aux.getH();
191 q = aux.getQ();
192 p = aux.getP();
193
194 // Retrograde factor
195 I = aux.getRetrogradeFactor();
196
197 // Eccentricity
198 ecc = aux.getEcc();
199
200 // Equinoctial coefficients
201 A = aux.getA();
202 B = aux.getB();
203 C = aux.getC();
204
205 // Equinoctial frame vectors
206 f = aux.getVectorF();
207 g = aux.getVectorG();
208 w = aux.getVectorW();
209
210 // Kepler mean motion
211 n = aux.getMeanMotion();
212
213 // Mean longitude
214 lm = aux.getLM();
215
216 // 1 / A
217 ooA = 1. / A;
218 // 1 / AB
219 ooAB = ooA / B;
220 // C / 2AB
221 co2AB = C * ooAB / 2.;
222 // 1 / (1 + B)
223 ooBpo = 1. / (1. + B);
224 // 2 / (n² * a)
225 ton2a = 2. / (n * n * a);
226 // mu
227 mu = aux.getMu();
228 // 1 / mu
229 ooMu = 1. / mu;
230 }
231
232 /** {@inheritDoc} */
233 @Override
234 public double[] getMeanElementRate(final SpacecraftState state) throws OrekitException {
235
236 double[] meanElementRate = new double[6];
237 // Computes the limits for the integral
238 final double[] ll = getLLimits(state);
239 // Computes integrated mean element rates if Llow < Lhigh
240 if (ll[0] < ll[1]) {
241 meanElementRate = getMeanElementRate(state, integrator, ll[0], ll[1]);
242 if (isDirty) {
243 boolean next = true;
244 for (int i = 0; i < MAX_ORDER_RANK && next; i++) {
245 final double[] meanRates = getMeanElementRate(state, new GaussQuadrature(GAUSS_ORDER[i]), ll[0], ll[1]);
246 if (getRatesDiff(meanElementRate, meanRates) < threshold) {
247 integrator = new GaussQuadrature(GAUSS_ORDER[i]);
248 next = false;
249 }
250 }
251 isDirty = false;
252 }
253 }
254 return meanElementRate;
255 }
256
257 /** Compute the acceleration due to the non conservative perturbing force.
258 *
259 * @param state current state information: date, kinematics, attitude
260 * @return the perturbing acceleration
261 * @exception OrekitException if some specific error occurs
262 */
263 protected Vector3D getAcceleration(final SpacecraftState state)
264 throws OrekitException {
265 final AccelerationRetriever retriever = new AccelerationRetriever(state);
266 contribution.addContribution(state, retriever);
267
268 return retriever.getAcceleration();
269 }
270
271 /** Compute the limits in L, the true longitude, for integration.
272 *
273 * @param state current state information: date, kinematics, attitude
274 * @return the integration limits in L
275 * @exception OrekitException if some specific error occurs
276 */
277 protected abstract double[] getLLimits(final SpacecraftState state) throws OrekitException;
278
279 /** Computes the mean equinoctial elements rates da<sub>i</sub> / dt.
280 *
281 * @param state current state
282 * @param gauss Gauss quadrature
283 * @param low lower bound of the integral interval
284 * @param high upper bound of the integral interval
285 * @return the mean element rates
286 * @throws OrekitException if some specific error occurs
287 */
288 private double[] getMeanElementRate(final SpacecraftState state,
289 final GaussQuadrature gauss,
290 final double low,
291 final double high) throws OrekitException {
292 final double[] meanElementRate = gauss.integrate(new IntegrableFunction(state, true, 0), low, high);
293 // Constant multiplier for integral
294 final double coef = 1. / (2. * FastMath.PI * B);
295 // Corrects mean element rates
296 for (int i = 0; i < 6; i++) {
297 meanElementRate[i] *= coef;
298 }
299 return meanElementRate;
300 }
301
302 /** Estimates the weighted magnitude of the difference between 2 sets of equinoctial elements rates.
303 *
304 * @param meanRef reference rates
305 * @param meanCur current rates
306 * @return estimated magnitude of weighted differences
307 */
308 private double getRatesDiff(final double[] meanRef, final double[] meanCur) {
309 double maxDiff = FastMath.abs(meanRef[0] - meanCur[0]) / a;
310 // Corrects mean element rates
311 for (int i = 1; i < meanRef.length; i++) {
312 final double diff = FastMath.abs(meanRef[i] - meanCur[i]);
313 if (maxDiff < diff) maxDiff = diff;
314 }
315 return maxDiff;
316 }
317
318 /** {@inheritDoc} */
319 @Override
320 public void registerAttitudeProvider(final AttitudeProvider provider) {
321 this.attitudeProvider = provider;
322 }
323
324 /** {@inheritDoc} */
325 @Override
326 public double[] getShortPeriodicVariations(final AbsoluteDate date,
327 final double[] meanElements) throws OrekitException {
328
329 // Build an Orbit object from the mean elements
330 final Orbit meanOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(
331 meanElements, PositionAngle.MEAN, date, this.mu, this.frame);
332
333 // Get the True longitude L
334 final double L = meanOrbit.getLv();
335
336 // Compute the center (l - λ)
337 final double center = L - meanElements[5];
338 // Compute (l - λ)²
339 final double center2 = center * center;
340
341 // Initialize short periodic variations
342 final double[] shortPeriodicVariation = new double[6];
343 for (int i = 0; i < 6; i++) {
344 shortPeriodicVariation[i] = gaussianSPCoefs.getCij(i, 0, date) +
345 center * gaussianSPCoefs.getDij(i, 1, date);
346 if (i == 5) {
347 shortPeriodicVariation[i] += center2 * gaussianSPCoefs.getDij(i, 2, date);
348 }
349 }
350
351 for (int j = 1; j <= JMAX; j++) {
352 for (int i = 0; i < 6; i++) {
353 // Get Cij and Sij
354 final double cij = gaussianSPCoefs.getCij(i, j, date);
355 final double sij = gaussianSPCoefs.getSij(i, j, date);
356
357 // add corresponding term to the short periodic variation
358 shortPeriodicVariation[i] += cij * FastMath.cos(j * L);
359 shortPeriodicVariation[i] += sij * FastMath.sin(j * L);
360 }
361 }
362
363 return shortPeriodicVariation;
364
365 }
366
367 /** {@inheritDoc} */
368 @Override
369 public void computeShortPeriodicsCoefficients(final SpacecraftState state)
370 throws OrekitException {
371
372 //Compute the coefficients
373 gaussianSPCoefs.computeCoefficients(state);
374 }
375
376 /** {@inheritDoc} */
377 @Override
378 public void resetShortPeriodicsCoefficients() {
379 if (gaussianSPCoefs != null) {
380 // reset the coefficients
381 gaussianSPCoefs.resetCoefficients();
382 }
383 }
384
385 /** Internal class for retrieving acceleration from a {@link ForceModel}. */
386 private static class AccelerationRetriever implements TimeDerivativesEquations {
387
388 /** acceleration vector. */
389 private Vector3D acceleration;
390
391 /** state. */
392 private final SpacecraftState state;
393
394 /** Simple constructor.
395 * @param state input state
396 */
397 public AccelerationRetriever(final SpacecraftState state) {
398 this.acceleration = Vector3D.ZERO;
399 this.state = state;
400 }
401
402 /** {@inheritDoc} */
403 @Override
404 public void addKeplerContribution(final double mu) {
405 }
406
407 /** {@inheritDoc} */
408 @Override
409 public void addXYZAcceleration(final double x, final double y,
410 final double z) {
411 //TODO How to be sure we are in the good frame ???
412 acceleration = new Vector3D(x, y, z);
413 }
414
415 /** {@inheritDoc} */
416 @Override
417 public void addAcceleration(final Vector3D gamma, final Frame frame)
418 throws OrekitException {
419 acceleration = frame.getTransformTo(state.getFrame(),
420 state.getDate()).transformVector(gamma);
421 }
422
423 /** {@inheritDoc} */
424 @Override
425 public void addMassDerivative(final double q) {
426 }
427
428 /** Get the acceleration vector.
429 * @return acceleration vector
430 */
431 public Vector3D getAcceleration() {
432 return acceleration;
433 }
434
435 }
436
437 /** Internal class for numerical quadrature. */
438 private class IntegrableFunction implements UnivariateVectorFunction {
439
440 /** Current state. */
441 private final SpacecraftState state;
442
443 /** Signal that this class is used to compute the values required by the mean element variations
444 * or by the short periodic element variations. */
445 private final boolean meanMode;
446
447 /** The j index.
448 * <p>
449 * Used only for short periodic variation. Ignored for mean elements variation.
450 * </p> */
451 private final int j;
452
453 /** Build a new instance.
454 * @param state current state information: date, kinematics, attitude
455 * @param meanMode if true return the value associated to the mean elements variation,
456 * if false return the values associated to the short periodic elements variation
457 * @param j the j index. used only for short periodic variation. Ignored for mean elements variation.
458 */
459 public IntegrableFunction(final SpacecraftState state, final boolean meanMode, final int j) {
460 this.state = state;
461 this.meanMode = meanMode;
462 this.j = j;
463 }
464
465 /** {@inheritDoc} */
466 @Override
467 public double[] value(final double x) {
468
469 //Compute the time difference from the true longitude difference
470 final double shiftedLm = trueToMean(x);
471 final double dLm = shiftedLm - lm;
472 final double dt = dLm / n;
473
474 final double cosL = FastMath.cos(x);
475 final double sinL = FastMath.sin(x);
476 final double roa = B * B / (1. + h * sinL + k * cosL);
477 final double roa2 = roa * roa;
478 final double r = a * roa;
479 final double X = r * cosL;
480 final double Y = r * sinL;
481 final double naob = n * a / B;
482 final double Xdot = -naob * (h + sinL);
483 final double Ydot = naob * (k + cosL);
484 final Vector3D vel = new Vector3D(Xdot, f, Ydot, g);
485
486 // Compute acceleration
487 Vector3D acc = Vector3D.ZERO;
488 try {
489
490 // shift the orbit to dt
491 final Orbit shiftedOrbit = state.getOrbit().shiftedBy(dt);
492
493 // Recompose an orbit with time held fixed to be compliant with DSST theory
494 final Orbit recomposedOrbit =
495 new EquinoctialOrbit(shiftedOrbit.getA(),
496 shiftedOrbit.getEquinoctialEx(),
497 shiftedOrbit.getEquinoctialEy(),
498 shiftedOrbit.getHx(),
499 shiftedOrbit.getHy(),
500 shiftedOrbit.getLv(),
501 PositionAngle.TRUE,
502 shiftedOrbit.getFrame(),
503 state.getDate(),
504 shiftedOrbit.getMu());
505
506 // Get the corresponding attitude
507 final Attitude recomposedAttitude =
508 attitudeProvider.getAttitude(recomposedOrbit,
509 recomposedOrbit.getDate(),
510 recomposedOrbit.getFrame());
511
512 // create shifted SpacecraftState with attitude at specified time
513 final SpacecraftState shiftedState =
514 new SpacecraftState(recomposedOrbit, recomposedAttitude, state.getMass());
515
516 acc = getAcceleration(shiftedState);
517
518 } catch (OrekitException oe) {
519 throw new OrekitExceptionWrapper(oe);
520 }
521 //Compute the derivatives of the elements by the speed
522 final double[] deriv = new double[6];
523 // da/dv
524 deriv[0] = getAoV(vel).dotProduct(acc);
525 // dex/dv
526 deriv[1] = getKoV(X, Y, Xdot, Ydot).dotProduct(acc);
527 // dey/dv
528 deriv[2] = getHoV(X, Y, Xdot, Ydot).dotProduct(acc);
529 // dhx/dv
530 deriv[3] = getQoV(X).dotProduct(acc);
531 // dhy/dv
532 deriv[4] = getPoV(Y).dotProduct(acc);
533 // dλ/dv
534 deriv[5] = getLoV(X, Y, Xdot, Ydot).dotProduct(acc);
535
536 // Compute mean elements rates
537 double[] val = null;
538 if (meanMode) {
539 val = new double[6];
540 for (int i = 0; i < 6; i++) {
541 // da<sub>i</sub>/dt
542 val[i] = roa2 * deriv[i];
543 }
544 } else {
545 val = new double [12];
546 //Compute cos(j*L) and sin(j*L);
547 final double cosjL = j == 1 ? cosL : FastMath.cos(j * x);
548 final double sinjL = j == 1 ? sinL : FastMath.sin(j * x);
549
550 for (int i = 0; i < 6; i++) {
551 // da<sub>i</sub>/dv * cos(jL)
552 val[i] = cosjL * deriv[i];
553 // da<sub>i</sub>/dv * sin(jL)
554 val[i + 6] = sinjL * deriv[i];
555 }
556 }
557 return val;
558 }
559
560 /** Converts true longitude to eccentric longitude.
561 * @param lv True longitude
562 * @return Eccentric longitude
563 */
564 private double trueToEccentric (final double lv) {
565 final double cosLv = FastMath.cos(lv);
566 final double sinLv = FastMath.sin(lv);
567 final double num = h * cosLv - k * sinLv;
568 final double den = B + 1 + k * cosLv + h * sinLv;
569 return lv + 2 * FastMath.atan(num / den);
570 }
571
572 /** Converts eccentric longitude to mean longitude.
573 * @param le Eccentric longitude
574 * @return Mean longitude
575 */
576 private double eccentricToMean (final double le) {
577 return le - k * FastMath.sin(le) + h * FastMath.cos(le);
578 }
579
580 /** Converts true longitude to mean longitude.
581 * @param lv True longitude
582 * @return Eccentric longitude
583 */
584 private double trueToMean (final double lv) {
585 return eccentricToMean(trueToEccentric(lv));
586 }
587
588 /** Compute δa/δv.
589 * @param vel satellite velocity
590 * @return δa/δv
591 */
592 private Vector3D getAoV(final Vector3D vel) {
593 return new Vector3D(ton2a, vel);
594 }
595
596 /** Compute δh/δv.
597 * @param X satellite position component along f, equinoctial reference frame 1st vector
598 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
599 * @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
600 * @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
601 * @return δh/δv
602 */
603 private Vector3D getHoV(final double X, final double Y, final double Xdot, final double Ydot) {
604 final double kf = (2. * Xdot * Y - X * Ydot) * ooMu;
605 final double kg = X * Xdot * ooMu;
606 final double kw = k * (I * q * Y - p * X) * ooAB;
607 return new Vector3D(kf, f, -kg, g, kw, w);
608 }
609
610 /** Compute δk/δv.
611 * @param X satellite position component along f, equinoctial reference frame 1st vector
612 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
613 * @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
614 * @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
615 * @return δk/δv
616 */
617 private Vector3D getKoV(final double X, final double Y, final double Xdot, final double Ydot) {
618 final double kf = Y * Ydot * ooMu;
619 final double kg = (2. * X * Ydot - Xdot * Y) * ooMu;
620 final double kw = h * (I * q * Y - p * X) * ooAB;
621 return new Vector3D(-kf, f, kg, g, -kw, w);
622 }
623
624 /** Compute δp/δv.
625 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
626 * @return δp/δv
627 */
628 private Vector3D getPoV(final double Y) {
629 return new Vector3D(co2AB * Y, w);
630 }
631
632 /** Compute δq/δv.
633 * @param X satellite position component along f, equinoctial reference frame 1st vector
634 * @return δq/δv
635 */
636 private Vector3D getQoV(final double X) {
637 return new Vector3D(I * co2AB * X, w);
638 }
639
640 /** Compute δλ/δv.
641 * @param X satellite position component along f, equinoctial reference frame 1st vector
642 * @param Y satellite position component along g, equinoctial reference frame 2nd vector
643 * @param Xdot satellite velocity component along f, equinoctial reference frame 1st vector
644 * @param Ydot satellite velocity component along g, equinoctial reference frame 2nd vector
645 * @return δλ/δv
646 */
647 private Vector3D getLoV(final double X, final double Y, final double Xdot, final double Ydot) {
648 final Vector3D pos = new Vector3D(X, f, Y, g);
649 final Vector3D v2 = new Vector3D(k, getHoV(X, Y, Xdot, Ydot), -h, getKoV(X, Y, Xdot, Ydot));
650 return new Vector3D(-2. * ooA, pos, ooBpo, v2, (I * q * Y - p * X) * ooA, w);
651 }
652
653 }
654
655 /** Class used to {@link #integrate(UnivariateVectorFunction, double, double) integrate}
656 * a {@link org.apache.commons.math3.analysis.UnivariateVectorFunction function}
657 * of the orbital elements using the Gaussian quadrature rule to get the acceleration.
658 */
659 private static class GaussQuadrature {
660
661 // CHECKSTYLE: stop NoWhitespaceAfter
662
663 // Points and weights for the available quadrature orders
664
665 /** Points for quadrature of order 12. */
666 private static final double[] P_12 = {
667 -0.98156063424671910000,
668 -0.90411725637047490000,
669 -0.76990267419430470000,
670 -0.58731795428661740000,
671 -0.36783149899818024000,
672 -0.12523340851146890000,
673 0.12523340851146890000,
674 0.36783149899818024000,
675 0.58731795428661740000,
676 0.76990267419430470000,
677 0.90411725637047490000,
678 0.98156063424671910000
679 };
680
681 /** Weights for quadrature of order 12. */
682 private static final double[] W_12 = {
683 0.04717533638651220000,
684 0.10693932599531830000,
685 0.16007832854334633000,
686 0.20316742672306584000,
687 0.23349253653835478000,
688 0.24914704581340286000,
689 0.24914704581340286000,
690 0.23349253653835478000,
691 0.20316742672306584000,
692 0.16007832854334633000,
693 0.10693932599531830000,
694 0.04717533638651220000
695 };
696
697 /** Points for quadrature of order 16. */
698 private static final double[] P_16 = {
699 -0.98940093499164990000,
700 -0.94457502307323260000,
701 -0.86563120238783160000,
702 -0.75540440835500310000,
703 -0.61787624440264380000,
704 -0.45801677765722737000,
705 -0.28160355077925890000,
706 -0.09501250983763745000,
707 0.09501250983763745000,
708 0.28160355077925890000,
709 0.45801677765722737000,
710 0.61787624440264380000,
711 0.75540440835500310000,
712 0.86563120238783160000,
713 0.94457502307323260000,
714 0.98940093499164990000
715 };
716
717 /** Weights for quadrature of order 16. */
718 private static final double[] W_16 = {
719 0.02715245941175405800,
720 0.06225352393864777000,
721 0.09515851168249283000,
722 0.12462897125553388000,
723 0.14959598881657685000,
724 0.16915651939500256000,
725 0.18260341504492360000,
726 0.18945061045506847000,
727 0.18945061045506847000,
728 0.18260341504492360000,
729 0.16915651939500256000,
730 0.14959598881657685000,
731 0.12462897125553388000,
732 0.09515851168249283000,
733 0.06225352393864777000,
734 0.02715245941175405800
735 };
736
737 /** Points for quadrature of order 20. */
738 private static final double[] P_20 = {
739 -0.99312859918509490000,
740 -0.96397192727791390000,
741 -0.91223442825132600000,
742 -0.83911697182221890000,
743 -0.74633190646015080000,
744 -0.63605368072651510000,
745 -0.51086700195082700000,
746 -0.37370608871541955000,
747 -0.22778585114164507000,
748 -0.07652652113349734000,
749 0.07652652113349734000,
750 0.22778585114164507000,
751 0.37370608871541955000,
752 0.51086700195082700000,
753 0.63605368072651510000,
754 0.74633190646015080000,
755 0.83911697182221890000,
756 0.91223442825132600000,
757 0.96397192727791390000,
758 0.99312859918509490000
759 };
760
761 /** Weights for quadrature of order 20. */
762 private static final double[] W_20 = {
763 0.01761400713915226400,
764 0.04060142980038684000,
765 0.06267204833410904000,
766 0.08327674157670477000,
767 0.10193011981724048000,
768 0.11819453196151844000,
769 0.13168863844917678000,
770 0.14209610931838212000,
771 0.14917298647260380000,
772 0.15275338713072600000,
773 0.15275338713072600000,
774 0.14917298647260380000,
775 0.14209610931838212000,
776 0.13168863844917678000,
777 0.11819453196151844000,
778 0.10193011981724048000,
779 0.08327674157670477000,
780 0.06267204833410904000,
781 0.04060142980038684000,
782 0.01761400713915226400
783 };
784
785 /** Points for quadrature of order 24. */
786 private static final double[] P_24 = {
787 -0.99518721999702130000,
788 -0.97472855597130950000,
789 -0.93827455200273270000,
790 -0.88641552700440100000,
791 -0.82000198597390300000,
792 -0.74012419157855440000,
793 -0.64809365193697550000,
794 -0.54542147138883950000,
795 -0.43379350762604520000,
796 -0.31504267969616340000,
797 -0.19111886747361634000,
798 -0.06405689286260563000,
799 0.06405689286260563000,
800 0.19111886747361634000,
801 0.31504267969616340000,
802 0.43379350762604520000,
803 0.54542147138883950000,
804 0.64809365193697550000,
805 0.74012419157855440000,
806 0.82000198597390300000,
807 0.88641552700440100000,
808 0.93827455200273270000,
809 0.97472855597130950000,
810 0.99518721999702130000
811 };
812
813 /** Weights for quadrature of order 24. */
814 private static final double[] W_24 = {
815 0.01234122979998733500,
816 0.02853138862893380600,
817 0.04427743881741981000,
818 0.05929858491543691500,
819 0.07334648141108027000,
820 0.08619016153195320000,
821 0.09761865210411391000,
822 0.10744427011596558000,
823 0.11550566805372553000,
824 0.12167047292780335000,
825 0.12583745634682825000,
826 0.12793819534675221000,
827 0.12793819534675221000,
828 0.12583745634682825000,
829 0.12167047292780335000,
830 0.11550566805372553000,
831 0.10744427011596558000,
832 0.09761865210411391000,
833 0.08619016153195320000,
834 0.07334648141108027000,
835 0.05929858491543691500,
836 0.04427743881741981000,
837 0.02853138862893380600,
838 0.01234122979998733500
839 };
840
841 /** Points for quadrature of order 32. */
842 private static final double[] P_32 = {
843 -0.99726386184948160000,
844 -0.98561151154526840000,
845 -0.96476225558750640000,
846 -0.93490607593773970000,
847 -0.89632115576605220000,
848 -0.84936761373256990000,
849 -0.79448379596794250000,
850 -0.73218211874028970000,
851 -0.66304426693021520000,
852 -0.58771575724076230000,
853 -0.50689990893222950000,
854 -0.42135127613063540000,
855 -0.33186860228212767000,
856 -0.23928736225213710000,
857 -0.14447196158279646000,
858 -0.04830766568773831000,
859 0.04830766568773831000,
860 0.14447196158279646000,
861 0.23928736225213710000,
862 0.33186860228212767000,
863 0.42135127613063540000,
864 0.50689990893222950000,
865 0.58771575724076230000,
866 0.66304426693021520000,
867 0.73218211874028970000,
868 0.79448379596794250000,
869 0.84936761373256990000,
870 0.89632115576605220000,
871 0.93490607593773970000,
872 0.96476225558750640000,
873 0.98561151154526840000,
874 0.99726386184948160000
875 };
876
877 /** Weights for quadrature of order 32. */
878 private static final double[] W_32 = {
879 0.00701861000947013600,
880 0.01627439473090571200,
881 0.02539206530926214200,
882 0.03427386291302141000,
883 0.04283589802222658600,
884 0.05099805926237621600,
885 0.05868409347853559000,
886 0.06582222277636193000,
887 0.07234579410884862000,
888 0.07819389578707042000,
889 0.08331192422694673000,
890 0.08765209300440380000,
891 0.09117387869576390000,
892 0.09384439908080441000,
893 0.09563872007927487000,
894 0.09654008851472784000,
895 0.09654008851472784000,
896 0.09563872007927487000,
897 0.09384439908080441000,
898 0.09117387869576390000,
899 0.08765209300440380000,
900 0.08331192422694673000,
901 0.07819389578707042000,
902 0.07234579410884862000,
903 0.06582222277636193000,
904 0.05868409347853559000,
905 0.05099805926237621600,
906 0.04283589802222658600,
907 0.03427386291302141000,
908 0.02539206530926214200,
909 0.01627439473090571200,
910 0.00701861000947013600
911 };
912
913 /** Points for quadrature of order 40. */
914 private static final double[] P_40 = {
915 -0.99823770971055930000,
916 -0.99072623869945710000,
917 -0.97725994998377420000,
918 -0.95791681921379170000,
919 -0.93281280827867660000,
920 -0.90209880696887420000,
921 -0.86595950321225960000,
922 -0.82461223083331170000,
923 -0.77830565142651940000,
924 -0.72731825518992710000,
925 -0.67195668461417960000,
926 -0.61255388966798030000,
927 -0.54946712509512820000,
928 -0.48307580168617870000,
929 -0.41377920437160500000,
930 -0.34199409082575850000,
931 -0.26815218500725370000,
932 -0.19269758070137110000,
933 -0.11608407067525522000,
934 -0.03877241750605081600,
935 0.03877241750605081600,
936 0.11608407067525522000,
937 0.19269758070137110000,
938 0.26815218500725370000,
939 0.34199409082575850000,
940 0.41377920437160500000,
941 0.48307580168617870000,
942 0.54946712509512820000,
943 0.61255388966798030000,
944 0.67195668461417960000,
945 0.72731825518992710000,
946 0.77830565142651940000,
947 0.82461223083331170000,
948 0.86595950321225960000,
949 0.90209880696887420000,
950 0.93281280827867660000,
951 0.95791681921379170000,
952 0.97725994998377420000,
953 0.99072623869945710000,
954 0.99823770971055930000
955 };
956
957 /** Weights for quadrature of order 40. */
958 private static final double[] W_40 = {
959 0.00452127709853309800,
960 0.01049828453115270400,
961 0.01642105838190797300,
962 0.02224584919416689000,
963 0.02793700698002338000,
964 0.03346019528254786500,
965 0.03878216797447199000,
966 0.04387090818567333000,
967 0.04869580763507221000,
968 0.05322784698393679000,
969 0.05743976909939157000,
970 0.06130624249292891000,
971 0.06480401345660108000,
972 0.06791204581523394000,
973 0.07061164739128681000,
974 0.07288658239580408000,
975 0.07472316905796833000,
976 0.07611036190062619000,
977 0.07703981816424793000,
978 0.07750594797842482000,
979 0.07750594797842482000,
980 0.07703981816424793000,
981 0.07611036190062619000,
982 0.07472316905796833000,
983 0.07288658239580408000,
984 0.07061164739128681000,
985 0.06791204581523394000,
986 0.06480401345660108000,
987 0.06130624249292891000,
988 0.05743976909939157000,
989 0.05322784698393679000,
990 0.04869580763507221000,
991 0.04387090818567333000,
992 0.03878216797447199000,
993 0.03346019528254786500,
994 0.02793700698002338000,
995 0.02224584919416689000,
996 0.01642105838190797300,
997 0.01049828453115270400,
998 0.00452127709853309800
999 };
1000
1001 /** Points for quadrature of order 48. */
1002 private static final double[] P_48 = {
1003 -0.99877100725242610000,
1004 -0.99353017226635080000,
1005 -0.98412458372282700000,
1006 -0.97059159254624720000,
1007 -0.95298770316043080000,
1008 -0.93138669070655440000,
1009 -0.90587913671556960000,
1010 -0.87657202027424800000,
1011 -0.84358826162439350000,
1012 -0.80706620402944250000,
1013 -0.76715903251574020000,
1014 -0.72403413092381470000,
1015 -0.67787237963266400000,
1016 -0.62886739677651370000,
1017 -0.57722472608397270000,
1018 -0.52316097472223300000,
1019 -0.46690290475095840000,
1020 -0.40868648199071680000,
1021 -0.34875588629216070000,
1022 -0.28736248735545555000,
1023 -0.22476379039468908000,
1024 -0.16122235606889174000,
1025 -0.09700469920946270000,
1026 -0.03238017096286937000,
1027 0.03238017096286937000,
1028 0.09700469920946270000,
1029 0.16122235606889174000,
1030 0.22476379039468908000,
1031 0.28736248735545555000,
1032 0.34875588629216070000,
1033 0.40868648199071680000,
1034 0.46690290475095840000,
1035 0.52316097472223300000,
1036 0.57722472608397270000,
1037 0.62886739677651370000,
1038 0.67787237963266400000,
1039 0.72403413092381470000,
1040 0.76715903251574020000,
1041 0.80706620402944250000,
1042 0.84358826162439350000,
1043 0.87657202027424800000,
1044 0.90587913671556960000,
1045 0.93138669070655440000,
1046 0.95298770316043080000,
1047 0.97059159254624720000,
1048 0.98412458372282700000,
1049 0.99353017226635080000,
1050 0.99877100725242610000
1051 };
1052
1053 /** Weights for quadrature of order 48. */
1054 private static final double[] W_48 = {
1055 0.00315334605230596250,
1056 0.00732755390127620800,
1057 0.01147723457923446900,
1058 0.01557931572294386600,
1059 0.01961616045735556700,
1060 0.02357076083932435600,
1061 0.02742650970835688000,
1062 0.03116722783279807000,
1063 0.03477722256477045000,
1064 0.03824135106583080600,
1065 0.04154508294346483000,
1066 0.04467456085669424000,
1067 0.04761665849249054000,
1068 0.05035903555385448000,
1069 0.05289018948519365000,
1070 0.05519950369998416500,
1071 0.05727729210040315000,
1072 0.05911483969839566000,
1073 0.06070443916589384000,
1074 0.06203942315989268000,
1075 0.06311419228625403000,
1076 0.06392423858464817000,
1077 0.06446616443595010000,
1078 0.06473769681268386000,
1079 0.06473769681268386000,
1080 0.06446616443595010000,
1081 0.06392423858464817000,
1082 0.06311419228625403000,
1083 0.06203942315989268000,
1084 0.06070443916589384000,
1085 0.05911483969839566000,
1086 0.05727729210040315000,
1087 0.05519950369998416500,
1088 0.05289018948519365000,
1089 0.05035903555385448000,
1090 0.04761665849249054000,
1091 0.04467456085669424000,
1092 0.04154508294346483000,
1093 0.03824135106583080600,
1094 0.03477722256477045000,
1095 0.03116722783279807000,
1096 0.02742650970835688000,
1097 0.02357076083932435600,
1098 0.01961616045735556700,
1099 0.01557931572294386600,
1100 0.01147723457923446900,
1101 0.00732755390127620800,
1102 0.00315334605230596250
1103 };
1104 // CHECKSTYLE: resume NoWhitespaceAfter
1105
1106 /** Node points. */
1107 private final double[] nodePoints;
1108
1109 /** Node weights. */
1110 private final double[] nodeWeights;
1111
1112 /** Creates a Gauss integrator of the given order.
1113 *
1114 * @param numberOfPoints Order of the integration rule.
1115 */
1116 public GaussQuadrature(final int numberOfPoints) {
1117
1118 switch(numberOfPoints) {
1119 case 12 :
1120 this.nodePoints = P_12.clone();
1121 this.nodeWeights = W_12.clone();
1122 break;
1123 case 16 :
1124 this.nodePoints = P_16.clone();
1125 this.nodeWeights = W_16.clone();
1126 break;
1127 case 20 :
1128 this.nodePoints = P_20.clone();
1129 this.nodeWeights = W_20.clone();
1130 break;
1131 case 24 :
1132 this.nodePoints = P_24.clone();
1133 this.nodeWeights = W_24.clone();
1134 break;
1135 case 32 :
1136 this.nodePoints = P_32.clone();
1137 this.nodeWeights = W_32.clone();
1138 break;
1139 case 40 :
1140 this.nodePoints = P_40.clone();
1141 this.nodeWeights = W_40.clone();
1142 break;
1143 case 48 :
1144 default :
1145 this.nodePoints = P_48.clone();
1146 this.nodeWeights = W_48.clone();
1147 break;
1148 }
1149
1150 }
1151
1152 /** Integrates a given function on the given interval.
1153 *
1154 * @param f Function to integrate.
1155 * @param lowerBound Lower bound of the integration interval.
1156 * @param upperBound Upper bound of the integration interval.
1157 * @return the integral of the weighted function.
1158 */
1159 public double[] integrate(final UnivariateVectorFunction f,
1160 final double lowerBound, final double upperBound) {
1161
1162 final double[] adaptedPoints = nodePoints.clone();
1163 final double[] adaptedWeights = nodeWeights.clone();
1164 transform(adaptedPoints, adaptedWeights, lowerBound, upperBound);
1165 return basicIntegrate(f, adaptedPoints, adaptedWeights);
1166 }
1167
1168 /** Performs a change of variable so that the integration
1169 * can be performed on an arbitrary interval {@code [a, b]}.
1170 * <p>
1171 * It is assumed that the natural interval is {@code [-1, 1]}.
1172 * </p>
1173 *
1174 * @param points Points to adapt to the new interval.
1175 * @param weights Weights to adapt to the new interval.
1176 * @param a Lower bound of the integration interval.
1177 * @param b Lower bound of the integration interval.
1178 */
1179 private void transform(final double[] points, final double[] weights,
1180 final double a, final double b) {
1181 // Scaling
1182 final double scale = (b - a) / 2;
1183 final double shift = a + scale;
1184 for (int i = 0; i < points.length; i++) {
1185 points[i] = points[i] * scale + shift;
1186 weights[i] *= scale;
1187 }
1188 }
1189
1190 /** Returns an estimate of the integral of {@code f(x) * w(x)},
1191 * where {@code w} is a weight function that depends on the actual
1192 * flavor of the Gauss integration scheme.
1193 *
1194 * @param f Function to integrate.
1195 * @param points Nodes.
1196 * @param weights Nodes weights.
1197 * @return the integral of the weighted function.
1198 */
1199 private double[] basicIntegrate(final UnivariateVectorFunction f,
1200 final double[] points,
1201 final double[] weights) {
1202 double x = points[0];
1203 double w = weights[0];
1204 double[] v = f.value(x);
1205 final double[] y = new double[v.length];
1206 for (int j = 0; j < v.length; j++) {
1207 y[j] = w * v[j];
1208 }
1209 final double[] t = y.clone();
1210 final double[] c = new double[v.length];
1211 final double[] s = t.clone();
1212 for (int i = 1; i < points.length; i++) {
1213 x = points[i];
1214 w = weights[i];
1215 v = f.value(x);
1216 for (int j = 0; j < v.length; j++) {
1217 y[j] = w * v[j] - c[j];
1218 t[j] = s[j] + y[j];
1219 c[j] = (t[j] - s[j]) - y[j];
1220 s[j] = t[j];
1221 }
1222 }
1223 return s;
1224 }
1225
1226 }
1227
1228 /** Compute the C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> coefficients.
1229 * <p>
1230 * Those coefficients are given in Danielson paper by expression 4.4-(6)
1231 * </p>
1232 * @author Petre Bazavan
1233 * @author Lucian Barbulescu
1234 */
1235 private class FourierCjSjCoefficients {
1236
1237 /** Maximum possible value for j. */
1238 private final int jMax;
1239
1240 /** The C<sub>i</sub><sup>j</sup> coefficients.
1241 * <p>
1242 * the index i corresponds to the following elements: <br/>
1243 * - 0 for a <br>
1244 * - 1 for k <br>
1245 * - 2 for h <br>
1246 * - 3 for q <br>
1247 * - 4 for p <br>
1248 * - 5 for λ <br>
1249 * </p>
1250 */
1251 private final double[][] cCoef;
1252
1253 /** The C<sub>i</sub><sup>j</sup> coefficients.
1254 * <p>
1255 * the index i corresponds to the following elements: <br/>
1256 * - 0 for a <br>
1257 * - 1 for k <br>
1258 * - 2 for h <br>
1259 * - 3 for q <br>
1260 * - 4 for p <br>
1261 * - 5 for λ <br>
1262 * </p>
1263 */
1264 private final double[][] sCoef;
1265
1266 /** Standard constructor.
1267 * @param state the current state
1268 * @param jMax maximum value for j
1269 * @throws OrekitException in case of an error
1270 */
1271 public FourierCjSjCoefficients(final SpacecraftState state, final int jMax) throws OrekitException {
1272 //Initialise the fields
1273 this.jMax = jMax;
1274
1275 //Allocate the arrays
1276 cCoef = new double[jMax + 1][6];
1277 sCoef = new double[jMax + 1][6];
1278
1279 //Compute the coefficients
1280 computeCoefficients(state);
1281 }
1282
1283 /**
1284 * Compute the Fourrier coefficients.
1285 * <p>
1286 * Only the C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup> coefficients need to be computed
1287 * as D<sub>i</sub><sup>m</sup> is always 0.
1288 * </p>
1289 * @param state the current state
1290 * @throws OrekitException in case of an error
1291 */
1292 private void computeCoefficients(final SpacecraftState state) throws OrekitException {
1293 // Computes the limits for the integral
1294 final double[] ll = getLLimits(state);
1295 // Computes integrated mean element rates if Llow < Lhigh
1296 if (ll[0] < ll[1]) {
1297 //Compute 1 / PI
1298 final double ooPI = 1 / FastMath.PI;
1299
1300 // loop through all values of j
1301 for (int j = 0; j <= jMax; j++) {
1302 final double[] curentCoefficients = integrator.integrate(new IntegrableFunction(state, false, j), ll[0], ll[1]);
1303
1304 //divide by PI and set the values for the coefficients
1305 for (int i = 0; i < 6; i++) {
1306 cCoef[j][i] = ooPI * curentCoefficients[i];
1307 sCoef[j][i] = ooPI * curentCoefficients[i + 6];
1308 }
1309 }
1310 }
1311 }
1312
1313 /** Get the coefficient C<sub>i</sub><sup>j</sup>.
1314 * @param i i index - corresponds to the required variation
1315 * @param j j index
1316 * @return the coefficient C<sub>i</sub><sup>j</sup>
1317 */
1318 public double getCij(final int i, final int j) {
1319 return cCoef[j][i];
1320 }
1321
1322 /** Get the coefficient S<sub>i</sub><sup>j</sup>.
1323 * @param i i index - corresponds to the required variation
1324 * @param j j index
1325 * @return the coefficient S<sub>i</sub><sup>j</sup>
1326 */
1327 public double getSij(final int i, final int j) {
1328 return sCoef[j][i];
1329 }
1330 }
1331
1332 /** This class handles the short periodic coefficients described in Danielson 2.5.3-26.
1333 *
1334 * <p>
1335 * The value of M is 0. Also, since the values of the Fourier coefficient D<sub>i</sub><sup>m</sup> is 0
1336 * then the values of the coefficients D<sub>i</sub><sup>m</sup> for m > 2 are also 0.
1337 * </p>
1338 * @author Petre Bazavan
1339 * @author Lucian Barbulescu
1340 *
1341 */
1342 private class GaussianShortPeriodicCoefficients {
1343
1344 /** Maximum value for j index. */
1345 private final int jMax;
1346
1347 /**The coefficients D<sub>i</sub><sup>j</sup>.
1348 * <p>
1349 * Only for j = 1 and j = 2 the coefficients are not 0. <br>
1350 * i corresponds to the equinoctial element, as follows:
1351 * - i=0 for a <br/>
1352 * - i=1 for k <br/>
1353 * - i=2 for h <br/>
1354 * - i=3 for q <br/>
1355 * - i=4 for p <br/>
1356 * - i=5 for λ <br/>
1357 * </p>
1358 */
1359 private final ShortPeriodicsInterpolatedCoefficient[][] dij;
1360
1361 /** The coefficients C<sub>i</sub><sup>j</sup>.
1362 * <p>
1363 * The index order is cij[j][i] <br/>
1364 * i corresponds to the equinoctial element, as follows: <br/>
1365 * - i=0 for a <br/>
1366 * - i=1 for k <br/>
1367 * - i=2 for h <br/>
1368 * - i=3 for q <br/>
1369 * - i=4 for p <br/>
1370 * - i=5 for λ <br/>
1371 * </p>
1372 */
1373 private final ShortPeriodicsInterpolatedCoefficient[][] cij;
1374
1375 /** The coefficients S<sub>i</sub><sup>j</sup>.
1376 * <p>
1377 * The index order is sij[j][i] <br/>
1378 * i corresponds to the equinoctial element, as follows: <br/>
1379 * - i=0 for a <br/>
1380 * - i=1 for k <br/>
1381 * - i=2 for h <br/>
1382 * - i=3 for q <br/>
1383 * - i=4 for p <br/>
1384 * - i=5 for λ <br/>
1385 * </p>
1386 */
1387 private final ShortPeriodicsInterpolatedCoefficient[][] sij;
1388
1389 /** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients.
1390 * <p>
1391 * Index 0 corresponds to ρ, index 1 corresponds to σ
1392 * Used to compute the U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients.
1393 * </p>
1394 */
1395 private final double[][] currentRhoSigmaj;
1396
1397 /** Constructor.
1398 * @param jMax maximum value for j index
1399 * @param interpolationPoints number of points used in the interpolation process
1400 */
1401 public GaussianShortPeriodicCoefficients(final int jMax, final int interpolationPoints) {
1402 //Initialise fields
1403 this.jMax = jMax;
1404
1405 this.dij = new ShortPeriodicsInterpolatedCoefficient[3][6];
1406
1407 this.cij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1][6];
1408 this.sij = new ShortPeriodicsInterpolatedCoefficient[jMax + 1][6];
1409
1410 this.currentRhoSigmaj = new double[2][3 * jMax + 1];
1411
1412 // Initialise the C<sub>i</sub><sup>j</sup>, S<sub>i</sub><sup>j</sup> and D<sub>i</sub><sup>j</sup> coefficients
1413 for (int j = 0; j <= jMax; j++) {
1414 for (int i = 0; i < 6; i++) {
1415 this.cij[j][i] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
1416 if (j > 0) {
1417 this.sij[j][i] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
1418 }
1419 // Initialise only the non-zero D<sub>i</sub><sup>j</sup> coefficients
1420 if (j == 1 || (j == 2 && i == 5)) {
1421 this.dij[j][i] = new ShortPeriodicsInterpolatedCoefficient(interpolationPoints);
1422 }
1423 }
1424 }
1425 }
1426
1427 /** Compute the short periodic coefficients.
1428 *
1429 * @param state current state information: date, kinematics, attitude
1430 * @throws OrekitException if an error occurs
1431 */
1432 public void computeCoefficients(final SpacecraftState state)
1433 throws OrekitException {
1434
1435 // get the current date
1436 final AbsoluteDate date = state.getDate();
1437
1438 // Compute ρ<sub>j</sub> and σ<sub>j</sub>
1439 computeRhoSigmaCoefficients(date);
1440
1441 // Compute the Fourier coefficients
1442 final FourierCjSjCoefficients fourierCjSj = new FourierCjSjCoefficients(state, jMax);
1443
1444 // Compute the required U and V coefficients
1445 final UijVijCoefficients uijvij = new UijVijCoefficients(currentRhoSigmaj, fourierCjSj, jMax);
1446
1447 // compute the k₂⁰ coefficient
1448 final double k20 = computeK20(jMax);
1449
1450 // 1. / n
1451 final double oon = 1. / n;
1452 // 3. / (2 * a * n)
1453 final double to2an = 1.5 * oon / a;
1454 // 3. / (4 * a * n)
1455 final double to4an = to2an / 2;
1456
1457 // Compute the coefficients for each element
1458 for (int i = 0; i < 6; i++) {
1459
1460 // compute D<sub>i</sub>¹ and D<sub>i</sub>² (all others are 0)
1461 double di1 = -oon * fourierCjSj.getCij(i, 0);
1462 if (i == 5) {
1463 di1 += to2an * uijvij.getU1(0, 0);
1464 }
1465 double di2 = 0.;
1466 if (i == 5) {
1467 di2 += -to4an * fourierCjSj.getCij(0, 0);
1468 }
1469
1470 //the C<sub>i</sub>⁰ is computed based on all others
1471 double currentCi0 = -di2 * k20;
1472
1473 for (int j = 1; j <= jMax; j++) {
1474 // compute the current C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup>
1475 double currentCij = oon * uijvij.getU1(j, i);
1476 if (i == 5) {
1477 currentCij += -to2an * uijvij.getU2(j);
1478 }
1479 double currentSij = oon * uijvij.getV1(j, i);
1480 if (i == 5) {
1481 currentSij += -to2an * uijvij.getV2(j);
1482 }
1483
1484 // add the computed coefficients to C<sub>i</sub>⁰
1485 currentCi0 += -(currentCij * currentRhoSigmaj[0][j] + currentSij * currentRhoSigmaj[1][j]);
1486
1487 // add the values to the interpolators
1488 cij[j][i].addGridPoint(date, currentCij);
1489 sij[j][i].addGridPoint(date, currentSij);
1490 }
1491
1492 // add the computed values to the interpolators
1493 cij[0][i].addGridPoint(date, currentCi0);
1494 dij[1][i].addGridPoint(date, di1);
1495 if (i == 5) {
1496 dij[2][i].addGridPoint(date, di2);
1497 }
1498 }
1499 }
1500 /** Reset the coefficients.
1501 * <p>
1502 * For each coefficient, clear history of computed points
1503 * </p>
1504 */
1505 public void resetCoefficients() {
1506
1507 for (int j = 0; j <= jMax; j++) {
1508 for (int i = 0; i < 6; i++) {
1509 this.cij[j][i].clearHistory();
1510 if (j > 0) {
1511 this.sij[j][i].clearHistory();
1512 }
1513 if (j == 1 || (j == 2 && i == 5)) {
1514 this.dij[j][i].clearHistory();
1515 }
1516 }
1517 }
1518 }
1519
1520 /**
1521 * Compute the auxiliary quantities ρ<sub>j</sub> and σ<sub>j</sub>.
1522 * <p>
1523 * The expressions used are equations 2.5.3-(4) from the Danielson paper. <br/>
1524 * ρ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>C<sub>j</sub>(k, h) <br/>
1525 * σ<sub>j</sub> = (1+jB)(-b)<sup>j</sup>S<sub>j</sub>(k, h) <br/>
1526 * </p>
1527 * @param date current date
1528 */
1529 private void computeRhoSigmaCoefficients(final AbsoluteDate date) {
1530 final CjSjCoefficient cjsjKH = new CjSjCoefficient(k, h);
1531 final double b = 1. / (1 + B);
1532
1533 // (-b)<sup>j</sup>
1534 double mbtj = 1;
1535
1536 for (int j = 1; j <= 3 * jMax; j++) {
1537
1538 //Compute current rho and sigma;
1539 mbtj *= -b;
1540 final double coef = (1 + j * B) * mbtj;
1541 currentRhoSigmaj[0][j] = coef * cjsjKH.getCj(j);
1542 currentRhoSigmaj[1][j] = coef * cjsjKH.getSj(j);
1543 }
1544 }
1545
1546 /** Compute the coefficient k₂⁰ by using the equation
1547 * 2.5.3-(9a) from Danielson.
1548 * <p>
1549 * After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:<br>
1550 * k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
1551 * </p>
1552 * @param kMax max value fot k index
1553 * @return the coefficient k₂⁰
1554 */
1555 private double computeK20(final int kMax) {
1556 double k20 = 0.;
1557
1558 for (int kIndex = 1; kIndex <= kMax; kIndex++) {
1559 // After inserting 2.5.3-(8) into 2.5.3-(9a) the result becomes:
1560 //k₂⁰ = Σ<sub>k=1</sub><sup>kMax</sup>[(2 / k²) * (σ<sub>k</sub>² + ρ<sub>k</sub>²)]
1561 double currentTerm = currentRhoSigmaj[1][kIndex] * currentRhoSigmaj[1][kIndex] +
1562 currentRhoSigmaj[0][kIndex] * currentRhoSigmaj[0][kIndex];
1563
1564 //multiply by 2 / k²
1565 currentTerm *= 2. / (kIndex * kIndex);
1566
1567 // add the term to the result
1568 k20 += currentTerm;
1569 }
1570
1571 return k20;
1572 }
1573
1574 /** Get C<sub>i</sub><sup>j</sup>.
1575 *
1576 * @param i i index
1577 * @param j j index
1578 * @param date the date
1579 * @return C<sub>i</sub><sup>j</sup>
1580 */
1581 public double getCij(final int i, final int j, final AbsoluteDate date) {
1582 return cij[j][i].value(date);
1583 }
1584
1585 /** Get S<sub>i</sub><sup>j</sup>.
1586 *
1587 * @param i i index
1588 * @param j j index
1589 * @param date the date
1590 * @return S<sub>i</sub><sup>j</sup>
1591 */
1592 public double getSij(final int i, final int j, final AbsoluteDate date) {
1593 return sij[j][i].value(date);
1594 }
1595
1596 /** Get D<sub>i</sub><sup>j</sup>.
1597 * @param i i index
1598 * @param j j index
1599 * @param date target date
1600 * @return D<sub>i</sub><sup>j</sup>
1601 */
1602 public double getDij(final int i, final int j, final AbsoluteDate date) {
1603 return dij[j][i].value(date);
1604 }
1605 }
1606
1607 /** The U<sub>i</sub><sup>j</sup> and V<sub>i</sub><sup>j</sup> coefficients described by
1608 * equations 2.5.3-(21) and 2.5.3-(22) from Danielson.
1609 * <p>
1610 * The index i takes only the values 1 and 2<br>
1611 * For U only the index 0 for j is used.
1612 * </p>
1613 *
1614 * @author Petre Bazavan
1615 * @author Lucian Barbulescu
1616 */
1617 private static class UijVijCoefficients {
1618
1619 /** The U₁<sup>j</sup> coefficients.
1620 * <p>
1621 * The first index identifies the Fourier coefficients used<br>
1622 * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
1623 * The only exception is when j = 0 when only the coefficient for fourier index = 1 (i == 0) is needed.<br>
1624 * Also, for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
1625 * to compute the coefficients U₂<sup>j</sup>
1626 * </p>
1627 */
1628 private final double[][] u1ij;
1629
1630 /** The V₁<sup>j</sup> coefficients.
1631 * <p>
1632 * The first index identifies the Fourier coefficients used<br>
1633 * Those coefficients are computed for all Fourier C<sub>i</sub><sup>j</sup> and S<sub>i</sub><sup>j</sup><br>
1634 * for fourier index = 1 (i == 0), the coefficients up to 2 * jMax are computed, because are required
1635 * to compute the coefficients V₂<sup>j</sup>
1636 * </p>
1637 */
1638 private final double[][] v1ij;
1639
1640 /** The U₂<sup>j</sup> coefficients.
1641 * <p>
1642 * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
1643 * </p>
1644 */
1645 private final double[] u2ij;
1646
1647 /** The V₂<sup>j</sup> coefficients.
1648 * <p>
1649 * Only the coefficients that use the Fourier index = 1 (i == 0) are computed as they are the only ones required.
1650 * </p>
1651 */
1652 private final double[] v2ij;
1653
1654 /** The current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients. */
1655 private final double[][] currentRhoSigmaj;
1656
1657 /** The C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup> Fourier coefficients. */
1658 private final FourierCjSjCoefficients fourierCjSj;
1659
1660 /** The maximum value for j index. */
1661 private final int jMax;
1662
1663 /** Constructor.
1664 * @param currentRhoSigmaj the current computed values for the ρ<sub>j</sub> and σ<sub>j</sub> coefficients
1665 * @param fourierCjSj the fourier coefficients C<sub>i</sub><sup>j</sup> and the S<sub>i</sub><sup>j</sup>
1666 * @param jMax maximum value for j index
1667 */
1668 public UijVijCoefficients(final double[][] currentRhoSigmaj, final FourierCjSjCoefficients fourierCjSj, final int jMax) {
1669 this.currentRhoSigmaj = currentRhoSigmaj;
1670 this.fourierCjSj = fourierCjSj;
1671 this.jMax = jMax;
1672
1673 // initialize the internal arrays.
1674 this.u1ij = new double[6][2 * jMax + 1];
1675 this.v1ij = new double[6][2 * jMax + 1];
1676 this.u2ij = new double[jMax + 1];
1677 this.v2ij = new double[jMax + 1];
1678
1679 //compute the coefficients
1680 computeU1V1Coefficients();
1681 computeU2V2Coefficients();
1682 }
1683
1684 /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients. */
1685 private void computeU1V1Coefficients() {
1686 // generate the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients
1687 // for j >= 1
1688 // also the U₁⁰ for Fourier index = 1 (i == 0) coefficient will be computed
1689 u1ij[0][0] = 0;
1690 for (int j = 1; j <= jMax; j++) {
1691 // compute 1 / j
1692 final double ooj = 1. / j;
1693
1694 for (int i = 0; i < 6; i++) {
1695 //j is aready between 1 and J
1696 u1ij[i][j] = fourierCjSj.getSij(i, j);
1697 v1ij[i][j] = fourierCjSj.getCij(i, j);
1698
1699 // 1 - δ<sub>1j</sub> is 1 for all j > 1
1700 if (j > 1) {
1701 // k starts with 1 because j-J is less than or equal to 0
1702 for (int kIndex = 1; kIndex <= j - 1; kIndex++) {
1703 // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
1704 // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
1705 u1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[1][kIndex] +
1706 fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[0][kIndex];
1707
1708 // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
1709 // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
1710 v1ij[i][j] += fourierCjSj.getCij(i, j - kIndex) * currentRhoSigmaj[0][kIndex] -
1711 fourierCjSj.getSij(i, j - kIndex) * currentRhoSigmaj[1][kIndex];
1712 }
1713 }
1714
1715 // since j must be between 1 and J-1 and is already between 1 and J
1716 // the following sum is skiped only for j = jMax
1717 if (j != jMax) {
1718 for (int kIndex = 1; kIndex <= jMax - j; kIndex++) {
1719 // -C<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub> +
1720 // S<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub>
1721 u1ij[i][j] += -fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[1][kIndex] +
1722 fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[0][kIndex];
1723
1724 // C<sub>i</sub><sup>j+k</sup> * ρ<sub>k</sub> +
1725 // S<sub>i</sub><sup>j+k</sup> * σ<sub>k</sub>
1726 v1ij[i][j] += fourierCjSj.getCij(i, j + kIndex) * currentRhoSigmaj[0][kIndex] +
1727 fourierCjSj.getSij(i, j + kIndex) * currentRhoSigmaj[1][kIndex];
1728 }
1729 }
1730
1731 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
1732 // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
1733 // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
1734 u1ij[i][j] += -fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[1][j + kIndex] -
1735 fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[0][j + kIndex];
1736
1737 // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
1738 // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
1739 v1ij[i][j] += fourierCjSj.getCij(i, kIndex) * currentRhoSigmaj[0][j + kIndex] +
1740 fourierCjSj.getSij(i, kIndex) * currentRhoSigmaj[1][j + kIndex];
1741 }
1742
1743 // divide by 1 / j
1744 u1ij[i][j] *= -ooj;
1745 v1ij[i][j] *= ooj;
1746
1747 // if index = 1 (i == 0) add the computed terms to U₁⁰
1748 if (i == 0) {
1749 //- (U₁<sup>j</sup> * ρ<sub>j</sub> + V₁<sup>j</sup> * σ<sub>j</sub>
1750 u1ij[0][0] += -u1ij[0][j] * currentRhoSigmaj[0][j] - v1ij[0][j] * currentRhoSigmaj[1][j];
1751 }
1752 }
1753 }
1754
1755 // Terms with j > jMax are required only when computing the coefficients
1756 // U₂<sup>j</sup> and V₂<sup>j</sup>
1757 // and those coefficients are only required for Fourier index = 1 (i == 0).
1758 for (int j = jMax + 1; j <= 2 * jMax; j++) {
1759 // compute 1 / j
1760 final double ooj = 1. / j;
1761 //the value of i is 0
1762 u1ij[0][j] = 0.;
1763 v1ij[0][j] = 0.;
1764
1765 //k starts from j-J as it is always greater than or equal to 1
1766 for (int kIndex = j - jMax; kIndex <= j - 1; kIndex++) {
1767 // C<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub> +
1768 // S<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub>
1769 u1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[1][kIndex] +
1770 fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[0][kIndex];
1771
1772 // C<sub>i</sub><sup>j-k</sup> * ρ<sub>k</sub> -
1773 // S<sub>i</sub><sup>j-k</sup> * σ<sub>k</sub>
1774 v1ij[0][j] += fourierCjSj.getCij(0, j - kIndex) * currentRhoSigmaj[0][kIndex] -
1775 fourierCjSj.getSij(0, j - kIndex) * currentRhoSigmaj[1][kIndex];
1776 }
1777 for (int kIndex = 1; kIndex <= jMax; kIndex++) {
1778 // C<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub> -
1779 // S<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub>
1780 u1ij[0][j] += -fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[1][j + kIndex] -
1781 fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[0][j + kIndex];
1782
1783 // C<sub>i</sub><sup>k</sup> * ρ<sub>j+k</sub> +
1784 // S<sub>i</sub><sup>k</sup> * σ<sub>j+k</sub>
1785 v1ij[0][j] += fourierCjSj.getCij(0, kIndex) * currentRhoSigmaj[0][j + kIndex] +
1786 fourierCjSj.getSij(0, kIndex) * currentRhoSigmaj[1][j + kIndex];
1787 }
1788
1789 // divide by 1 / j
1790 u1ij[0][j] *= -ooj;
1791 v1ij[0][j] *= ooj;
1792 }
1793 }
1794
1795 /** Build the U₁<sup>j</sup> and V₁<sup>j</sup> coefficients.
1796 * <p>
1797 * Only the coefficients for Fourier index = 1 (i == 0) are required.
1798 * </p>
1799 */
1800 private void computeU2V2Coefficients() {
1801 for (int j = 1; j <= jMax; j++) {
1802 // compute 1 / j
1803 final double ooj = 1. / j;
1804
1805 // only the values for i == 0 are computed
1806 u2ij[j] = v1ij[0][j];
1807 v2ij[j] = u1ij[0][j];
1808
1809 // 1 - δ<sub>1j</sub> is 1 for all j > 1
1810 if (j > 1) {
1811 for (int l = 1; l <= j - 1; l++) {
1812 // U₁<sup>j-l</sup> * σ<sub>l</sub> +
1813 // V₁<sup>j-l</sup> * ρ<sub>l</sub>
1814 u2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[1][l] +
1815 v1ij[0][j - l] * currentRhoSigmaj[0][l];
1816
1817 // U₁<sup>j-l</sup> * ρ<sub>l</sub> -
1818 // V₁<sup>j-l</sup> * σ<sub>l</sub>
1819 v2ij[j] += u1ij[0][j - l] * currentRhoSigmaj[0][l] -
1820 v1ij[0][j - l] * currentRhoSigmaj[1][l];
1821 }
1822 }
1823
1824 for (int l = 1; l <= jMax; l++) {
1825 // -U₁<sup>j+l</sup> * σ<sub>l</sub> +
1826 // U₁<sup>l</sup> * σ<sub>j+l</sub> +
1827 // V₁<sup>j+l</sup> * ρ<sub>l</sub> -
1828 // V₁<sup>l</sup> * ρ<sub>j+l</sub>
1829 u2ij[j] += -u1ij[0][j + l] * currentRhoSigmaj[1][l] +
1830 u1ij[0][l] * currentRhoSigmaj[1][j + l] +
1831 v1ij[0][j + l] * currentRhoSigmaj[0][l] -
1832 v1ij[0][l] * currentRhoSigmaj[0][j + l];
1833
1834 // U₁<sup>j+l</sup> * ρ<sub>l</sub> +
1835 // U₁<sup>l</sup> * ρ<sub>j+l</sub> +
1836 // V₁<sup>j+l</sup> * σ<sub>l</sub> +
1837 // V₁<sup>l</sup> * σ<sub>j+l</sub>
1838 u2ij[j] += u1ij[0][j + l] * currentRhoSigmaj[0][l] +
1839 u1ij[0][l] * currentRhoSigmaj[0][j + l] +
1840 v1ij[0][j + l] * currentRhoSigmaj[1][l] +
1841 v1ij[0][l] * currentRhoSigmaj[1][j + l];
1842 }
1843
1844 // divide by 1 / j
1845 u2ij[j] *= -ooj;
1846 v2ij[j] *= ooj;
1847 }
1848 }
1849
1850 /** Get the coefficient U₁<sup>j</sup> for Fourier index i.
1851 *
1852 * @param j j index
1853 * @param i Fourier index (starts at 0)
1854 * @return the coefficient U₁<sup>j</sup> for the given Fourier index i
1855 */
1856 public double getU1(final int j, final int i) {
1857 return u1ij[i][j];
1858 }
1859
1860 /** Get the coefficient V₁<sup>j</sup> for Fourier index i.
1861 *
1862 * @param j j index
1863 * @param i Fourier index (starts at 0)
1864 * @return the coefficient V₁<sup>j</sup> for the given Fourier index i
1865 */
1866 public double getV1(final int j, final int i) {
1867 return v1ij[i][j];
1868 }
1869
1870 /** Get the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0).
1871 *
1872 * @param j j index
1873 * @return the coefficient U₂<sup>j</sup> for Fourier index = 1 (i == 0)
1874 */
1875 public double getU2(final int j) {
1876 return u2ij[j];
1877 }
1878
1879 /** Get the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0).
1880 *
1881 * @param j j index
1882 * @return the coefficient V₂<sup>j</sup> for Fourier index = 1 (i == 0)
1883 */
1884 public double getV2(final int j) {
1885 return v2ij[j];
1886 }
1887 }
1888 }