ICGEMFormatReader.java

  1. /* Copyright 2002-2025 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 org.hipparchus.util.FastMath;
  20. import org.hipparchus.util.MathUtils;
  21. import org.hipparchus.util.Precision;
  22. import org.orekit.annotation.DefaultDataContext;
  23. import org.orekit.data.DataContext;
  24. import org.orekit.errors.OrekitException;
  25. import org.orekit.errors.OrekitMessages;
  26. import org.orekit.time.AbsoluteDate;
  27. import org.orekit.time.DateComponents;
  28. import org.orekit.time.TimeComponents;
  29. import org.orekit.time.TimeScale;
  30. import org.orekit.utils.Constants;
  31. import org.orekit.utils.TimeSpanMap;

  33. import java.io.BufferedReader;
  34. import java.io.IOException;
  35. import java.io.InputStream;
  36. import java.io.InputStreamReader;
  37. import java.nio.charset.StandardCharsets;
  38. import java.text.ParseException;
  39. import java.util.ArrayList;
  40. import java.util.List;
  41. import java.util.Locale;
  42. import java.util.regex.Matcher;
  43. import java.util.regex.Pattern;

  45. /** Reader for the ICGEM gravity field format.
  46.  *
  47.  * <p>This format is used to describe the gravity field of EIGEN models
  48.  * published by the GFZ Potsdam since 2004. It is described in Franz
  49.  * Barthelmes and Christoph F&ouml;rste paper: "the ICGEM-format".
  50.  * The 2006-02-28 version of this paper can be found <a
  51.  * href="http://op.gfz-potsdam.de/grace/results/grav/g005_ICGEM-Format.pdf">here</a>
  52.  * and the 2011-06-07 version of this paper can be found <a
  53.  * href="http://icgem.gfz-potsdam.de/ICGEM-Format-2011.pdf">here</a>.
  54.  * These versions differ in time-dependent coefficients, which are linear-only prior
  55.  * to 2011 (up to eigen-5 model) and have also harmonic effects after that date
  56.  * (starting with eigen-6 model). A third (undocumented as of 2018-05-14) version
  57.  * of the file format also adds a time-span for time-dependent coefficients, allowing
  58.  * for piecewise models. All three versions are supported by the class.</p>
  59.  * <p>
  60.  * This reader uses relaxed check on the gravity constant key so any key ending
  61.  * in gravity_constant is accepted and not only earth_gravity_constant as specified
  62.  * in the previous documents. This allows to read also non Earth gravity fields
  63.  * as found in <a href="http://icgem.gfz-potsdam.de/tom_celestial">ICGEM
  64.  * - Gravity Field Models of other Celestial Bodies</a> page to be read.
  65.  * </p>
  66.  *
  67.  * <p> The proper way to use this class is to call the {@link GravityFieldFactory}
  68.  *  which will determine which reader to use with the selected gravity field file.</p>
  69.  *
  70.  * @see GravityFields
  71.  * @author Luc Maisonobe
  72.  */
  73. public class ICGEMFormatReader extends PotentialCoefficientsReader {

  75.     /** Format. */
  76.     private static final String FORMAT                  = "format";

  78.     /** Maximum supported formats. */
  79.     private static final double MAX_FORMAT              = 2.0;

  81.     /** Product type. */
  82.     private static final String PRODUCT_TYPE            = "product_type";

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

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

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

  93.     /** Max degree. */
  94.     private static final String MAX_DEGREE              = "max_degree";

  96.     /** Errors indicator. */
  97.     private static final String ERRORS                  = "errors";

  99.     /** Tide system indicator. */
  100.     private static final String TIDE_SYSTEM_INDICATOR   = "tide_system";

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

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

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

  111.     /** Normalization indicator. */
  112.     private static final String NORMALIZATION_INDICATOR = "norm";

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

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

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

  123.     /** Gravity field coefficient. */
  124.     private static final String GFC                     = "gfc";

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

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

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

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

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

  141.     /** Name of base coefficients. */
  142.     private static final String BASE_NAMES              = "C/S";

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

  147.     /** Pattern for supported formats. */
  148.     private static final Pattern SUPPORTED_FORMAT = Pattern.compile("icgem(\\d+\\.\\d+)");

  150.     /** Flag for Gravitational coefficient. */
  151.     private static final int MU = 0x1;

  153.     /** Flag for scaling radius. */
  154.     private static final int AE = 0x2;

  156.     /** Flag for degree/order. */
  157.     private static final int LIMITS = 0x4;

  159.     /** Flag for errors. */
  160.     private static final int ERR = 0x8;

  162.     /** Flag for coefficients. */
  163.     private static final int COEFFS = 0x10;

  165.     /** Indicator for normalized coefficients. */
  166.     private boolean normalized;

  168.     /** Reference dates. */
  169.     private List<AbsoluteDate> referenceDates;

  171.     /** Pulsations. */
  172.     private List<Double>       pulsations;

  174.     /** Time map of the harmonics. */
  175.     private TimeSpanMap<Container> containers;

  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.     }

  206.     /** {@inheritDoc} */
  207.     public void loadData(final InputStream input, final String name)
  208.         throws IOException, ParseException, OrekitException {

  210.         // reset the indicator before loading any data
  211.         setReadComplete(false);
  212.         containers     = null;
  213.         referenceDates = new ArrayList<>();
  214.         pulsations     = new ArrayList<>();

  216.         // by default, the field is normalized with unknown tide system
  217.         // (will be overridden later if non-default)
  218.         normalized            = true;
  219.         TideSystem tideSystem = TideSystem.UNKNOWN;
  220.         Errors     errors     = Errors.NO;

  222.         double    version        = 1.0;
  223.         boolean   inHeader       = true;
  224.         Flattener flattener      = null;
  225.         int       flags          = 0;
  226.         double[]  c0             = null;
  227.         double[]  s0             = null;
  228.         int       lineNumber     = 0;
  229.         String    line           = null;
  230.         try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
  231.             for (line = r.readLine(); line != null; line = r.readLine()) {
  232.                 boolean parseError = false;
  233.                 ++lineNumber;
  234.                 line = line.trim();
  235.                 if (line.isEmpty()) {
  236.                     continue;
  237.                 }
  238.                 final String[] tab = SEPARATOR.split(line);
  239.                 if (inHeader) {
  240.                     if (tab.length >= 2 && FORMAT.equals(tab[0])) {
  241.                         final Matcher matcher = SUPPORTED_FORMAT.matcher(tab[1]);
  242.                         if (matcher.matches()) {
  243.                             version = Double.parseDouble(matcher.group(1));
  244.                             if (version > MAX_FORMAT) {
  245.                                 parseError = true;
  246.                             }
  247.                         } else {
  248.                             parseError = true;
  249.                         }
  250.                     } else if (tab.length >= 2 && PRODUCT_TYPE.equals(tab[0])) {
  251.                         parseError = !GRAVITY_FIELD.equals(tab[1]);
  252.                     } else if (tab.length >= 2 && tab[0].endsWith(GRAVITY_CONSTANT)) {
  253.                         setMu(parseDouble(tab[1]));
  254.                         flags |= MU;
  255.                     } else if (tab.length >= 2 && REFERENCE_RADIUS.equals(tab[0])) {
  256.                         setAe(parseDouble(tab[1]));
  257.                         flags |= AE;
  258.                     } else if (tab.length >= 2 && MAX_DEGREE.equals(tab[0])) {

  260.                         final int degree = FastMath.min(getMaxParseDegree(), Integer.parseInt(tab[1]));
  261.                         final int order  = FastMath.min(getMaxParseOrder(), degree);
  262.                         flattener  = new Flattener(degree, order);
  263.                         c0         = buildFlatArray(flattener, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
  264.                         s0         = buildFlatArray(flattener, missingCoefficientsAllowed() ? 0.0 : Double.NaN);
  265.                         flags     |= LIMITS;

  267.                     } else if (tab.length >= 2 && ERRORS.equals(tab[0])) {
  268.                         try {
  269.                             errors = Errors.valueOf(tab[1].toUpperCase(Locale.US));
  270.                             flags |= ERR;
  271.                         } catch (IllegalArgumentException iae) {
  272.                             parseError = true;
  273.                         }
  274.                     } else if (tab.length >= 2 && TIDE_SYSTEM_INDICATOR.equals(tab[0])) {
  275.                         if (ZERO_TIDE.equals(tab[1])) {
  276.                             tideSystem = TideSystem.ZERO_TIDE;
  277.                         } else if (TIDE_FREE.equals(tab[1])) {
  278.                             tideSystem = TideSystem.TIDE_FREE;
  279.                         } else if (TIDE_UNKNOWN.equals(tab[1])) {
  280.                             tideSystem = TideSystem.UNKNOWN;
  281.                         } else {
  282.                             parseError = true;
  283.                         }
  284.                     } else if (tab.length >= 2 && NORMALIZATION_INDICATOR.equals(tab[0])) {
  285.                         if (NORMALIZED.equals(tab[1])) {
  286.                             normalized = true;
  287.                         } else if (UNNORMALIZED.equals(tab[1])) {
  288.                             normalized = false;
  289.                         } else {
  290.                             parseError = true;
  291.                         }
  292.                     } else if (tab.length >= 1 && END_OF_HEADER.equals(tab[0])) {
  293.                         inHeader   = false;
  294.                     }
  295.                     if (parseError) {
  296.                         throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  297.                                                   lineNumber, name, line);
  298.                     }
  299.                 } else if (dataKeyKnown(tab) && tab.length > 2) {

  301.                     final int n = Integer.parseInt(tab[1]);
  302.                     final int m = Integer.parseInt(tab[2]);
  303.                     flags |= COEFFS;
  304.                     if (!flattener.withinRange(n, m)) {
  305.                         // just ignore coefficients we don't need
  306.                         continue;
  307.                     }

  309.                     if (tab.length > 4 && GFC.equals(tab[0])) {
  310.                         // fixed coefficient

  312.                         parseCoefficient(tab[3], flattener, c0, n, m, "C", name);
  313.                         parseCoefficient(tab[4], flattener, s0, n, m, "S", name);

  315.                     } else if (version < 2.0 && tab.length > 5 + errors.fields && GFCT.equals(tab[0])) {
  316.                         // base of linear coefficient with infinite time range

  318.                         if (containers == null) {
  319.                             // prepare the single container (it will be populated when next lines are parsed)
  320.                             containers = new TimeSpanMap<>(new Container(flattener));
  321.                         }

  323.                         // set the constant coefficients to 0 as they will be managed by the piecewise model
  324.                         final int globalIndex = flattener.index(n, m);
  325.                         c0[globalIndex]       = 0.0;
  326.                         s0[globalIndex]       = 0.0;

  328.                         // store the single reference date valid for the whole field
  329.                         final AbsoluteDate lineDate = parseDate(tab[5 + errors.fields]);
  330.                         final int referenceIndex    = referenceDateIndex(referenceDates, lineDate);
  331.                         if (referenceIndex != 0) {
  332.                             // we already know the reference date, check this lines does not define a new one
  333.                             throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
  334.                                                       referenceDates.get(0), lineDate, name,
  335.                                                       lineDate.durationFrom(referenceDates.get(0)));
  336.                         }

  338.                         final Container single = containers.getFirstSpan().getData();
  339.                         final int       index  = single.flattener.index(n, m);
  340.                         if (single.components[index] != null) {
  341.                             throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  342.                                                       BASE_NAMES, n, m, name);
  343.                         }
  344.                         single.components[index] = new TimeDependentHarmonic(referenceIndex, parseDouble(tab[3]), parseDouble(tab[4]));


  347.                     } else if (version >= 2.0 && tab.length > 6 + errors.fields && GFCT.equals(tab[0])) {
  348.                         // base of linear coefficient with limited time range

  350.                         if (containers == null) {
  351.                             // prepare empty map to hold containers as they are parsed
  352.                             containers = new TimeSpanMap<>(null);
  353.                         }

  355.                         // set the constant coefficients to 0 as they will be managed by the piecewise model
  356.                         final int globalIndex = flattener.index(n, m);
  357.                         c0[globalIndex]       = 0.0;
  358.                         s0[globalIndex]       = 0.0;

  360.                         final AbsoluteDate t0 = parseDate(tab[5 + errors.fields]);
  361.                         final AbsoluteDate t1 = parseDate(tab[6 + errors.fields]);

  363.                         // get the containers active for the specified time range
  364.                         final List<TimeSpanMap.Span<Container>> active = getActive(flattener, t0, t1);
  365.                         for (final TimeSpanMap.Span<Container> span : active) {
  366.                             final Container             container = span.getData();
  367.                             final int                   index     = container.flattener.index(n, m);
  368.                             if (container.components[index] != null) {
  369.                                 throw new OrekitException(OrekitMessages.DUPLICATED_GRAVITY_FIELD_COEFFICIENT_IN_FILE,
  370.                                                           BASE_NAMES, n, m, name);
  371.                             }
  372.                             container.components[index] = new TimeDependentHarmonic(referenceDateIndex(referenceDates, t0),
  373.                                                                                     parseDouble(tab[3]), parseDouble(tab[4]));
  374.                         }

  376.                     } else if (version < 2.0 && tab.length > 4 && (DOT.equals(tab[0]) || TRND.equals(tab[0]))) {
  377.                         // slope of linear coefficient with infinite range

  379.                         // store the secular trend coefficients
  380.                         final Container single = containers.getFirstSpan().getData();
  381.                         final TimeDependentHarmonic harmonic = single.components[single.flattener.index(n, m)];
  382.                         if (harmonic == null) {
  383.                             parseError = true;
  384.                         } else {
  385.                             harmonic.setTrend(parseDouble(tab[3]) / Constants.JULIAN_YEAR,
  386.                                               parseDouble(tab[4]) / Constants.JULIAN_YEAR);
  387.                         }

  389.                     } else if (version >= 2.0 && tab.length > 6 + errors.fields && TRND.equals(tab[0])) {
  390.                         // slope of linear coefficient with limited range

  392.                         final AbsoluteDate t0 = parseDate(tab[5 + errors.fields]);
  393.                         final AbsoluteDate t1 = parseDate(tab[6 + errors.fields]);

  395.                         // get the containers active for the specified time range
  396.                         final List<TimeSpanMap.Span<Container>> active = getActive(flattener, t0, t1);
  397.                         for (final TimeSpanMap.Span<Container> span : active) {
  398.                             final Container             container = span.getData();
  399.                             final int                   index     = container.flattener.index(n, m);
  400.                             if (container.components[index] == null) {
  401.                                 parseError = true;
  402.                                 break;
  403.                             } else {
  404.                                 container.components[index].setTrend(parseDouble(tab[3]) / Constants.JULIAN_YEAR,
  405.                                                                      parseDouble(tab[4]) / Constants.JULIAN_YEAR);
  406.                             }
  407.                         }

  409.                     } else if (version < 2.0 && tab.length > 5 + errors.fields && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {
  410.                         // periodic coefficient with infinite range

  412.                         final double period = parseDouble(tab[5 + errors.fields]) * Constants.JULIAN_YEAR;
  413.                         final int    pIndex = pulsationIndex(pulsations, MathUtils.TWO_PI / period);

  415.                         // store the periodic effects coefficients
  416.                         final Container single = containers.getFirstSpan().getData();
  417.                         final TimeDependentHarmonic harmonic = single.components[single.flattener.index(n, m)];
  418.                         if (harmonic == null) {
  419.                             parseError = true;
  420.                         } else {
  421.                             if (ACOS.equals(tab[0])) {
  422.                                 harmonic.addCosine(-1, pIndex, parseDouble(tab[3]), parseDouble(tab[4]));
  423.                             } else {
  424.                                 harmonic.addSine(-1, pIndex, parseDouble(tab[3]), parseDouble(tab[4]));
  425.                             }
  426.                         }

  428.                     } else if (version >= 2.0 && tab.length > 7 + errors.fields && (ASIN.equals(tab[0]) || ACOS.equals(tab[0]))) {
  429.                         // periodic coefficient with limited range

  431.                         final AbsoluteDate t0      = parseDate(tab[5 + errors.fields]);
  432.                         final AbsoluteDate t1      = parseDate(tab[6 + errors.fields]);
  433.                         final int          tIndex  = referenceDateIndex(referenceDates, t0);
  434.                         final double       period  = parseDouble(tab[7 + errors.fields]) * Constants.JULIAN_YEAR;
  435.                         final int          pIndex  = pulsationIndex(pulsations, MathUtils.TWO_PI / period);

  437.                         // get the containers active for the specified time range
  438.                         final List<TimeSpanMap.Span<Container>> active = getActive(flattener, t0, t1);
  439.                         for (final TimeSpanMap.Span<Container> span : active) {
  440.                             final Container             container = span.getData();
  441.                             final int                   index     = container.flattener.index(n, m);
  442.                             if (container.components[index] == null) {
  443.                                 parseError = true;
  444.                                 break;
  445.                             } else {
  446.                                 if (ASIN.equals(tab[0])) {
  447.                                     container.components[index].addSine(tIndex, pIndex,
  448.                                                                         parseDouble(tab[3]), parseDouble(tab[4]));
  449.                                 } else {
  450.                                     container.components[index].addCosine(tIndex, pIndex,
  451.                                                                           parseDouble(tab[3]), parseDouble(tab[4]));
  452.                                 }
  453.                             }
  454.                         }

  456.                     } else {
  457.                         parseError = true;
  458.                     }

  460.                 } else if (dataKeyKnown(tab)) {
  461.                     // this was an expected data key, but the line is truncated
  462.                     parseError = true;
  463.                 }

  465.                 if (parseError) {
  466.                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  467.                                               lineNumber, name, line);
  468.                 }

  470.             }

  472.         } catch (NumberFormatException nfe) {
  473.             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  474.                                       lineNumber, name, line);
  475.         }

  477.         if (flags != (MU | AE | LIMITS | ERR | COEFFS)) {
  478.             String loaderName = getClass().getName();
  479.             loaderName = loaderName.substring(loaderName.lastIndexOf('.') + 1);
  480.             throw new OrekitException(OrekitMessages.UNEXPECTED_FILE_FORMAT_ERROR_FOR_LOADER,
  481.                                       name, loaderName);
  482.         }

  484.         if (missingCoefficientsAllowed()) {
  485.             // ensure at least the (0, 0) element is properly set
  486.             if (Precision.equals(c0[flattener.index(0, 0)], 0.0, 0)) {
  487.                 c0[flattener.index(0, 0)] = 1.0;
  488.             }
  489.         }

  491.         setRawCoefficients(normalized, flattener, c0, s0, name);
  492.         setTideSystem(tideSystem);
  493.         setReadComplete(true);

  495.     }

  497.     /** Check if a line starts with a known data key.
  498.      * @param tab line fields
  499.      * @return true if line starts with a known data key
  500.      * @since 11.1
  501.      */
  502.     private boolean dataKeyKnown(final String[] tab) {
  503.         return tab.length > 0 &&
  504.                (GFC.equals(tab[0])  || GFCT.equals(tab[0]) ||
  505.                 DOT.equals(tab[0])  || TRND.equals(tab[0]) ||
  506.                 ASIN.equals(tab[0]) || ACOS.equals(tab[0]));
  507.     }

  509.     /** Get the spans with containers active over a time range.
  510.      * @param flattener converter from triangular form to flat form
  511.      * @param t0 start of time span
  512.      * @param t1 end of time span
  513.      * @return span active from {@code t0} to {@code t1}
  514.      * @since 11.1
  515.      */
  516.     private List<TimeSpanMap.Span<Container>> getActive(final Flattener flattener,
  517.                                                         final AbsoluteDate t0, final AbsoluteDate t1) {

  519.         final List<TimeSpanMap.Span<Container>> active = new ArrayList<>();

  521.         TimeSpanMap.Span<Container> span = containers.getSpan(t0);
  522.         if (span.getStart().isBefore(t0)) {
  523.             if (span.getEnd().isAfterOrEqualTo(t1)) {
  524.                 if (span.getData() == null) {
  525.                     // the specified time range lies on an empty range
  526.                     span = containers.addValidBetween(new Container(flattener), t0, t1);
  527.                 } else {
  528.                     // the specified time range splits an existing container in three parts
  529.                     containers.addValidAfter(copyContainer(span.getData(), flattener), t1, false);
  530.                     span = containers.addValidAfter(copyContainer(span.getData(), flattener), t0, false);
  531.                 }
  532.             } else {
  533.                 span = containers.addValidAfter(copyContainer(span.getData(), flattener), t0, false);
  534.             }
  535.         }

  537.         while (span.getStart().isBefore(t1)) {
  538.             if (span.getEnd().isAfter(t1)) {
  539.                 // this span extends past t1, we must split it
  540.                 span = containers.addValidBefore(copyContainer(span.getData(), flattener), t1, false);
  541.             }
  542.             active.add(span);
  543.             span = span.next();
  544.         }

  546.         return active;

  548.     }

  550.     /** Copy a container.
  551.      * @param original time span to copy (may be null)
  552.      * @param flattener converter between triangular and flat forms
  553.      * @return fresh copy
  554.      */
  555.     private Container copyContainer(final Container original, final Flattener flattener) {
  556.         return original == null ?
  557.                new Container(flattener) :
  558.                original.resize(flattener.getDegree(), flattener.getOrder());
  559.     }

  561.     /** Get the index of a reference date, adding it if needed.
  562.      * @param known known reference dates
  563.      * @param referenceDate reference date to select
  564.      * @return index of the reference date in the {@code known} list
  565.      * @since 11.1
  566.      */
  567.     private int referenceDateIndex(final List<AbsoluteDate> known, final AbsoluteDate referenceDate) {
  568.         for (int i = 0; i < known.size(); ++i) {
  569.             if (known.get(i).equals(referenceDate)) {
  570.                 return i;
  571.             }
  572.         }
  573.         known.add(referenceDate);
  574.         return known.size() - 1;
  575.     }

  577.     /** Get the index of a pulsation, adding it if needed.
  578.      * @param known known pulsations
  579.      * @param pulsation pulsation to select
  580.      * @return index of the pulsation in the {@code known} list
  581.      * @since 11.1
  582.      */
  583.     private int pulsationIndex(final List<Double> known, final double pulsation) {
  584.         for (int i = 0; i < known.size(); ++i) {
  585.             if (Precision.equals(known.get(i), pulsation, 1)) {
  586.                 return i;
  587.             }
  588.         }
  589.         known.add(pulsation);
  590.         return known.size() - 1;
  591.     }

  593.     /** {@inheritDoc} */
  594.     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
  595.                                                      final int degree, final int order) {

  597.         // get the constant part of the field
  598.         final ConstantSphericalHarmonics constant = getBaseProvider(wantNormalized, degree, order);
  599.         if (containers == null) {
  600.             // there are no time-dependent parts in the field
  601.             return constant;
  602.         }

  604.         // create the shared parts of the model
  605.         final AbsoluteDate[] dates = new AbsoluteDate[referenceDates.size()];
  606.         for (int i = 0; i < dates.length; ++i) {
  607.             dates[i] = referenceDates.get(i);
  608.         }
  609.         final double[] puls = new double[pulsations.size()];
  610.         for (int i = 0; i < puls.length; ++i) {
  611.             puls[i] = pulsations.get(i);
  612.         }

  614.         // convert the mutable containers to piecewise parts with desired normalization
  615.         final TimeSpanMap<PiecewisePart> pieces = new TimeSpanMap<>(null);
  616.         for (TimeSpanMap.Span<Container> span = containers.getFirstSpan(); span != null; span = span.next()) {
  617.             if (span.getData() != null) {
  618.                 final Flattener spanFlattener = span.getData().flattener;
  619.                 final Flattener rescaledFlattener = new Flattener(FastMath.min(degree, spanFlattener.getDegree()),
  620.                                                                   FastMath.min(order, spanFlattener.getOrder()));
  621.                 pieces.addValidBetween(new PiecewisePart(rescaledFlattener,
  622.                                                          rescale(wantNormalized, rescaledFlattener, span.getData().flattener,
  623.                                                                  span.getData().components)),
  624.                                        span.getStart(), span.getEnd());
  625.             }
  626.         }

  628.         return new PiecewiseSphericalHarmonics(constant, dates, puls, pieces);

  630.     }

  632.     /** Parse a reference date.
  633.      * <p>
  634.      * The reference dates have either the format yyyymmdd (for 2011 format)
  635.      * or the format yyyymmdd.xxxx (for format version 2.0).
  636.      * </p>
  637.      * <p>
  638.      * The documentation for 2011 format does not specify the time scales,
  639.      * but on of the example reads "The reference time t0 is: t0 = 2005.0 y"
  640.      * when the corresponding field in the data section reads "20050101",
  641.      * so we assume the dates are consistent with astronomical conventions
  642.      * and 2005.0 is 2005-01-01T12:00:00 (i.e. noon).
  643.      * </p>
  644.      * <p>
  645.      * The 2.0 format is not described anywhere (at least I did not find any
  646.      * references), but the .xxxx fractional part seems to be hours and minutes chosen
  647.      * close to some strong earthquakes looking at the dates in Eigen 6S4 file
  648.      * with non-zero fractional part and the corresponding earthquakes hours
  649.      * (19850109.1751 vs. 1985-01-09T19:28:21, but it was not really a big quake,
  650.      * maybe there is a confusion with the 1985 Mexico earthquake at 1985-09-19T13:17:47,
  651.      * 20020815.0817 vs 2002-08-15:05:30:26, 20041226.0060 vs. 2004-12-26T00:58:53,
  652.      * 20100227.0735 vs. 2010-02-27T06:34:11, and 20110311.0515 vs. 2011-03-11T05:46:24).
  653.      * We guess the .0060 fractional part for the 2004 Sumatra-Andaman islands
  654.      * earthquake results from sloppy rounding when writing the file.
  655.      * </p>
  656.      * @param field text field containing the date
  657.      * @return parsed date
  658.      * @since 11.1
  659.      */
  660.     private AbsoluteDate parseDate(final String field) {

  662.         // check the date part (format yyyymmdd)
  663.         final DateComponents dc = new DateComponents(Integer.parseInt(field.substring(0, 4)),
  664.                                                      Integer.parseInt(field.substring(4, 6)),
  665.                                                      Integer.parseInt(field.substring(6, 8)));

  667.         // check the hour part (format .hhmm, working around checks on minutes)
  668.         final TimeComponents tc;
  669.         if (field.length() > 8) {
  670.             // we convert from hours and minutes here in order to allow
  671.             // the strange special case found in Eigen 6S4 file with date 20041226.0060
  672.             tc = new TimeComponents(Integer.parseInt(field.substring(9, 11)) * 3600 +
  673.                                     Integer.parseInt(field.substring(11, 13)) * 60);
  674.         } else {
  675.             // assume astronomical convention for hour
  676.             tc = TimeComponents.H12;
  677.         }

  679.         return toDate(dc, tc);

  681.     }

  683.     /** Temporary container for reading coefficients.
  684.      * @since 11.1
  685.      */
  686.     private static class Container {

  688.         /** Converter between (degree, order) indices and flatten array. */
  689.         private final Flattener flattener;

  691.         /** Components of the spherical harmonics. */
  692.         private final TimeDependentHarmonic[] components;

  694.         /** Build a container with given degree and order.
  695.          * @param flattener converter between (degree, order) indices and flatten array
  696.          */
  697.         Container(final Flattener flattener) {
  698.             this.flattener  = flattener;
  699.             this.components = new TimeDependentHarmonic[flattener.arraySize()];
  700.         }

  702.         /** Build a resized container.
  703.          * @param degree degree of the container
  704.          * @param order order of the container
  705.          * @return resized container
  706.          */
  707.         Container resize(final int degree, final int order) {

  709.             // create new instance
  710.             final Container resized = new Container(new Flattener(degree, order));

  712.             // copy harmonics
  713.             for (int n = 0; n <= degree; ++n) {
  714.                 for (int m = 0; m <= FastMath.min(n, order); ++m) {
  715.                     resized.components[resized.flattener.index(n, m)] = components[flattener.index(n, m)];
  716.                 }
  717.             }

  719.             return resized;

  721.         }

  723.     }

  725.     /** Errors present in the gravity field.
  726.      * @since 11.1
  727.      */
  728.     private enum Errors {

  730.         /** No errors. */
  731.         NO(0),

  733.         /** Calibrated errors. */
  734.         CALIBRATED(2),

  736.         /** Formal errors. */
  737.         FORMAL(2),

  739.         /** Both calibrated and formal. */
  740.         CALIBRATED_AND_FORMAL(4);

  742.         /** Number of error fields in data lines. */
  743.         private final int fields;

  745.         /** Simple constructor.
  746.          * @param fields umber of error fields in data lines
  747.          */
  748.         Errors(final int fields) {
  749.             this.fields = fields;
  750.         }

  752.     }

  754. }
Created with JaCoCo 0.8.12.202403310830