org.orekit.data.PoissonSeries Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of orekit Show documentation
Show all versions of orekit Show documentation
OREKIT (ORbits Extrapolation KIT) is a low level space dynamics library.
It provides basic elements (orbits, dates, attitude, frames ...) and
various algorithms to handle them (conversions, analytical and numerical
propagation, pointing ...).
/* Copyright 2002-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.data;
import java.util.HashMap;
import java.util.Map;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.MathUtils.SumAndResidual;
/**
* Class representing a Poisson series for nutation or ephemeris computations.
*
* A Poisson series is composed of a time polynomial part and a non-polynomial
* part which consist in summation series. The {@link SeriesTerm series terms}
* are harmonic functions (combination of sines and cosines) of polynomial
* arguments. The polynomial arguments are combinations of luni-solar or
* planetary {@link BodiesElements elements}.
*
* @author Luc Maisonobe
* @see PoissonSeriesParser
* @see SeriesTerm
* @see PolynomialNutation
*/
public class PoissonSeries {
/** Polynomial part. */
private final PolynomialNutation polynomial;
/** Non-polynomial series. */
private final Map series;
/** Build a Poisson series from an IERS table file.
* @param polynomial polynomial part (may be null)
* @param series non-polynomial part
*/
public PoissonSeries(final PolynomialNutation polynomial, final Map series) {
this.polynomial = polynomial;
this.series = series;
}
/** Get the polynomial part of the series.
* @return polynomial part of the series.
*/
public PolynomialNutation getPolynomial() {
return polynomial;
}
/** Get the number of different terms in the non-polynomial part.
* @return number of different terms in the non-polynomial part
*/
public int getNonPolynomialSize() {
return series.size();
}
/** Evaluate the value of the series.
* @param elements bodies elements for nutation
* @return value of the series
*/
public double value(final BodiesElements elements) {
// polynomial part
final double p = polynomial.value(elements.getTC());
// non-polynomial part
double npHigh = 0;
double npLow = 0;
for (final Map.Entry entry : series.entrySet()) {
final double v = entry.getValue().value(elements)[0];
// Use 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh, v);
npHigh = sumAndResidual.getSum();
npLow += sumAndResidual.getResidual();
}
// add the polynomial and the non-polynomial parts
return p + (npHigh + npLow);
}
/** Evaluate the value of the series.
* @param elements bodies elements for nutation
* @param type fo the field elements
* @return value of the series
*/
public > T value(final FieldBodiesElements elements) {
// polynomial part
final T tc = elements.getTC();
final T p = polynomial.value(tc);
// non-polynomial part
T sum = tc.getField().getZero();
for (final Map.Entry entry : series.entrySet()) {
sum = sum.add(entry.getValue().value(elements)[0]);
}
// add the polynomial and the non-polynomial parts
return p.add(sum);
}
/** This interface represents a fast evaluator for Poisson series.
* @see PoissonSeries#compile(PoissonSeries...)
* @since 6.1
*/
public interface CompiledSeries {
/** Evaluate a set of Poisson series.
* @param elements bodies elements for nutation
* @return value of the series
*/
double[] value(BodiesElements elements);
/** Evaluate time derivative of a set of Poisson series.
* @param elements bodies elements for nutation
* @return time derivative of the series
*/
double[] derivative(BodiesElements elements);
/** Evaluate a set of Poisson series.
* @param elements bodies elements for nutation
* @param the type of the field elements
* @return value of the series
*/
> S[] value(FieldBodiesElements elements);
/** Evaluate time derivative of a set of Poisson series.
* @param elements bodies elements for nutation
* @param the type of the field elements
* @return time derivative of the series
*/
> S[] derivative(FieldBodiesElements elements);
}
/** Join several nutation series, for fast simultaneous evaluation.
* @param poissonSeries Poisson series to join
* @return a single function that evaluates all series together
* @since 6.1
*/
@SafeVarargs
public static CompiledSeries compile(final PoissonSeries... poissonSeries) {
// store all polynomials
final PolynomialNutation[] polynomials = new PolynomialNutation[poissonSeries.length];
for (int i = 0; i < polynomials.length; ++i) {
polynomials[i] = poissonSeries[i].polynomial;
}
// gather all series terms
final Map joinedMap = new HashMap();
for (final PoissonSeries ps : poissonSeries) {
for (Map.Entry entry : ps.series.entrySet()) {
final long key = entry.getKey();
if (!joinedMap.containsKey(key)) {
// retrieve all Delaunay and planetary multipliers from the key
final int[] m = NutationCodec.decode(key);
// prepare a new term, ready to handle the required dimension
final SeriesTerm term =
SeriesTerm.buildTerm(m[0],
m[1], m[2], m[3], m[4], m[5],
m[6], m[7], m[8], m[9], m[10], m[11], m[12], m[13], m[14]);
term.add(poissonSeries.length - 1, -1, Double.NaN, Double.NaN);
// store it
joinedMap.put(key, term);
}
}
}
// join series by sharing terms, in order to speed up evaluation
// which is dominated by the computation of sine/cosine in each term
for (int i = 0; i < poissonSeries.length; ++i) {
for (final Map.Entry entry : poissonSeries[i].series.entrySet()) {
final SeriesTerm singleTerm = entry.getValue();
final SeriesTerm joinedTerm = joinedMap.get(entry.getKey());
for (int degree = 0; degree <= singleTerm.getDegree(0); ++degree) {
joinedTerm.add(i, degree,
singleTerm.getSinCoeff(0, degree),
singleTerm.getCosCoeff(0, degree));
}
}
}
// use a single array for faster access
final SeriesTerm[] joinedTerms = new SeriesTerm[joinedMap.size()];
int index = 0;
for (final Map.Entry entry : joinedMap.entrySet()) {
joinedTerms[index++] = entry.getValue();
}
return new CompiledSeries() {
/** {@inheritDoc} */
@Override
public double[] value(final BodiesElements elements) {
// non-polynomial part
final double[] npHigh = new double[polynomials.length];
final double[] npLow = new double[polynomials.length];
for (final SeriesTerm term : joinedTerms) {
final double[] termValue = term.value(elements);
for (int i = 0; i < termValue.length; ++i) {
// Use 2Sum for high precision.
final SumAndResidual sumAndResidual = MathUtils.twoSum(npHigh[i], termValue[i]);
npHigh[i] = sumAndResidual.getSum();
npLow[i] += sumAndResidual.getResidual();
}
}
// add residual and polynomial part
for (int i = 0; i < npHigh.length; ++i) {
npHigh[i] += npLow[i] + polynomials[i].value(elements.getTC());
}
return npHigh;
}
/** {@inheritDoc} */
@Override
public double[] derivative(final BodiesElements elements) {
// non-polynomial part
final double[] v = new double[polynomials.length];
for (final SeriesTerm term : joinedTerms) {
final double[] termDerivative = term.derivative(elements);
for (int i = 0; i < termDerivative.length; ++i) {
v[i] += termDerivative[i];
}
}
// add polynomial part
for (int i = 0; i < v.length; ++i) {
v[i] += polynomials[i].derivative(elements.getTC());
}
return v;
}
/** {@inheritDoc} */
@Override
public > S[] value(final FieldBodiesElements elements) {
// non-polynomial part
final S[] v = MathArrays.buildArray(elements.getTC().getField(), polynomials.length);
for (final SeriesTerm term : joinedTerms) {
final S[] termValue = term.value(elements);
for (int i = 0; i < termValue.length; ++i) {
v[i] = v[i].add(termValue[i]);
}
}
// add polynomial part
final S tc = elements.getTC();
for (int i = 0; i < v.length; ++i) {
v[i] = v[i].add(polynomials[i].value(tc));
}
return v;
}
/** {@inheritDoc} */
@Override
public > S[] derivative(final FieldBodiesElements elements) {
// non-polynomial part
final S[] v = MathArrays.buildArray(elements.getTC().getField(), polynomials.length);
for (final SeriesTerm term : joinedTerms) {
final S[] termDerivative = term.derivative(elements);
for (int i = 0; i < termDerivative.length; ++i) {
v[i] = v[i].add(termDerivative[i]);
}
}
// add polynomial part
final S tc = elements.getTC();
for (int i = 0; i < v.length; ++i) {
v[i] = v[i].add(polynomials[i].derivative(tc));
}
return v;
}
};
}
}