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

  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.Precision;
  27. import org.orekit.annotation.DefaultDataContext;
  28. import org.orekit.data.DataContext;
  29. import org.orekit.data.DataLoader;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitMessages;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.time.DateComponents;
  34. import org.orekit.time.TimeComponents;
  35. import org.orekit.time.TimeScale;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  89.     /** Simple constructor.
  90.      * <p>Build an uninitialized reader.</p>
  91.      *
  92.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  93.      *
  94.      * @param supportedNames regular expression for supported files names
  95.      * @param missingCoefficientsAllowed allow missing coefficients in the input data
  96.      * @see #PotentialCoefficientsReader(String, boolean, TimeScale)
  97.      */
  98.     @DefaultDataContext
  99.     protected PotentialCoefficientsReader(final String supportedNames,
  100.                                           final boolean missingCoefficientsAllowed) {
  101.         this(supportedNames, missingCoefficientsAllowed,
  102.                 DataContext.getDefault().getTimeScales().getTT());
  103.     }

  105.     /** Simple constructor.
  106.      * <p>Build an uninitialized reader.</p>
  107.      * @param supportedNames regular expression for supported files names
  108.      * @param missingCoefficientsAllowed allow missing coefficients in the input data
  109.      * @param timeScale to use when parsing dates.
  110.      * @since 10.1
  111.      */
  112.     protected PotentialCoefficientsReader(final String supportedNames,
  113.                                           final boolean missingCoefficientsAllowed,
  114.                                           final TimeScale timeScale) {
  115.         this.supportedNames             = supportedNames;
  116.         this.missingCoefficientsAllowed = missingCoefficientsAllowed;
  117.         this.maxParseDegree             = Integer.MAX_VALUE;
  118.         this.maxParseOrder              = Integer.MAX_VALUE;
  119.         this.readComplete               = false;
  120.         this.ae                         = Double.NaN;
  121.         this.mu                         = Double.NaN;
  122.         this.flattener                  = null;
  123.         this.rawC                       = null;
  124.         this.rawS                       = null;
  125.         this.normalized                 = false;
  126.         this.tideSystem                 = TideSystem.UNKNOWN;
  127.         this.timeScale                  = timeScale;
  128.     }

  130.     /** Get the regular expression for supported files names.
  131.      * @return regular expression for supported files names
  132.      */
  133.     public String getSupportedNames() {
  134.         return supportedNames;
  135.     }

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

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

  153.     /** Get the degree limit for the next file parsing.
  154.      * @return degree limit for the next file parsing
  155.      * @since 6.0
  156.      */
  157.     public int getMaxParseDegree() {
  158.         return maxParseDegree;
  159.     }

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

  170.     /** Get the order limit for the next file parsing.
  171.      * @return order limit for the next file parsing
  172.      * @since 6.0
  173.      */
  174.     public int getMaxParseOrder() {
  175.         return maxParseOrder;
  176.     }

  178.     /** {@inheritDoc} */
  179.     public boolean stillAcceptsData() {
  180.         return !(readComplete &&
  181.                  getMaxAvailableDegree() >= getMaxParseDegree() &&
  182.                  getMaxAvailableOrder()  >= getMaxParseOrder());
  183.     }

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

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

  199.     /** Get the central body reference radius.
  200.      * @return central body reference radius
  201.      */
  202.     protected double getAe() {
  203.         return ae;
  204.     }

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

  213.     /** Get the central body attraction coefficient.
  214.      * @return central body attraction coefficient
  215.      */
  216.     protected double getMu() {
  217.         return mu;
  218.     }

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

  227.     /** Get the {@link TideSystem} used in the gravity field.
  228.      * @return tide system used in the gravity field
  229.      */
  230.     protected TideSystem getTideSystem() {
  231.         return tideSystem;
  232.     }

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

  246.         this.flattener = f;

  248.         // normalization indicator
  249.         normalized = rawNormalized;

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

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

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

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

  292.     }

  294.     /**
  295.      * Set a coefficient if it has not been set already.
  296.      * <p>
  297.      * If {@code array[i]} is 0 or NaN this method sets it to {@code value} and returns
  298.      * {@code true}. Otherwise the original value of {@code array[i]} is preserved and
  299.      * this method return {@code false}.
  300.      * <p>
  301.      * If {@code array[i]} does not exist then this method returns {@code false}.
  302.      *
  303.      * @param array the coefficient array.
  304.      * @param i     index in array.
  305.      * @param value the new value to set.
  306.      * @return {@code true} if the coefficient was set to {@code value}, {@code false} if
  307.      * the coefficient was not set to {@code value}. A {@code false} return indicates the
  308.      * coefficient has previously been set to a non-NaN, non-zero value.
  309.      */
  310.     private boolean setIfUnset(final double[] array, final int i, final double value) {
  311.         if (array.length > i && (Double.isNaN(array[i]) || Precision.equals(array[i], 0.0, 0))) {
  312.             // the coefficient was not already initialized
  313.             array[i] = value;
  314.             return true;
  315.         } else {
  316.             return false;
  317.         }
  318.     }

  320.     /** Get the maximal degree available in the last file parsed.
  321.      * @return maximal degree available in the last file parsed
  322.      * @since 6.0
  323.      */
  324.     public int getMaxAvailableDegree() {
  325.         return flattener.getDegree();
  326.     }

  328.     /** Get the maximal order available in the last file parsed.
  329.      * @return maximal order available in the last file parsed
  330.      * @since 6.0
  331.      */
  332.     public int getMaxAvailableOrder() {
  333.         return flattener.getOrder();
  334.     }

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

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

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

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

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

  375.     }

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

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

  400.     /** Build a coefficients row.
  401.      * @param degree row degree
  402.      * @param order row order
  403.      * @param value initial value to put in array elements
  404.      * @return built row
  405.      */
  406.     protected static double[] buildRow(final int degree, final int order, final double value) {
  407.         final double[] row = new double[FastMath.min(order, degree) + 1];
  408.         Arrays.fill(row, value);
  409.         return row;
  410.     }

  412.     /** Parse a coefficient.
  413.      * @param field text field to parse
  414.      * @param f converter from triangular to flat form
  415.      * @param array array where to put the coefficient
  416.      * @param i first index in the list
  417.      * @param j second index in the list
  418.      * @param cName name of the coefficient
  419.      * @param name name of the file
  420.      * @since 11.1
  421.      */
  422.     protected void parseCoefficient(final String field, final Flattener f,
  423.                                     final double[] array, final int i, final int j,
  424.                                     final String cName, final String name) {
  425.         final int    index    = f.index(i, j);
  426.         final double value    = parseDouble(field);
  427.         final double oldValue = array[index];
  428.         if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 0)) {
  429.             // the coefficient was not already initialized
  430.             array[index] = value;
  431.         } else {
  432.             throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  433.                                       name, i, j, name);
  434.         }
  435.     }

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

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

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

  463.         // scaling and normalization factors
  464.         final FactorsGenerator generator;
  465.         if (wantNormalized == normalized) {
  466.             // the parsed coefficients already match the specified normalization
  467.             generator = (n, m) -> scale;
  468.         } else {
  469.             // we need to normalize/unnormalize parsed coefficients
  470.             final double[][] unnormalizationFactors =
  471.                             GravityFieldFactory.getUnnormalizationFactors(rescaledFlattener.getDegree(),
  472.                                                                           rescaledFlattener.getOrder());
  473.             generator = wantNormalized ?
  474.                 (n, m) -> scale / unnormalizationFactors[n][m] :
  475.                 (n, m) -> scale * unnormalizationFactors[n][m];
  476.         }

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

  488.         return rescaled;

  490.     }

  492.     /** Rescale coefficients arrays.
  493.      * <p>
  494.      * The normalized/unnormalized nature of original coefficients is inherited from previous parsing.
  495.      * </p>
  496.      * @param wantNormalized if true, the rescaled coefficients must be normalized,
  497.      * otherwise they must be un-normalized
  498.      * @param rescaledFlattener converter from triangular to flat form
  499.      * @param originalFlattener converter from triangular to flat form
  500.      * @param original original coefficients
  501.      * @return rescaled coefficients
  502.      * @since 11.1
  503.      */
  504.     protected TimeDependentHarmonic[] rescale(final boolean wantNormalized, final Flattener rescaledFlattener,
  505.                                               final Flattener originalFlattener, final TimeDependentHarmonic[] original) {

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

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

  517.         // scaling and normalization factors
  518.         final FactorsGenerator generator;
  519.         if (wantNormalized == normalized) {
  520.             // the parsed coefficients already match the specified normalization
  521.             generator = (n, m) -> 1.0;
  522.         } else {
  523.             // we need to normalize/unnormalize parsed coefficients
  524.             final double[][] unnormalizationFactors =
  525.                             GravityFieldFactory.getUnnormalizationFactors(rescaledFlattener.getDegree(),
  526.                                                                           rescaledFlattener.getOrder());
  527.             generator = wantNormalized ?
  528.                 (n, m) -> 1.0 / unnormalizationFactors[n][m] :
  529.                 (n, m) -> unnormalizationFactors[n][m];
  530.         }

  532.         // perform rescaling
  533.         final TimeDependentHarmonic[] rescaled = new TimeDependentHarmonic[rescaledFlattener.arraySize()];
  534.         for (int n = 0; n <= rescaledFlattener.getDegree(); ++n) {
  535.             for (int m = 0; m <= FastMath.min(n, rescaledFlattener.getOrder()); ++m) {
  536.                 final int originalndex = originalFlattener.index(n, m);
  537.                 if (original[originalndex] != null) {
  538.                     final int    rescaledndex = rescaledFlattener.index(n, m);
  539.                     final double factor       = generator.factor(n, m);
  540.                     rescaled[rescaledndex]    = new TimeDependentHarmonic(factor, original[originalndex]);
  541.                 }
  542.             }
  543.         }

  545.         return rescaled;

  547.     }

  549.     /**
  550.      * Create a date from components. Assumes the time part is noon.
  551.      *
  552.      * @param components year, month, day.
  553.      * @return date.
  554.      */
  555.     protected AbsoluteDate toDate(final DateComponents components) {
  556.         return toDate(components, TimeComponents.H12);
  557.     }

  559.     /**
  560.      * Create a date from components.
  561.      *
  562.      * @param dc dates components.
  563.      * @param tc time components
  564.      * @return date.
  565.      * @since 11.1
  566.      */
  567.     protected AbsoluteDate toDate(final DateComponents dc, final TimeComponents tc) {
  568.         return new AbsoluteDate(dc, tc, timeScale);
  569.     }

  571.     /** Generator for normalization/unnormalization factors.
  572.      * @since 11.1
  573.      */
  574.     private interface FactorsGenerator {
  575.         /** Generator the normalization/unnormalization factors.
  576.          * @param n degree of the gravity field component
  577.          * @param m order of the gravity field component
  578.          * @return factor to apply to term
  579.          */
  580.         double factor(int n, int m);
  581.     }

  583. }
Created with JaCoCo 0.8.12.202403310830