GravityFieldFactory.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.potential;

  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. import org.orekit.time.AbsoluteDate;

  27. import java.util.List;

  29. /** Factory used to read gravity field files in several supported formats.
  30.  * @author Fabien Maussion
  31.  * @author Pascal Parraud
  32.  * @author Luc Maisonobe
  33.  */
  34. public class GravityFieldFactory {

  36.     /* These constants were left here instead of being moved to LazyLoadedGravityFields
  37.      * because they are public.
  38.      */

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

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

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

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

  52.     /** Default regular expression for SHA files. */
  53.     public static final String SHA_FILENAME = "^sha\\..*$";

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

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

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

  64.     /** Private constructor.
  65.      * <p>This class is a utility class, it should neither have a public
  66.      * nor a default constructor. This private constructor prevents
  67.      * the compiler from generating one automatically.</p>
  68.      */
  69.     private GravityFieldFactory() {
  70.     }

  72.     /* Data loading methods. */

  74.     /**
  75.      * Get the instance of {@link GravityFields} that is called by the static methods of
  76.      * this class.
  77.      *
  78.      * @return the gravity fields used by this factory.
  79.      * @since 10.1
  80.      */
  81.     @DefaultDataContext
  82.     public static LazyLoadedGravityFields getGravityFields() {
  83.         return DataContext.getDefault().getGravityFields();
  84.     }

  86.     /** Add a reader for gravity fields.
  87.      * @param reader custom reader to add for the gravity field
  88.      * @see #addDefaultPotentialCoefficientsReaders()
  89.      * @see #clearPotentialCoefficientsReaders()
  90.      */
  91.     @DefaultDataContext
  92.     public static void addPotentialCoefficientsReader(final PotentialCoefficientsReader reader) {
  93.         getGravityFields().addPotentialCoefficientsReader(reader);
  94.     }

  96.     /** Add the default readers for gravity fields.
  97.      * <p>
  98.      * The default READERS supports ICGEM, SHM, EGM, GRGS and SHA formats with the
  99.      * default names {@link #ICGEM_FILENAME}, {@link #SHM_FILENAME}, {@link
  100.      * #EGM_FILENAME}, {@link #GRGS_FILENAME}, {@link #SHA_FILENAME}
  101.      * and don't allow missing coefficients.
  102.      * </p>
  103.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  104.      * @see #clearPotentialCoefficientsReaders()
  105.      */
  106.     @DefaultDataContext
  107.     public static void addDefaultPotentialCoefficientsReaders() {
  108.         getGravityFields().addDefaultPotentialCoefficientsReaders();
  109.     }

  111.     /** Clear gravity field readers.
  112.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  113.      * @see #addDefaultPotentialCoefficientsReaders()
  114.      */
  115.     @DefaultDataContext
  116.     public static void clearPotentialCoefficientsReaders() {
  117.         getGravityFields().clearPotentialCoefficientsReaders();
  118.     }

  120.     /** Add a reader for ocean tides.
  121.      * @param reader custom reader to add for the gravity field
  122.      * @see #addDefaultPotentialCoefficientsReaders()
  123.      * @see #clearPotentialCoefficientsReaders()
  124.      */
  125.     @DefaultDataContext
  126.     public static void addOceanTidesReader(final OceanTidesReader reader) {
  127.         getGravityFields().addOceanTidesReader(reader);
  128.     }

  130.     /** Configure ocean load deformation coefficients.
  131.      * @param oldc ocean load deformation coefficients
  132.      * @see #getOceanLoadDeformationCoefficients()
  133.      */
  134.     @DefaultDataContext
  135.     public static void configureOceanLoadDeformationCoefficients(final OceanLoadDeformationCoefficients oldc) {
  136.         getGravityFields().configureOceanLoadDeformationCoefficients(oldc);
  137.     }

  139.     /** Get the configured ocean load deformation coefficients.
  140.      * <p>
  141.      * If {@link #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  142.      * configureOceanLoadDeformationCoefficients} has never been called, the default
  143.      * value will be the {@link OceanLoadDeformationCoefficients#IERS_2010 IERS 2010}
  144.      * coefficients.
  145.      * </p>
  146.      * @return ocean load deformation coefficients
  147.      * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  148.      */
  149.     @DefaultDataContext
  150.     public static OceanLoadDeformationCoefficients getOceanLoadDeformationCoefficients() {
  151.         return getGravityFields().getOceanLoadDeformationCoefficients();
  152.     }

  154.     /** Add the default READERS for ocean tides.
  155.      * <p>
  156.      * The default READERS supports files similar to the fes2004_Cnm-Snm.dat and
  157.      * fes2004.dat as published by IERS, using the {@link
  158.      * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  159.      * configured} ocean load deformation coefficients, which by default are the
  160.      * IERS 2010 coefficients, which are limited to degree 6. If higher degree
  161.      * coefficients are needed, the {@link
  162.      * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  163.      * configureOceanLoadDeformationCoefficients} method can be called prior to
  164.      * loading the ocean tides model with the {@link
  165.      * OceanLoadDeformationCoefficients#GEGOUT high degree coefficients} computed
  166.      * by Pascal Gégout.
  167.      * </p>
  168.      * <p>
  169.      * WARNING: the files referenced in the published conventions have some errors.
  170.      * These errors have been corrected and the updated files can be found here:
  171.      * <a href="http://tai.bipm.org/iers/convupdt/convupdt_c6.html">
  172.      * http://tai.bipm.org/iers/convupdt/convupdt_c6.html</a>.
  173.      * </p>
  174.           * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  175.      * @see #clearPotentialCoefficientsReaders()
  176.      * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
  177.      * @see #getOceanLoadDeformationCoefficients()
  178.      */
  179.     @DefaultDataContext
  180.     public static void addDefaultOceanTidesReaders() {
  181.         getGravityFields().addDefaultOceanTidesReaders();
  182.     }

  184.     /** Clear ocean tides readers.
  185.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  186.      * @see #addDefaultPotentialCoefficientsReaders()
  187.      */
  188.     @DefaultDataContext
  189.     public static void clearOceanTidesReaders() {
  190.         getGravityFields().clearOceanTidesReaders();
  191.     }

  193.     /** Get the constant gravity field coefficients provider from the first supported file.
  194.      * <p>
  195.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  196.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  197.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  198.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  199.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  200.      * method will be called automatically.
  201.      * </p>
  202.      * @param degree maximal degree
  203.      * @param order maximal order
  204.      * @param freezingDate freezing epoch
  205.      * @return a gravity field coefficients provider containing already loaded data
  206.      * @since 12.0
  207.      * @see #getNormalizedProvider(int, int)
  208.      */
  209.     @DefaultDataContext
  210.     public static NormalizedSphericalHarmonicsProvider getConstantNormalizedProvider(final int degree, final int order,
  211.                                                                                      final AbsoluteDate freezingDate) {
  212.         return getGravityFields().getConstantNormalizedProvider(degree, order, freezingDate);
  213.     }

  215.     /** Get the gravity field coefficients provider from the first supported file.
  216.      * <p>
  217.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  218.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  219.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  220.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  221.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  222.      * method will be called automatically.
  223.      * </p>
  224.      * @param degree maximal degree
  225.      * @param order maximal order
  226.      * @return a gravity field coefficients provider containing already loaded data
  227.      * @since 6.0
  228.      * @see #getConstantNormalizedProvider(int, int, AbsoluteDate)
  229.      */
  230.     @DefaultDataContext
  231.     public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final int degree,
  232.                                                                              final int order) {
  233.         return getGravityFields().getNormalizedProvider(degree, order);
  234.     }

  236.     /** Get the constant gravity field coefficients provider from the first supported file.
  237.      * <p>
  238.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  239.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  240.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  241.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  242.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  243.      * method will be called automatically.
  244.      * </p>
  245.      * @param degree maximal degree
  246.      * @param order maximal order
  247.      * @param freezingDate freezing epoch
  248.      * @return a gravity field coefficients provider containing already loaded data
  249.      * @since 6.0
  250.      * @see #getUnnormalizedProvider(int, int)
  251.      */
  252.     @DefaultDataContext
  253.     public static UnnormalizedSphericalHarmonicsProvider getConstantUnnormalizedProvider(final int degree, final int order,
  254.                                                                                          final AbsoluteDate freezingDate) {
  255.         return getGravityFields().getConstantUnnormalizedProvider(degree, order, freezingDate);
  256.     }

  258.     /** Get the gravity field coefficients provider from the first supported file.
  259.      * <p>
  260.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  261.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  262.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  263.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  264.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  265.      * method will be called automatically.
  266.      * </p>
  267.      * @param degree maximal degree
  268.      * @param order maximal order
  269.      * @return a gravity field coefficients provider containing already loaded data
  270.      * @since 6.0
  271.      * @see #getConstantUnnormalizedProvider(int, int, AbsoluteDate)
  272.      */
  273.     @DefaultDataContext
  274.     public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final int degree,
  275.                                                                                  final int order) {
  276.         return getGravityFields().getUnnormalizedProvider(degree, order);
  277.     }

  279.     /** Read a gravity field coefficients provider from the first supported file.
  280.      * <p>
  281.      * If no {@link PotentialCoefficientsReader} has been added by calling {@link
  282.      * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  283.      * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
  284.      * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
  285.      * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
  286.      * method will be called automatically.
  287.      * </p>
  288.      * @param maxParseDegree maximal degree to parse
  289.      * @param maxParseOrder maximal order to parse
  290.      * @return a reader containing already loaded data
  291.      * @since 6.0
  292.      */
  293.     @DefaultDataContext
  294.     public static PotentialCoefficientsReader readGravityField(final int maxParseDegree,
  295.                                                                final int maxParseOrder) {
  296.         return getGravityFields().readGravityField(maxParseDegree, maxParseOrder);
  297.     }

  299.     /** Get the ocean tides waves from the first supported file.
  300.      * <p>
  301.      * If no {@link OceanTidesReader} has been added by calling {@link
  302.      * #addOceanTidesReader(OceanTidesReader)
  303.      * addOceanTidesReader} or if {@link #clearOceanTidesReaders()
  304.      * clearOceanTidesReaders} has been called afterwards, the {@link
  305.      * #addDefaultOceanTidesReaders() addDefaultOceanTidesReaders}
  306.      * method will be called automatically.
  307.      * </p>
  308.      * <p><span style="color:red">
  309.      * WARNING: as of 2013-11-17, there seem to be an inconsistency when loading
  310.      * one or the other file, for wave Sa (Doodson number 56.554) and P1 (Doodson
  311.      * number 163.555). The sign of the coefficients are different. We think the
  312.      * problem lies in the input files from IERS and not in the conversion (which
  313.      * works for all other waves), but cannot be sure. For this reason, ocean
  314.      * tides are still considered experimental at this date.
  315.      * </span></p>
  316.      * @param degree maximal degree
  317.      * @param order maximal order
  318.      * @return list of tides waves containing already loaded data
  319.      * @since 6.1
  320.      */
  321.     @DefaultDataContext
  322.     public static List<OceanTidesWave> getOceanTidesWaves(final int degree, final int order) {
  323.         return getGravityFields().getOceanTidesWaves(degree, order);
  324.     }

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

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

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

  365.     /** Create a time-independent {@link UnnormalizedSphericalHarmonicsProvider} from canonical coefficients.
  366.      * <p>
  367.      * Note that contrary to the other factory method, this one does not read any data, it simply uses
  368.      * the provided data
  369.      * </p>
  370.      * @param ae central body reference radius
  371.      * @param mu central body attraction coefficient
  372.      * @param tideSystem tide system
  373.      * @param unnormalizedC un-normalized tesseral-sectorial coefficients (cosine part)
  374.      * @param unnormalizedS un-normalized tesseral-sectorial coefficients (sine part)
  375.      * @return provider for un-normalized coefficients
  376.      * @since 6.0
  377.      */
  378.     public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final double ae, final double mu,
  379.                                                                                  final TideSystem tideSystem,
  380.                                                                                  final double[][] unnormalizedC,
  381.                                                                                  final double[][] unnormalizedS) {
  382.         final Flattener flattener = new Flattener(unnormalizedC.length - 1, unnormalizedC[unnormalizedC.length - 1].length - 1);
  383.         final RawSphericalHarmonicsProvider constant =
  384.                         new ConstantSphericalHarmonics(ae, mu, tideSystem, flattener,
  385.                                                        flattener.flatten(unnormalizedC), flattener.flatten(unnormalizedS));
  386.         return new WrappingUnnormalizedProvider(constant);
  387.     }

  389.     /** Create an {@link UnnormalizedSphericalHarmonicsProvider} from a {@link NormalizedSphericalHarmonicsProvider}.
  390.      * <p>
  391.      * Note that contrary to the other factory method, this one does not read any data, it simply uses
  392.      * the provided data.
  393.      * </p>
  394.      * @param normalized provider to un-normalize
  395.      * @return provider for un-normalized coefficients
  396.      * @since 6.0
  397.      */
  398.     public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final NormalizedSphericalHarmonicsProvider normalized) {
  399.         return new Unnormalizer(normalized);
  400.     }

  402.     /** Get a un-normalization factors array.
  403.      * <p>
  404.      * Un-normalized coefficients are obtained by multiplying normalized
  405.      * coefficients by the factors array elements.
  406.      * </p>
  407.      * @param degree maximal degree
  408.      * @param order maximal order
  409.      * @return triangular un-normalization factors array
  410.      * @since 6.0
  411.      */
  412.     public static double[][] getUnnormalizationFactors(final int degree, final int order) {

  414.         // allocate a triangular array
  415.         final int rows = degree + 1;
  416.         final double[][] factor = new double[rows][];
  417.         factor[0] = new double[] {1.0};

  419.         // compute the factors
  420.         for (int n = 1; n <= degree; n++) {
  421.             final double[] row = new double[FastMath.min(n, order) + 1];
  422.             row[0] = FastMath.sqrt(2 * n + 1);
  423.             double coeff = 2.0 * (2 * n + 1);
  424.             for (int m = 1; m < row.length; m++) {
  425.                 coeff /= (n - m + 1) * (n + m);
  426.                 row[m] = FastMath.sqrt(coeff);
  427.                 if (row[m] < Precision.SAFE_MIN) {
  428.                     throw new OrekitException(OrekitMessages.GRAVITY_FIELD_NORMALIZATION_UNDERFLOW,
  429.                             n, m);
  430.                 }
  431.             }
  432.             factor[n] = row;
  433.         }

  435.         return factor;

  437.     }

  439. }
Created with JaCoCo 0.8.12.202403310830