SolidTidesField.java

  1. /* Copyright 2002-2017 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.forces.gravity;

  18. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  19. import org.hipparchus.util.FastMath;
  20. import org.orekit.bodies.CelestialBody;
  21. import org.orekit.errors.OrekitException;
  22. import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
  23. import org.orekit.forces.gravity.potential.TideSystem;
  24. import org.orekit.frames.Frame;
  25. import org.orekit.time.AbsoluteDate;
  26. import org.orekit.time.TimeVectorFunction;
  27. import org.orekit.utils.LoveNumbers;

  28. /** Gravity field corresponding to solid tides.
  29.  * <p>
  30.  * This solid tides force model implementation corresponds to the method described
  31.  * in <a href="http://www.iers.org/nn_11216/IERS/EN/Publications/TechnicalNotes/tn36.html">
  32.  * IERS conventions (2010)</a>, chapter 6, section 6.2.
  33.  * </p>
  34.  * <p>
  35.  * The computation of the spherical harmonics part is done using the algorithm
  36.  * designed by S. A. Holmes and W. E. Featherstone from Department of Spatial Sciences,
  37.  * Curtin University of Technology, Perth, Australia and described in their 2002 paper:
  38.  * <a href="http://cct.gfy.ku.dk/publ_others/ccta1870.pdf">A unified approach to
  39.  * the Clenshaw summation and the recursive computation of very high degree and
  40.  * order normalised associated Legendre functions</a> (Journal of Geodesy (2002)
  41.  * 76: 279–299).
  42.  * </p>
  43.  * <p>
  44.  * Note that this class is <em>not</em> thread-safe, and that tides computation
  45.  * are computer intensive if repeated. So this class is really expected to
  46.  * be wrapped within a {@link
  47.  * org.orekit.forces.gravity.potential.CachedNormalizedSphericalHarmonicsProvider}
  48.  * that will both ensure thread safety and improve performances using caching.
  49.  * </p>
  50.  * @see SolidTides
  51.  * @author Luc Maisonobe
  52.  * @since 6.1
  53.  */
  54. class SolidTidesField implements NormalizedSphericalHarmonicsProvider {

  55.     /** Maximum degree for normalized Legendre functions. */
  56.     private static final int MAX_LEGENDRE_DEGREE = 4;

  57.     /** Love numbers. */
  58.     private final LoveNumbers love;

  59.     /** Function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂). */
  60.     private final TimeVectorFunction deltaCSFunction;

  61.     /** Permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used. */
  62.     private final double deltaC20PermanentTide;

  63.     /** Function computing pole tide terms (ΔC₂₁, ΔS₂₁). */
  64.     private final TimeVectorFunction poleTideFunction;

  65.     /** Rotating body frame. */
  66.     private final Frame centralBodyFrame;

  67.     /** Central body reference radius. */
  68.     private final double ae;

  69.     /** Central body attraction coefficient. */
  70.     private final double mu;

  71.     /** Tide system used in the central attraction model. */
  72.     private final TideSystem centralTideSystem;

  73.     /** Tide generating bodies (typically Sun and Moon). Read only after construction. */
  74.     private final CelestialBody[] bodies;

  75.     /** First recursion coefficients for tesseral terms. Read only after construction. */
  76.     private final double[][] anm;

  77.     /** Second recursion coefficients for tesseral terms. Read only after construction. */
  78.     private final double[][] bnm;

  79.     /** Third recursion coefficients for sectorial terms. Read only after construction. */
  80.     private final double[] dmm;

  81.     /** Simple constructor.
  82.      * @param love Love numbers
  83.      * @param deltaCSFunction function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂)
  84.      * @param deltaC20PermanentTide permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used
  85.      * @param poleTideFunction function computing pole tide terms (ΔC₂₁, ΔS₂₁), may be null
  86.      * @param centralBodyFrame rotating body frame
  87.      * @param ae central body reference radius
  88.      * @param mu central body attraction coefficient
  89.      * @param centralTideSystem tide system used in the central attraction model
  90.      * @param bodies tide generating bodies (typically Sun and Moon)
  91.      */
  92.     SolidTidesField(final LoveNumbers love, final TimeVectorFunction deltaCSFunction,
  93.                            final double deltaC20PermanentTide, final TimeVectorFunction poleTideFunction,
  94.                            final Frame centralBodyFrame, final double ae, final double mu,
  95.                            final TideSystem centralTideSystem, final CelestialBody... bodies) {

  96.         // store mode parameters
  97.         this.centralBodyFrame  = centralBodyFrame;
  98.         this.ae                = ae;
  99.         this.mu                = mu;
  100.         this.centralTideSystem = centralTideSystem;
  101.         this.bodies            = bodies;

  102.         // compute recursion coefficients for Legendre functions
  103.         this.anm               = buildTriangularArray(5, false);
  104.         this.bnm               = buildTriangularArray(5, false);
  105.         this.dmm               = new double[love.getSize()];
  106.         recursionCoefficients();

  107.         // Love numbers
  108.         this.love = love;

  109.         // frequency dependent terms
  110.         this.deltaCSFunction = deltaCSFunction;

  111.         // permanent tide
  112.         this.deltaC20PermanentTide = deltaC20PermanentTide;

  113.         // pole tide
  114.         this.poleTideFunction = poleTideFunction;

  115.     }

  116.     /** {@inheritDoc} */
  117.     @Override
  118.     public int getMaxDegree() {
  119.         return MAX_LEGENDRE_DEGREE;
  120.     }

  121.     /** {@inheritDoc} */
  122.     @Override
  123.     public int getMaxOrder() {
  124.         return MAX_LEGENDRE_DEGREE;
  125.     }

  126.     /** {@inheritDoc} */
  127.     @Override
  128.     public double getMu() {
  129.         return mu;
  130.     }

  131.     /** {@inheritDoc} */
  132.     @Override
  133.     public double getAe() {
  134.         return ae;
  135.     }

  136.     /** {@inheritDoc} */
  137.     @Override
  138.     public AbsoluteDate getReferenceDate() {
  139.         return AbsoluteDate.J2000_EPOCH;
  140.     }

  141.     /** {@inheritDoc} */
  142.     @Override
  143.     public double getOffset(final AbsoluteDate date) {
  144.         return date.durationFrom(AbsoluteDate.J2000_EPOCH);
  145.     }

  146.     /** {@inheritDoc} */
  147.     @Override
  148.     public TideSystem getTideSystem() {
  149.         // not really used here, but for consistency we can state that either
  150.         // we add the permanent tide or it was already in the central attraction
  151.         return TideSystem.ZERO_TIDE;
  152.     }

  153.     /** {@inheritDoc} */
  154.     @Override
  155.     public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) throws OrekitException {

  156.         // computed Cnm and Snm coefficients
  157.         final double[][] cnm = buildTriangularArray(5, true);
  158.         final double[][] snm = buildTriangularArray(5, true);

  159.         // work array to hold Legendre coefficients
  160.         final double[][] pnm = buildTriangularArray(5, true);

  161.         // step 1: frequency independent part
  162.         // equations 6.6 (for degrees 2 and 3) and 6.7 (for degree 4) in IERS conventions 2010
  163.         for (final CelestialBody body : bodies) {

  164.             // compute tide generating body state
  165.             final Vector3D position = body.getPVCoordinates(date, centralBodyFrame).getPosition();

  166.             // compute polar coordinates
  167.             final double x    = position.getX();
  168.             final double y    = position.getY();
  169.             final double z    = position.getZ();
  170.             final double x2   = x * x;
  171.             final double y2   = y * y;
  172.             final double z2   = z * z;
  173.             final double r2   = x2 + y2 + z2;
  174.             final double r    = FastMath.sqrt (r2);
  175.             final double rho2 = x2 + y2;
  176.             final double rho  = FastMath.sqrt(rho2);

  177.             // evaluate Pnm
  178.             evaluateLegendre(z / r, rho / r, pnm);

  179.             // update spherical harmonic coefficients
  180.             frequencyIndependentPart(r, body.getGM(), x / rho, y / rho, pnm, cnm, snm);

  181.         }

  182.         // step 2: frequency dependent corrections
  183.         frequencyDependentPart(date, cnm, snm);

  184.         if (centralTideSystem == TideSystem.ZERO_TIDE) {
  185.             // step 3: remove permanent tide which is already considered
  186.             // in the central body gravity field
  187.             removePermanentTide(cnm);
  188.         }

  189.         if (poleTideFunction != null) {
  190.             // add pole tide
  191.             poleTide(date, cnm, snm);
  192.         }

  193.         return new TideHarmonics(date, cnm, snm);

  194.     }

  195.     /** Compute recursion coefficients. */
  196.     private void recursionCoefficients() {

  197.         // pre-compute the recursion coefficients
  198.         // (see equations 11 and 12 from Holmes and Featherstone paper)
  199.         for (int n = 0; n < anm.length; ++n) {
  200.             for (int m = 0; m < n; ++m) {
  201.                 final double f = (n - m) * (n + m );
  202.                 anm[n][m] = FastMath.sqrt((2 * n - 1) * (2 * n + 1) / f);
  203.                 bnm[n][m] = FastMath.sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) / (f * (2 * n - 3)));
  204.             }
  205.         }

  206.         // sectorial terms corresponding to equation 13 in Holmes and Featherstone paper
  207.         dmm[0] = Double.NaN; // dummy initialization for unused term
  208.         dmm[1] = Double.NaN; // dummy initialization for unused term
  209.         for (int m = 2; m < dmm.length; ++m) {
  210.             dmm[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m));
  211.         }

  212.     }

  213.     /** Evaluate Legendre functions.
  214.      * @param t cos(θ), where θ is the polar angle
  215.      * @param u sin(θ), where θ is the polar angle
  216.      * @param pnm the computed coefficients. Modified in place.
  217.      */
  218.     private void evaluateLegendre(final double t, final double u, final double[][] pnm) {

  219.         // as the degree is very low, we use the standard forward column method
  220.         // and store everything (see equations 11 and 13 from Holmes and Featherstone paper)
  221.         pnm[0][0] = 1;
  222.         pnm[1][0] = anm[1][0] * t;
  223.         pnm[1][1] = FastMath.sqrt(3) * u;
  224.         for (int m = 2; m < dmm.length; ++m) {
  225.             pnm[m][m - 1] = anm[m][m - 1] * t * pnm[m - 1][m - 1];
  226.             pnm[m][m]     = dmm[m] * u * pnm[m - 1][m - 1];
  227.         }
  228.         for (int m = 0; m < dmm.length; ++m) {
  229.             for (int n = m + 2; n < pnm.length; ++n) {
  230.                 pnm[n][m] = anm[n][m] * t * pnm[n - 1][m] - bnm[n][m] * pnm[n - 2][m];
  231.             }
  232.         }

  233.     }

  234.     /** Update coefficients applying frequency independent step, for one tide generating body.
  235.      * @param r distance to tide generating body
  236.      * @param gm tide generating body attraction coefficient
  237.      * @param cosLambda cosine of the tide generating body longitude
  238.      * @param sinLambda sine of the tide generating body longitude
  239.      * @param pnm the Legendre coefficients. See {@link #evaluateLegendre(double, double, double[][])}.
  240.      * @param cnm the computed Cnm coefficients. Modified in place.
  241.      * @param snm the computed Snm coefficients. Modified in place.
  242.      */
  243.     private void frequencyIndependentPart(final double r,
  244.                                           final double gm,
  245.                                           final double cosLambda,
  246.                                           final double sinLambda,
  247.                                           final double[][] pnm,
  248.                                           final double[][] cnm,
  249.                                           final double[][] snm) {

  250.         final double rRatio = ae / r;
  251.         double fM           = gm / mu;
  252.         double cosMLambda   = 1;
  253.         double sinMLambda   = 0;
  254.         for (int m = 0; m < love.getSize(); ++m) {

  255.             double fNPlus1 = fM;
  256.             for (int n = m; n < love.getSize(); ++n) {
  257.                 fNPlus1 *= rRatio;
  258.                 final double coeff = (fNPlus1 / (2 * n + 1)) * pnm[n][m];
  259.                 final double cCos  = coeff * cosMLambda;
  260.                 final double cSin  = coeff * sinMLambda;

  261.                 // direct effect of degree n tides on degree n coefficients
  262.                 // equation 6.6 from IERS conventions (2010)
  263.                 final double kR = love.getReal(n, m);
  264.                 final double kI = love.getImaginary(n, m);
  265.                 cnm[n][m] += kR * cCos + kI * cSin;
  266.                 snm[n][m] += kR * cSin - kI * cCos;

  267.                 if (n == 2) {
  268.                     // indirect effect of degree 2 tides on degree 4 coefficients
  269.                     // equation 6.7 from IERS conventions (2010)
  270.                     final double kP = love.getPlus(n, m);
  271.                     cnm[4][m] += kP * cCos;
  272.                     snm[4][m] += kP * cSin;
  273.                 }

  274.             }

  275.             // prepare next iteration on order
  276.             final double tmp = cosMLambda * cosLambda - sinMLambda * sinLambda;
  277.             sinMLambda = sinMLambda * cosLambda + cosMLambda * sinLambda;
  278.             cosMLambda = tmp;
  279.             fM        *= rRatio;

  280.         }

  281.     }

  282.     /** Update coefficients applying frequency dependent step.
  283.      * @param date current date
  284.      * @param cnm the Cnm coefficients. Modified in place.
  285.      * @param snm the Snm coefficients. Modified in place.
  286.      */
  287.     private void frequencyDependentPart(final AbsoluteDate date,
  288.                                         final double[][] cnm,
  289.                                         final double[][] snm) {
  290.         final double[] deltaCS = deltaCSFunction.value(date);
  291.         cnm[2][0] += deltaCS[0]; // ΔC₂₀
  292.         cnm[2][1] += deltaCS[1]; // ΔC₂₁
  293.         snm[2][1] += deltaCS[2]; // ΔS₂₁
  294.         cnm[2][2] += deltaCS[3]; // ΔC₂₂
  295.         snm[2][2] += deltaCS[4]; // ΔS₂₂
  296.     }

  297.     /** Remove the permanent tide already counted in zero-tide central gravity fields.
  298.      * @param cnm the Cnm coefficients. Modified in place.
  299.      */
  300.     private void removePermanentTide(final double[][] cnm) {
  301.         cnm[2][0] -= deltaC20PermanentTide;
  302.     }

  303.     /** Update coefficients applying pole tide.
  304.      * @param date current date
  305.      * @param cnm the Cnm coefficients. Modified in place.
  306.      * @param snm the Snm coefficients. Modified in place.
  307.      */
  308.     private void poleTide(final AbsoluteDate date,
  309.                           final double[][] cnm,
  310.                           final double[][] snm) {
  311.         final double[] deltaCS = poleTideFunction.value(date);
  312.         cnm[2][1] += deltaCS[0]; // ΔC₂₁
  313.         snm[2][1] += deltaCS[1]; // ΔS₂₁
  314.     }

  315.     /** Create a triangular array.
  316.      * @param size array size
  317.      * @param withDiagonal if true, the array contains a[p][p] terms, otherwise each
  318.      * row ends up at a[p][p-1]
  319.      * @return new triangular array
  320.      */
  321.     private double[][] buildTriangularArray(final int size, final boolean withDiagonal) {
  322.         final double[][] array = new double[size][];
  323.         for (int i = 0; i < array.length; ++i) {
  324.             array[i] = new double[withDiagonal ? i + 1 : i];
  325.         }
  326.         return array;
  327.     }

  328.     /** The Tidal geopotential evaluated on a specific date. */
  329.     private static class TideHarmonics implements NormalizedSphericalHarmonics {

  330.         /** evaluation date. */
  331.         private final AbsoluteDate date;

  332.         /** Cached cnm. */
  333.         private final double[][] cnm;

  334.         /** Cached snm. */
  335.         private final double[][] snm;

  336.         /** Construct the tidal harmonics on the given date.
  337.          *
  338.          * @param date of evaluation
  339.          * @param cnm the Cnm coefficients. Not copied.
  340.          * @param snm the Snm coeffiecients. Not copied.
  341.          */
  342.         private TideHarmonics(final AbsoluteDate date,
  343.                               final double[][] cnm,
  344.                               final double[][] snm) {
  345.             this.date = date;
  346.             this.cnm = cnm;
  347.             this.snm = snm;
  348.         }

  349.         /** {@inheritDoc} */
  350.         @Override
  351.         public AbsoluteDate getDate() {
  352.             return date;
  353.         }

  354.         /** {@inheritDoc} */
  355.         @Override
  356.         public double getNormalizedCnm(final int n, final int m)
  357.             throws OrekitException {
  358.             return cnm[n][m];
  359.         }

  360.         /** {@inheritDoc} */
  361.         @Override
  362.         public double getNormalizedSnm(final int n, final int m)
  363.             throws OrekitException {
  364.             return snm[n][m];
  365.         }

  366.     }

  367. }