Work In Progress on ICGEM V2.0 format.

src/main/java/org/orekit/forces/gravity/potential/ICGEMFormatReader.java

@@ -38,8 +38,10 @@ import org.orekit.errors.OrekitMessages;
 import org.orekit.errors.OrekitParseException;
 import org.orekit.time.AbsoluteDate;
 import org.orekit.time.DateComponents;
+import org.orekit.time.TimeComponents;
 import org.orekit.time.TimeScale;
 import org.orekit.utils.Constants;
+import org.orekit.utils.TimeSpanMap;
 
 /** Reader for the ICGEM gravity field format.
  *
@@ -143,26 +145,8 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
     /** Indicator for normalized coefficients. */
     private boolean normalized;
 
-    /** Reference date. */
-    private AbsoluteDate referenceDate;
-
-    /** Secular trend of the cosine coefficients. */
-    private final List<List<Double>> cTrend;
-
-    /** Secular trend of the sine coefficients. */
-    private final List<List<Double>> sTrend;
-
-    /** Cosine part of the cosine coefficients pulsation. */
-    private final Map<Double, List<List<Double>>> cCos;
-
-    /** Sine part of the cosine coefficients pulsation. */
-    private final Map<Double, List<List<Double>>> cSin;
-
-    /** Cosine part of the sine coefficients pulsation. */
-    private final Map<Double, List<List<Double>>> sCos;
-
-    /** Sine part of the sine coefficients pulsation. */
-    private final Map<Double, List<List<Double>>> sSin;
+    /** Time map of the harmonics. */
+    private TimeSpanMap<PiecewiseSphericalHarmonics> map;
 
     /** Simple constructor.
      *
@@ -191,13 +175,6 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
                              final boolean missingCoefficientsAllowed,
                              final TimeScale timeScale) {
         super(supportedNames, missingCoefficientsAllowed, timeScale);
-        referenceDate = null;
-        cTrend = new ArrayList<>();
-        sTrend = new ArrayList<>();
-        cCos   = new HashMap<>();
-        cSin   = new HashMap<>();
-        sCos   = new HashMap<>();
-        sSin   = new HashMap<>();
     }
 
     /** {@inheritDoc} */
@@ -206,27 +183,21 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
 
         // reset the indicator before loading any data
         setReadComplete(false);
-        referenceDate = null;
-        cTrend.clear();
-        sTrend.clear();
-        cCos.clear();
-        cSin.clear();
-        sCos.clear();
-        sSin.clear();
+        map = new TimeSpanMap<>(null);
 
         // by default, the field is normalized with unknown tide system
         // (will be overridden later if non-default)
         normalized = true;
         TideSystem tideSystem = TideSystem.UNKNOWN;
 
+        double version   = 1.0;
         boolean inHeader = true;
-        double[][] c     = null;
-        double[][] s     = null;
-        boolean okCoeffs = false;
+        TimeSpanMap<PiecewiseSphericalHarmonics> map = new TimeSpanMap<>(null);
         int lineNumber   = 0;
         String line      = null;
         try (BufferedReader r = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
             for (line = r.readLine(); line != null; line = r.readLine()) {
+                boolean parseError = false;
                 ++lineNumber;
                 line = line.trim();
                 if (line.length() == 0) {
@@ -234,10 +205,16 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
                 }
                 final String[] tab = SEPARATOR.split(line);
                 if (inHeader) {
-                    boolean parseError = false;
                     if (tab.length == 2 && FORMAT.equals(tab[0])) {
                         final Matcher matcher = Pattern.compile(SUPPORTED_FORMAT).matcher(tab[1]);
-                        parseError = !matcher.matches() || Double.parseDouble(matcher.group(1)) > MAX_FORMAT;
+                        if (matcher.matches()) {
+                            version = Double.parseDouble(matcher.group(1));
+                            if (version > MAX_FORMAT) {
+                                parseError = true;
+                            }
+                        } else {
+                            parseError = true;
+                        }
                     } else if (tab.length == 2 && PRODUCT_TYPE.equals(tab[0])) {
                         parseError = !GRAVITY_FIELD.equals(tab[1]);
                     } else if (tab.length == 2 && tab[0].endsWith(GRAVITY_CONSTANT)) {
@@ -277,7 +254,7 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
                                                        lineNumber, name, line);
                     }
                 } else {
-                    if (tab.length == 7 && GFC.equals(tab[0]) || tab.length == 8 && GFCT.equals(tab[0])) {
+                    if (tab.length == 7 && GFC.equals(tab[0]) || (tab.length == 8 || tab.length == 9) && GFCT.equals(tab[0])) {
 
                         final int i = Integer.parseInt(tab[1]);
                         final int j = Integer.parseInt(tab[2]);
@@ -287,19 +264,30 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
                             parseCoefficient(tab[4], s, i, j, "S", name);
                             okCoeffs = true;
 
-                            if (tab.length == 8) {
-                                // check the reference date (format yyyymmdd)
-                                final DateComponents localRef = new DateComponents(Integer.parseInt(tab[7].substring(0, 4)),
-                                                                                   Integer.parseInt(tab[7].substring(4, 6)),
-                                                                                   Integer.parseInt(tab[7].substring(6, 8)));
-                                if (referenceDate == null) {
-                                    // first reference found, store it
-                                    referenceDate = toDate(localRef);
-                                } else if (!referenceDate.equals(toDate(localRef))) {
-                                    final AbsoluteDate localDate = toDate(localRef);
-                                    throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
-                                                              referenceDate, localDate, name,
-                                                              localDate.durationFrom(referenceDate));
+                            if (tab.length > 8) {
+                                if (version < 2.0) {
+                                    // before version 2.0, there is only one reference date
+                                    if (tab.length != 8) {
+                                        parseError = true;
+                                    } else {
+                                        final AbsoluteDate lineDate = parseDate(tab[7]);
+                                        if (referenceDate == null) {
+                                            // first reference found, store it
+                                            referenceDate = lineDate;
+                                        } else if (!referenceDate.equals(lineDate)) {
+                                            // we already know the reference date, check this lines does not define a new one
+                                            throw new OrekitException(OrekitMessages.SEVERAL_REFERENCE_DATES_IN_GRAVITY_FIELD,
+                                                                      referenceDate, lineDate, name,
+                                                                      lineDate.durationFrom(referenceDate));
+                                        }
+                                    }
+                                } else {
+                                    // starting with version 2.0, two reference dates define validity intervals
+                                    if (tab.length != 9) {
+                                        parseError = true;
+                                    } else {
+                                        // TODO
+                                    }
                                 }
                             }
 
@@ -353,9 +341,14 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
                         }
 
                     } else {
+                        parseError = true;
+                    }
+
+                    if (parseError) {
                         throw new OrekitParseException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
                                                        lineNumber, name, line);
                     }
+
                 }
 
             }
@@ -400,39 +393,98 @@ public class ICGEMFormatReader extends PotentialCoefficientsReader {
     public RawSphericalHarmonicsProvider getProvider(final boolean wantNormalized,
                                                      final int degree, final int order) {
 
-        RawSphericalHarmonicsProvider provider = getConstantProvider(wantNormalized, degree, order);
-        if (cTrend.isEmpty() && cCos.isEmpty()) {
-            // there are no time-dependent coefficients
-            return provider;
-        }
+        return new RawSphericalHarmonicsProvider() {
+
+            /** {@inheritDoc} */
+            @Override
+            public int getMaxDegree() {
+                return map.getFirstSpan().getData().getConstant().getMaxDegree();
+            }
 
-        if (!cTrend.isEmpty()) {
+            /** {@inheritDoc} */
+            @Override
+            public int getMaxOrder() {
+                return map.getFirstSpan().getData().getConstant().getMaxOrder();
+            }
 
-            // add the secular trend layer
-            final double[][] cArrayTrend = toArray(cTrend);
-            final double[][] sArrayTrend = toArray(sTrend);
-            rescale(1.0 / Constants.JULIAN_YEAR, normalized, cArrayTrend, sArrayTrend, wantNormalized, cArrayTrend, sArrayTrend);
-            provider = new SecularTrendSphericalHarmonics(provider, referenceDate, cArrayTrend, sArrayTrend);
+            /** {@inheritDoc} */
+            @Override
+            public double getMu() {
+                return map.getFirstSpan().getData().getConstant().getMu();
+            }
 
-        }
+            /** {@inheritDoc} */
+            @Override
+            public double getAe() {
+                return map.getFirstSpan().getData().getConstant().getAe();
+            }
+
+            /** {@inheritDoc} */
+            @Override
+            public AbsoluteDate getReferenceDate() {
+                return map.getFirstSpan().getData().getConstant().getReferenceDate();
+            }
+
+            /** {@inheritDoc} */
+            @Override
+            public TideSystem getTideSystem() {
+                return map.getFirstSpan().getData().getConstant().getTideSystem();
+            }
+
+            /** {@inheritDoc} */
+            @Override
+            public double getOffset(final AbsoluteDate date) {
+                return map.get(date).getConstant().getOffset(date);
+            }
 
-        for (final Map.Entry<Double, List<List<Double>>> entry : cCos.entrySet()) {
+            /** {@inheritDoc} */
+            @Override
+            public RawSphericalHarmonics onDate(final AbsoluteDate date) {
+                return map.get(date).onDate(date);
+            }
 
-            final double period = entry.getKey();
+        };
 
-            // add the pulsating layer for the current period
-            final double[][] cArrayCos = toArray(cCos.get(period));
-            final double[][] sArrayCos = toArray(sCos.get(period));
-            final double[][] cArraySin = toArray(cSin.get(period));
-            final double[][] sArraySin = toArray(sSin.get(period));
-            rescale(1.0, normalized, cArrayCos, sArrayCos, wantNormalized, cArrayCos, sArrayCos);
-            rescale(1.0, normalized, cArraySin, sArraySin, wantNormalized, cArraySin, sArraySin);
-            provider = new PulsatingSphericalHarmonics(provider, period * Constants.JULIAN_YEAR,
-                                                       cArrayCos, cArraySin, sArrayCos, sArraySin);
+    }
 
+    /** Parse a reference date.
+     * <p>
+     * The reference dates have either the format yyyymmdd (for 2011 format)
+     * or the format yyyymmdd.xxxx (for format version 2.0). The 2.0 format
+     * is not described anywhere (at least I did not find any references),
+     * but the .xxxx fractional part seems to be hours and minutes chosen
+     * close to some strong earthquakes looking at the dates in Eigen 6S4 file
+     * with non-zero fractional part and the corresponding earthquakes hours
+     * (19850109.1751 vs. 1985-01-09T19:28:21, but it was not really a big quake,
+     * maybe there is a confusion with the 1985 Mexico earthquake at 1985-09-19T13:17:47,
+     * 20020815.0817 vs 2002-08-15:05:30:26, 20041226.0060 vs. 2004-12-26T00:58:53,
+     * 20100227.0735 vs. 2010-02-27T06:34:11, and 20110311.0515 vs. 2011-03-11T05:46:24).
+     * We guess the .0060 fractional part for the 2004 Sumatra-Andaman islands
+     * earthquake results from sloppy rounding when writing the file.
+     * </p>
+     * @param field text field containing the date
+     * @return parsed date
+     * @since 11.1
+     */
+    private AbsoluteDate parseDate(final String field) {
+
+        // check the date part (format yyyymmdd)
+        final DateComponents dc = new DateComponents(Integer.parseInt(field.substring(0, 4)),
+                                                     Integer.parseInt(field.substring(4, 6)),
+                                                     Integer.parseInt(field.substring(6, 8)));
+
+        // check the hour part (format .hhmm, working around checks on minutes)
+        final TimeComponents tc;
+        if (field.length() > 8) {
+            // we convert from hours and minutes here in order to allow
+            // the strange special case found in Eigen 6S4 file with date 20041226.0060
+            tc = new TimeComponents(Integer.parseInt(field.substring(9, 11)) / 24.0 +
+                                    Integer.parseInt(field.substring(11, 13)) / 3600.0);
+        } else {
+            tc = TimeComponents.H00;
         }
 
-        return provider;
+        return toDate(dc, tc);
 
     }
 

src/main/java/org/orekit/forces/gravity/potential/PiecewiseSphericalHarmonics.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.forces.gravity.potential;
+
+import org.hipparchus.util.FastMath;
+import org.hipparchus.util.SinCos;
+import org.orekit.forces.gravity.potential.RawSphericalHarmonicsProvider.RawSphericalHarmonics;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.utils.TimeSpanMap;
+
+/** Part of a piecewise gravity fields valid for one time interval.
+ * @author Luc Maisonobe
+ * @since 11.1
+ */
+class PiecewiseSphericalHarmonics {
+
+    /** Constant part of the field. */
+    private final ConstantSphericalHarmonics constant;
+
+    /** Reference dates. */
+    private final AbsoluteDate[] references;
+
+    /** Pulsations (rad/s). */
+    private final double[] pulsations;
+
+    /** Harmonics. */
+    private final TimeSpanMap<TimeDependentHarmonic[]> harmonics;
+
+    /** Simple constructor.
+     * @param constant constant part of the field
+     * @param references references dates
+     * @param pulsations pulsations (rad/s)
+     * @param harmonics map fo harmonics components
+     */
+    PiecewiseSphericalHarmonics(final ConstantSphericalHarmonics constant,
+                                final AbsoluteDate[] references, final double[] pulsations,
+                                final TimeSpanMap<TimeDependentHarmonic[]> harmonics) {
+        this.constant   = constant;
+        this.references = references.clone();
+        this.pulsations = pulsations.clone();
+        this.harmonics  = harmonics;
+    }
+
+    /** Get the constant part of the field.
+     * @return constant part of the field
+     */
+    public ConstantSphericalHarmonics getConstant() {
+        return constant;
+    }
+
+    /** Convert (degree, order) indices to one-dimensional flatten index.
+     * <p>
+     * As the outer loop in {@link org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel}
+     * if on decreasing order and the inner loop is in increasing degree (starting
+     * from the diagonal), the flatten index increases as these loops are performed.
+     * </p>
+     * @param n degree index
+     * @param m order index
+     * @return one-dimensional flatten index
+     */
+    private int index(final int n, final int m) {
+        // TODO
+    }
+
+    /** Get the raw spherical harmonic coefficients on a specific date.
+     * @param date to evaluate the spherical harmonics
+     * @return the raw spherical harmonics on {@code date}.
+     */
+    public RawSphericalHarmonics onDate(final AbsoluteDate date) {
+
+        // raw (constant) harmonics
+        final RawSphericalHarmonics raw = constant.onDate(date);
+
+        // select part of the piecewise model that is active at specified date
+        final TimeDependentHarmonic[] active = harmonics.get(date);
+
+        // pre-compute canonical functions
+        final double[]   offsets = new double[references.length];
+        final SinCos[][] sinCos  = new SinCos[references.length][pulsations.length];
+        for (int i = 0; i < references.length; ++i) {
+            final double offset = date.durationFrom(references[i]);
+            offsets[i] = offset;
+            for (int j = 0; j < pulsations.length; ++j) {
+                sinCos[i][j] = FastMath.sinCos(offset * pulsations[j]);
+            }
+        }
+
+        return new RawSphericalHarmonics() {
+
+            @Override
+            public AbsoluteDate getDate() {
+                return date;
+            }
+
+            /** {@inheritDoc} */
+            public double getRawCnm(final int n, final int m) {
+                return active[index(n, m)].getRawCnm(raw.getRawCnm(n, m), offsets, sinCos);
+            }
+
+            /** {@inheritDoc} */
+            public double getRawSnm(final int n, final int m) {
+                return active[index(n, m)].getRawSnm(raw.getRawSnm(n, m), offsets, sinCos);
+            }
+
+        };
+
+    }
+
+}

src/main/java/org/orekit/forces/gravity/potential/PotentialCoefficientsReader.java

@@ -122,7 +122,7 @@ public abstract class PotentialCoefficientsReader implements DataLoader {
         this.rawS                       = null;
         this.normalized                 = false;
         this.tideSystem                 = TideSystem.UNKNOWN;
-        this.timeScale = timeScale;
+        this.timeScale                  = timeScale;
     }
 
     /** Get the regular expression for supported files names.
@@ -569,7 +569,19 @@ public abstract class PotentialCoefficientsReader implements DataLoader {
      * @return date.
      */
     protected AbsoluteDate toDate(final DateComponents components) {
-        return new AbsoluteDate(components, TimeComponents.H12, timeScale);
+        return toDate(components, TimeComponents.H12);
+    }
+
+    /**
+     * Create a date from components.
+     *
+     * @param dc dates components.
+     * @param tc time components
+     * @return date.
+     * @since 11.1
+     */
+    protected AbsoluteDate toDate(final DateComponents dc, final TimeComponents tc) {
+        return new AbsoluteDate(dc, tc, timeScale);
     }
 
 }

src/main/java/org/orekit/forces/gravity/potential/TimeDependentHarmonic.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.forces.gravity.potential;
+
+import org.hipparchus.util.SinCos;
+
+/** Time-dependent part of a single component of spherical harmonics.
+ * @author Luc Maisonobe
+ * @since 11.1
+ */
+class TimeDependentHarmonic {
+
+    /** Index of the trend reference in the gravity field. */
+    private final int trendReferenceIndex;
+
+    /** Secular trend of the cosine coefficient. */
+    private final double cTrend;
+
+    /** Secular trend of the sine coefficient. */
+    private final double sTrend;
+
+    /** Indices of the reference dates in the gravity field. */
+    private int[] referenceIndices;
+
+    /** Indices of the harmonic pulsations in the gravity field. */
+    private int[] pulsationIndices;
+
+    /** Cosine component of the cosine coefficient. */
+    private double[] cosC;
+
+    /** Sine component of the cosine coefficient. */
+    private double[] sinC;
+
+    /** Cosine component of the sine coefficient. */
+    private double[] cosS;
+
+    /** Sine component of the sine coefficient. */
+    private double[] sinS;
+
+    /** Build a part with only trend.
+     * @param trendReferenceIndex index of the trend reference in the gravity field
+     * @param cTrend secular trend of the cosine coefficient (s⁻¹)
+     * @param sTrend secular trend of the sine coefficient (s⁻¹)
+     */
+    TimeDependentHarmonic(final int trendReferenceIndex, final double cTrend, final double sTrend) {
+
+        // linear part
+        this.trendReferenceIndex = trendReferenceIndex;
+        this.cTrend              = cTrend;
+        this.sTrend              = sTrend;
+
+        // empty harmonic part
+        this.referenceIndices    = new int[0];
+        this.pulsationIndices    = new int[0];
+        this.cosC                = new double[0];
+        this.sinC                = new double[0];
+        this.cosS                = new double[0];
+        this.sinS                = new double[0];
+
+    }
+
+    /** Add an harmonic component.
+     * @param referenceIndex index of the reference date in the gravity field
+     * @param pulsationIndex index of the harmonic pulsation in the gravity field
+     * @param cosineC cosine component of the cosine coefficient
+     * @param sineC sine component of the cosine coefficient
+     * @param cosineS cosine component of the sine coefficient
+     * @param sineS sine component of the sine coefficient
+     */
+    public void addHarmonic(final int referenceIndex, final int pulsationIndex,
+                            final double cosineC, final double sineC,
+                            final double cosineS, final double sineS) {
+        this.referenceIndices = addInt(referenceIndex, this.referenceIndices);
+        this.pulsationIndices = addInt(pulsationIndex, this.pulsationIndices);
+        this.cosC             = addDouble(cosineC,     this.cosC);
+        this.sinC             = addDouble(sineC,       this.sinC);
+        this.cosS             = addDouble(cosineS,     this.cosS);
+        this.sinS             = addDouble(sineS,       this.sinS);
+    }
+
+    /** Add an integer to an array.
+     * <p>
+     * Expanding the array one element at a time may seem a waste of time,
+     * but we expect the array to be 0, 1 or 2 elements long only, and this
+     * if performed only when reading gravity field, so its is worth doing
+     * it this way.
+     * </p>
+     * @param n integer to add
+     * @param array array where to add the integer
+     * @return new array
+     */
+    private static int[] addInt(final int n, final int[] array) {
+        final int[] newArray = new int[array.length];
+        System.arraycopy(array, 0, newArray, 0, array.length);
+        newArray[array.length] = n;
+        return newArray;
+    }
+
+    /** Add a double number to an array.
+     * <p>
+     * Expanding the array one element at a time may seem a waste of time,
+     * but we expect the array to be 0, 1 or 2 elements long only, and this
+     * if performed only when reading gravity field, so its is worth doing
+     * it this way.
+     * </p>
+     * @param d double number to add
+     * @param array array where to add the double number
+     * @return new array
+     */
+    private static double[] addDouble(final double d, final double[] array) {
+        final double[] newArray = new double[array.length];
+        System.arraycopy(array, 0, newArray, 0, array.length);
+        newArray[array.length] = d;
+        return newArray;
+    }
+
+    /** Get a spherical harmonic cosine coefficient.
+     * @param constantCnm constant part of the cosine coefficient
+     * @param offsets offsets to reference dates in the gravity field
+     * @param pulsations angular pulsations in the gravity field
+     * @return raw coefficient Cnm
+     */
+    public double getRawCnm(final double constantCnm, final double[] offsets, final SinCos[][] pulsations) {
+
+        // linear part
+        double cnm = constantCnm + offsets[trendReferenceIndex] * cTrend;
+
+        for (int i = 0; i < pulsationIndices.length; ++i) {
+            // add pulsation
+            final SinCos pulsation = pulsations[referenceIndices[i]][pulsationIndices[i]];
+            cnm += cosC[i] * pulsation.cos() + sinC[i] * pulsation.sin();
+        }
+
+        return cnm;
+
+    }
+
+    /** Get a spherical harmonic sine coefficient.
+     * @param constantSnm constant part of the sine coefficient
+     * @param offsets offsets to reference dates in the gravity field
+     * @param pulsations angular pulsations in the gravity field
+     * @return raw coefficient Snm
+     */
+    public double getRawSnm(final double constantSnm, final double[] offsets, final SinCos[][] pulsations) {
+
+        // linear part
+        double snm = constantSnm + offsets[trendReferenceIndex] * sTrend;
+
+        for (int i = 0; i < pulsationIndices.length; ++i) {
+            // add pulsation
+            final SinCos pulsation = pulsations[referenceIndices[i]][pulsationIndices[i]];
+            snm += cosS[i] * pulsation.cos() + sinS[i] * pulsation.sin();
+        }
+
+        return snm;
+
+    }
+
+}