GravityFieldFactory.java

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

  18. import java.util.List;

  19. import org.hipparchus.util.FastMath;
  20. import org.hipparchus.util.Precision;
  21. import org.orekit.annotation.DefaultDataContext;
  22. import org.orekit.data.DataContext;
  23. import org.orekit.errors.OrekitException;
  24. import org.orekit.errors.OrekitMessages;

  25. /** Factory used to read gravity field files in several supported formats.
  26.  * @author Fabien Maussion
  27.  * @author Pascal Parraud
  28.  * @author Luc Maisonobe
  29.  */
  30. public class GravityFieldFactory {

  31.     /* These constants were left here instead of being moved to LazyLoadedGravityFields
  32.      * because they are public.
  33.      */

  34.     /** Default regular expression for ICGEM files. */
  35.     public static final String ICGEM_FILENAME = "^(.*\\.gfc)|(g(\\d)+_eigen[-_](\\w)+_coef)$";

  36.     /** Default regular expression for SHM files. */
  37.     public static final String SHM_FILENAME = "^eigen[-_](\\w)+_coef$";

  38.     /** Default regular expression for EGM files. */
  39.     public static final String EGM_FILENAME = "^egm\\d\\d_to\\d.*$";

  40.     /** Default regular expression for GRGS files. */
  41.     public static final String GRGS_FILENAME = "^grim\\d_.*$";

  42.     /** Default regular expression for FES Cnm, Snm tides files. */
  43.     public static final String FES_CNM_SNM_FILENAME = "^fes(\\d)+_Cnm-Snm.dat$";

  44.     /** Default regular expression for FES C hat and epsilon tides files. */
  45.     public static final String FES_CHAT_EPSILON_FILENAME = "^fes(\\d)+.dat$";

  46.     /** Default regular expression for FES Hf tides files. */
  47.     public static final String FES_HF_FILENAME = "^hf-fes(\\d)+.dat$";

  48.     /** Private constructor.
  49.      * <p>This class is a utility class, it should neither have a public
  50.      * nor a default constructor. This private constructor prevents
  51.      * the compiler from generating one automatically.</p>
  52.      */
  53.     private GravityFieldFactory() {
  54.     }

  55.     /* Data loading methods. */

  56.     /**
  57.      * Get the instance of {@link GravityFields} that is called by the static methods of
  58.      * this class.
  59.      *
  60.      * @return the gravity fields used by this factory.
  61.      * @since 10.1
  62.      */
  63.     @DefaultDataContext
  64.     public static LazyLoadedGravityFields getGravityFields() {
  65.         return DataContext.getDefault().getGravityFields();
  66.     }

  67.     /** Add a reader for gravity fields.
  68.      * @param reader custom reader to add for the gravity field
  69.      * @see #addDefaultPotentialCoefficientsReaders()
  70.      * @see #clearPotentialCoefficientsReaders()
  71.      */
  72.     @DefaultDataContext
  73.     public static void addPotentialCoefficientsReader(final PotentialCoefficientsReader reader) {
  74.         getGravityFields().addPotentialCoefficientsReader(reader);
  75.     }

  76.     /** Add the default readers for gravity fields.
  77.      * <p>
  78.      * The default READERS supports ICGEM, SHM, EGM and GRGS formats with the
  79.      * default names {@link #ICGEM_FILENAME}, {@link #SHM_FILENAME}, {@link
  80.      * #EGM_FILENAME}, {@link #GRGS_FILENAME} and don't allow missing coefficients.
  81.      * </p>
  82.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  83.      * @see #clearPotentialCoefficientsReaders()
  84.      */
  85.     @DefaultDataContext
  86.     public static void addDefaultPotentialCoefficientsReaders() {
  87.         getGravityFields().addDefaultPotentialCoefficientsReaders();
  88.     }

  89.     /** Clear gravity field readers.
  90.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  91.      * @see #addDefaultPotentialCoefficientsReaders()
  92.      */
  93.     @DefaultDataContext
  94.     public static void clearPotentialCoefficientsReaders() {
  95.         getGravityFields().clearPotentialCoefficientsReaders();
  96.     }

  97.     /** Add a reader for ocean tides.
  98.      * @param reader custom reader to add for the gravity field
  99.      * @see #addDefaultPotentialCoefficientsReaders()
  100.      * @see #clearPotentialCoefficientsReaders()
  101.      */
  102.     @DefaultDataContext
  103.     public static void addOceanTidesReader(final OceanTidesReader reader) {
  104.         getGravityFields().addOceanTidesReader(reader);
  105.     }

  106.     /** Configure ocean load deformation coefficients.
  107.      * @param oldc ocean load deformation coefficients
  108.      * @see #getOceanLoadDeformationCoefficients()
  109.      */
  110.     @DefaultDataContext
  111.     public static void configureOceanLoadDeformationCoefficients(final OceanLoadDeformationCoefficients oldc) {
  112.         getGravityFields().configureOceanLoadDeformationCoefficients(oldc);
  113.     }

  114.     /** Get the configured ocean load deformation coefficients.
  115.      * <p>
  116.      * If {@link #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  117.      * configureOceanLoadDeformationCoefficients} has never been called, the default
  118.      * value will be the {@link OceanLoadDeformationCoefficients#IERS_2010 IERS 2010}
  119.      * coefficients.
  120.      * </p>
  121.      * @return ocean load deformation coefficients
  122.      * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  123.      */
  124.     @DefaultDataContext
  125.     public static OceanLoadDeformationCoefficients getOceanLoadDeformationCoefficients() {
  126.         return getGravityFields().getOceanLoadDeformationCoefficients();
  127.     }

  128.     /** Add the default READERS for ocean tides.
  129.      * <p>
  130.      * The default READERS supports files similar to the fes2004_Cnm-Snm.dat and
  131.      * fes2004.dat as published by IERS, using the {@link
  132.      * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  133.      * configured} ocean load deformation coefficients, which by default are the
  134.      * IERS 2010 coefficients, which are limited to degree 6. If higher degree
  135.      * coefficients are needed, the {@link
  136.      * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  137.      * configureOceanLoadDeformationCoefficients} method can be called prior to
  138.      * loading the ocean tides model with the {@link
  139.      * OceanLoadDeformationCoefficients#GEGOUT high degree coefficients} computed
  140.      * by Pascal Gégout.
  141.      * </p>
  142.      * <p>
  143.      * WARNING: the files referenced in the published conventions have some errors.
  144.      * These errors have been corrected and the updated files can be found here:
  145.      * <a href="http://tai.bipm.org/iers/convupdt/convupdt_c6.html">
  146.      * http://tai.bipm.org/iers/convupdt/convupdt_c6.html</a>.
  147.      * </p>
  148.           * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  149.      * @see #clearPotentialCoefficientsReaders()
  150.      * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  151.      * @see #getOceanLoadDeformationCoefficients()
  152.      */
  153.     @DefaultDataContext
  154.     public static void addDefaultOceanTidesReaders() {
  155.         getGravityFields().addDefaultOceanTidesReaders();
  156.     }

  157.     /** Clear ocean tides readers.
  158.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  159.      * @see #addDefaultPotentialCoefficientsReaders()
  160.      */
  161.     @DefaultDataContext
  162.     public static void clearOceanTidesReaders() {
  163.         getGravityFields().clearOceanTidesReaders();
  164.     }

  165.     /** Get the constant gravity field coefficients provider from the first supported file.
  166.      * <p>
  167.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  168.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  169.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  170.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  171.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  172.      * method will be called automatically.
  173.      * </p>
  174.      * @param degree maximal degree
  175.      * @param order maximal order
  176.      * @return a gravity field coefficients provider containing already loaded data
  177.      * @since 6.0
  178.      * @see #getNormalizedProvider(int, int)
  179.      */
  180.     @DefaultDataContext
  181.     public static NormalizedSphericalHarmonicsProvider getConstantNormalizedProvider(final int degree,
  182.                                                                                      final int order) {
  183.         return getGravityFields().getConstantNormalizedProvider(degree, order);
  184.     }

  185.     /** Get the gravity field coefficients provider from the first supported file.
  186.      * <p>
  187.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  188.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  189.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  190.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  191.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  192.      * method will be called automatically.
  193.      * </p>
  194.      * @param degree maximal degree
  195.      * @param order maximal order
  196.      * @return a gravity field coefficients provider containing already loaded data
  197.      * @since 6.0
  198.      * @see #getConstantNormalizedProvider(int, int)
  199.      */
  200.     @DefaultDataContext
  201.     public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final int degree,
  202.                                                                              final int order) {
  203.         return getGravityFields().getNormalizedProvider(degree, order);
  204.     }

  205.     /** Get the constant gravity field coefficients provider from the first supported file.
  206.      * <p>
  207.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  208.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  209.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  210.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  211.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  212.      * method will be called automatically.
  213.      * </p>
  214.      * @param degree maximal degree
  215.      * @param order maximal order
  216.      * @return a gravity field coefficients provider containing already loaded data
  217.      * @since 6.0
  218.      * @see #getUnnormalizedProvider(int, int)
  219.      */
  220.     @DefaultDataContext
  221.     public static UnnormalizedSphericalHarmonicsProvider getConstantUnnormalizedProvider(final int degree,
  222.                                                                                          final int order) {
  223.         return getGravityFields().getConstantUnnormalizedProvider(degree, order);
  224.     }

  225.     /** Get the gravity field coefficients provider from the first supported file.
  226.      * <p>
  227.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  228.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  229.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  230.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  231.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  232.      * method will be called automatically.
  233.      * </p>
  234.      * @param degree maximal degree
  235.      * @param order maximal order
  236.      * @return a gravity field coefficients provider containing already loaded data
  237.      * @since 6.0
  238.      * @see #getConstantUnnormalizedProvider(int, int)
  239.      */
  240.     @DefaultDataContext
  241.     public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final int degree,
  242.                                                                                  final int order) {
  243.         return getGravityFields().getUnnormalizedProvider(degree, order);
  244.     }

  245.     /** Read a gravity field coefficients provider from the first supported file.
  246.      * <p>
  247.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  248.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  249.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  250.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  251.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  252.      * method will be called automatically.
  253.      * </p>
  254.      * @param maxParseDegree maximal degree to parse
  255.      * @param maxParseOrder maximal order to parse
  256.      * @return a reader containing already loaded data
  257.      * @since 6.0
  258.      */
  259.     @DefaultDataContext
  260.     public static PotentialCoefficientsReader readGravityField(final int maxParseDegree,
  261.                                                                final int maxParseOrder) {
  262.         return getGravityFields().readGravityField(maxParseDegree, maxParseOrder);
  263.     }

  264.     /** Get the ocean tides waves from the first supported file.
  265.      * <p>
  266.      * If no {@link OceanTidesReader} has been added by calling {@link
  267.      * #addOceanTidesReader(OceanTidesReader)
  268.      * addOceanTidesReader} or if {@link #clearOceanTidesReaders()
  269.      * clearOceanTidesReaders} has been called afterwards, the {@link
  270.      * #addDefaultOceanTidesReaders() addDefaultOceanTidesReaders}
  271.      * method will be called automatically.
  272.      * </p>
  273.      * <p><span style="color:red">
  274.      * WARNING: as of 2013-11-17, there seem to be an inconsistency when loading
  275.      * one or the other file, for wave Sa (Doodson number 56.554) and P1 (Doodson
  276.      * number 163.555). The sign of the coefficients are different. We think the
  277.      * problem lies in the input files from IERS and not in the conversion (which
  278.      * works for all other waves), but cannot be sure. For this reason, ocean
  279.      * tides are still considered experimental at this date.
  280.      * </span></p>
  281.      * @param degree maximal degree
  282.      * @param order maximal order
  283.      * @return list of tides waves containing already loaded data
  284.      * @since 6.1
  285.      */
  286.     @DefaultDataContext
  287.     public static List<OceanTidesWave> getOceanTidesWaves(final int degree, final int order) {
  288.         return getGravityFields().getOceanTidesWaves(degree, order);
  289.     }

  290.     /* static helper methods that don't load data. */

  291.     /** Create a time-independent {@link NormalizedSphericalHarmonicsProvider} from canonical coefficients.
  292.      * <p>
  293.      * Note that contrary to the other factory method, this one does not read any data, it simply uses
  294.      * the provided data
  295.      * </p>
  296.      * @param ae central body reference radius
  297.      * @param mu central body attraction coefficient
  298.      * @param tideSystem tide system
  299.      * @param normalizedC normalized tesseral-sectorial coefficients (cosine part)
  300.      * @param normalizedS normalized tesseral-sectorial coefficients (sine part)
  301.      * @return provider for normalized coefficients
  302.      * @since 6.0
  303.      */
  304.     public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final double ae, final double mu,
  305.                                                                              final TideSystem tideSystem,
  306.                                                                              final double[][] normalizedC,
  307.                                                                              final double[][] normalizedS) {
  308.         final Flattener flattener = new Flattener(normalizedC.length - 1, normalizedC[normalizedC.length - 1].length - 1);
  309.         final RawSphericalHarmonicsProvider constant =
  310.                         new ConstantSphericalHarmonics(ae, mu, tideSystem, flattener,
  311.                                                        flattener.flatten(normalizedC), flattener.flatten(normalizedS));
  312.         return new WrappingNormalizedProvider(constant);
  313.     }

  314.     /** Create a {@link NormalizedSphericalHarmonicsProvider} from an {@link UnnormalizedSphericalHarmonicsProvider}.
  315.      * <p>
  316.      * Note that contrary to the other factory method, this one does not read any data, it simply uses
  317.      * the provided data.
  318.      * </p>
  319.      * @param unnormalized provider to normalize
  320.      * @return provider for normalized coefficients
  321.      * @since 6.0
  322.      */
  323.     public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final UnnormalizedSphericalHarmonicsProvider unnormalized) {
  324.         return new Normalizer(unnormalized);
  325.     }

  326.     /** Create a time-independent {@link UnnormalizedSphericalHarmonicsProvider} from canonical coefficients.
  327.      * <p>
  328.      * Note that contrary to the other factory method, this one does not read any data, it simply uses
  329.      * the provided data
  330.      * </p>
  331.      * @param ae central body reference radius
  332.      * @param mu central body attraction coefficient
  333.      * @param tideSystem tide system
  334.      * @param unnormalizedC un-normalized tesseral-sectorial coefficients (cosine part)
  335.      * @param unnormalizedS un-normalized tesseral-sectorial coefficients (sine part)
  336.      * @return provider for un-normalized coefficients
  337.      * @since 6.0
  338.      */
  339.     public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final double ae, final double mu,
  340.                                                                                  final TideSystem tideSystem,
  341.                                                                                  final double[][] unnormalizedC,
  342.                                                                                  final double[][] unnormalizedS) {
  343.         final Flattener flattener = new Flattener(unnormalizedC.length - 1, unnormalizedC[unnormalizedC.length - 1].length - 1);
  344.         final RawSphericalHarmonicsProvider constant =
  345.                         new ConstantSphericalHarmonics(ae, mu, tideSystem, flattener,
  346.                                                        flattener.flatten(unnormalizedC), flattener.flatten(unnormalizedS));
  347.         return new WrappingUnnormalizedProvider(constant);
  348.     }

  349.     /** Create an {@link UnnormalizedSphericalHarmonicsProvider} from a {@link NormalizedSphericalHarmonicsProvider}.
  350.      * <p>
  351.      * Note that contrary to the other factory method, this one does not read any data, it simply uses
  352.      * the provided data.
  353.      * </p>
  354.      * @param normalized provider to un-normalize
  355.      * @return provider for un-normalized coefficients
  356.      * @since 6.0
  357.      */
  358.     public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final NormalizedSphericalHarmonicsProvider normalized) {
  359.         return new Unnormalizer(normalized);
  360.     }

  361.     /** Get a un-normalization factors array.
  362.      * <p>
  363.      * Un-normalized coefficients are obtained by multiplying normalized
  364.      * coefficients by the factors array elements.
  365.      * </p>
  366.      * @param degree maximal degree
  367.      * @param order maximal order
  368.      * @return triangular un-normalization factors array
  369.      * @since 6.0
  370.      */
  371.     public static double[][] getUnnormalizationFactors(final int degree, final int order) {

  372.         // allocate a triangular array
  373.         final int rows = degree + 1;
  374.         final double[][] factor = new double[rows][];
  375.         factor[0] = new double[] {1.0};

  376.         // compute the factors
  377.         for (int n = 1; n <= degree; n++) {
  378.             final double[] row = new double[FastMath.min(n, order) + 1];
  379.             row[0] = FastMath.sqrt(2 * n + 1);
  380.             double coeff = 2.0 * (2 * n + 1);
  381.             for (int m = 1; m < row.length; m++) {
  382.                 coeff /= (n - m + 1) * (n + m);
  383.                 row[m] = FastMath.sqrt(coeff);
  384.                 if (row[m] < Precision.SAFE_MIN) {
  385.                     throw new OrekitException(OrekitMessages.GRAVITY_FIELD_NORMALIZATION_UNDERFLOW,
  386.                             n, m);
  387.                 }
  388.             }
  389.             factor[n] = row;
  390.         }

  391.         return factor;

  392.     }

  393. }