GravityFieldFactory.java

/* Copyright 2002-2023 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.forces.gravity.potential;

import java.util.List;

import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.time.AbsoluteDate;

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

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

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

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

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

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

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

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

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

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

    /* Data loading methods. */

    /**
     * Get the instance of {@link GravityFields} that is called by the static methods of
     * this class.
     *
     * @return the gravity fields used by this factory.
     * @since 10.1
     */
    @DefaultDataContext
    public static LazyLoadedGravityFields getGravityFields() {
        return DataContext.getDefault().getGravityFields();
    }

    /** Add a reader for gravity fields.
     * @param reader custom reader to add for the gravity field
     * @see #addDefaultPotentialCoefficientsReaders()
     * @see #clearPotentialCoefficientsReaders()
     */
    @DefaultDataContext
    public static void addPotentialCoefficientsReader(final PotentialCoefficientsReader reader) {
        getGravityFields().addPotentialCoefficientsReader(reader);
    }

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

    /** Clear gravity field readers.
     * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * @see #addDefaultPotentialCoefficientsReaders()
     */
    @DefaultDataContext
    public static void clearPotentialCoefficientsReaders() {
        getGravityFields().clearPotentialCoefficientsReaders();
    }

    /** Add a reader for ocean tides.
     * @param reader custom reader to add for the gravity field
     * @see #addDefaultPotentialCoefficientsReaders()
     * @see #clearPotentialCoefficientsReaders()
     */
    @DefaultDataContext
    public static void addOceanTidesReader(final OceanTidesReader reader) {
        getGravityFields().addOceanTidesReader(reader);
    }

    /** Configure ocean load deformation coefficients.
     * @param oldc ocean load deformation coefficients
     * @see #getOceanLoadDeformationCoefficients()
     */
    @DefaultDataContext
    public static void configureOceanLoadDeformationCoefficients(final OceanLoadDeformationCoefficients oldc) {
        getGravityFields().configureOceanLoadDeformationCoefficients(oldc);
    }

    /** Get the configured ocean load deformation coefficients.
     * <p>
     * If {@link #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
     * configureOceanLoadDeformationCoefficients} has never been called, the default
     * value will be the {@link OceanLoadDeformationCoefficients#IERS_2010 IERS 2010}
     * coefficients.
     * </p>
     * @return ocean load deformation coefficients
     * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
     */
    @DefaultDataContext
    public static OceanLoadDeformationCoefficients getOceanLoadDeformationCoefficients() {
        return getGravityFields().getOceanLoadDeformationCoefficients();
    }

    /** Add the default READERS for ocean tides.
     * <p>
     * The default READERS supports files similar to the fes2004_Cnm-Snm.dat and
     * fes2004.dat as published by IERS, using the {@link
     * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
     * configured} ocean load deformation coefficients, which by default are the
     * IERS 2010 coefficients, which are limited to degree 6. If higher degree
     * coefficients are needed, the {@link
     * #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
     * configureOceanLoadDeformationCoefficients} method can be called prior to
     * loading the ocean tides model with the {@link
     * OceanLoadDeformationCoefficients#GEGOUT high degree coefficients} computed
     * by Pascal Gégout.
     * </p>
     * <p>
     * WARNING: the files referenced in the published conventions have some errors.
     * These errors have been corrected and the updated files can be found here:
     * <a href="http://tai.bipm.org/iers/convupdt/convupdt_c6.html">
     * http://tai.bipm.org/iers/convupdt/convupdt_c6.html</a>.
     * </p>
          * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * @see #clearPotentialCoefficientsReaders()
     * @see #configureOceanLoadDeformationCoefficients(OceanLoadDeformationCoefficients)
     * @see #getOceanLoadDeformationCoefficients()
     */
    @DefaultDataContext
    public static void addDefaultOceanTidesReaders() {
        getGravityFields().addDefaultOceanTidesReaders();
    }

    /** Clear ocean tides readers.
     * @see #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * @see #addDefaultPotentialCoefficientsReaders()
     */
    @DefaultDataContext
    public static void clearOceanTidesReaders() {
        getGravityFields().clearOceanTidesReaders();
    }

    /** Get the constant gravity field coefficients provider from the first supported file.
     * <p>
     * If no {@link PotentialCoefficientsReader} has been added by calling {@link
     * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
     * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
     * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
     * method will be called automatically.
     * </p>
     * @param degree maximal degree
     * @param order maximal order
     * @param freezingDate freezing epoch
     * @return a gravity field coefficients provider containing already loaded data
     * @since 12.0
     * @see #getNormalizedProvider(int, int)
     */
    @DefaultDataContext
    public static NormalizedSphericalHarmonicsProvider getConstantNormalizedProvider(final int degree, final int order,
                                                                                     final AbsoluteDate freezingDate) {
        return getGravityFields().getConstantNormalizedProvider(degree, order, freezingDate);
    }

    /** Get the gravity field coefficients provider from the first supported file.
     * <p>
     * If no {@link PotentialCoefficientsReader} has been added by calling {@link
     * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
     * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
     * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
     * method will be called automatically.
     * </p>
     * @param degree maximal degree
     * @param order maximal order
     * @return a gravity field coefficients provider containing already loaded data
     * @since 6.0
     * @see #getConstantNormalizedProvider(int, int, AbsoluteDate)
     */
    @DefaultDataContext
    public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final int degree,
                                                                             final int order) {
        return getGravityFields().getNormalizedProvider(degree, order);
    }

    /** Get the constant gravity field coefficients provider from the first supported file.
     * <p>
     * If no {@link PotentialCoefficientsReader} has been added by calling {@link
     * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
     * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
     * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
     * method will be called automatically.
     * </p>
     * @param degree maximal degree
     * @param order maximal order
     * @param freezingDate freezing epoch
     * @return a gravity field coefficients provider containing already loaded data
     * @since 6.0
     * @see #getUnnormalizedProvider(int, int)
     */
    @DefaultDataContext
    public static UnnormalizedSphericalHarmonicsProvider getConstantUnnormalizedProvider(final int degree, final int order,
                                                                                         final AbsoluteDate freezingDate) {
        return getGravityFields().getConstantUnnormalizedProvider(degree, order, freezingDate);
    }

    /** Get the gravity field coefficients provider from the first supported file.
     * <p>
     * If no {@link PotentialCoefficientsReader} has been added by calling {@link
     * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
     * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
     * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
     * method will be called automatically.
     * </p>
     * @param degree maximal degree
     * @param order maximal order
     * @return a gravity field coefficients provider containing already loaded data
     * @since 6.0
     * @see #getConstantUnnormalizedProvider(int, int, AbsoluteDate)
     */
    @DefaultDataContext
    public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final int degree,
                                                                                 final int order) {
        return getGravityFields().getUnnormalizedProvider(degree, order);
    }

    /** Read a gravity field coefficients provider from the first supported file.
     * <p>
     * If no {@link PotentialCoefficientsReader} has been added by calling {@link
     * #addPotentialCoefficientsReader(PotentialCoefficientsReader)
     * addPotentialCoefficientsReader} or if {@link #clearPotentialCoefficientsReaders()
     * clearPotentialCoefficientsReaders} has been called afterwards, the {@link
     * #addDefaultPotentialCoefficientsReaders() addDefaultPotentialCoefficientsReaders}
     * method will be called automatically.
     * </p>
     * @param maxParseDegree maximal degree to parse
     * @param maxParseOrder maximal order to parse
     * @return a reader containing already loaded data
     * @since 6.0
     */
    @DefaultDataContext
    public static PotentialCoefficientsReader readGravityField(final int maxParseDegree,
                                                               final int maxParseOrder) {
        return getGravityFields().readGravityField(maxParseDegree, maxParseOrder);
    }

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

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

    /** Create a time-independent {@link NormalizedSphericalHarmonicsProvider} from canonical coefficients.
     * <p>
     * Note that contrary to the other factory method, this one does not read any data, it simply uses
     * the provided data
     * </p>
     * @param ae central body reference radius
     * @param mu central body attraction coefficient
     * @param tideSystem tide system
     * @param normalizedC normalized tesseral-sectorial coefficients (cosine part)
     * @param normalizedS normalized tesseral-sectorial coefficients (sine part)
     * @return provider for normalized coefficients
     * @since 6.0
     */
    public static NormalizedSphericalHarmonicsProvider getNormalizedProvider(final double ae, final double mu,
                                                                             final TideSystem tideSystem,
                                                                             final double[][] normalizedC,
                                                                             final double[][] normalizedS) {
        final Flattener flattener = new Flattener(normalizedC.length - 1, normalizedC[normalizedC.length - 1].length - 1);
        final RawSphericalHarmonicsProvider constant =
                        new ConstantSphericalHarmonics(ae, mu, tideSystem, flattener,
                                                       flattener.flatten(normalizedC), flattener.flatten(normalizedS));
        return new WrappingNormalizedProvider(constant);
    }

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

    /** Create a time-independent {@link UnnormalizedSphericalHarmonicsProvider} from canonical coefficients.
     * <p>
     * Note that contrary to the other factory method, this one does not read any data, it simply uses
     * the provided data
     * </p>
     * @param ae central body reference radius
     * @param mu central body attraction coefficient
     * @param tideSystem tide system
     * @param unnormalizedC un-normalized tesseral-sectorial coefficients (cosine part)
     * @param unnormalizedS un-normalized tesseral-sectorial coefficients (sine part)
     * @return provider for un-normalized coefficients
     * @since 6.0
     */
    public static UnnormalizedSphericalHarmonicsProvider getUnnormalizedProvider(final double ae, final double mu,
                                                                                 final TideSystem tideSystem,
                                                                                 final double[][] unnormalizedC,
                                                                                 final double[][] unnormalizedS) {
        final Flattener flattener = new Flattener(unnormalizedC.length - 1, unnormalizedC[unnormalizedC.length - 1].length - 1);
        final RawSphericalHarmonicsProvider constant =
                        new ConstantSphericalHarmonics(ae, mu, tideSystem, flattener,
                                                       flattener.flatten(unnormalizedC), flattener.flatten(unnormalizedS));
        return new WrappingUnnormalizedProvider(constant);
    }

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

    /** Get a un-normalization factors array.
     * <p>
     * Un-normalized coefficients are obtained by multiplying normalized
     * coefficients by the factors array elements.
     * </p>
     * @param degree maximal degree
     * @param order maximal order
     * @return triangular un-normalization factors array
     * @since 6.0
     */
    public static double[][] getUnnormalizationFactors(final int degree, final int order) {

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

        // compute the factors
        for (int n = 1; n <= degree; n++) {
            final double[] row = new double[FastMath.min(n, order) + 1];
            row[0] = FastMath.sqrt(2 * n + 1);
            double coeff = 2.0 * (2 * n + 1);
            for (int m = 1; m < row.length; m++) {
                coeff /= (n - m + 1) * (n + m);
                row[m] = FastMath.sqrt(coeff);
                if (row[m] < Precision.SAFE_MIN) {
                    throw new OrekitException(OrekitMessages.GRAVITY_FIELD_NORMALIZATION_UNDERFLOW,
                            n, m);
                }
            }
            factor[n] = row;
        }

        return factor;

    }

}