PotentialCoefficientsReader.java

  1. /* Copyright 2002-2015 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.ArrayList;
  23. import java.util.Arrays;
  24. import java.util.List;
  25. import java.util.Locale;

  27. import org.apache.commons.math3.util.FastMath;
  28. import org.apache.commons.math3.util.Precision;
  29. import org.orekit.data.DataLoader;
  30. import org.orekit.errors.OrekitException;
  31. import org.orekit.errors.OrekitMessages;

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

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

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

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

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

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

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

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

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

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

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

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

  79.     /** Simple constructor.
  80.      * <p>Build an uninitialized reader.</p>
  81.      * @param supportedNames regular expression for supported files names
  82.      * @param missingCoefficientsAllowed allow missing coefficients in the input data
  83.      */
  84.     protected PotentialCoefficientsReader(final String supportedNames,
  85.                                           final boolean missingCoefficientsAllowed) {
  86.         this.supportedNames             = supportedNames;
  87.         this.missingCoefficientsAllowed = missingCoefficientsAllowed;
  88.         this.maxParseDegree             = Integer.MAX_VALUE;
  89.         this.maxParseOrder              = Integer.MAX_VALUE;
  90.         this.readComplete               = false;
  91.         this.ae                         = Double.NaN;
  92.         this.mu                         = Double.NaN;
  93.         this.rawC                       = null;
  94.         this.rawS                       = null;
  95.         this.normalized                 = false;
  96.         this.tideSystem                 = TideSystem.UNKNOWN;
  97.     }

  99.     /** Get the regular expression for supported files names.
  100.      * @return regular expression for supported files names
  101.      */
  102.     public String getSupportedNames() {
  103.         return supportedNames;
  104.     }

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

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

  122.     /** Get the degree limit for the next file parsing.
  123.      * @return degree limit for the next file parsing
  124.      * @since 6.0
  125.      */
  126.     public int getMaxParseDegree() {
  127.         return maxParseDegree;
  128.     }

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

  139.     /** Get the order limit for the next file parsing.
  140.      * @return order limit for the next file parsing
  141.      * @since 6.0
  142.      */
  143.     public int getMaxParseOrder() {
  144.         return maxParseOrder;
  145.     }

  147.     /** {@inheritDoc} */
  148.     public boolean stillAcceptsData() {
  149.         return !(readComplete &&
  150.                  getMaxAvailableDegree() >= getMaxParseDegree() &&
  151.                  getMaxAvailableOrder()  >= getMaxParseOrder());
  152.     }

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

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

  168.     /** Get the central body reference radius.
  169.      * @return central body reference radius
  170.      */
  171.     protected double getAe() {
  172.         return ae;
  173.     }

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

  182.     /** Get the central body attraction coefficient.
  183.      * @return central body attraction coefficient
  184.      */
  185.     protected double getMu() {
  186.         return mu;
  187.     }

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

  196.     /** Get the {@link TideSystem} used in the gravity field.
  197.      * @return tide system used in the gravity field
  198.      */
  199.     protected TideSystem getTideSystem() {
  200.         return tideSystem;
  201.     }

  203.     /** Set the tesseral-sectorial coefficients matrix.
  204.      * @param rawNormalized if true, raw coefficients are normalized
  205.      * @param c raw tesseral-sectorial coefficients matrix
  206.      * (a reference to the array will be stored)
  207.      * @param s raw tesseral-sectorial coefficients matrix
  208.      * (a reference to the array will be stored)
  209.      * @param name name of the file (or zip entry)
  210.      * @exception OrekitException if a coefficient is missing
  211.      */
  212.     protected void setRawCoefficients(final boolean rawNormalized,
  213.                                       final double[][] c, final double[][] s,
  214.                                       final String name)
  215.         throws OrekitException {

  217.         // normalization indicator
  218.         normalized = rawNormalized;

  220.         // set known constant values, if they were not defined in the file.
  221.         // See Hofmann-Wellenhof and Moritz, "Physical Geodesy",
  222.         // section 2.6 Harmonics of Lower Degree.
  223.         // All S_i,0 are irrelevant because they are multiplied by zero.
  224.         // C0,0 is 1, the central part, since all coefficients are normalized by GM.
  225.         setIfUnset(c, 0, 0, 1);
  226.         setIfUnset(s, 0, 0, 0);
  227.         // C1,0, C1,1, and S1,1 are the x,y,z coordinates of the center of mass,
  228.         // which are 0 since all coefficients are given in an Earth centered frame
  229.         setIfUnset(c, 1, 0, 0);
  230.         setIfUnset(s, 1, 0, 0);
  231.         setIfUnset(c, 1, 1, 0);
  232.         setIfUnset(s, 1, 1, 0);

  234.         // cosine part
  235.         for (int i = 0; i < c.length; ++i) {
  236.             for (int j = 0; j < c[i].length; ++j) {
  237.                 if (Double.isNaN(c[i][j])) {
  238.                     throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  239.                                               'C', i, j, name);
  240.                 }
  241.             }
  242.         }
  243.         rawC = c;

  245.         // sine part
  246.         for (int i = 0; i < s.length; ++i) {
  247.             for (int j = 0; j < s[i].length; ++j) {
  248.                 if (Double.isNaN(s[i][j])) {
  249.                     throw new OrekitException(OrekitMessages.MISSING_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  250.                                               'S', i, j, name);
  251.                 }
  252.             }
  253.         }
  254.         rawS = s;

  256.     }

  258.     /**
  259.      * Set a coefficient if it has not been set already.
  260.      * <p>
  261.      * If {@code array[i][j]} is 0 or NaN this method sets it to {@code value} and returns
  262.      * {@code true}. Otherwise the original value of {@code array[i][j]} is preserved and
  263.      * this method return {@code false}.
  264.      * <p>
  265.      * If {@code array[i][j]} does not exist then this method returns {@code false}.
  266.      *
  267.      * @param array the coefficient array.
  268.      * @param i     degree, the first index to {@code array}.
  269.      * @param j     order, the second index to {@code array}.
  270.      * @param value the new value to set.
  271.      * @return {@code true} if the coefficient was set to {@code value}, {@code false} if
  272.      * the coefficient was not set to {@code value}. A {@code false} return indicates the
  273.      * coefficient has previously been set to a non-NaN, non-zero value.
  274.      */
  275.     private boolean setIfUnset(final double[][] array,
  276.                                final int i,
  277.                                final int j,
  278.                                final double value) {
  279.         if (array.length > i && array[i].length > j &&
  280.                 (Double.isNaN(array[i][j]) || Precision.equals(array[i][j], 0.0, 1))) {
  281.             // the coefficient was not already initialized
  282.             array[i][j] = value;
  283.             return true;
  284.         } else {
  285.             return false;
  286.         }
  287.     }

  289.     /** Get the maximal degree available in the last file parsed.
  290.      * @return maximal degree available in the last file parsed
  291.      * @since 6.0
  292.      */
  293.     public int getMaxAvailableDegree() {
  294.         return rawC.length - 1;
  295.     }

  297.     /** Get the maximal order available in the last file parsed.
  298.      * @return maximal order available in the last file parsed
  299.      * @since 6.0
  300.      */
  301.     public int getMaxAvailableOrder() {
  302.         return rawC[rawC.length - 1].length - 1;
  303.     }

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

  309.     /** Get a provider for read spherical harmonics coefficients.
  310.      * @param wantNormalized if true, the provider will provide normalized coefficients,
  311.      * otherwise it will provide un-normalized coefficients
  312.      * @param degree maximal degree
  313.      * @param order maximal order
  314.      * @return a new provider
  315.      * @exception OrekitException if the requested maximal degree or order exceeds the
  316.      * available degree or order or if no gravity field has read yet
  317.      * @see #getConstantProvider(boolean, int, int)
  318.      * @since 6.0
  319.      */
  320.     public abstract RawSphericalHarmonicsProvider getProvider(boolean wantNormalized, int degree, int order)
  321.         throws OrekitException;

  323.     /** Get a time-independent provider for read spherical harmonics coefficients.
  324.      * @param wantNormalized if true, the raw provider must provide normalized coefficients,
  325.      * otherwise it will provide un-normalized coefficients
  326.      * @param degree maximal degree
  327.      * @param order maximal order
  328.      * @return a new provider, with no time-dependent parts
  329.      * @exception OrekitException if the requested maximal degree or order exceeds the
  330.      * available degree or order or if no gravity field has read yet
  331.      * @see #getProvider(boolean, int, int)
  332.      * @since 6.0
  333.      */
  334.     protected ConstantSphericalHarmonics getConstantProvider(final boolean wantNormalized,
  335.                                                              final int degree, final int order)
  336.         throws OrekitException {

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

  342.         if (degree >= rawC.length) {
  343.             throw new OrekitException(OrekitMessages.TOO_LARGE_DEGREE_FOR_GRAVITY_FIELD,
  344.                                       degree, rawC.length - 1);
  345.         }

  347.         if (order >= rawC[rawC.length - 1].length) {
  348.             throw new OrekitException(OrekitMessages.TOO_LARGE_ORDER_FOR_GRAVITY_FIELD,
  349.                                       order, rawC[rawC.length - 1].length);
  350.         }

  352.         // fix normalization
  353.         final double[][] truncatedC = buildTriangularArray(degree, order, 0.0);
  354.         final double[][] truncatedS = buildTriangularArray(degree, order, 0.0);
  355.         rescale(1.0, normalized, rawC, rawS, wantNormalized, truncatedC, truncatedS);

  357.         return new ConstantSphericalHarmonics(ae, mu, tideSystem, truncatedC, truncatedS);

  359.     }

  361.     /** Build a coefficients triangular array.
  362.      * @param degree array degree
  363.      * @param order array order
  364.      * @param value initial value to put in array elements
  365.      * @return built array
  366.      */
  367.     protected static double[][] buildTriangularArray(final int degree, final int order, final double value) {
  368.         final double[][] array = new double[degree + 1][];
  369.         for (int k = 0; k < array.length; ++k) {
  370.             array[k] = buildRow(k, order, value);
  371.         }
  372.         return array;
  373.     }

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

  386.     /** Build a coefficients row.
  387.      * @param degree row degree
  388.      * @param order row order
  389.      * @param value initial value to put in array elements
  390.      * @return built row
  391.      */
  392.     protected static double[] buildRow(final int degree, final int order, final double value) {
  393.         final double[] row = new double[FastMath.min(order, degree) + 1];
  394.         Arrays.fill(row, value);
  395.         return row;
  396.     }

  398.     /** Extend a list of lists of coefficients if needed.
  399.      * @param list list of lists of coefficients
  400.      * @param degree degree required to be present
  401.      * @param order order required to be present
  402.      * @param value initial value to put in list elements
  403.      */
  404.     protected void extendListOfLists(final List<List<Double>> list, final int degree, final int order,
  405.                                      final double value) {
  406.         for (int i = list.size(); i <= degree; ++i) {
  407.             list.add(new ArrayList<Double>());
  408.         }
  409.         final List<Double> listN = list.get(degree);
  410.         final Double v = Double.valueOf(value);
  411.         for (int j = listN.size(); j <= order; ++j) {
  412.             listN.add(v);
  413.         }
  414.     }

  416.     /** Convert a list of list into an array.
  417.      * @param list list of lists of coefficients
  418.      * @return a new array
  419.      */
  420.     protected double[][] toArray(final List<List<Double>> list) {
  421.         final double[][] array = new double[list.size()][];
  422.         for (int i = 0; i < array.length; ++i) {
  423.             array[i] = new double[list.get(i).size()];
  424.             for (int j = 0; j < array[i].length; ++j) {
  425.                 array[i][j] = list.get(i).get(j);
  426.             }
  427.         }
  428.         return array;
  429.     }

  431.     /** Parse a coefficient.
  432.      * @param field text field to parse
  433.      * @param list list where to put the coefficient
  434.      * @param i first index in the list
  435.      * @param j second index in the list
  436.      * @param cName name of the coefficient
  437.      * @param name name of the file
  438.      * @exception OrekitException if the coefficient is already set
  439.      */
  440.     protected void parseCoefficient(final String field, final List<List<Double>> list,
  441.                                     final int i, final int j,
  442.                                     final String cName, final String name)
  443.         throws OrekitException {
  444.         final double value    = parseDouble(field);
  445.         final double oldValue = list.get(i).get(j);
  446.         if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 1)) {
  447.             // the coefficient was not already initialized
  448.             list.get(i).set(j, value);
  449.         } else {
  450.             throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  451.                                       name, i, j, name);
  452.         }
  453.     }

  455.     /** Parse a coefficient.
  456.      * @param field text field to parse
  457.      * @param array array where to put the coefficient
  458.      * @param i first index in the list
  459.      * @param j second index in the list
  460.      * @param cName name of the coefficient
  461.      * @param name name of the file
  462.      * @exception OrekitException if the coefficient is already set
  463.      */
  464.     protected void parseCoefficient(final String field, final double[][] array,
  465.                                     final int i, final int j,
  466.                                     final String cName, final String name)
  467.         throws OrekitException {
  468.         final double value    = parseDouble(field);
  469.         final double oldValue = array[i][j];
  470.         if (Double.isNaN(oldValue) || Precision.equals(oldValue, 0.0, 1)) {
  471.             // the coefficient was not already initialized
  472.             array[i][j] = value;
  473.         } else {
  474.             throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  475.                                       name, i, j, name);
  476.         }
  477.     }

  479.     /** Rescale coefficients arrays.
  480.      * @param scale general scaling factor to apply to all elements
  481.      * @param normalizedOrigin if true, the origin coefficients are normalized
  482.      * @param originC cosine part of the origina coefficients
  483.      * @param originS sine part of the origin coefficients
  484.      * @param wantNormalized if true, the rescaled coefficients must be normalized
  485.      * @param rescaledC cosine part of the rescaled coefficients to fill in (may be the originC array)
  486.      * @param rescaledS sine part of the rescaled coefficients to fill in (may be the originS array)
  487.      * @exception OrekitException if normalization/unnormalization fails because of an underflow
  488.      * due to too high degree/order
  489.      */
  490.     protected static void rescale(final double scale,
  491.                                   final boolean normalizedOrigin, final double[][] originC,
  492.                                   final double[][] originS, final boolean wantNormalized,
  493.                                   final double[][] rescaledC, final double[][] rescaledS)
  494.         throws OrekitException {

  496.         if (wantNormalized == normalizedOrigin) {
  497.             // apply only the general scaling factor
  498.             for (int i = 0; i < rescaledC.length; ++i) {
  499.                 final double[] rCi = rescaledC[i];
  500.                 final double[] rSi = rescaledS[i];
  501.                 final double[] oCi = originC[i];
  502.                 final double[] oSi = originS[i];
  503.                 for (int j = 0; j < rCi.length; ++j) {
  504.                     rCi[j] = oCi[j] * scale;
  505.                     rSi[j] = oSi[j] * scale;
  506.                 }
  507.             }
  508.         } else {

  510.             // we have to re-scale the coefficients
  511.             // (we use rescaledC.length - 1 for the order instead of rescaledC[rescaledC.length - 1].length - 1
  512.             //  because typically trend or pulsation arrays are irregular, some test cases have
  513.             //  order 2 elements at degree 2, but only order 1 elements for higher degrees for example)
  514.             final double[][] factors = GravityFieldFactory.getUnnormalizationFactors(rescaledC.length - 1,
  515.                                                                                      rescaledC.length - 1);

  517.             if (wantNormalized) {
  518.                 // normalize the coefficients
  519.                 for (int i = 0; i < rescaledC.length; ++i) {
  520.                     final double[] rCi = rescaledC[i];
  521.                     final double[] rSi = rescaledS[i];
  522.                     final double[] oCi = originC[i];
  523.                     final double[] oSi = originS[i];
  524.                     final double[] fi  = factors[i];
  525.                     for (int j = 0; j < rCi.length; ++j) {
  526.                         final double factor = scale / fi[j];
  527.                         rCi[j] = oCi[j] * factor;
  528.                         rSi[j] = oSi[j] * factor;
  529.                     }
  530.                 }
  531.             } else {
  532.                 // un-normalize the coefficients
  533.                 for (int i = 0; i < rescaledC.length; ++i) {
  534.                     final double[] rCi = rescaledC[i];
  535.                     final double[] rSi = rescaledS[i];
  536.                     final double[] oCi = originC[i];
  537.                     final double[] oSi = originS[i];
  538.                     final double[] fi  = factors[i];
  539.                     for (int j = 0; j < rCi.length; ++j) {
  540.                         final double factor = scale * fi[j];
  541.                         rCi[j] = oCi[j] * factor;
  542.                         rSi[j] = oSi[j] * factor;
  543.                     }
  544.                 }
  545.             }

  547.         }
  548.     }

  550. }
Created with JaCoCo 0.7.2.201409121644