SolidTidesField.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (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.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.util.FastMath;
  21. import org.orekit.bodies.CelestialBody;
  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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  137.     }

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

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

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

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

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

  169.     /** {@inheritDoc} */
  170.     @Override
  171.     public TideSystem getTideSystem() {
  172.         // not really used here, but for consistency we can state that either
  173.         // we add the permanent tide or it was already in the central attraction
  174.         return TideSystem.ZERO_TIDE;
  175.     }

  177.     /** {@inheritDoc} */
  178.     @Override
  179.     public NormalizedSphericalHarmonics onDate(final AbsoluteDate date) {

  181.         // computed Cnm and Snm coefficients
  182.         final double[][] cnm = buildTriangularArray(5, true);
  183.         final double[][] snm = buildTriangularArray(5, true);

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

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

  192.             // compute tide generating body state
  193.             final Vector3D position = body.getPosition(date, centralBodyFrame);

  195.             // compute polar coordinates
  196.             final double x    = position.getX();
  197.             final double y    = position.getY();
  198.             final double z    = position.getZ();
  199.             final double x2   = x * x;
  200.             final double y2   = y * y;
  201.             final double z2   = z * z;
  202.             final double r2   = x2 + y2 + z2;
  203.             final double r    = FastMath.sqrt (r2);
  204.             final double rho2 = x2 + y2;
  205.             final double rho  = FastMath.sqrt(rho2);

  207.             // evaluate Pnm
  208.             evaluateLegendre(z / r, rho / r, pnm);

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

  213.         }

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

  218.         if (centralTideSystem == TideSystem.ZERO_TIDE) {
  219.             // step 3: remove permanent tide which is already considered
  220.             // in the central body gravity field
  221.             removePermanentTide(cnm);
  222.         }

  224.         if (poleTideFunction != null) {
  225.             // add pole tide
  226.             poleTide(date, cnm, snm);
  227.         }

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

  231.     }

  233.     /** Compute recursion coefficients. */
  234.     private void recursionCoefficients() {

  236.         // pre-compute the recursion coefficients
  237.         // (see equations 11 and 12 from Holmes and Featherstone paper)
  238.         for (int n = 0; n < anm.length; ++n) {
  239.             for (int m = 0; m < n; ++m) {
  240.                 final double f = (n - m) * (n + m );
  241.                 anm[n][m] = FastMath.sqrt((2 * n - 1) * (2 * n + 1) / f);
  242.                 bnm[n][m] = FastMath.sqrt((2 * n + 1) * (n + m - 1) * (n - m - 1) / (f * (2 * n - 3)));
  243.             }
  244.         }

  246.         // sectorial terms corresponding to equation 13 in Holmes and Featherstone paper
  247.         dmm[0] = Double.NaN; // dummy initialization for unused term
  248.         dmm[1] = Double.NaN; // dummy initialization for unused term
  249.         for (int m = 2; m < dmm.length; ++m) {
  250.             dmm[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m));
  251.         }

  253.     }

  255.     /** Evaluate Legendre functions.
  256.      * @param t cos(θ), where θ is the polar angle
  257.      * @param u sin(θ), where θ is the polar angle
  258.      * @param pnm the computed coefficients. Modified in place.
  259.      */
  260.     private void evaluateLegendre(final double t, final double u, final double[][] pnm) {

  262.         // as the degree is very low, we use the standard forward column method
  263.         // and store everything (see equations 11 and 13 from Holmes and Featherstone paper)
  264.         pnm[0][0] = 1;
  265.         pnm[1][0] = anm[1][0] * t;
  266.         pnm[1][1] = FastMath.sqrt(3) * u;
  267.         for (int m = 2; m < dmm.length; ++m) {
  268.             pnm[m][m - 1] = anm[m][m - 1] * t * pnm[m - 1][m - 1];
  269.             pnm[m][m]     = dmm[m] * u * pnm[m - 1][m - 1];
  270.         }
  271.         for (int m = 0; m < dmm.length; ++m) {
  272.             for (int n = m + 2; n < pnm.length; ++n) {
  273.                 pnm[n][m] = anm[n][m] * t * pnm[n - 1][m] - bnm[n][m] * pnm[n - 2][m];
  274.             }
  275.         }

  277.     }

  279.     /** Update coefficients applying frequency independent step, for one tide generating body.
  280.      * @param r distance to tide generating body
  281.      * @param gm tide generating body attraction coefficient
  282.      * @param cosLambda cosine of the tide generating body longitude
  283.      * @param sinLambda sine of the tide generating body longitude
  284.      * @param pnm the Legendre coefficients. See {@link #evaluateLegendre(double, double, double[][])}.
  285.      * @param cnm the computed Cnm coefficients. Modified in place.
  286.      * @param snm the computed Snm coefficients. Modified in place.
  287.      */
  288.     private void frequencyIndependentPart(final double r,
  289.                                           final double gm,
  290.                                           final double cosLambda,
  291.                                           final double sinLambda,
  292.                                           final double[][] pnm,
  293.                                           final double[][] cnm,
  294.                                           final double[][] snm) {

  296.         final double rRatio = ae / r;
  297.         double fM           = gm / mu;
  298.         double cosMLambda   = 1;
  299.         double sinMLambda   = 0;
  300.         for (int m = 0; m < love.getSize(); ++m) {

  302.             double fNPlus1 = fM;
  303.             for (int n = m; n < love.getSize(); ++n) {
  304.                 fNPlus1 *= rRatio;
  305.                 final double coeff = (fNPlus1 / (2 * n + 1)) * pnm[n][m];
  306.                 final double cCos  = coeff * cosMLambda;
  307.                 final double cSin  = coeff * sinMLambda;

  309.                 // direct effect of degree n tides on degree n coefficients
  310.                 // equation 6.6 from IERS conventions (2010)
  311.                 final double kR = love.getReal(n, m);
  312.                 final double kI = love.getImaginary(n, m);
  313.                 cnm[n][m] += kR * cCos + kI * cSin;
  314.                 snm[n][m] += kR * cSin - kI * cCos;

  316.                 if (n == 2) {
  317.                     // indirect effect of degree 2 tides on degree 4 coefficients
  318.                     // equation 6.7 from IERS conventions (2010)
  319.                     final double kP = love.getPlus(n, m);
  320.                     cnm[4][m] += kP * cCos;
  321.                     snm[4][m] += kP * cSin;
  322.                 }

  324.             }

  326.             // prepare next iteration on order
  327.             final double tmp = cosMLambda * cosLambda - sinMLambda * sinLambda;
  328.             sinMLambda = sinMLambda * cosLambda + cosMLambda * sinLambda;
  329.             cosMLambda = tmp;
  330.             fM        *= rRatio;

  332.         }

  334.     }

  336.     /** Update coefficients applying frequency dependent step.
  337.      * @param date current date
  338.      * @param cnm the Cnm coefficients. Modified in place.
  339.      * @param snm the Snm coefficients. Modified in place.
  340.      */
  341.     private void frequencyDependentPart(final AbsoluteDate date,
  342.                                         final double[][] cnm,
  343.                                         final double[][] snm) {
  344.         final double[] deltaCS = deltaCSFunction.value(date);
  345.         cnm[2][0] += deltaCS[0]; // ΔC₂₀
  346.         cnm[2][1] += deltaCS[1]; // ΔC₂₁
  347.         snm[2][1] += deltaCS[2]; // ΔS₂₁
  348.         cnm[2][2] += deltaCS[3]; // ΔC₂₂
  349.         snm[2][2] += deltaCS[4]; // ΔS₂₂
  350.     }

  352.     /** Remove the permanent tide already counted in zero-tide central gravity fields.
  353.      * @param cnm the Cnm coefficients. Modified in place.
  354.      */
  355.     private void removePermanentTide(final double[][] cnm) {
  356.         cnm[2][0] -= deltaC20PermanentTide;
  357.     }

  359.     /** Update coefficients applying pole tide.
  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 poleTide(final AbsoluteDate date,
  365.                           final double[][] cnm,
  366.                           final double[][] snm) {
  367.         final double[] deltaCS = poleTideFunction.value(date);
  368.         cnm[2][1] += deltaCS[0]; // ΔC₂₁
  369.         snm[2][1] += deltaCS[1]; // ΔS₂₁
  370.     }

  372.     /** Create a triangular array.
  373.      * @param size array size
  374.      * @param withDiagonal if true, the array contains a[p][p] terms, otherwise each
  375.      * row ends up at a[p][p-1]
  376.      * @return new triangular array
  377.      */
  378.     private double[][] buildTriangularArray(final int size, final boolean withDiagonal) {
  379.         final double[][] array = new double[size][];
  380.         for (int i = 0; i < array.length; ++i) {
  381.             array[i] = new double[withDiagonal ? i + 1 : i];
  382.         }
  383.         return array;
  384.     }

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

  389.         /** evaluation date. */
  390.         private final AbsoluteDate date;

  392.         /** Cached cnm. */
  393.         private final double[][] cnm;

  395.         /** Cached snm. */
  396.         private final double[][] snm;

  398.         /** Construct the tidal harmonics on the given date.
  399.          *
  400.          * @param date of evaluation
  401.          * @param cnm the Cnm coefficients. Not copied.
  402.          * @param snm the Snm coeffiecients. Not copied.
  403.          */
  404.         private TideHarmonics(final AbsoluteDate date,
  405.                               final double[][] cnm,
  406.                               final double[][] snm) {
  407.             this.date = date;
  408.             this.cnm = cnm;
  409.             this.snm = snm;
  410.         }

  412.         /** {@inheritDoc} */
  413.         @Override
  414.         public AbsoluteDate getDate() {
  415.             return date;
  416.         }

  418.         /** {@inheritDoc} */
  419.         @Override
  420.         public double getNormalizedCnm(final int n, final int m) {
  421.             return cnm[n][m];
  422.         }

  424.         /** {@inheritDoc} */
  425.         @Override
  426.         public double getNormalizedSnm(final int n, final int m) {
  427.             return snm[n][m];
  428.         }

  430.     }

  432. }
Created with JaCoCo 0.8.12.202403310830