SolidTidesField.java

  1. /* Copyright 2002-2013 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;

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

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

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

  61.     /** Love numbers. */
  62.     private final LoveNumbers love;

  64.     /** Function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂). */
  65.     private final TimeFunction<double []> deltaCSFunction;

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

  70.     /** Function computing pole tide terms (ΔC₂₁, ΔS₂₁). */
  71.     private final TimeFunction<double []> poleTideFunction;

  73.     /** Rotating body frame. */
  74.     private final Frame centralBodyFrame;

  76.     /** Central body reference radius. */
  77.     private final double ae;

  79.     /** Central body attraction coefficient. */
  80.     private final double mu;

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

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

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

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

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

  97.     /** Simple constructor.
  98.      * @param love Love numbers
  99.      * @param deltaCSFunction function computing frequency dependent terms (ΔC₂₀, ΔC₂₁, ΔS₂₁, ΔC₂₂, ΔS₂₂)
  100.      * @param deltaC20PermanentTide permanent tide to be <em>removed</em> from ΔC₂₀ when zero-tide potentials are used
  101.      * @param poleTideFunction function computing pole tide terms (ΔC₂₁, ΔS₂₁), may be null
  102.      * @param centralBodyFrame rotating body frame
  103.      * @param ae central body reference radius
  104.      * @param mu central body attraction coefficient
  105.      * @param centralTideSystem tide system used in the central attraction model
  106.      * @param bodies tide generating bodies (typically Sun and Moon)
  107.      */
  108.     public SolidTidesField(final LoveNumbers love, final TimeFunction<double []> deltaCSFunction,
  109.                            final double deltaC20PermanentTide, final TimeFunction<double []> poleTideFunction,
  110.                            final Frame centralBodyFrame, final double ae, final double mu,
  111.                            final TideSystem centralTideSystem, final CelestialBody ... bodies) {

  113.         // store mode parameters
  114.         this.centralBodyFrame  = centralBodyFrame;
  115.         this.ae                = ae;
  116.         this.mu                = mu;
  117.         this.centralTideSystem = centralTideSystem;
  118.         this.bodies            = bodies;

  120.         // compute recursion coefficients for Legendre functions
  121.         this.anm               = buildTriangularArray(5, false);
  122.         this.bnm               = buildTriangularArray(5, false);
  123.         this.dmm               = new double[love.getSize()];
  124.         recursionCoefficients();

  126.         // Love numbers
  127.         this.love = love;

  129.         // frequency dependent terms
  130.         this.deltaCSFunction = deltaCSFunction;

  132.         // permanent tide
  133.         this.deltaC20PermanentTide = deltaC20PermanentTide;

  135.         // pole tide
  136.         this.poleTideFunction = poleTideFunction;

  138.     }

  140.     /** {@inheritDoc} */
  141.     @Override
  142.     public int getMaxDegree() {
  143.         return MAX_LEGENDRE_DEGREE;
  144.     }

  146.     /** {@inheritDoc} */
  147.     @Override
  148.     public int getMaxOrder() {
  149.         return MAX_LEGENDRE_DEGREE;
  150.     }

  152.     /** {@inheritDoc} */
  153.     @Override
  154.     public double getMu() {
  155.         return mu;
  156.     }

  158.     /** {@inheritDoc} */
  159.     @Override
  160.     public double getAe() {
  161.         return ae;
  162.     }

  164.     /** {@inheritDoc} */
  165.     @Override
  166.     public AbsoluteDate getReferenceDate() {
  167.         return AbsoluteDate.J2000_EPOCH;
  168.     }

  170.     /** {@inheritDoc} */
  171.     @Override
  172.     public double getOffset(final AbsoluteDate date) {
  173.         return date.durationFrom(AbsoluteDate.J2000_EPOCH);
  174.     }

  176.     /** {@inheritDoc} */
  177.     @Override
  178.     public TideSystem getTideSystem() {
  179.         // not really used here, but for consistency we can state that either
  180.         // we add the permanent tide or it was already in the central attraction
  181.         return TideSystem.ZERO_TIDE;
  182.     }

  184.     /** {@inheritDoc} */
  185.     @Override
  186.     @Deprecated
  187.     public double getNormalizedCnm(final double dateOffset, final int n, final int m)
  188.         throws OrekitException {
  189.         return onDate(getReferenceDate().shiftedBy(dateOffset)).getNormalizedCnm(n, m);
  190.     }

  192.     /** {@inheritDoc} */
  193.     @Override
  194.     @Deprecated
  195.     public double getNormalizedSnm(final double dateOffset, final int n, final int m)
  196.         throws OrekitException {
  197.         return onDate(getReferenceDate().shiftedBy(dateOffset)).getNormalizedSnm(n, m);
  198.     }

  200.     /** {@inheritDoc} */
  201.     @Override
  202.     public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) throws OrekitException {

  204.         // computed Cnm and Snm coefficients
  205.         final double[][] cnm = buildTriangularArray(5, true);
  206.         final double[][] snm = buildTriangularArray(5, true);

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

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

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

  218.             // compute polar coordinates
  219.             final double x    = position.getX();
  220.             final double y    = position.getY();
  221.             final double z    = position.getZ();
  222.             final double x2   = x * x;
  223.             final double y2   = y * y;
  224.             final double z2   = z * z;
  225.             final double r2   = x2 + y2 + z2;
  226.             final double r    = FastMath.sqrt (r2);
  227.             final double rho2 = x2 + y2;
  228.             final double rho  = FastMath.sqrt(rho2);

  230.             // evaluate Pnm
  231.             evaluateLegendre(z / r, rho / r, pnm);

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

  236.         }

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

  241.         if (centralTideSystem == TideSystem.ZERO_TIDE) {
  242.             // step 3: remove permanent tide which is already considered
  243.             // in the central body gravity field
  244.             removePermanentTide(cnm);
  245.         }

  247.         if (poleTideFunction != null) {
  248.             // add pole tide
  249.             poleTide(date, cnm, snm);
  250.         }

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

  254.     }

  256.     /** Compute recursion coefficients. */
  257.     private void recursionCoefficients() {

  259.         // pre-compute the recursion coefficients
  260.         // (see equations 11 and 12 from Holmes and Featherstone paper)
  261.         for (int n = 0; n < anm.length; ++n) {
  262.             for (int m = 0; m < n; ++m) {
  263.                 final double f = (n - m) * (n + m );
  264.                 anm[n][m] = FastMath.sqrt((2 * n - 1) * (2 * n + 1) / f);
  265.                 bnm[n][m] = FastMath.sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) / (f * (2 * n - 3)));
  266.             }
  267.         }

  269.         // sectorial terms corresponding to equation 13 in Holmes and Featherstone paper
  270.         dmm[0] = Double.NaN; // dummy initialization for unused term
  271.         dmm[1] = Double.NaN; // dummy initialization for unused term
  272.         for (int m = 2; m < dmm.length; ++m) {
  273.             dmm[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m));
  274.         }

  276.     }

  278.     /** Evaluate Legendre functions.
  279.      * @param t cos(&theta;), where &theta; is the polar angle
  280.      * @param u sin(&theta;), where &theta; is the polar angle
  281.      * @param pnm the computed coefficients. Modified in place.
  282.      */
  283.     private void evaluateLegendre(final double t, final double u, final double[][] pnm) {

  285.         // as the degree is very low, we use the standard forward column method
  286.         // and store everything (see equations 11 and 13 from Holmes and Featherstone paper)
  287.         pnm[0][0] = 1;
  288.         pnm[1][0] = anm[1][0] * t;
  289.         pnm[1][1] = FastMath.sqrt(3) * u;
  290.         for (int m = 2; m < dmm.length; ++m) {
  291.             pnm[m][m - 1] = anm[m][m - 1] * t * pnm[m - 1][m - 1];
  292.             pnm[m][m]     = dmm[m] * u * pnm[m - 1][m - 1];
  293.         }
  294.         for (int m = 0; m < dmm.length; ++m) {
  295.             for (int n = m + 2; n < pnm.length; ++n) {
  296.                 pnm[n][m] = anm[n][m] * t * pnm[n - 1][m] - bnm[n][m] * pnm[n - 2][m];
  297.             }
  298.         }

  300.     }

  302.     /** Update coefficients applying frequency independent step, for one tide generating body.
  303.      * @param r distance to tide generating body
  304.      * @param gm tide generating body attraction coefficient
  305.      * @param cosLambda cosine of the tide generating body longitude
  306.      * @param sinLambda sine of the tide generating body longitude
  307.      * @param pnm the Legendre coefficients. See {@link #evaluateLegendre(double, double, double[][])}.
  308.      * @param cnm the computed Cnm coefficients. Modified in place.
  309.      * @param snm the computed Snm coefficients. Modified in place.
  310.      */
  311.     private void frequencyIndependentPart(final double r,
  312.                                           final double gm,
  313.                                           final double cosLambda,
  314.                                           final double sinLambda,
  315.                                           final double[][] pnm,
  316.                                           final double[][] cnm,
  317.                                           final double[][] snm) {

  319.         final double rRatio = ae / r;
  320.         double fM           = gm / mu;
  321.         double cosMLambda   = 1;
  322.         double sinMLambda   = 0;
  323.         for (int m = 0; m < love.getSize(); ++m) {

  325.             double fNPlus1 = fM;
  326.             for (int n = m; n < love.getSize(); ++n) {
  327.                 fNPlus1 *= rRatio;
  328.                 final double coeff = (fNPlus1 / (2 * n + 1)) * pnm[n][m];
  329.                 final double cCos  = coeff * cosMLambda;
  330.                 final double cSin  = coeff * sinMLambda;

  332.                 // direct effect of degree n tides on degree n coefficients
  333.                 // equation 6.6 from IERS conventions (2010)
  334.                 final double kR = love.getReal(n, m);
  335.                 final double kI = love.getImaginary(n, m);
  336.                 cnm[n][m] += kR * cCos + kI * cSin;
  337.                 snm[n][m] += kR * cSin - kI * cCos;

  339.                 if (n == 2) {
  340.                     // indirect effect of degree 2 tides on degree 4 coefficients
  341.                     // equation 6.7 from IERS conventions (2010)
  342.                     final double kP = love.getPlus(n, m);
  343.                     cnm[4][m] += kP * cCos;
  344.                     snm[4][m] += kP * cSin;
  345.                 }

  347.             }

  349.             // prepare next iteration on order
  350.             final double tmp = cosMLambda * cosLambda - sinMLambda * sinLambda;
  351.             sinMLambda = sinMLambda * cosLambda + cosMLambda * sinLambda;
  352.             cosMLambda = tmp;
  353.             fM        *= rRatio;

  355.         }

  357.     }

  359.     /** Update coefficients applying frequency dependent step.
  360.      * @param date current date
  361.      * @param cnm the Cnm coefficients. Modified in place.
  362.      * @param snm the Snm coefficients. Modified in place.
  363.      */
  364.     private void frequencyDependentPart(final AbsoluteDate date,
  365.                                         final double[][] cnm,
  366.                                         final double[][] snm) {
  367.         final double[] deltaCS = deltaCSFunction.value(date);
  368.         cnm[2][0] += deltaCS[0]; // ΔC₂₀
  369.         cnm[2][1] += deltaCS[1]; // ΔC₂₁
  370.         snm[2][1] += deltaCS[2]; // ΔS₂₁
  371.         cnm[2][2] += deltaCS[3]; // ΔC₂₂
  372.         snm[2][2] += deltaCS[4]; // ΔS₂₂
  373.     }

  375.     /** Remove the permanent tide already counted in zero-tide central gravity fields.
  376.      * @param cnm the Cnm coefficients. Modified in place.
  377.      */
  378.     private void removePermanentTide(final double[][] cnm) {
  379.         cnm[2][0] -= deltaC20PermanentTide;
  380.     }

  382.     /** Update coefficients applying pole tide.
  383.      * @param date current date
  384.      * @param cnm the Cnm coefficients. Modified in place.
  385.      * @param snm the Snm coefficients. Modified in place.
  386.      */
  387.     private void poleTide(final AbsoluteDate date,
  388.                           final double[][] cnm,
  389.                           final double[][] snm) {
  390.         final double[] deltaCS = poleTideFunction.value(date);
  391.         cnm[2][1] += deltaCS[0]; // ΔC₂₁
  392.         snm[2][1] += deltaCS[1]; // ΔS₂₁
  393.     }

  395.     /** Create a triangular array.
  396.      * @param size array size
  397.      * @param withDiagonal if true, the array contains a[p][p] terms, otherwise each
  398.      * row ends up at a[p][p-1]
  399.      * @return new triangular array
  400.      */
  401.     private double[][] buildTriangularArray(final int size, final boolean withDiagonal) {
  402.         final double[][] array = new double[size][];
  403.         for (int i = 0; i < array.length; ++i) {
  404.             array[i] = new double[withDiagonal ? i + 1 : i];
  405.         }
  406.         return array;
  407.     }

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

  412.         /** evaluation date. */
  413.         private final AbsoluteDate date;

  415.         /** Cached cnm. */
  416.         private final double[][] cnm;

  418.         /** Cached snm. */
  419.         private final double[][] snm;

  421.         /** Construct the tidal harmonics on the given date.
  422.          *
  423.          * @param date of evaluation
  424.          * @param cnm the Cnm coefficients. Not copied.
  425.          * @param snm the Snm coeffiecients. Not copied.
  426.          */
  427.         private TideHarmonics(final AbsoluteDate date,
  428.                               final double[][] cnm,
  429.                               final double[][] snm) {
  430.             this.date = date;
  431.             this.cnm = cnm;
  432.             this.snm = snm;
  433.         }

  435.         /** {@inheritDoc} */
  436.         @Override
  437.         public AbsoluteDate getDate() {
  438.             return date;
  439.         }

  441.         /** {@inheritDoc} */
  442.         @Override
  443.         public double getNormalizedCnm(final int n, final int m)
  444.             throws OrekitException {
  445.             return cnm[n][m];
  446.         }

  448.         /** {@inheritDoc} */
  449.         @Override
  450.         public double getNormalizedSnm(final int n, final int m)
  451.             throws OrekitException {
  452.             return snm[n][m];
  453.         }

  455.     }

  457. }
Created with JaCoCo 0.6.3.201306030806