ICGEMFormatReader.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.BufferedReader;
  20. import java.io.IOException;
  21. import java.io.InputStream;
  22. import java.io.InputStreamReader;
  23. import java.text.ParseException;
  24. import java.util.ArrayList;
  25. import java.util.HashMap;
  26. import java.util.List;
  27. import java.util.Map;

  29. import org.apache.commons.math3.util.FastMath;
  30. import org.apache.commons.math3.util.Precision;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitMessages;
  33. import org.orekit.time.DateComponents;
  34. import org.orekit.utils.Constants;

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

  79.     /** Product type. */
  80.     private static final String PRODUCT_TYPE            = "product_type";

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

  85.     /** Gravity constant marker. */
  86.     private static final String GRAVITY_CONSTANT        = "earth_gravity_constant";

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

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

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

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

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

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

  106.     /** Normalization indicator. */
  107.     private static final String NORMALIZATION_INDICATOR = "norm";

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

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

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

  118.     /** Gravity field coefficient. */
  119.     private static final String GFC                     = "gfc";

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

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

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

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

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

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

  139.     /** Indicator for normalized coeffcients. */
  140.     private boolean normalized;

  142.     /** Reference date. */
  143.     private DateComponents referenceDate;

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

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

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

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

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

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

  163.     /** Simple constructor.
  164.      * @param supportedNames regular expression for supported files names
  165.      * @param missingCoefficientsAllowed if true, allows missing coefficients in the input data
  166.      */
  167.     public ICGEMFormatReader(final String supportedNames, final boolean missingCoefficientsAllowed) {
  168.         super(supportedNames, missingCoefficientsAllowed);
  169.         referenceDate = null;
  170.         cTrend = new ArrayList<List<Double>>();
  171.         sTrend = new ArrayList<List<Double>>();
  172.         cCos   = new HashMap<Double, List<List<Double>>>();
  173.         cSin   = new HashMap<Double, List<List<Double>>>();
  174.         sCos   = new HashMap<Double, List<List<Double>>>();
  175.         sSin   = new HashMap<Double, List<List<Double>>>();
  176.     }

  178.     /** {@inheritDoc} */
  179.     public void loadData(final InputStream input, final String name)
  180.         throws IOException, ParseException, OrekitException {

  182.         // reset the indicator before loading any data
  183.         setReadComplete(false);
  184.         referenceDate = null;
  185.         cTrend.clear();
  186.         sTrend.clear();
  187.         cCos.clear();
  188.         cSin.clear();
  189.         sCos.clear();
  190.         sSin.clear();

  192.         // by default, the field is normalized with unknown tide system
  193.         // (will be overridden later if non-default)
  194.         normalized = true;
  195.         tideSystem = TideSystem.UNKNOWN;

  197.         final BufferedReader r = new BufferedReader(new InputStreamReader(input, "UTF-8"));
  198.         boolean inHeader = true;
  199.         double[][] c               = null;
  200.         double[][] s               = null;
  201.         boolean okCoeffs           = false;
  202.         int lineNumber   = 0;
  203.         for (String line = r.readLine(); line != null; line = r.readLine()) {
  204.             try {
  205.                 ++lineNumber;
  206.                 if (line.trim().length() == 0) {
  207.                     continue;
  208.                 }
  209.                 final String[] tab = line.split("\\s+");
  210.                 if (inHeader) {
  211.                     if ((tab.length == 2) && PRODUCT_TYPE.equals(tab[0])) {
  212.                         if (!GRAVITY_FIELD.equals(tab[1])) {
  213.                             throw OrekitException.createParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  214.                                                                        lineNumber, name, line);
  215.                         }
  216.                     } else if ((tab.length == 2) && GRAVITY_CONSTANT.equals(tab[0])) {
  217.                         setMu(parseDouble(tab[1]));
  218.                     } else if ((tab.length == 2) && REFERENCE_RADIUS.equals(tab[0])) {
  219.                         setAe(parseDouble(tab[1]));
  220.                     } else if ((tab.length == 2) && MAX_DEGREE.equals(tab[0])) {

  222.                         final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
  223.                         final int order  = FastMath.min(getMaxParseOrder(), degree);
  224.                         c = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
  225.                         s = buildTriangularArray(degree, order, missingCoefficientsAllowed() ? 0.0 : Double.NaN);

  227.                     } else if ((tab.length == 2) && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
  228.                         if (ZERO_TIDE.equals(tab[1])) {
  229.                             tideSystem = TideSystem.ZERO_TIDE;
  230.                         } else if (TIDE_FREE.equals(tab[1])) {
  231.                             tideSystem = TideSystem.TIDE_FREE;
  232.                         } else if (TIDE_UNKNOWN.equals(tab[1])) {
  233.                             tideSystem = TideSystem.UNKNOWN;
  234.                         } else {
  235.                             throw OrekitException.createParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  236.                                                                        lineNumber, name, line);
  237.                         }
  238.                     } else if ((tab.length == 2) && NORMALIZATION_INDICATOR.equals(tab[0])) {
  239.                         if (NORMALIZED.equals(tab[1])) {
  240.                             normalized = true;
  241.                         } else if (UNNORMALIZED.equals(tab[1])) {
  242.                             normalized = false;
  243.                         } else {
  244.                             throw OrekitException.createParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  245.                                                                        lineNumber, name, line);
  246.                         }
  247.                     } else if ((tab.length == 2) && END_OF_HEADER.equals(tab[0])) {
  248.                         inHeader = false;
  249.                     }
  250.                 } else {
  251.                     if ((tab.length == 7 && GFC.equals(tab[0])) || (tab.length == 8 && GFCT.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) {

  257.                             parseCoefficient(tab[3], c, i, j, "C", name);
  258.                             parseCoefficient(tab[4], s, i, j, "S", name);
  259.                             okCoeffs = true;

  261.                             if (tab.length == 8) {
  262.                                 // check the reference date (format yyyymmdd)
  263.                                 final DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
  264.                                                                                    Integer.parseInt(tab[7].substring(4, 6)),
  265.                                                                                    Integer.parseInt(tab[7].substring(6, 8)));
  266.                                 if (referenceDate == null) {
  267.                                     // first reference found, store it
  268.                                     referenceDate = localRef;
  269.                                 } else if (!referenceDate.equals(localRef)) {
  270.                                     throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
  271.                                                               referenceDate, localRef, name);
  272.                                 }
  273.                             }

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

  278.                         final int i = Integer.parseInt(tab[1]);
  279.                         final int j = Integer.parseInt(tab[2]);
  280.                         if (i < c.length && j < c[i].length) {

  282.                             // store the secular trend coefficients
  283.                             extendListOfLists(cTrend, i, j, 0.0);
  284.                             extendListOfLists(sTrend, i, j, 0.0);
  285.                             parseCoefficient(tab[3], cTrend, i, j, "Ctrend", name);
  286.                             parseCoefficient(tab[4], sTrend, i, j, "Strend", name);

  288.                         }

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

  292.                         final int i = Integer.parseInt(tab[1]);
  293.                         final int j = Integer.parseInt(tab[2]);
  294.                         if (i < c.length && j < c[i].length) {

  296.                             // extract arrays corresponding to period
  297.                             final Double period = Double.valueOf(tab[7]);
  298.                             if (!cCos.containsKey(period)) {
  299.                                 cCos.put(period, new ArrayList<List<Double>>());
  300.                                 cSin.put(period, new ArrayList<List<Double>>());
  301.                                 sCos.put(period, new ArrayList<List<Double>>());
  302.                                 sSin.put(period, new ArrayList<List<Double>>());
  303.                             }
  304.                             final List<List<Double>> cCosPeriod = cCos.get(period);
  305.                             final List<List<Double>> cSinPeriod = cSin.get(period);
  306.                             final List<List<Double>> sCosPeriod = sCos.get(period);
  307.                             final List<List<Double>> sSinPeriod = sSin.get(period);

  309.                             // store the pulsation coefficients
  310.                             extendListOfLists(cCosPeriod, i, j, 0.0);
  311.                             extendListOfLists(cSinPeriod, i, j, 0.0);
  312.                             extendListOfLists(sCosPeriod, i, j, 0.0);
  313.                             extendListOfLists(sSinPeriod, i, j, 0.0);
  314.                             if (ACOS.equals(tab[0])) {
  315.                                 parseCoefficient(tab[3], cCosPeriod, i, j, "Ccos", name);
  316.                                 parseCoefficient(tab[4], sCosPeriod, i, j, "SCos", name);
  317.                             } else {
  318.                                 parseCoefficient(tab[3], cSinPeriod, i, j, "Csin", name);
  319.                                 parseCoefficient(tab[4], sSinPeriod, i, j, "Ssin", name);
  320.                             }

  322.                         }

  324.                     } else {
  325.                         throw OrekitException.createParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  326.                                                                    lineNumber, name, line);
  327.                     }
  328.                 }
  329.             } catch (NumberFormatException nfe) {
  330.                 final ParseException pe = OrekitException.createParseException(
  331.                         OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, lineNumber, name, line);
  332.                 pe.initCause(nfe);
  333.                 throw pe;
  334.             }
  335.         }

  337.         if (missingCoefficientsAllowed() && c.length > 0 && c[0].length > 0) {
  338.             // ensure at least the (0, 0) element is properly set
  339.             if (Precision.equals(c[0][0], 0.0, 1)) {
  340.                 c[0][0] = 1.0;
  341.             }
  342.         }

  344.         if (Double.isNaN(getAe()) || Double.isNaN(getMu()) || !okCoeffs) {
  345.             String loaderName = getClass().getName();
  346.             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
  347.             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
  348.                                       name, loaderName);
  349.         }

  351.         setRawCoefficients(normalized, c, s, name);
  352.         setTideSystem(tideSystem);
  353.         setReadComplete(true);

  355.     }

  357.     /** Get a provider for read spherical harmonics coefficients.
  358.      * <p>
  359.      * ICGEM fields do include time-dependent parts which are taken into account
  360.      * in the returned provider.
  361.      * </p>
  362.      * @param wantNormalized if true, the provider will provide normalized coefficients,
  363.      * otherwise it will provide un-normalized coefficients
  364.      * @param degree maximal degree
  365.      * @param order maximal order
  366.      * @return a new provider
  367.      * @exception OrekitException if the requested maximal degree or order exceeds the
  368.      * available degree or order or if no gravity field has read yet
  369.      * @since 6.0
  370.      */
  371.     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
  372.                                                      final int degree, final int order)
  373.         throws OrekitException {

  375.         RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
  376.         if (cTrend.isEmpty() && cCos.isEmpty()) {
  377.             // there are no time-dependent coefficients
  378.             return provider;
  379.         }

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

  383.             // add the secular trend layer
  384.             final double[][] cArrayTrend = toArray(cTrend);
  385.             final double[][] sArrayTrend = toArray(sTrend);
  386.             rescale(1.0 / Constants.JULIAN_YEAR, normalized, cArrayTrend, sArrayTrend, wantNormalized, cArrayTrend, sArrayTrend);
  387.             provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArrayTrend, sArrayTrend);

  389.         }

  391.         for (final Double period : cCos.keySet()) {

  393.             // add the pulsating layer for the current period
  394.             final double[][] cArrayCos = toArray(cCos.get(period));
  395.             final double[][] sArrayCos = toArray(sCos.get(period));
  396.             final double[][] cArraySin = toArray(cSin.get(period));
  397.             final double[][] sArraySin = toArray(sSin.get(period));
  398.             rescale(1.0, normalized, cArrayCos, sArrayCos, wantNormalized, cArrayCos, sArrayCos);
  399.             rescale(1.0, normalized, cArraySin, sArraySin, wantNormalized, cArraySin, sArraySin);
  400.             provider = new PulsatingSphericalHarmonics(provider, period * Constants.JULIAN_YEAR,
  401.                                                        cArrayCos, cArraySin, sArrayCos, sArraySin);

  403.         }

  405.         return provider;

  407.     }

  409. }
Created with JaCoCo 0.7.2.201409121644