/* Copyright 2002-2017 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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.forces.gravity;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.CelestialBody;
import org.orekit.errors.OrekitException;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.TideSystem;
import org.orekit.frames.Frame;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeVectorFunction;
import org.orekit.utils.LoveNumbers;
/** Gravity field corresponding to solid tides.
* <p>
* This solid tides force model implementation corresponds to the method described
* in <a href="http://www.iers.org/nn_11216/IERS/EN/Publications/TechnicalNotes/tn36.html">
* IERS conventions (2010)</a>, chapter 6, section 6.2.
* </p>
* <p>
* The computation of the spherical harmonics part is done using the algorithm
* designed by S. A. Holmes and W. E. Featherstone from Department of Spatial Sciences,
* Curtin University of Technology, Perth, Australia and described in their 2002 paper:
* <a href="http://cct.gfy.ku.dk/publ_others/ccta1870.pdf">A unified approach to
* the Clenshaw summation and the recursive computation of very high degree and
* order normalised associated Legendre functions</a> (Journal of Geodesy (2002)
* 76: 279–299).
* </p>
* <p>
* Note that this class is <em>not</em> thread-safe, and that tides computation
* are computer intensive if repeated. So this class is really expected to
* be wrapped within a {@link
* org.orekit.forces.gravity.potential.CachedNormalizedSphericalHarmonicsProvider}
* that will both ensure thread safety and improve performances using caching.
* </p>
* @see SolidTides
* @author Luc Maisonobe
* @since 6.1
*/
class SolidTidesField implements NormalizedSphericalHarmonicsProvider {
/** Maximum degree for normalized Legendre functions. */
private static final int MAX_LEGENDRE_DEGREE = 4;
/** Love numbers. */
private final LoveNumbers love;
/** Function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂). */
private final TimeVectorFunction deltaCSFunction;
/** Permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used. */
private final double deltaC20PermanentTide;
/** Function computing pole tide terms (ΔC₂₁, ΔS₂₁). */
private final TimeVectorFunction poleTideFunction;
/** Rotating body frame. */
private final Frame centralBodyFrame;
/** Central body reference radius. */
private final double ae;
/** Central body attraction coefficient. */
private final double mu;
/** Tide system used in the central attraction model. */
private final TideSystem centralTideSystem;
/** Tide generating bodies (typically Sun and Moon). Read only after construction. */
private final CelestialBody[] bodies;
/** First recursion coefficients for tesseral terms. Read only after construction. */
private final double[][] anm;
/** Second recursion coefficients for tesseral terms. Read only after construction. */
private final double[][] bnm;
/** Third recursion coefficients for sectorial terms. Read only after construction. */
private final double[] dmm;
/** Simple constructor.
* @param love Love numbers
* @param deltaCSFunction function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂)
* @param deltaC20PermanentTide permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used
* @param poleTideFunction function computing pole tide terms (ΔC₂₁, ΔS₂₁), may be null
* @param centralBodyFrame rotating body frame
* @param ae central body reference radius
* @param mu central body attraction coefficient
* @param centralTideSystem tide system used in the central attraction model
* @param bodies tide generating bodies (typically Sun and Moon)
*/
SolidTidesField(final LoveNumbers love, final TimeVectorFunction deltaCSFunction,
final double deltaC20PermanentTide, final TimeVectorFunction poleTideFunction,
final Frame centralBodyFrame, final double ae, final double mu,
final TideSystem centralTideSystem, final CelestialBody... bodies) {
// store mode parameters
this.centralBodyFrame = centralBodyFrame;
this.ae = ae;
this.mu = mu;
this.centralTideSystem = centralTideSystem;
this.bodies = bodies;
// compute recursion coefficients for Legendre functions
this.anm = buildTriangularArray(5, false);
this.bnm = buildTriangularArray(5, false);
this.dmm = new double[love.getSize()];
recursionCoefficients();
// Love numbers
this.love = love;
// frequency dependent terms
this.deltaCSFunction = deltaCSFunction;
// permanent tide
this.deltaC20PermanentTide = deltaC20PermanentTide;
// pole tide
this.poleTideFunction = poleTideFunction;
}
/** {@inheritDoc} */
@Override
public int getMaxDegree() {
return MAX_LEGENDRE_DEGREE;
}
/** {@inheritDoc} */
@Override
public int getMaxOrder() {
return MAX_LEGENDRE_DEGREE;
}
/** {@inheritDoc} */
@Override
public double getMu() {
return mu;
}
/** {@inheritDoc} */
@Override
public double getAe() {
return ae;
}
/** {@inheritDoc} */
@Override
public AbsoluteDate getReferenceDate() {
return AbsoluteDate.J2000_EPOCH;
}
/** {@inheritDoc} */
@Override
public double getOffset(final AbsoluteDate date) {
return date.durationFrom(AbsoluteDate.J2000_EPOCH);
}
/** {@inheritDoc} */
@Override
public TideSystem getTideSystem() {
// not really used here, but for consistency we can state that either
// we add the permanent tide or it was already in the central attraction
return TideSystem.ZERO_TIDE;
}
/** {@inheritDoc} */
@Override
public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) throws OrekitException {
// computed Cnm and Snm coefficients
final double[][] cnm = buildTriangularArray(5, true);
final double[][] snm = buildTriangularArray(5, true);
// work array to hold Legendre coefficients
final double[][] pnm = buildTriangularArray(5, true);
// step 1: frequency independent part
// equations 6.6 (for degrees 2 and 3) and 6.7 (for degree 4) in IERS conventions 2010
for (final CelestialBody body : bodies) {
// compute tide generating body state
final Vector3D position = body.getPVCoordinates(date, centralBodyFrame).getPosition();
// compute polar coordinates
final double x = position.getX();
final double y = position.getY();
final double z = position.getZ();
final double x2 = x * x;
final double y2 = y * y;
final double z2 = z * z;
final double r2 = x2 + y2 + z2;
final double r = FastMath.sqrt (r2);
final double rho2 = x2 + y2;
final double rho = FastMath.sqrt(rho2);
// evaluate Pnm
evaluateLegendre(z / r, rho / r, pnm);
// update spherical harmonic coefficients
frequencyIndependentPart(r, body.getGM(), x / rho, y / rho, pnm, cnm, snm);
}
// step 2: frequency dependent corrections
frequencyDependentPart(date, cnm, snm);
if (centralTideSystem == TideSystem.ZERO_TIDE) {
// step 3: remove permanent tide which is already considered
// in the central body gravity field
removePermanentTide(cnm);
}
if (poleTideFunction != null) {
// add pole tide
poleTide(date, cnm, snm);
}
return new TideHarmonics(date, cnm, snm);
}
/** Compute recursion coefficients. */
private void recursionCoefficients() {
// pre-compute the recursion coefficients
// (see equations 11 and 12 from Holmes and Featherstone paper)
for (int n = 0; n < anm.length; ++n) {
for (int m = 0; m < n; ++m) {
final double f = (n - m) * (n + m );
anm[n][m] = FastMath.sqrt((2 * n - 1) * (2 * n + 1) / f);
bnm[n][m] = FastMath.sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) / (f * (2 * n - 3)));
}
}
// sectorial terms corresponding to equation 13 in Holmes and Featherstone paper
dmm[0] = Double.NaN; // dummy initialization for unused term
dmm[1] = Double.NaN; // dummy initialization for unused term
for (int m = 2; m < dmm.length; ++m) {
dmm[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m));
}
}
/** Evaluate Legendre functions.
* @param t cos(θ), where θ is the polar angle
* @param u sin(θ), where θ is the polar angle
* @param pnm the computed coefficients. Modified in place.
*/
private void evaluateLegendre(final double t, final double u, final double[][] pnm) {
// as the degree is very low, we use the standard forward column method
// and store everything (see equations 11 and 13 from Holmes and Featherstone paper)
pnm[0][0] = 1;
pnm[1][0] = anm[1][0] * t;
pnm[1][1] = FastMath.sqrt(3) * u;
for (int m = 2; m < dmm.length; ++m) {
pnm[m][m - 1] = anm[m][m - 1] * t * pnm[m - 1][m - 1];
pnm[m][m] = dmm[m] * u * pnm[m - 1][m - 1];
}
for (int m = 0; m < dmm.length; ++m) {
for (int n = m + 2; n < pnm.length; ++n) {
pnm[n][m] = anm[n][m] * t * pnm[n - 1][m] - bnm[n][m] * pnm[n - 2][m];
}
}
}
/** Update coefficients applying frequency independent step, for one tide generating body.
* @param r distance to tide generating body
* @param gm tide generating body attraction coefficient
* @param cosLambda cosine of the tide generating body longitude
* @param sinLambda sine of the tide generating body longitude
* @param pnm the Legendre coefficients. See {@link #evaluateLegendre(double, double, double[][])}.
* @param cnm the computed Cnm coefficients. Modified in place.
* @param snm the computed Snm coefficients. Modified in place.
*/
private void frequencyIndependentPart(final double r,
final double gm,
final double cosLambda,
final double sinLambda,
final double[][] pnm,
final double[][] cnm,
final double[][] snm) {
final double rRatio = ae / r;
double fM = gm / mu;
double cosMLambda = 1;
double sinMLambda = 0;
for (int m = 0; m < love.getSize(); ++m) {
double fNPlus1 = fM;
for (int n = m; n < love.getSize(); ++n) {
fNPlus1 *= rRatio;
final double coeff = (fNPlus1 / (2 * n + 1)) * pnm[n][m];
final double cCos = coeff * cosMLambda;
final double cSin = coeff * sinMLambda;
// direct effect of degree n tides on degree n coefficients
// equation 6.6 from IERS conventions (2010)
final double kR = love.getReal(n, m);
final double kI = love.getImaginary(n, m);
cnm[n][m] += kR * cCos + kI * cSin;
snm[n][m] += kR * cSin - kI * cCos;
if (n == 2) {
// indirect effect of degree 2 tides on degree 4 coefficients
// equation 6.7 from IERS conventions (2010)
final double kP = love.getPlus(n, m);
cnm[4][m] += kP * cCos;
snm[4][m] += kP * cSin;
}
}
// prepare next iteration on order
final double tmp = cosMLambda * cosLambda - sinMLambda * sinLambda;
sinMLambda = sinMLambda * cosLambda + cosMLambda * sinLambda;
cosMLambda = tmp;
fM *= rRatio;
}
}
/** Update coefficients applying frequency dependent step.
* @param date current date
* @param cnm the Cnm coefficients. Modified in place.
* @param snm the Snm coefficients. Modified in place.
*/
private void frequencyDependentPart(final AbsoluteDate date,
final double[][] cnm,
final double[][] snm) {
final double[] deltaCS = deltaCSFunction.value(date);
cnm[2][0] += deltaCS[0]; // ΔC₂₀
cnm[2][1] += deltaCS[1]; // ΔC₂₁
snm[2][1] += deltaCS[2]; // ΔS₂₁
cnm[2][2] += deltaCS[3]; // ΔC₂₂
snm[2][2] += deltaCS[4]; // ΔS₂₂
}
/** Remove the permanent tide already counted in zero-tide central gravity fields.
* @param cnm the Cnm coefficients. Modified in place.
*/
private void removePermanentTide(final double[][] cnm) {
cnm[2][0] -= deltaC20PermanentTide;
}
/** Update coefficients applying pole tide.
* @param date current date
* @param cnm the Cnm coefficients. Modified in place.
* @param snm the Snm coefficients. Modified in place.
*/
private void poleTide(final AbsoluteDate date,
final double[][] cnm,
final double[][] snm) {
final double[] deltaCS = poleTideFunction.value(date);
cnm[2][1] += deltaCS[0]; // ΔC₂₁
snm[2][1] += deltaCS[1]; // ΔS₂₁
}
/** Create a triangular array.
* @param size array size
* @param withDiagonal if true, the array contains a[p][p] terms, otherwise each
* row ends up at a[p][p-1]
* @return new triangular array
*/
private double[][] buildTriangularArray(final int size, final boolean withDiagonal) {
final double[][] array = new double[size][];
for (int i = 0; i < array.length; ++i) {
array[i] = new double[withDiagonal ? i + 1 : i];
}
return array;
}
/** The Tidal geopotential evaluated on a specific date. */
private static class TideHarmonics implements NormalizedSphericalHarmonics {
/** evaluation date. */
private final AbsoluteDate date;
/** Cached cnm. */
private final double[][] cnm;
/** Cached snm. */
private final double[][] snm;
/** Construct the tidal harmonics on the given date.
*
* @param date of evaluation
* @param cnm the Cnm coefficients. Not copied.
* @param snm the Snm coeffiecients. Not copied.
*/
private TideHarmonics(final AbsoluteDate date,
final double[][] cnm,
final double[][] snm) {
this.date = date;
this.cnm = cnm;
this.snm = snm;
}
/** {@inheritDoc} */
@Override
public AbsoluteDate getDate() {
return date;
}
/** {@inheritDoc} */
@Override
public double getNormalizedCnm(final int n, final int m)
throws OrekitException {
return cnm[n][m];
}
/** {@inheritDoc} */
@Override
public double getNormalizedSnm(final int n, final int m)
throws OrekitException {
return snm[n][m];
}
}
}