ICGEMFormatReader.java

  1. /* Copyright 2002-2020 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.BufferedReader;
  20. import java.io.IOException;
  21. import java.io.InputStream;
  22. import java.io.InputStreamReader;
  23. import java.nio.charset.StandardCharsets;
  24. import java.text.ParseException;
  25. import java.util.ArrayList;
  26. import java.util.HashMap;
  27. import java.util.List;
  28. import java.util.Map;
  29. import java.util.regex.Pattern;

  31. import org.hipparchus.util.FastMath;
  32. import org.hipparchus.util.Precision;
  33. import org.orekit.annotation.DefaultDataContext;
  34. import org.orekit.data.DataContext;
  35. import org.orekit.errors.OrekitException;
  36. import org.orekit.errors.OrekitMessages;
  37. import org.orekit.errors.OrekitParseException;
  38. import org.orekit.time.AbsoluteDate;
  39. import org.orekit.time.DateComponents;
  40. import org.orekit.time.TimeScale;
  41. import org.orekit.utils.Constants;

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

  93.     /** Product type. */
  94.     private static final String PRODUCT_TYPE            = "product_type";

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

  99.     /** Gravity constant marker. */
  100.     private static final String GRAVITY_CONSTANT        = "gravity_constant";

  102.     /** Reference radius. */
  103.     private static final String REFERENCE_RADIUS        = "radius";

  105.     /** Max degree. */
  106.     private static final String MAX_DEGREE              = "max_degree";

  108.     /** Tide system indicator. */
  109.     private static final String TIDE_SYSTEM_INDICATOR   = "tide_system";

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

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

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

  120.     /** Normalization indicator. */
  121.     private static final String NORMALIZATION_INDICATOR = "norm";

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

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

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

  132.     /** Gravity field coefficient. */
  133.     private static final String GFC                     = "gfc";

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

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

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

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

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

  150.     /** Pattern for delimiting regular expressions. */
  151.     private static final Pattern SEPARATOR = Pattern.compile("\\s+");

  153.     /** Indicator for normalized coefficients. */
  154.     private boolean normalized;

  156.     /** Reference date. */
  157.     private AbsoluteDate referenceDate;

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

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

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

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

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

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

  177.     /** Simple constructor.
  178.      *
  179.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  180.      *
  181.      * @param supportedNames regular expression for supported files names
  182.      * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
  183.      * @see #ICGEMFormatReader(String, boolean, TimeScale)
  184.      */
  185.     @DefaultDataContext
  186.     public ICGEMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
  187.         this(supportedNames, missingCoefficientsAllowed,
  188.                 DataContext.getDefault().getTimeScales().getTT());
  189.     }

  191.     /**
  192.      * Simple constructor.
  193.      *
  194.      * @param supportedNames             regular expression for supported files names
  195.      * @param missingCoefficientsAllowed if true, allows missing coefficients in the input
  196.      *                                   data
  197.      * @param timeScale                  to use when parsing dates.
  198.      * @since 10.1
  199.      */
  200.     public ICGEMFormatReader(final String supportedNames,
  201.                              final boolean missingCoefficientsAllowed,
  202.                              final TimeScale timeScale) {
  203.         super(supportedNames, missingCoefficientsAllowed, timeScale);
  204.         referenceDate = null;
  205.         cTrend = new ArrayList<>();
  206.         sTrend = new ArrayList<>();
  207.         cCos   = new HashMap<>();
  208.         cSin   = new HashMap<>();
  209.         sCos   = new HashMap<>();
  210.         sSin   = new HashMap<>();
  211.     }

  213.     /** {@inheritDoc} */
  214.     public void loadData(final InputStream input, final String name)
  215.         throws IOException, ParseException, OrekitException {

  217.         // reset the indicator before loading any data
  218.         setReadComplete(false);
  219.         referenceDate = null;
  220.         cTrend.clear();
  221.         sTrend.clear();
  222.         cCos.clear();
  223.         cSin.clear();
  224.         sCos.clear();
  225.         sSin.clear();

  227.         // by default, the field is normalized with unknown tide system
  228.         // (will be overridden later if non-default)
  229.         normalized = true;
  230.         TideSystem tideSystem = TideSystem.UNKNOWN;

  232.         boolean inHeader = true;
  233.         double[][] c     = null;
  234.         double[][] s     = null;
  235.         boolean okCoeffs = false;
  236.         int lineNumber   = 0;
  237.         String line      = null;
  238.         try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
  239.             for (line = r.readLine(); line != null; line = r.readLine()) {
  240.                 ++lineNumber;
  241.                 line = line.trim();
  242.                 if (line.length() == 0) {
  243.                     continue;
  244.                 }
  245.                 final String[] tab = SEPARATOR.split(line);
  246.                 if (inHeader) {
  247.                     if ((tab.length == 2) && PRODUCT_TYPE.equals(tab[0])) {
  248.                         if (!GRAVITY_FIELD.equals(tab[1])) {
  249.                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  250.                                                            lineNumber, name, line);
  251.                         }
  252.                     } else if ((tab.length == 2) && tab[0].endsWith(GRAVITY_CONSTANT)) {
  253.                         setMu(parseDouble(tab[1]));
  254.                     } else if ((tab.length == 2) && REFERENCE_RADIUS.equals(tab[0])) {
  255.                         setAe(parseDouble(tab[1]));
  256.                     } else if ((tab.length == 2) && MAX_DEGREE.equals(tab[0])) {

  258.                         final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
  259.                         final int order  = FastMath.min(getMaxParseOrder(), degree);
  260.                         c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
  261.                         s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);

  263.                     } else if ((tab.length == 2) && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
  264.                         if (ZERO_TIDE.equals(tab[1])) {
  265.                             tideSystem = TideSystem.ZERO_TIDE;
  266.                         } else if (TIDE_FREE.equals(tab[1])) {
  267.                             tideSystem = TideSystem.TIDE_FREE;
  268.                         } else if (TIDE_UNKNOWN.equals(tab[1])) {
  269.                             tideSystem = TideSystem.UNKNOWN;
  270.                         } else {
  271.                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  272.                                                            lineNumber, name, line);
  273.                         }
  274.                     } else if ((tab.length == 2) && NORMALIZATION_INDICATOR.equals(tab[0])) {
  275.                         if (NORMALIZED.equals(tab[1])) {
  276.                             normalized = true;
  277.                         } else if (UNNORMALIZED.equals(tab[1])) {
  278.                             normalized = false;
  279.                         } else {
  280.                             throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  281.                                                            lineNumber, name, line);
  282.                         }
  283.                     } else if ((tab.length == 2) && END_OF_HEADER.equals(tab[0])) {
  284.                         inHeader = false;
  285.                     }
  286.                 } else {
  287.                     if ((tab.length == 7 && GFC.equals(tab[0])) || (tab.length == 8 && GFCT.equals(tab[0]))) {

  289.                         final int i = Integer.parseInt(tab[1]);
  290.                         final int j = Integer.parseInt(tab[2]);
  291.                         if (i < c.length && j < c[i].length) {

  293.                             parseCoefficient(tab[3], c, i, j, "C", name);
  294.                             parseCoefficient(tab[4], s, i, j, "S", name);
  295.                             okCoeffs = true;

  297.                             if (tab.length == 8) {
  298.                                 // check the reference date (format yyyymmdd)
  299.                                 final DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
  300.                                                                                    Integer.parseInt(tab[7].substring(4, 6)),
  301.                                                                                    Integer.parseInt(tab[7].substring(6, 8)));
  302.                                 if (referenceDate == null) {
  303.                                     // first reference found, store it
  304.                                     referenceDate = toDate(localRef);
  305.                                 } else if (!referenceDate.equals(toDate(localRef))) {
  306.                                     throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
  307.                                                               referenceDate, toDate(localRef), name);
  308.                                 }
  309.                             }

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

  314.                         final int i = Integer.parseInt(tab[1]);
  315.                         final int j = Integer.parseInt(tab[2]);
  316.                         if (i < c.length && j < c[i].length) {

  318.                             // store the secular trend coefficients
  319.                             extendListOfLists(cTrend, i, j, 0.0);
  320.                             extendListOfLists(sTrend, i, j, 0.0);
  321.                             parseCoefficient(tab[3], cTrend, i, j, "Ctrend", name);
  322.                             parseCoefficient(tab[4], sTrend, i, j, "Strend", name);

  324.                         }

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

  328.                         final int i = Integer.parseInt(tab[1]);
  329.                         final int j = Integer.parseInt(tab[2]);
  330.                         if (i < c.length && j < c[i].length) {

  332.                             // extract arrays corresponding to period
  333.                             final Double period = Double.valueOf(tab[7]);
  334.                             if (!cCos.containsKey(period)) {
  335.                                 cCos.put(period, new ArrayList<>());
  336.                                 cSin.put(period, new ArrayList<>());
  337.                                 sCos.put(period, new ArrayList<>());
  338.                                 sSin.put(period, new ArrayList<>());
  339.                             }
  340.                             final List<List<Double>> cCosPeriod = cCos.get(period);
  341.                             final List<List<Double>> cSinPeriod = cSin.get(period);
  342.                             final List<List<Double>> sCosPeriod = sCos.get(period);
  343.                             final List<List<Double>> sSinPeriod = sSin.get(period);

  345.                             // store the pulsation coefficients
  346.                             extendListOfLists(cCosPeriod, i, j, 0.0);
  347.                             extendListOfLists(cSinPeriod, i, j, 0.0);
  348.                             extendListOfLists(sCosPeriod, i, j, 0.0);
  349.                             extendListOfLists(sSinPeriod, i, j, 0.0);
  350.                             if (ACOS.equals(tab[0])) {
  351.                                 parseCoefficient(tab[3], cCosPeriod, i, j, "Ccos", name);
  352.                                 parseCoefficient(tab[4], sCosPeriod, i, j, "SCos", name);
  353.                             } else {
  354.                                 parseCoefficient(tab[3], cSinPeriod, i, j, "Csin", name);
  355.                                 parseCoefficient(tab[4], sSinPeriod, i, j, "Ssin", name);
  356.                             }

  358.                         }

  360.                     } else {
  361.                         throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  362.                                                        lineNumber, name, line);
  363.                     }
  364.                 }

  366.             }
  367.         } catch (NumberFormatException nfe) {
  368.             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  369.                                       lineNumber, name, line);
  370.         }


  373.         if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
  374.             // ensure at least the (0, 0) element is properly set
  375.             if (Precision.equals(c[0][0], 0.0, 0)) {
  376.                 c[0][0] = 1.0;
  377.             }
  378.         }

  380.         if (Double.isNaN(getAe()) || Double.isNaN(getMu()) || !okCoeffs) {
  381.             String loaderName = getClass().getName();
  382.             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
  383.             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
  384.                                       name, loaderName);
  385.         }

  387.         setRawCoefficients(normalized, c, s, name);
  388.         setTideSystem(tideSystem);
  389.         setReadComplete(true);

  391.     }

  393.     /** Get a provider for read spherical harmonics coefficients.
  394.      * <p>
  395.      * ICGEM fields do include time-dependent parts which are taken into account
  396.      * in the returned provider.
  397.      * </p>
  398.      * @param wantNormalized if true, the provider will provide normalized coefficients,
  399.      * otherwise it will provide un-normalized coefficients
  400.      * @param degree maximal degree
  401.      * @param order maximal order
  402.      * @return a new provider
  403.      * @since 6.0
  404.      */
  405.     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
  406.                                                      final int degree, final int order) {

  408.         RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
  409.         if (cTrend.isEmpty() && cCos.isEmpty()) {
  410.             // there are no time-dependent coefficients
  411.             return provider;
  412.         }

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

  416.             // add the secular trend layer
  417.             final double[][] cArrayTrend = toArray(cTrend);
  418.             final double[][] sArrayTrend = toArray(sTrend);
  419.             rescale(1.0 / Constants.JULIAN_YEAR, normalized, cArrayTrend, sArrayTrend, wantNormalized, cArrayTrend, sArrayTrend);
  420.             provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArrayTrend, sArrayTrend);

  422.         }

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

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

  428.             // add the pulsating layer for the current period
  429.             final double[][] cArrayCos = toArray(cCos.get(period));
  430.             final double[][] sArrayCos = toArray(sCos.get(period));
  431.             final double[][] cArraySin = toArray(cSin.get(period));
  432.             final double[][] sArraySin = toArray(sSin.get(period));
  433.             rescale(1.0, normalized, cArrayCos, sArrayCos, wantNormalized, cArrayCos, sArrayCos);
  434.             rescale(1.0, normalized, cArraySin, sArraySin, wantNormalized, cArraySin, sArraySin);
  435.             provider = new PulsatingSphericalHarmonics(provider, period * Constants.JULIAN_YEAR,
  436.                                                        cArrayCos, cArraySin, sArrayCos, sArraySin);

  438.         }

  440.         return provider;

  442.     }

  444. }
Created with JaCoCo 0.8.6.202009150832