OceanLoading.java

  1. /* Copyright 2002-2024 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.models.earth.displacement;

  19. import java.util.HashMap;
  20. import java.util.Map;
  21. import java.util.function.Function;

  23. import org.hipparchus.analysis.UnivariateFunction;
  24. import org.hipparchus.analysis.interpolation.SplineInterpolator;
  25. import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.SinCos;
  29. import org.orekit.bodies.GeodeticPoint;
  30. import org.orekit.bodies.OneAxisEllipsoid;
  31. import org.orekit.data.BodiesElements;
  32. import org.orekit.frames.Frame;

  34. /**
  35.  * Modeling of displacement of reference points due to ocean loading.
  36.  * <p>
  37.  * This class implements the same model as IERS HARDIP.F program. For a
  38.  * given site, this model uses a set of amplitudes and phases for the 11
  39.  * main tides (M₂, S₂, N₂, K₂, K₁, O₁, P₁, Q₁, Mf, Mm, and Ssa) in BLQ
  40.  * format as provided by the <a
  41.  * href="http://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>
  42.  * at Onsala Space Observatory. From these elements, additional admittances
  43.  * are derived using spline interpolation based on tides frequencies for
  44.  * a total of 342 tides, including the 11 main tides.
  45.  * </p>
  46.  * <p>
  47.  * This implementation is a complete rewrite of the original HARDISP.F program
  48.  * developed by Duncan Agnew and copyright 2008 IERS Conventions center. This
  49.  * derived work is not endorsed by the IERS conventions center. What remains
  50.  * from the original program is the model (spline interpolation and coefficients).
  51.  * The code by itself is completely different, using the underlying mathematical
  52.  * library for spline interpolation and the existing Orekit features for nutation
  53.  * arguments, time and time scales handling, tides modeling...
  54.  * </p>
  55.  * <p>
  56.  * Instances of this class are guaranteed to be immutable
  57.  * </p>
  58.  * <p>
  59.  * The original HARDISP.F program is distributed with the following notice:
  60.  * </p>
  61.  * <pre>
  62.  *  Copyright (C) 2008
  63.  *  IERS Conventions Center
  64.  *
  65.  *  ==================================
  66.  *  IERS Conventions Software License
  67.  *  ==================================
  68.  *
  69.  *  NOTICE TO USER:
  70.  *
  71.  *  BY USING THIS SOFTWARE YOU ACCEPT THE FOLLOWING TERMS AND CONDITIONS
  72.  *  WHICH APPLY TO ITS USE.
  73.  *
  74.  *  1. The Software is provided by the IERS Conventions Center ("the
  75.  *     Center").
  76.  *
  77.  *  2. Permission is granted to anyone to use the Software for any
  78.  *     purpose, including commercial applications, free of charge,
  79.  *     subject to the conditions and restrictions listed below.
  80.  *
  81.  *  3. You (the user) may adapt the Software and its algorithms for your
  82.  *     own purposes and you may distribute the resulting "derived work"
  83.  *     to others, provided that the derived work complies with the
  84.  *     following requirements:
  85.  *
  86.  *     a) Your work shall be clearly identified so that it cannot be
  87.  *        mistaken for IERS Conventions software and that it has been
  88.  *        neither distributed by nor endorsed by the Center.
  89.  *
  90.  *     b) Your work (including source code) must contain descriptions of
  91.  *        how the derived work is based upon and/or differs from the
  92.  *        original Software.
  93.  *
  94.  *     c) The name(s) of all modified routine(s) that you distribute
  95.  *        shall be changed.
  96.  *
  97.  *     d) The origin of the IERS Conventions components of your derived
  98.  *        work must not be misrepresented; you must not claim that you
  99.  *        wrote the original Software.
  100.  *
  101.  *     e) The source code must be included for all routine(s) that you
  102.  *        distribute.  This notice must be reproduced intact in any
  103.  *        source distribution.
  104.  *
  105.  *  4. In any published work produced by the user and which includes
  106.  *     results achieved by using the Software, you shall acknowledge
  107.  *     that the Software was used in obtaining those results.
  108.  *
  109.  *  5. The Software is provided to the user "as is" and the Center makes
  110.  *     no warranty as to its use or performance.   The Center does not
  111.  *     and cannot warrant the performance or results which the user may
  112.  *     obtain by using the Software.  The Center makes no warranties,
  113.  *     express or implied, as to non-infringement of third party rights,
  114.  *     merchantability, or fitness for any particular purpose.  In no
  115.  *     event will the Center be liable to the user for any consequential,
  116.  *     incidental, or special damages, including any lost profits or lost
  117.  *     savings, even if a Center representative has been advised of such
  118.  *     damages, or for any claim by any third party.
  119.  *
  120.  *  Correspondence concerning IERS Conventions software should be
  121.  *  addressed as follows:
  122.  *
  123.  *                     Gerard Petit
  124.  *     Internet email: gpetit[at]bipm.org
  125.  *     Postal address: IERS Conventions Center
  126.  *                     Time, frequency and gravimetry section, BIPM
  127.  *                     Pavillon de Breteuil
  128.  *                     92312 Sevres  FRANCE
  129.  *
  130.  *     or
  131.  *
  132.  *                     Brian Luzum
  133.  *     Internet email: brian.luzum[at]usno.navy.mil
  134.  *     Postal address: IERS Conventions Center
  135.  *                     Earth Orientation Department
  136.  *                     3450 Massachusetts Ave, NW
  137.  *                     Washington, DC 20392
  138.  * </pre>
  139.  * @see org.orekit.estimation.measurements.GroundStation
  140.  * @since 9.1
  141.  * @author Luc Maisonobe
  142.  */
  143. public class OceanLoading implements StationDisplacement {

  145.     // CHECKSTYLE: stop Indentation check
  146.     /** Amplitudes of all tides used. */
  147.     private static final double[] CARTWRIGHT_EDDEN_AMPLITUDE = {
  148.         0.632208,  0.294107,  0.121046,  0.079915,  0.023818, -0.023589,  0.022994,
  149.         0.019333, -0.017871,  0.017192,  0.016018,  0.004671, -0.004662, -0.004519,
  150.         0.004470,  0.004467,  0.002589, -0.002455, -0.002172,  0.001972,  0.001947,
  151.         0.001914, -0.001898,  0.001802,  0.001304,  0.001170,  0.001130,  0.001061,
  152.        -0.001022, -0.001017,  0.001014,  0.000901, -0.000857,  0.000855,  0.000855,
  153.         0.000772,  0.000741,  0.000741, -0.000721,  0.000698,  0.000658,  0.000654,
  154.        -0.000653,  0.000633,  0.000626, -0.000598,  0.000590,  0.000544,  0.000479,
  155.        -0.000464,  0.000413, -0.000390,  0.000373,  0.000366,  0.000366, -0.000360,
  156.        -0.000355,  0.000354,  0.000329,  0.000328,  0.000319,  0.000302,  0.000279,
  157.        -0.000274, -0.000272,  0.000248, -0.000225,  0.000224, -0.000223, -0.000216,
  158.         0.000211,  0.000209,  0.000194,  0.000185, -0.000174, -0.000171,  0.000159,
  159.         0.000131,  0.000127,  0.000120,  0.000118,  0.000117,  0.000108,  0.000107,
  160.         0.000105, -0.000102,  0.000102,  0.000099, -0.000096,  0.000095, -0.000089,
  161.        -0.000085, -0.000084, -0.000081, -0.000077, -0.000072, -0.000067,  0.000066,
  162.         0.000064,  0.000063,  0.000063,  0.000063,  0.000062,  0.000062, -0.000060,
  163.         0.000056,  0.000053,  0.000051,  0.000050,  0.368645, -0.262232, -0.121995,
  164.        -0.050208,  0.050031, -0.049470,  0.020620,  0.020613,  0.011279, -0.009530,
  165.        -0.009469, -0.008012,  0.007414, -0.007300,  0.007227, -0.007131, -0.006644,
  166.         0.005249,  0.004137,  0.004087,  0.003944,  0.003943,  0.003420,  0.003418,
  167.         0.002885,  0.002884,  0.002160, -0.001936,  0.001934, -0.001798,  0.001690,
  168.         0.001689,  0.001516,  0.001514, -0.001511,  0.001383,  0.001372,  0.001371,
  169.        -0.001253, -0.001075,  0.001020,  0.000901,  0.000865, -0.000794,  0.000788,
  170.         0.000782, -0.000747, -0.000745,  0.000670, -0.000603, -0.000597,  0.000542,
  171.         0.000542, -0.000541, -0.000469, -0.000440,  0.000438,  0.000422,  0.000410,
  172.        -0.000374, -0.000365,  0.000345,  0.000335, -0.000321, -0.000319,  0.000307,
  173.         0.000291,  0.000290, -0.000289,  0.000286,  0.000275,  0.000271,  0.000263,
  174.        -0.000245,  0.000225,  0.000225,  0.000221, -0.000202, -0.000200, -0.000199,
  175.         0.000192,  0.000183,  0.000183,  0.000183, -0.000170,  0.000169,  0.000168,
  176.         0.000162,  0.000149, -0.000147, -0.000141,  0.000138,  0.000136,  0.000136,
  177.         0.000127,  0.000127, -0.000126, -0.000121, -0.000121,  0.000117, -0.000116,
  178.        -0.000114, -0.000114, -0.000114,  0.000114,  0.000113,  0.000109,  0.000108,
  179.         0.000106, -0.000106, -0.000106,  0.000105,  0.000104, -0.000103, -0.000100,
  180.        -0.000100, -0.000100,  0.000099, -0.000098,  0.000093,  0.000093,  0.000090,
  181.        -0.000088,  0.000083, -0.000083, -0.000082, -0.000081, -0.000079, -0.000077,
  182.        -0.000075, -0.000075, -0.000075,  0.000071,  0.000071, -0.000071,  0.000068,
  183.         0.000068,  0.000065,  0.000065,  0.000064,  0.000064,  0.000064, -0.000064,
  184.        -0.000060,  0.000056,  0.000056,  0.000053,  0.000053,  0.000053, -0.000053,
  185.         0.000053,  0.000053,  0.000052,  0.000050, -0.066607, -0.035184, -0.030988,
  186.         0.027929, -0.027616, -0.012753, -0.006728, -0.005837, -0.005286, -0.004921,
  187.        -0.002884, -0.002583, -0.002422,  0.002310,  0.002283, -0.002037,  0.001883,
  188.        -0.001811, -0.001687, -0.001004, -0.000925, -0.000844,  0.000766,  0.000766,
  189.        -0.000700, -0.000495, -0.000492,  0.000491,  0.000483,  0.000437, -0.000416,
  190.        -0.000384,  0.000374, -0.000312, -0.000288, -0.000273,  0.000259,  0.000245,
  191.        -0.000232,  0.000229, -0.000216,  0.000206, -0.000204, -0.000202,  0.000200,
  192.         0.000195, -0.000190,  0.000187,  0.000180, -0.000179,  0.000170,  0.000153,
  193.        -0.000137, -0.000119, -0.000119, -0.000112, -0.000110, -0.000110,  0.000107,
  194.        -0.000095, -0.000095, -0.000091, -0.000090, -0.000081, -0.000079, -0.000079,
  195.         0.000077, -0.000073,  0.000069, -0.000067, -0.000066,  0.000065,  0.000064,
  196.        -0.000062,  0.000060,  0.000059, -0.000056,  0.000055, -0.000051
  197.     };
  198.     // CHECKSTYLE: resume Indentation check

  200.     // CHECKSTYLE: stop NoWhitespaceAfter check
  201.     /** Doodson arguments for all tides used. */
  202.     private static final int[][] DOODSON_ARGUMENTS = {
  203.         { 2,  0,  0,  0,  0,  0 }, { 2,  2, -2,  0,  0,  0 }, { 2, -1,  0,  1,  0,  0 },
  204.         { 2,  2,  0,  0,  0,  0 }, { 2,  2,  0,  0,  1,  0 }, { 2,  0,  0,  0, -1,  0 },
  205.         { 2, -1,  2, -1,  0,  0 }, { 2, -2,  2,  0,  0,  0 }, { 2,  1,  0, -1,  0,  0 },
  206.         { 2,  2, -3,  0,  0,  1 }, { 2, -2,  0,  2,  0,  0 }, { 2, -3,  2,  1,  0,  0 },
  207.         { 2,  1, -2,  1,  0,  0 }, { 2, -1,  0,  1, -1,  0 }, { 2,  3,  0, -1,  0,  0 },
  208.         { 2,  1,  0,  1,  0,  0 }, { 2,  2,  0,  0,  2,  0 }, { 2,  2, -1,  0,  0, -1 },
  209.         { 2,  0, -1,  0,  0,  1 }, { 2,  1,  0,  1,  1,  0 }, { 2,  3,  0, -1,  1,  0 },
  210.         { 2,  0,  1,  0,  0, -1 }, { 2,  0, -2,  2,  0,  0 }, { 2, -3,  0,  3,  0,  0 },
  211.         { 2, -2,  3,  0,  0, -1 }, { 2,  4,  0,  0,  0,  0 }, { 2, -1,  1,  1,  0, -1 },
  212.         { 2, -1,  3, -1,  0, -1 }, { 2,  2,  0,  0, -1,  0 }, { 2, -1, -1,  1,  0,  1 },
  213.         { 2,  4,  0,  0,  1,  0 }, { 2, -3,  4, -1,  0,  0 }, { 2, -1,  2, -1, -1,  0 },
  214.         { 2,  3, -2,  1,  0,  0 }, { 2,  1,  2, -1,  0,  0 }, { 2, -4,  2,  2,  0,  0 },
  215.         { 2,  4, -2,  0,  0,  0 }, { 2,  0,  2,  0,  0,  0 }, { 2, -2,  2,  0, -1,  0 },
  216.         { 2,  2, -4,  0,  0,  2 }, { 2,  2, -2,  0, -1,  0 }, { 2,  1,  0, -1, -1,  0 },
  217.         { 2, -1,  1,  0,  0,  0 }, { 2,  2, -1,  0,  0,  1 }, { 2,  2,  1,  0,  0, -1 },
  218.         { 2, -2,  0,  2, -1,  0 }, { 2, -2,  4, -2,  0,  0 }, { 2,  2,  2,  0,  0,  0 },
  219.         { 2, -4,  4,  0,  0,  0 }, { 2, -1,  0, -1, -2,  0 }, { 2,  1,  2, -1,  1,  0 },
  220.         { 2, -1, -2,  3,  0,  0 }, { 2,  3, -2,  1,  1,  0 }, { 2,  4,  0, -2,  0,  0 },
  221.         { 2,  0,  0,  2,  0,  0 }, { 2,  0,  2, -2,  0,  0 }, { 2,  0,  2,  0,  1,  0 },
  222.         { 2, -3,  3,  1,  0, -1 }, { 2,  0,  0,  0, -2,  0 }, { 2,  4,  0,  0,  2,  0 },
  223.         { 2,  4, -2,  0,  1,  0 }, { 2,  0,  0,  0,  0,  2 }, { 2,  1,  0,  1,  2,  0 },
  224.         { 2,  0, -2,  0, -2,  0 }, { 2, -2,  1,  0,  0,  1 }, { 2, -2,  1,  2,  0, -1 },
  225.         { 2, -1,  1, -1,  0,  1 }, { 2,  5,  0, -1,  0,  0 }, { 2,  1, -3,  1,  0,  1 },
  226.         { 2, -2, -1,  2,  0,  1 }, { 2,  3,  0, -1,  2,  0 }, { 2,  1, -2,  1, -1,  0 },
  227.         { 2,  5,  0, -1,  1,  0 }, { 2, -4,  0,  4,  0,  0 }, { 2, -3,  2,  1, -1,  0 },
  228.         { 2, -2,  1,  1,  0,  0 }, { 2,  4,  0, -2,  1,  0 }, { 2,  0,  0,  2,  1,  0 },
  229.         { 2, -5,  4,  1,  0,  0 }, { 2,  0,  2,  0,  2,  0 }, { 2, -1,  2,  1,  0,  0 },
  230.         { 2,  5, -2, -1,  0,  0 }, { 2,  1, -1,  0,  0,  0 }, { 2,  2, -2,  0,  0,  2 },
  231.         { 2, -5,  2,  3,  0,  0 }, { 2, -1, -2,  1, -2,  0 }, { 2, -3,  5, -1,  0, -1 },
  232.         { 2, -1,  0,  0,  0,  1 }, { 2, -2,  0,  0, -2,  0 }, { 2,  0, -1,  1,  0,  0 },
  233.         { 2, -3,  1,  1,  0,  1 }, { 2,  3,  0, -1, -1,  0 }, { 2,  1,  0,  1, -1,  0 },
  234.         { 2, -1,  2,  1,  1,  0 }, { 2,  0, -3,  2,  0,  1 }, { 2,  1, -1, -1,  0,  1 },
  235.         { 2, -3,  0,  3, -1,  0 }, { 2,  0, -2,  2, -1,  0 }, { 2, -4,  3,  2,  0, -1 },
  236.         { 2, -1,  0,  1, -2,  0 }, { 2,  5,  0, -1,  2,  0 }, { 2, -4,  5,  0,  0, -1 },
  237.         { 2, -2,  4,  0,  0, -2 }, { 2, -1,  0,  1,  0,  2 }, { 2, -2, -2,  4,  0,  0 },
  238.         { 2,  3, -2, -1, -1,  0 }, { 2, -2,  5, -2,  0, -1 }, { 2,  0, -1,  0, -1,  1 },
  239.         { 2,  5, -2, -1,  1,  0 }, { 1,  1,  0,  0,  0,  0 }, { 1, -1,  0,  0,  0,  0 },
  240.         { 1,  1, -2,  0,  0,  0 }, { 1, -2,  0,  1,  0,  0 }, { 1,  1,  0,  0,  1,  0 },
  241.         { 1, -1,  0,  0, -1,  0 }, { 1,  2,  0, -1,  0,  0 }, { 1,  0,  0,  1,  0,  0 },
  242.         { 1,  3,  0,  0,  0,  0 }, { 1, -2,  2, -1,  0,  0 }, { 1, -2,  0,  1, -1,  0 },
  243.         { 1, -3,  2,  0,  0,  0 }, { 1,  0,  0, -1,  0,  0 }, { 1,  1,  0,  0, -1,  0 },
  244.         { 1,  3,  0,  0,  1,  0 }, { 1,  1, -3,  0,  0,  1 }, { 1, -3,  0,  2,  0,  0 },
  245.         { 1,  1,  2,  0,  0,  0 }, { 1,  0,  0,  1,  1,  0 }, { 1,  2,  0, -1,  1,  0 },
  246.         { 1,  0,  2, -1,  0,  0 }, { 1,  2, -2,  1,  0,  0 }, { 1,  3, -2,  0,  0,  0 },
  247.         { 1, -1,  2,  0,  0,  0 }, { 1,  1,  1,  0,  0, -1 }, { 1,  1, -1,  0,  0,  1 },
  248.         { 1,  4,  0, -1,  0,  0 }, { 1, -4,  2,  1,  0,  0 }, { 1,  0, -2,  1,  0,  0 },
  249.         { 1, -2,  2, -1, -1,  0 }, { 1,  3,  0, -2,  0,  0 }, { 1, -1,  0,  2,  0,  0 },
  250.         { 1, -1,  0,  0, -2,  0 }, { 1,  3,  0,  0,  2,  0 }, { 1, -3,  2,  0, -1,  0 },
  251.         { 1,  4,  0, -1,  1,  0 }, { 1,  0,  0, -1, -1,  0 }, { 1,  1, -2,  0, -1,  0 },
  252.         { 1, -3,  0,  2, -1,  0 }, { 1,  1,  0,  0,  2,  0 }, { 1,  1, -1,  0,  0, -1 },
  253.         { 1, -1, -1,  0,  0,  1 }, { 1,  0,  2, -1,  1,  0 }, { 1, -1,  1,  0,  0, -1 },
  254.         { 1, -1, -2,  2,  0,  0 }, { 1,  2, -2,  1,  1,  0 }, { 1, -4,  0,  3,  0,  0 },
  255.         { 1, -1,  2,  0,  1,  0 }, { 1,  3, -2,  0,  1,  0 }, { 1,  2,  0, -1, -1,  0 },
  256.         { 1,  0,  0,  1, -1,  0 }, { 1, -2,  2,  1,  0,  0 }, { 1,  4, -2, -1,  0,  0 },
  257.         { 1, -3,  3,  0,  0, -1 }, { 1, -2,  1,  1,  0, -1 }, { 1, -2,  3, -1,  0, -1 },
  258.         { 1,  0, -2,  1, -1,  0 }, { 1, -2, -1,  1,  0,  1 }, { 1,  4, -2,  1,  0,  0 },
  259.         { 1, -4,  4, -1,  0,  0 }, { 1, -4,  2,  1, -1,  0 }, { 1,  5, -2,  0,  0,  0 },
  260.         { 1,  3,  0, -2,  1,  0 }, { 1, -5,  2,  2,  0,  0 }, { 1,  2,  0,  1,  0,  0 },
  261.         { 1,  1,  3,  0,  0, -1 }, { 1, -2,  0,  1, -2,  0 }, { 1,  4,  0, -1,  2,  0 },
  262.         { 1,  1, -4,  0,  0,  2 }, { 1,  5,  0, -2,  0,  0 }, { 1, -1,  0,  2,  1,  0 },
  263.         { 1, -2,  1,  0,  0,  0 }, { 1,  4, -2,  1,  1,  0 }, { 1, -3,  4, -2,  0,  0 },
  264.         { 1, -1,  3,  0,  0, -1 }, { 1,  3, -3,  0,  0,  1 }, { 1,  5, -2,  0,  1,  0 },
  265.         { 1,  1,  2,  0,  1,  0 }, { 1,  2,  0,  1,  1,  0 }, { 1, -5,  4,  0,  0,  0 },
  266.         { 1, -2,  0, -1, -2,  0 }, { 1,  5,  0, -2,  1,  0 }, { 1,  1,  2, -2,  0,  0 },
  267.         { 1,  1, -2,  2,  0,  0 }, { 1, -2,  2,  1,  1,  0 }, { 1,  0,  3, -1,  0, -1 },
  268.         { 1,  2, -3,  1,  0,  1 }, { 1, -2, -2,  3,  0,  0 }, { 1, -1,  2, -2,  0,  0 },
  269.         { 1, -4,  3,  1,  0, -1 }, { 1, -4,  0,  3, -1,  0 }, { 1, -1, -2,  2, -1,  0 },
  270.         { 1, -2,  0,  3,  0,  0 }, { 1,  4,  0, -3,  0,  0 }, { 1,  0,  1,  1,  0, -1 },
  271.         { 1,  2, -1, -1,  0,  1 }, { 1,  2, -2,  1, -1,  0 }, { 1,  0,  0, -1, -2,  0 },
  272.         { 1,  2,  0,  1,  2,  0 }, { 1,  2, -2, -1, -1,  0 }, { 1,  0,  0,  1,  2,  0 },
  273.         { 1,  0,  1,  0,  0,  0 }, { 1,  2, -1,  0,  0,  0 }, { 1,  0,  2, -1, -1,  0 },
  274.         { 1, -1, -2,  0, -2,  0 }, { 1, -3,  1,  0,  0,  1 }, { 1,  3, -2,  0, -1,  0 },
  275.         { 1, -1, -1,  0, -1,  1 }, { 1,  4, -2, -1,  1,  0 }, { 1,  2,  1, -1,  0, -1 },
  276.         { 1,  0, -1,  1,  0,  1 }, { 1, -2,  4, -1,  0,  0 }, { 1,  4, -4,  1,  0,  0 },
  277.         { 1, -3,  1,  2,  0, -1 }, { 1, -3,  3,  0, -1, -1 }, { 1,  1,  2,  0,  2,  0 },
  278.         { 1,  1, -2,  0, -2,  0 }, { 1,  3,  0,  0,  3,  0 }, { 1, -1,  2,  0, -1,  0 },
  279.         { 1, -2,  1, -1,  0,  1 }, { 1,  0, -3,  1,  0,  1 }, { 1, -3, -1,  2,  0,  1 },
  280.         { 1,  2,  0, -1,  2,  0 }, { 1,  6, -2, -1,  0,  0 }, { 1,  2,  2, -1,  0,  0 },
  281.         { 1, -1,  1,  0, -1, -1 }, { 1, -2,  3, -1, -1, -1 }, { 1, -1,  0,  0,  0,  2 },
  282.         { 1, -5,  0,  4,  0,  0 }, { 1,  1,  0,  0,  0, -2 }, { 1, -2,  1,  1, -1, -1 },
  283.         { 1,  1, -1,  0,  1,  1 }, { 1,  1,  2,  0,  0, -2 }, { 1, -3,  1,  1,  0,  0 },
  284.         { 1, -4,  4, -1, -1,  0 }, { 1,  1,  0, -2, -1,  0 }, { 1, -2, -1,  1, -1,  1 },
  285.         { 1, -3,  2,  2,  0,  0 }, { 1,  5, -2, -2,  0,  0 }, { 1,  3, -4,  2,  0,  0 },
  286.         { 1,  1, -2,  0,  0,  2 }, { 1, -1,  4, -2,  0,  0 }, { 1,  2,  2, -1,  1,  0 },
  287.         { 1, -5,  2,  2, -1,  0 }, { 1,  1, -3,  0, -1,  1 }, { 1,  1,  1,  0,  1, -1 },
  288.         { 1,  6, -2, -1,  1,  0 }, { 1, -2,  2, -1, -2,  0 }, { 1,  4, -2,  1,  2,  0 },
  289.         { 1, -6,  4,  1,  0,  0 }, { 1,  5, -4,  0,  0,  0 }, { 1, -3,  4,  0,  0,  0 },
  290.         { 1,  1,  2, -2,  1,  0 }, { 1, -2,  1,  0, -1,  0 }, { 0,  2,  0,  0,  0,  0 },
  291.         { 0,  1,  0, -1,  0,  0 }, { 0,  0,  2,  0,  0,  0 }, { 0,  0,  0,  0,  1,  0 },
  292.         { 0,  2,  0,  0,  1,  0 }, { 0,  3,  0, -1,  0,  0 }, { 0,  1, -2,  1,  0,  0 },
  293.         { 0,  2, -2,  0,  0,  0 }, { 0,  3,  0, -1,  1,  0 }, { 0,  0,  1,  0,  0, -1 },
  294.         { 0,  2,  0, -2,  0,  0 }, { 0,  2,  0,  0,  2,  0 }, { 0,  3, -2,  1,  0,  0 },
  295.         { 0,  1,  0, -1, -1,  0 }, { 0,  1,  0, -1,  1,  0 }, { 0,  4, -2,  0,  0,  0 },
  296.         { 0,  1,  0,  1,  0,  0 }, { 0,  0,  3,  0,  0, -1 }, { 0,  4,  0, -2,  0,  0 },
  297.         { 0,  3, -2,  1,  1,  0 }, { 0,  3, -2, -1,  0,  0 }, { 0,  4, -2,  0,  1,  0 },
  298.         { 0,  0,  2,  0,  1,  0 }, { 0,  1,  0,  1,  1,  0 }, { 0,  4,  0, -2,  1,  0 },
  299.         { 0,  3,  0, -1,  2,  0 }, { 0,  5, -2, -1,  0,  0 }, { 0,  1,  2, -1,  0,  0 },
  300.         { 0,  1, -2,  1, -1,  0 }, { 0,  1, -2,  1,  1,  0 }, { 0,  2, -2,  0, -1,  0 },
  301.         { 0,  2, -3,  0,  0,  1 }, { 0,  2, -2,  0,  1,  0 }, { 0,  0,  2, -2,  0,  0 },
  302.         { 0,  1, -3,  1,  0,  1 }, { 0,  0,  0,  0,  2,  0 }, { 0,  0,  1,  0,  0,  1 },
  303.         { 0,  1,  2, -1,  1,  0 }, { 0,  3,  0, -3,  0,  0 }, { 0,  2,  1,  0,  0, -1 },
  304.         { 0,  1, -1, -1,  0,  1 }, { 0,  1,  0,  1,  2,  0 }, { 0,  5, -2, -1,  1,  0 },
  305.         { 0,  2, -1,  0,  0,  1 }, { 0,  2,  2, -2,  0,  0 }, { 0,  1, -1,  0,  0,  0 },
  306.         { 0,  5,  0, -3,  0,  0 }, { 0,  2,  0, -2,  1,  0 }, { 0,  1,  1, -1,  0, -1 },
  307.         { 0,  3, -4,  1,  0,  0 }, { 0,  0,  2,  0,  2,  0 }, { 0,  2,  0, -2, -1,  0 },
  308.         { 0,  4, -3,  0,  0,  1 }, { 0,  3, -1, -1,  0,  1 }, { 0,  0,  2,  0,  0, -2 },
  309.         { 0,  3, -3,  1,  0,  1 }, { 0,  2, -4,  2,  0,  0 }, { 0,  4, -2, -2,  0,  0 },
  310.         { 0,  3,  1, -1,  0, -1 }, { 0,  5, -4,  1,  0,  0 }, { 0,  3, -2, -1, -1,  0 },
  311.         { 0,  3, -2,  1,  2,  0 }, { 0,  4, -4,  0,  0,  0 }, { 0,  6, -2, -2,  0,  0 },
  312.         { 0,  5,  0, -3,  1,  0 }, { 0,  4, -2,  0,  2,  0 }, { 0,  2,  2, -2,  1,  0 },
  313.         { 0,  0,  4,  0,  0, -2 }, { 0,  3, -1,  0,  0,  0 }, { 0,  3, -3, -1,  0,  1 },
  314.         { 0,  4,  0, -2,  2,  0 }, { 0,  1, -2, -1, -1,  0 }, { 0,  2, -1,  0,  0, -1 },
  315.         { 0,  4, -4,  2,  0,  0 }, { 0,  2,  1,  0,  1, -1 }, { 0,  3, -2, -1,  1,  0 },
  316.         { 0,  4, -3,  0,  1,  1 }, { 0,  2,  0,  0,  3,  0 }, { 0,  6, -4,  0,  0,  0 }
  317.     };
  318.     // CHECKSTYLE: resume NoWhitespaceAfter check

  320.     /** Cartwright-Edden amplitudes for all tides. */
  321.     private static final Map<Tide, Double> CARTWRIGHT_EDDEN_AMPLITUDE_MAP;

  323.     static {
  324.         CARTWRIGHT_EDDEN_AMPLITUDE_MAP = new HashMap<>(CARTWRIGHT_EDDEN_AMPLITUDE.length);
  325.         for (int i = 0; i < CARTWRIGHT_EDDEN_AMPLITUDE.length; ++i) {
  326.             CARTWRIGHT_EDDEN_AMPLITUDE_MAP.put(new Tide(DOODSON_ARGUMENTS[i][0], DOODSON_ARGUMENTS[i][1], DOODSON_ARGUMENTS[i][2],
  327.                                                         DOODSON_ARGUMENTS[i][3], DOODSON_ARGUMENTS[i][4], DOODSON_ARGUMENTS[i][5]),
  328.                                                CARTWRIGHT_EDDEN_AMPLITUDE[i]);
  329.         }
  330.     }

  332.     /** Earth shape. */
  333.     private final OneAxisEllipsoid earth;

  335.     /** Data for main tides, for which we have ocean loading coefficients. */
  336.     private final MainTideData[][] mainTides;

  338.     /** Simple constructor.
  339.      * @param earth Earth shape
  340.      * @param coefficients coefficients for the considered site
  341.      * @see OceanLoadingCoefficientsBLQFactory
  342.      */
  343.     public OceanLoading(final OneAxisEllipsoid earth, final OceanLoadingCoefficients coefficients) {

  345.         this.earth = earth;

  347.         // set up complex admittances, scaled to Cartwright-Edden amplitudes
  348.         // and grouped by species (0: long period, 1: diurnal, 2: semi-diurnal)
  349.         mainTides  = new MainTideData[coefficients.getNbSpecies()][];
  350.         for (int i = 0; i < mainTides.length; ++i) {
  351.             mainTides[i] = new MainTideData[coefficients.getNbTides(i)];
  352.             for (int j = 0; j < mainTides[i].length; ++j) {
  353.                 final double amplitude = CARTWRIGHT_EDDEN_AMPLITUDE_MAP.get(coefficients.getTide(i, j));
  354.                 mainTides[i][j] = new MainTideData(coefficients, i, j, FastMath.abs(amplitude));
  355.             }
  356.         }

  358.     }

  360.     /** {@inheritDoc} */
  361.     @Override
  362.     public Vector3D displacement(final BodiesElements elements, final Frame earthFrame,
  363.                                  final Vector3D referencePoint) {

  365.         // allocate arrays for each species splines
  366.         final UnivariateFunction[] realZSpline      = new UnivariateFunction[mainTides.length];
  367.         final UnivariateFunction[] imaginaryZSpline = new UnivariateFunction[mainTides.length];
  368.         final UnivariateFunction[] realWSpline      = new UnivariateFunction[mainTides.length];
  369.         final UnivariateFunction[] imaginaryWSpline = new UnivariateFunction[mainTides.length];
  370.         final UnivariateFunction[] realSSpline      = new UnivariateFunction[mainTides.length];
  371.         final UnivariateFunction[] imaginarySSpline = new UnivariateFunction[mainTides.length];

  373.         // prepare splines for each species
  374.         for (int i = 0; i < mainTides.length; ++i) {

  376.             // compute current rates
  377.             final double[] rates = new double[mainTides[i].length];
  378.             for (int j = 0; j < rates.length; ++j) {
  379.                 rates[j] = mainTides[i][j].tide.getRate(elements);
  380.             }

  382.             // set up splines for the current rates
  383.             realZSpline[i]      = spline(rates, mainTides[i], d -> d.realZ);
  384.             imaginaryZSpline[i] = spline(rates, mainTides[i], d -> d.imaginaryZ);
  385.             realWSpline[i]      = spline(rates, mainTides[i], d -> d.realW);
  386.             imaginaryWSpline[i] = spline(rates, mainTides[i], d -> d.imaginaryW);
  387.             realSSpline[i]      = spline(rates, mainTides[i], d -> d.realS);
  388.             imaginarySSpline[i] = spline(rates, mainTides[i], d -> d.imaginaryS);

  390.         }

  392.         // evaluate all harmonics using interpolated admittances
  393.         double dz = 0;
  394.         double dw = 0;
  395.         double ds = 0;
  396.         for (final Map.Entry<Tide, Double> entry : CARTWRIGHT_EDDEN_AMPLITUDE_MAP.entrySet()) {

  398.             final Tide   tide      = entry.getKey();
  399.             final double amplitude = entry.getValue();
  400.             final int    i         = tide.getTauMultiplier();
  401.             final double rate      = tide.getRate(elements);

  403.             // apply splines for the current rate of this tide
  404.             final double rZ = realZSpline[i].value(rate);
  405.             final double iZ = imaginaryZSpline[i].value(rate);
  406.             final double rW = realWSpline[i].value(rate);
  407.             final double iW = imaginaryWSpline[i].value(rate);
  408.             final double rS = realSSpline[i].value(rate);
  409.             final double iS = imaginarySSpline[i].value(rate);

  411.             // phase for the current tide, including corrections
  412.             final double correction;
  413.             if (tide.getTauMultiplier() == 0) {
  414.                 correction = FastMath.PI;
  415.             } else if (tide.getTauMultiplier() == 1) {
  416.                 correction = 0.5 * FastMath.PI;
  417.             } else {
  418.                 correction = 0.0;
  419.             }
  420.             final double phase = tide.getPhase(elements) + correction;

  422.             dz += amplitude * FastMath.hypot(rZ, iZ) * FastMath.cos(phase + FastMath.atan2(iZ, rZ));
  423.             dw += amplitude * FastMath.hypot(rW, iW) * FastMath.cos(phase + FastMath.atan2(iW, rW));
  424.             ds += amplitude * FastMath.hypot(rS, iS) * FastMath.cos(phase + FastMath.atan2(iS, rS));

  426.         }

  428.         // convert to proper frame
  429.         final GeodeticPoint gp = earth.transform(referencePoint, earthFrame, elements.getDate());
  430.         return new Vector3D(dz, gp.getZenith(),
  431.                             dw, gp.getWest(),
  432.                             ds, gp.getSouth());

  434.     }

  436.     /** Get a spline function for interpolating between main tide data.
  437.      * @param rates rates for the tides species
  438.      * @param data data for the tides species
  439.      * @param selector data selector
  440.      * @return spline function for interpolating the selected data
  441.      */
  442.     private UnivariateFunction spline(final double[] rates, final MainTideData[] data,
  443.                                       final Function<MainTideData, Double> selector) {
  444.         final double[] y = new double[data.length];
  445.         for (int i = 0; i < y.length; ++i) {
  446.             y[i] = selector.apply(data[i]);
  447.         }
  448.         final PolynomialSplineFunction psf = new SplineInterpolator().interpolate(rates, y);

  450.         // as per HARDISP program EVAL subroutine, if spline evaluation is outside of range,
  451.         // the closest value is used. This occurs for example for long period tides.
  452.         // The main tides have rates 0.0821°/h, 0.5444°/h and 1.0980°/h. However,
  453.         // tide 55565 has rate 0.0022°/h, which is below the min rate and tide 75565 has
  454.         // rate 1.1002°/h, which is above max rate
  455.         final double[] knots          = psf.getKnots();
  456.         final double   minRate        = knots[0];
  457.         final double   valueAtMinRate = psf.value(minRate);
  458.         final double   maxRate        = knots[knots.length - 1];
  459.         final double   valueAtMaxRate = psf.value(maxRate);
  460.         return t -> (t < minRate) ? valueAtMinRate : (t > maxRate) ? valueAtMaxRate : psf.value(t);

  462.     }

  464.     /** Container for main tide data. */
  465.     private static class MainTideData {

  467.         /** Tide for which we have ocean loading coefficients. */
  468.         private final Tide tide;

  470.         /** Scaled real part of admittance along zenith axis. */
  471.         private final double realZ;

  473.         /** Scaled imaginary part of admittance along zenith axis. */
  474.         private final double imaginaryZ;

  476.         /** Scaled real part of admittance along west axis. */
  477.         private final double realW;

  479.         /** Scaled imaginary part of admittance along west axis. */
  480.         private final double imaginaryW;

  482.         /** Scaled real part of admittance along south axis. */
  483.         private final double realS;

  485.         /** Scaled imaginary part of admittance south axis. */
  486.         private final double imaginaryS;

  488.         /** Simple constructor.
  489.          * @param coefficients coefficients for the considered site
  490.          * @param i tide species
  491.          * @param j tide index in the species
  492.          * @param absAmplitude absolute value of the Cartwright-Edden amplitude of the tide
  493.          */
  494.         MainTideData(final OceanLoadingCoefficients coefficients, final int i, final int j, final double absAmplitude) {
  495.             // Sine and Cosine of difference angles
  496.             final SinCos scZenith = FastMath.sinCos(coefficients.getZenithPhase(i, j));
  497.             final SinCos scWest   = FastMath.sinCos(coefficients.getWestPhase(i, j));
  498.             final SinCos scSouth  = FastMath.sinCos(coefficients.getSouthPhase(i, j));
  499.             // Initialize attributes
  500.             tide       = coefficients.getTide(i, j);
  501.             realZ      = coefficients.getZenithAmplitude(i, j) * scZenith.cos() / absAmplitude;
  502.             imaginaryZ = coefficients.getZenithAmplitude(i, j) * scZenith.sin() / absAmplitude;
  503.             realW      = coefficients.getWestAmplitude(i, j)   * scWest.cos()   / absAmplitude;
  504.             imaginaryW = coefficients.getWestAmplitude(i, j)   * scWest.sin()   / absAmplitude;
  505.             realS      = coefficients.getSouthAmplitude(i, j)  * scSouth.cos()  / absAmplitude;
  506.             imaginaryS = coefficients.getSouthAmplitude(i, j)  * scSouth.sin()  / absAmplitude;
  507.         }

  509.     }

  511. }
Created with JaCoCo 0.8.12.202403310830