ICGEMFormatReader.java

  1. /* Copyright 2002-2017 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;

  18. import java.io.BufferedReader;
  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.InputStreamReader;
  22. import java.text.ParseException;
  23. import java.util.ArrayList;
  24. import java.util.HashMap;
  25. import java.util.List;
  26. import java.util.Map;

  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.Precision;
  29. import org.orekit.errors.OrekitException;
  30. import org.orekit.errors.OrekitMessages;
  31. import org.orekit.errors.OrekitParseException;
  32. import org.orekit.time.DateComponents;
  33. import org.orekit.utils.Constants;

  34. /** Reader for the ICGEM gravity field format.
  35.  *
  36.  * <p>This format is used to describe the gravity field of EIGEN models
  37.  * published by the GFZ Potsdam since 2004. It is described in Franz
  38.  * Barthelmes and Christoph F&ouml;rste paper: "the ICGEM-format".
  39.  * The 2006-02-28 version of this paper can be found <a
  40.  * href="http://op.gfz-potsdam.de/grace/results/grav/g005_ICGEM-Format.pdf">here</a>
  41.  * and the 2011-06-07 version of this paper can be found <a
  42.  * href="http://icgem.gfz-potsdam.de/ICGEM/documents/ICGEM-Format-2011.pdf">here</a>.
  43.  * These versions differ in time-dependent coefficients, which are linear-only prior
  44.  * to 2011 (up to eigen-5 model) and have also harmonic effects after that date
  45.  * (starting with eigen-6 model). Both versions are supported by the class.</p>
  46.  * <p>
  47.  * This reader uses relaxed check on the gravity constant key so any key ending
  48.  * in gravity_constant is accepted and not only earth_gravity_constant as specified
  49.  * in the previous documents. This allows to read also non Earth gravity fields
  50.  * as found in <a href="http://icgem.gfz-potsdam.de/ICGEM/ModelstabBodies.html">ICGEM
  51.  * - Gravity Field Models of other Celestial Bodies</a> page to be read.
  52.  * </p>
  53.  * <p>
  54.  * In order to simplify implementation, some design choices have been made: the
  55.  * reference date and the periods of harmonic pulsations are stored globally and
  56.  * not on a per-coefficient basis. This has the following implications:
  57.  * </p>
  58.  * <ul>
  59.  *   <li>
  60.  *     all time-stamped coefficients must share the same reference date, otherwise
  61.  *     an error will be triggered during parsing,
  62.  *   </li>
  63.  *   <li>
  64.  *     in order to avoid too large memory and CPU consumption, only a few different
  65.  *     periods should appear in the file,
  66.  *   </li>
  67.  *   <li>
  68.  *     only one occurrence of each coefficient may appear in the file, otherwise
  69.  *     an error will be triggered during parsing. Multiple occurrences with different
  70.  *     time stamps are forbidden (both because they correspond to a duplicated entry
  71.  *     and because they define two different reference dates as per previous design
  72.  *     choice).
  73.  *   </li>
  74.  * </ul>
  75.  *
  76.  * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
  77.  *  which will determine which reader to use with the selected gravity field file.</p>
  78.  *
  79.  * @see GravityFieldFactory
  80.  * @author Luc Maisonobe
  81.  */
  82. public class ICGEMFormatReader extends PotentialCoefficientsReader {

  83.     /** Product type. */
  84.     private static final String PRODUCT_TYPE            = "product_type";

  85.     /** Gravity field product type. */
  86.     private static final String GRAVITY_FIELD           = "gravity_field";

  87.     /** Gravity constant marker. */
  88.     private static final String GRAVITY_CONSTANT        = "gravity_constant";

  89.     /** Reference radius. */
  90.     private static final String REFERENCE_RADIUS        = "radius";

  91.     /** Max degree. */
  92.     private static final String MAX_DEGREE              = "max_degree";

  93.     /** Tide system indicator. */
  94.     private static final String TIDE_SYSTEM_INDICATOR   = "tide_system";

  95.     /** Indicator value for zero-tide system. */
  96.     private static final String ZERO_TIDE               = "zero_tide";

  97.     /** Indicator value for tide-free system. */
  98.     private static final String TIDE_FREE               = "tide_free";

  99.     /** Indicator value for unknown tide system. */
  100.     private static final String TIDE_UNKNOWN            = "unknown";

  101.     /** Normalization indicator. */
  102.     private static final String NORMALIZATION_INDICATOR = "norm";

  103.     /** Indicator value for normalized coefficients. */
  104.     private static final String NORMALIZED              = "fully_normalized";

  105.     /** Indicator value for un-normalized coefficients. */
  106.     private static final String UNNORMALIZED            = "unnormalized";

  107.     /** End of header marker. */
  108.     private static final String END_OF_HEADER           = "end_of_head";

  109.     /** Gravity field coefficient. */
  110.     private static final String GFC                     = "gfc";

  111.     /** Time stamped gravity field coefficient. */
  112.     private static final String GFCT                    = "gfct";

  113.     /** Gravity field coefficient first time derivative. */
  114.     private static final String DOT                     = "dot";

  115.     /** Gravity field coefficient trend. */
  116.     private static final String TRND                    = "trnd";

  117.     /** Gravity field coefficient sine amplitude. */
  118.     private static final String ASIN                    = "asin";

  119.     /** Gravity field coefficient cosine amplitude. */
  120.     private static final String ACOS                    = "acos";

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

  123.     /** Indicator for normalized coefficients. */
  124.     private boolean normalized;

  125.     /** Reference date. */
  126.     private DateComponents referenceDate;

  127.     /** Secular trend of the cosine coefficients. */
  128.     private final List<List<Double>> cTrend;

  129.     /** Secular trend of the sine coefficients. */
  130.     private final List<List<Double>> sTrend;

  131.     /** Cosine part of the cosine coefficients pulsation. */
  132.     private final Map<Double, List<List<Double>>> cCos;

  133.     /** Sine part of the cosine coefficients pulsation. */
  134.     private final Map<Double, List<List<Double>>> cSin;

  135.     /** Cosine part of the sine coefficients pulsation. */
  136.     private final Map<Double, List<List<Double>>> sCos;

  137.     /** Sine part of the sine coefficients pulsation. */
  138.     private final Map<Double, List<List<Double>>> sSin;

  139.     /** Simple constructor.
  140.      * @param supportedNames regular expression for supported files names
  141.      * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
  142.      */
  143.     public ICGEMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
  144.         super(supportedNames, missingCoefficientsAllowed);
  145.         referenceDate = null;
  146.         cTrend = new ArrayList<List<Double>>();
  147.         sTrend = new ArrayList<List<Double>>();
  148.         cCos   = new HashMap<Double, List<List<Double>>>();
  149.         cSin   = new HashMap<Double, List<List<Double>>>();
  150.         sCos   = new HashMap<Double, List<List<Double>>>();
  151.         sSin   = new HashMap<Double, List<List<Double>>>();
  152.     }

  153.     /** {@inheritDoc} */
  154.     public void loadData(final InputStream input, final String name)
  155.         throws IOException, ParseException, OrekitException {

  156.         // reset the indicator before loading any data
  157.         setReadComplete(false);
  158.         referenceDate = null;
  159.         cTrend.clear();
  160.         sTrend.clear();
  161.         cCos.clear();
  162.         cSin.clear();
  163.         sCos.clear();
  164.         sSin.clear();

  165.         // by default, the field is normalized with unknown tide system
  166.         // (will be overridden later if non-default)
  167.         normalized = true;
  168.         tideSystem = TideSystem.UNKNOWN;

  169.         final BufferedReader r = new BufferedReader(new InputStreamReader(input, "UTF-8"));
  170.         boolean inHeader = true;
  171.         double[][] c               = null;
  172.         double[][] s               = null;
  173.         boolean okCoeffs           = false;
  174.         int lineNumber   = 0;
  175.         for (String line = r.readLine(); line != null; line = r.readLine()) {
  176.             try {
  177.                 ++lineNumber;
  178.                 if (line.trim().length() == 0) {
  179.                     continue;
  180.                 }
  181.                 final String[] tab = line.split("\\s+");
  182.                 if (inHeader) {
  183.                     if ((tab.length == 2) && PRODUCT_TYPE.equals(tab[0])) {
  184.                         if (!GRAVITY_FIELD.equals(tab[1])) {
  185.                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  186.                                                            lineNumber, name, line);
  187.                         }
  188.                     } else if ((tab.length == 2) && tab[0].endsWith(GRAVITY_CONSTANT)) {
  189.                         setMu(parseDouble(tab[1]));
  190.                     } else if ((tab.length == 2) && REFERENCE_RADIUS.equals(tab[0])) {
  191.                         setAe(parseDouble(tab[1]));
  192.                     } else if ((tab.length == 2) && MAX_DEGREE.equals(tab[0])) {

  193.                         final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
  194.                         final int order  = FastMath.min(getMaxParseOrder(), degree);
  195.                         c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
  196.                         s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);

  197.                     } else if ((tab.length == 2) && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
  198.                         if (ZERO_TIDE.equals(tab[1])) {
  199.                             tideSystem = TideSystem.ZERO_TIDE;
  200.                         } else if (TIDE_FREE.equals(tab[1])) {
  201.                             tideSystem = TideSystem.TIDE_FREE;
  202.                         } else if (TIDE_UNKNOWN.equals(tab[1])) {
  203.                             tideSystem = TideSystem.UNKNOWN;
  204.                         } else {
  205.                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  206.                                                            lineNumber, name, line);
  207.                         }
  208.                     } else if ((tab.length == 2) && NORMALIZATION_INDICATOR.equals(tab[0])) {
  209.                         if (NORMALIZED.equals(tab[1])) {
  210.                             normalized = true;
  211.                         } else if (UNNORMALIZED.equals(tab[1])) {
  212.                             normalized = false;
  213.                         } else {
  214.                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  215.                                                            lineNumber, name, line);
  216.                         }
  217.                     } else if ((tab.length == 2) && END_OF_HEADER.equals(tab[0])) {
  218.                         inHeader = false;
  219.                     }
  220.                 } else {
  221.                     if ((tab.length == 7 && GFC.equals(tab[0])) || (tab.length == 8 && GFCT.equals(tab[0]))) {

  222.                         final int i = Integer.parseInt(tab[1]);
  223.                         final int j = Integer.parseInt(tab[2]);
  224.                         if (i < c.length && j < c[i].length) {

  225.                             parseCoefficient(tab[3], c, i, j, "C", name);
  226.                             parseCoefficient(tab[4], s, i, j, "S", name);
  227.                             okCoeffs = true;

  228.                             if (tab.length == 8) {
  229.                                 // check the reference date (format yyyymmdd)
  230.                                 final DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
  231.                                                                                    Integer.parseInt(tab[7].substring(4, 6)),
  232.                                                                                    Integer.parseInt(tab[7].substring(6, 8)));
  233.                                 if (referenceDate == null) {
  234.                                     // first reference found, store it
  235.                                     referenceDate = localRef;
  236.                                 } else if (!referenceDate.equals(localRef)) {
  237.                                     throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
  238.                                                               referenceDate, localRef, name);
  239.                                 }
  240.                             }

  241.                         }
  242.                     } else if (tab.length == 7 && (DOT.equals(tab[0]) || TRND.equals(tab[0]))) {

  243.                         final int i = Integer.parseInt(tab[1]);
  244.                         final int j = Integer.parseInt(tab[2]);
  245.                         if (i < c.length && j < c[i].length) {

  246.                             // store the secular trend coefficients
  247.                             extendListOfLists(cTrend, i, j, 0.0);
  248.                             extendListOfLists(sTrend, i, j, 0.0);
  249.                             parseCoefficient(tab[3], cTrend, i, j, "Ctrend", name);
  250.                             parseCoefficient(tab[4], sTrend, i, j, "Strend", name);

  251.                         }

  252.                     } else if (tab.length == 8 && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {

  253.                         final int i = Integer.parseInt(tab[1]);
  254.                         final int j = Integer.parseInt(tab[2]);
  255.                         if (i < c.length && j < c[i].length) {

  256.                             // extract arrays corresponding to period
  257.                             final Double period = Double.valueOf(tab[7]);
  258.                             if (!cCos.containsKey(period)) {
  259.                                 cCos.put(period, new ArrayList<List<Double>>());
  260.                                 cSin.put(period, new ArrayList<List<Double>>());
  261.                                 sCos.put(period, new ArrayList<List<Double>>());
  262.                                 sSin.put(period, new ArrayList<List<Double>>());
  263.                             }
  264.                             final List<List<Double>> cCosPeriod = cCos.get(period);
  265.                             final List<List<Double>> cSinPeriod = cSin.get(period);
  266.                             final List<List<Double>> sCosPeriod = sCos.get(period);
  267.                             final List<List<Double>> sSinPeriod = sSin.get(period);

  268.                             // store the pulsation coefficients
  269.                             extendListOfLists(cCosPeriod, i, j, 0.0);
  270.                             extendListOfLists(cSinPeriod, i, j, 0.0);
  271.                             extendListOfLists(sCosPeriod, i, j, 0.0);
  272.                             extendListOfLists(sSinPeriod, i, j, 0.0);
  273.                             if (ACOS.equals(tab[0])) {
  274.                                 parseCoefficient(tab[3], cCosPeriod, i, j, "Ccos", name);
  275.                                 parseCoefficient(tab[4], sCosPeriod, i, j, "SCos", name);
  276.                             } else {
  277.                                 parseCoefficient(tab[3], cSinPeriod, i, j, "Csin", name);
  278.                                 parseCoefficient(tab[4], sSinPeriod, i, j, "Ssin", name);
  279.                             }

  280.                         }

  281.                     } else {
  282.                         throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  283.                                                        lineNumber, name, line);
  284.                     }
  285.                 }
  286.             } catch (NumberFormatException nfe) {
  287.                 final OrekitParseException pe = new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  288.                                                                          lineNumber, name, line);
  289.                 pe.initCause(nfe);
  290.                 throw pe;
  291.             }
  292.         }

  293.         if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
  294.             // ensure at least the (0, 0) element is properly set
  295.             if (Precision.equals(c[0][0], 0.0, 0)) {
  296.                 c[0][0] = 1.0;
  297.             }
  298.         }

  299.         if (Double.isNaN(getAe()) || Double.isNaN(getMu()) || !okCoeffs) {
  300.             String loaderName = getClass().getName();
  301.             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
  302.             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
  303.                                       name, loaderName);
  304.         }

  305.         setRawCoefficients(normalized, c, s, name);
  306.         setTideSystem(tideSystem);
  307.         setReadComplete(true);

  308.     }

  309.     /** Get a provider for read spherical harmonics coefficients.
  310.      * <p>
  311.      * ICGEM fields do include time-dependent parts which are taken into account
  312.      * in the returned provider.
  313.      * </p>
  314.      * @param wantNormalized if true, the provider will provide normalized coefficients,
  315.      * otherwise it will provide un-normalized coefficients
  316.      * @param degree maximal degree
  317.      * @param order maximal order
  318.      * @return a new provider
  319.      * @exception OrekitException if the requested maximal degree or order exceeds the
  320.      * available degree or order or if no gravity field has read yet
  321.      * @since 6.0
  322.      */
  323.     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
  324.                                                      final int degree, final int order)
  325.         throws OrekitException {

  326.         RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
  327.         if (cTrend.isEmpty() && cCos.isEmpty()) {
  328.             // there are no time-dependent coefficients
  329.             return provider;
  330.         }

  331.         if (!cTrend.isEmpty()) {

  332.             // add the secular trend layer
  333.             final double[][] cArrayTrend = toArray(cTrend);
  334.             final double[][] sArrayTrend = toArray(sTrend);
  335.             rescale(1.0 / Constants.JULIAN_YEAR, normalized, cArrayTrend, sArrayTrend, wantNormalized, cArrayTrend, sArrayTrend);
  336.             provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArrayTrend, sArrayTrend);

  337.         }

  338.         for (final Map.Entry<Double, List<List<Double>>> entry : cCos.entrySet()) {

  339.             final double period = entry.getKey();

  340.             // add the pulsating layer for the current period
  341.             final double[][] cArrayCos = toArray(cCos.get(period));
  342.             final double[][] sArrayCos = toArray(sCos.get(period));
  343.             final double[][] cArraySin = toArray(cSin.get(period));
  344.             final double[][] sArraySin = toArray(sSin.get(period));
  345.             rescale(1.0, normalized, cArrayCos, sArrayCos, wantNormalized, cArrayCos, sArrayCos);
  346.             rescale(1.0, normalized, cArraySin, sArraySin, wantNormalized, cArraySin, sArraySin);
  347.             provider = new PulsatingSphericalHarmonics(provider, period * Constants.JULIAN_YEAR,
  348.                                                        cArrayCos, cArraySin, sArrayCos, sArraySin);

  349.         }

  350.         return provider;

  351.     }

  352. }