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

  19. import java.util.List;

  21. import org.hipparchus.util.FastMath;
  22. import org.hipparchus.util.Precision;
  23. import org.orekit.annotation.DefaultDataContext;
  24. import org.orekit.data.DataContext;
  25. import org.orekit.errors.OrekitException;
  26. import org.orekit.errors.OrekitMessages;
  27. import org.orekit.time.AbsoluteDate;

  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 FES Cnm, Snm tides files. */
  53.     public static final String FES_CNM_SNM_FILENAME = "^fes(\\d)+_Cnm-Snm.dat$";

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

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

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

  69.     /* Data loading methods. */

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

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

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

  107.     /** Clear gravity field readers.
  108.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  109.      * @see #addDefaultPotentialCoefficientsReaders()
  110.      */
  111.     @DefaultDataContext
  112.     public static void clearPotentialCoefficientsReaders() {
  113.         getGravityFields().clearPotentialCoefficientsReaders();
  114.     }

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

  126.     /** Configure ocean load deformation coefficients.
  127.      * @param oldc ocean load deformation coefficients
  128.      * @see #getOceanLoadDeformationCoefficients()
  129.      */
  130.     @DefaultDataContext
  131.     public static void configureOceanLoadDeformationCoefficients(final OceanLoadDeformationCoefficients oldc) {
  132.         getGravityFields().configureOceanLoadDeformationCoefficients(oldc);
  133.     }

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

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

  180.     /** Clear ocean tides readers.
  181.      * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
  182.      * @see #addDefaultPotentialCoefficientsReaders()
  183.      */
  184.     @DefaultDataContext
  185.     public static void clearOceanTidesReaders() {
  186.         getGravityFields().clearOceanTidesReaders();
  187.     }

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

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

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

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

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

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

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

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

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

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

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

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

  410.         // allocate a triangular array
  411.         final int rows = degree + 1;
  412.         final double[][] factor = new double[rows][];
  413.         factor[0] = new double[] {1.0};

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

  431.         return factor;

  433.     }

  435. }
Created with JaCoCo 0.8.11.202310140853