PotentialCoefficientsReader.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.io.IOException;
import java.io.InputStream;
import java.text.ParseException;
import java.util.Arrays;
import java.util.Locale;

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

/**This abstract class represents a Gravitational Potential Coefficients file reader.
 *
 * <p> As it exits many different coefficients models and containers this
 *  interface represents all the methods that should be implemented by a reader.
 *  The proper way to use this interface is to call the {@link GravityFieldFactory}
 *  which will determine which reader to use with the selected potential
 *  coefficients file.</p>
 *
 * @see GravityFields
 * @author Fabien Maussion
 */
public abstract class PotentialCoefficientsReader implements DataLoader {

    /** Maximal degree to parse. */
    private int maxParseDegree;

    /** Maximal order to parse. */
    private int maxParseOrder;

    /** Regular expression for supported files names. */
    private final String supportedNames;

    /** Allow missing coefficients in the input data. */
    private final boolean missingCoefficientsAllowed;

    /** Time scale for parsing dates. */
    private final TimeScale timeScale;

    /** Indicator for complete read. */
    private boolean readComplete;

    /** Central body reference radius. */
    private double ae;

    /** Central body attraction coefficient. */
    private double mu;

    /** Converter from triangular to flat form. */
    private Flattener flattener;

    /** Raw tesseral-sectorial coefficients matrix. */
    private double[] rawC;

    /** Raw tesseral-sectorial coefficients matrix. */
    private double[] rawS;

    /** Indicator for normalized raw coefficients. */
    private boolean normalized;

    /** Tide system. */
    private TideSystem tideSystem;

    /** Simple constructor.
     * <p>Build an uninitialized reader.</p>
     *
     * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
     *
     * @param supportedNames regular expression for supported files names
     * @param missingCoefficientsAllowed allow missing coefficients in the input data
     * @see #PotentialCoefficientsReader(String, boolean, TimeScale)
     */
    @DefaultDataContext
    protected PotentialCoefficientsReader(final String supportedNames,
                                          final boolean missingCoefficientsAllowed) {
        this(supportedNames, missingCoefficientsAllowed,
                DataContext.getDefault().getTimeScales().getTT());
    }

    /** Simple constructor.
     * <p>Build an uninitialized reader.</p>
     * @param supportedNames regular expression for supported files names
     * @param missingCoefficientsAllowed allow missing coefficients in the input data
     * @param timeScale to use when parsing dates.
     * @since 10.1
     */
    protected PotentialCoefficientsReader(final String supportedNames,
                                          final boolean missingCoefficientsAllowed,
                                          final TimeScale timeScale) {
        this.supportedNames             = supportedNames;
        this.missingCoefficientsAllowed = missingCoefficientsAllowed;
        this.maxParseDegree             = Integer.MAX_VALUE;
        this.maxParseOrder              = Integer.MAX_VALUE;
        this.readComplete               = false;
        this.ae                         = Double.NaN;
        this.mu                         = Double.NaN;
        this.flattener                  = null;
        this.rawC                       = null;
        this.rawS                       = null;
        this.normalized                 = false;
        this.tideSystem                 = TideSystem.UNKNOWN;
        this.timeScale                  = timeScale;
    }

    /** Get the regular expression for supported files names.
     * @return regular expression for supported files names
     */
    public String getSupportedNames() {
        return supportedNames;
    }

    /** Check if missing coefficients are allowed in the input data.
     * @return true if missing coefficients are allowed in the input data
     */
    public boolean missingCoefficientsAllowed() {
        return missingCoefficientsAllowed;
    }

    /** Set the degree limit for the next file parsing.
     * @param maxParseDegree maximal degree to parse (may be safely
     * set to {@link Integer#MAX_VALUE} to parse all available coefficients)
     * @since 6.0
     */
    public void setMaxParseDegree(final int maxParseDegree) {
        this.maxParseDegree = maxParseDegree;
    }

    /** Get the degree limit for the next file parsing.
     * @return degree limit for the next file parsing
     * @since 6.0
     */
    public int getMaxParseDegree() {
        return maxParseDegree;
    }

    /** Set the order limit for the next file parsing.
     * @param maxParseOrder maximal order to parse (may be safely
     * set to {@link Integer#MAX_VALUE} to parse all available coefficients)
     * @since 6.0
     */
    public void setMaxParseOrder(final int maxParseOrder) {
        this.maxParseOrder = maxParseOrder;
    }

    /** Get the order limit for the next file parsing.
     * @return order limit for the next file parsing
     * @since 6.0
     */
    public int getMaxParseOrder() {
        return maxParseOrder;
    }

    /** {@inheritDoc} */
    public boolean stillAcceptsData() {
        return !(readComplete &&
                 getMaxAvailableDegree() >= getMaxParseDegree() &&
                 getMaxAvailableOrder()  >= getMaxParseOrder());
    }

    /** Set the indicator for completed read.
     * @param readComplete if true, a gravity field has been completely read
     */
    protected void setReadComplete(final boolean readComplete) {
        this.readComplete = readComplete;
    }

    /** Set the central body reference radius.
     * @param ae central body reference radius
     */
    protected void setAe(final double ae) {
        this.ae = ae;
    }

    /** Get the central body reference radius.
     * @return central body reference radius
     */
    protected double getAe() {
        return ae;
    }

    /** Set the central body attraction coefficient.
     * @param mu central body attraction coefficient
     */
    protected void setMu(final double mu) {
        this.mu = mu;
    }

    /** Get the central body attraction coefficient.
     * @return central body attraction coefficient
     */
    protected double getMu() {
        return mu;
    }

    /** Set the {@link TideSystem} used in the gravity field.
     * @param tideSystem tide system used in the gravity field
     */
    protected void setTideSystem(final TideSystem tideSystem) {
        this.tideSystem = tideSystem;
    }

    /** Get the {@link TideSystem} used in the gravity field.
     * @return tide system used in the gravity field
     */
    protected TideSystem getTideSystem() {
        return tideSystem;
    }

    /** Set the tesseral-sectorial coefficients matrix.
     * @param rawNormalized if true, raw coefficients are normalized
     * @param f converter from triangular to flat form
     * @param c raw tesseral-sectorial coefficients matrix
     * @param s raw tesseral-sectorial coefficients matrix
     * @param name name of the file (or zip entry)
     * @since 11.1
     */
    protected void setRawCoefficients(final boolean rawNormalized, final Flattener f,
                                      final double[] c, final double[] s,
                                      final String name) {

        this.flattener = f;

        // normalization indicator
        normalized = rawNormalized;

        // set known constant values, if they were not defined in the file.
        // See Hofmann-Wellenhof and Moritz, "Physical Geodesy",
        // section 2.6 Harmonics of Lower Degree.
        // All S_i,0 are irrelevant because they are multiplied by zero.
        // C0,0 is 1, the central part, since all coefficients are normalized by GM.
        setIfUnset(c, flattener.index(0, 0), 1);
        setIfUnset(s, flattener.index(0, 0), 0);

        if (flattener.getDegree() >= 1) {
            // C1,0, C1,1, and S1,1 are the x,y,z coordinates of the center of mass,
            // which are 0 since all coefficients are given in an Earth centered frame
            setIfUnset(c, flattener.index(1, 0), 0);
            setIfUnset(s, flattener.index(1, 0), 0);
            if (flattener.getOrder() >= 1) {
                setIfUnset(c, flattener.index(1, 1), 0);
                setIfUnset(s, flattener.index(1, 1), 0);
            }
        }

        // cosine part
        for (int i = 0; i <= flattener.getDegree(); ++i) {
            for (int j = 0; j <= FastMath.min(i, flattener.getOrder()); ++j) {
                if (Double.isNaN(c[flattener.index(i, j)])) {
                    throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
                                              'C', i, j, name);
                }
            }
        }
        rawC = c.clone();

        // sine part
        for (int i = 0; i <= flattener.getDegree(); ++i) {
            for (int j = 0; j <= FastMath.min(i, flattener.getOrder()); ++j) {
                if (Double.isNaN(s[flattener.index(i, j)])) {
                    throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
                                              'S', i, j, name);
                }
            }
        }
        rawS = s.clone();

    }

    /**
     * Set a coefficient if it has not been set already.
     * <p>
     * If {@code array[i]} is 0 or NaN this method sets it to {@code value} and returns
     * {@code true}. Otherwise the original value of {@code array[i]} is preserved and
     * this method return {@code false}.
     * <p>
     * If {@code array[i]} does not exist then this method returns {@code false}.
     *
     * @param array the coefficient array.
     * @param i     index in array.
     * @param value the new value to set.
     * @return {@code true} if the coefficient was set to {@code value}, {@code false} if
     * the coefficient was not set to {@code value}. A {@code false} return indicates the
     * coefficient has previously been set to a non-NaN, non-zero value.
     */
    private boolean setIfUnset(final double[] array, final int i, final double value) {
        if (array.length > i && (Double.isNaN(array[i]) || Precision.equals(array[i], 0.0, 0))) {
            // the coefficient was not already initialized
            array[i] = value;
            return true;
        } else {
            return false;
        }
    }

    /** Get the maximal degree available in the last file parsed.
     * @return maximal degree available in the last file parsed
     * @since 6.0
     */
    public int getMaxAvailableDegree() {
        return flattener.getDegree();
    }

    /** Get the maximal order available in the last file parsed.
     * @return maximal order available in the last file parsed
     * @since 6.0
     */
    public int getMaxAvailableOrder() {
        return flattener.getOrder();
    }

    /** {@inheritDoc} */
    public abstract void loadData(InputStream input, String name)
        throws IOException, ParseException, OrekitException;

    /** Get a provider for read spherical harmonics coefficients.
     * @param wantNormalized if true, the provider will provide normalized coefficients,
     * otherwise it will provide un-normalized coefficients
     * @param degree maximal degree
     * @param order maximal order
     * @return a new provider
     * @since 6.0
     */
    public abstract RawSphericalHarmonicsProvider getProvider(boolean wantNormalized, int degree, int order);

    /** Get a time-independent provider containing base harmonics coefficients.
     * <p>
     * Beware that some coeefficients may be missing here, if they are managed as time-dependent
     * piecewise models (as in ICGEM V2.0).
     * </p>
     * @param wantNormalized if true, the raw provider must provide normalized coefficients,
     * otherwise it will provide un-normalized coefficients
     * @param degree maximal degree
     * @param order maximal order
     * @return a new provider, with no time-dependent parts
     * @see #getProvider(boolean, int, int)
     * @since 11.1
     */
    protected ConstantSphericalHarmonics getBaseProvider(final boolean wantNormalized,
                                                         final int degree, final int order) {

        if (!readComplete) {
            throw new OrekitException(OrekitMessages.NO_GRAVITY_FIELD_DATA_LOADED);
        }

        final Flattener truncatedFlattener = new Flattener(degree, order);
        return new ConstantSphericalHarmonics(ae, mu, tideSystem, truncatedFlattener,
                                              rescale(1.0, wantNormalized, truncatedFlattener, flattener, rawC),
                                              rescale(1.0, wantNormalized, truncatedFlattener, flattener, rawS));

    }

    /** Build a coefficients array in flat form.
     * @param flattener converter from triangular to flat form
     * @param value initial value to put in array elements
     * @return built array
     * @since 11.1
     */
    protected static double[] buildFlatArray(final Flattener flattener, final double value) {
        final double[] array = new double[flattener.arraySize()];
        Arrays.fill(array, value);
        return array;
    }

    /**
     * Parse a double from a string. Accept the Fortran convention of using a 'D' or
     * 'd' instead of an 'E' or 'e'.
     *
     * @param string to be parsed.
     * @return the double value of {@code string}.
     */
    protected static double parseDouble(final String string) {
        return Double.parseDouble(string.toUpperCase(Locale.ENGLISH).replace('D', 'E'));
    }

    /** Build a coefficients row.
     * @param degree row degree
     * @param order row order
     * @param value initial value to put in array elements
     * @return built row
     */
    protected static double[] buildRow(final int degree, final int order, final double value) {
        final double[] row = new double[FastMath.min(order, degree) + 1];
        Arrays.fill(row, value);
        return row;
    }

    /** Parse a coefficient.
     * @param field text field to parse
     * @param f converter from triangular to flat form
     * @param array array where to put the coefficient
     * @param i first index in the list
     * @param j second index in the list
     * @param cName name of the coefficient
     * @param name name of the file
     * @since 11.1
     */
    protected void parseCoefficient(final String field, final Flattener f,
                                    final double[] array, final int i, final int j,
                                    final String cName, final String name) {
        final int    index    = f.index(i, j);
        final double value    = parseDouble(field);
        final double oldValue = array[index];
        if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 0)) {
            // the coefficient was not already initialized
            array[index] = value;
        } else {
            throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
                                      name, i, j, name);
        }
    }

    /** Rescale coefficients arrays.
     * <p>
     * The normalized/unnormalized nature of original coefficients is inherited from previous parsing.
     * </p>
     * @param scale general scaling factor to apply to all elements
     * @param wantNormalized if true, the rescaled coefficients must be normalized,
     * otherwise they must be un-normalized
     * @param rescaledFlattener converter from triangular to flat form
     * @param originalFlattener converter from triangular to flat form
     * @param original original coefficients
     * @return rescaled coefficients
     * @since 11.1
     */
    protected double[] rescale(final double scale, final boolean wantNormalized, final Flattener rescaledFlattener,
                               final Flattener originalFlattener, final double[] original) {

        if (rescaledFlattener.getDegree() > originalFlattener.getDegree()) {
            throw new OrekitException(OrekitMessages.TOO_LARGE_DEGREE_FOR_GRAVITY_FIELD,
                                      rescaledFlattener.getDegree(), flattener.getDegree());
        }

        if (rescaledFlattener.getOrder() > originalFlattener.getOrder()) {
            throw new OrekitException(OrekitMessages.TOO_LARGE_ORDER_FOR_GRAVITY_FIELD,
                                      rescaledFlattener.getOrder(), flattener.getOrder());
        }

        // scaling and normalization factors
        final FactorsGenerator generator;
        if (wantNormalized == normalized) {
            // the parsed coefficients already match the specified normalization
            generator = (n, m) -> scale;
        } else {
            // we need to normalize/unnormalize parsed coefficients
            final double[][] unnormalizationFactors =
                            GravityFieldFactory.getUnnormalizationFactors(rescaledFlattener.getDegree(),
                                                                          rescaledFlattener.getOrder());
            generator = wantNormalized ?
                (n, m) -> scale / unnormalizationFactors[n][m] :
                (n, m) -> scale * unnormalizationFactors[n][m];
        }

        // perform rescaling
        final double[] rescaled = buildFlatArray(rescaledFlattener, 0.0);
        for (int n = 0; n <= rescaledFlattener.getDegree(); ++n) {
            for (int m = 0; m <= FastMath.min(n, rescaledFlattener.getOrder()); ++m) {
                final int    rescaledndex  = rescaledFlattener.index(n, m);
                final int    originalndex  = originalFlattener.index(n, m);
                rescaled[rescaledndex] = original[originalndex] * generator.factor(n, m);
            }
        }

        return rescaled;

    }

    /** Rescale coefficients arrays.
     * <p>
     * The normalized/unnormalized nature of original coefficients is inherited from previous parsing.
     * </p>
     * @param wantNormalized if true, the rescaled coefficients must be normalized,
     * otherwise they must be un-normalized
     * @param rescaledFlattener converter from triangular to flat form
     * @param originalFlattener converter from triangular to flat form
     * @param original original coefficients
     * @return rescaled coefficients
     * @since 11.1
     */
    protected TimeDependentHarmonic[] rescale(final boolean wantNormalized, final Flattener rescaledFlattener,
                                              final Flattener originalFlattener, final TimeDependentHarmonic[] original) {

        if (rescaledFlattener.getDegree() > originalFlattener.getDegree()) {
            throw new OrekitException(OrekitMessages.TOO_LARGE_DEGREE_FOR_GRAVITY_FIELD,
                                      rescaledFlattener.getDegree(), flattener.getDegree());
        }

        if (rescaledFlattener.getOrder() > originalFlattener.getOrder()) {
            throw new OrekitException(OrekitMessages.TOO_LARGE_ORDER_FOR_GRAVITY_FIELD,
                                      rescaledFlattener.getOrder(), flattener.getOrder());
        }

        // scaling and normalization factors
        final FactorsGenerator generator;
        if (wantNormalized == normalized) {
            // the parsed coefficients already match the specified normalization
            generator = (n, m) -> 1.0;
        } else {
            // we need to normalize/unnormalize parsed coefficients
            final double[][] unnormalizationFactors =
                            GravityFieldFactory.getUnnormalizationFactors(rescaledFlattener.getDegree(),
                                                                          rescaledFlattener.getOrder());
            generator = wantNormalized ?
                (n, m) -> 1.0 / unnormalizationFactors[n][m] :
                (n, m) -> unnormalizationFactors[n][m];
        }

        // perform rescaling
        final TimeDependentHarmonic[] rescaled = new TimeDependentHarmonic[rescaledFlattener.arraySize()];
        for (int n = 0; n <= rescaledFlattener.getDegree(); ++n) {
            for (int m = 0; m <= FastMath.min(n, rescaledFlattener.getOrder()); ++m) {
                final int originalndex = originalFlattener.index(n, m);
                if (original[originalndex] != null) {
                    final int    rescaledndex = rescaledFlattener.index(n, m);
                    final double factor       = generator.factor(n, m);
                    rescaled[rescaledndex]    = new TimeDependentHarmonic(factor, original[originalndex]);
                }
            }
        }

        return rescaled;

    }

    /**
     * Create a date from components. Assumes the time part is noon.
     *
     * @param components year, month, day.
     * @return date.
     */
    protected AbsoluteDate toDate(final DateComponents components) {
        return toDate(components, TimeComponents.H12);
    }

    /**
     * Create a date from components.
     *
     * @param dc dates components.
     * @param tc time components
     * @return date.
     * @since 11.1
     */
    protected AbsoluteDate toDate(final DateComponents dc, final TimeComponents tc) {
        return new AbsoluteDate(dc, tc, timeScale);
    }

    /** Generator for normalization/unnormalization factors.
     * @since 11.1
     */
    private interface FactorsGenerator {
        /** Generator the normalization/unnormalization factors.
         * @param n degree of the gravity field component
         * @param m order of the gravity field component
         * @return factor to apply to term
         */
        double factor(int n, int m);
    }

}