GlobalIonosphereMapModel.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.models.earth.ionosphere;

  18. import java.io.BufferedInputStream;
  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.util.ArrayList;
  25. import java.util.Collections;
  26. import java.util.List;

  27. import org.hipparchus.CalculusFieldElement;
  28. import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
  29. import org.hipparchus.exception.DummyLocalizable;
  30. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  31. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  32. import org.hipparchus.util.FastMath;
  33. import org.orekit.annotation.DefaultDataContext;
  34. import org.orekit.bodies.BodyShape;
  35. import org.orekit.bodies.GeodeticPoint;
  36. import org.orekit.data.DataContext;
  37. import org.orekit.data.DataLoader;
  38. import org.orekit.data.DataProvidersManager;
  39. import org.orekit.data.DataSource;
  40. import org.orekit.errors.OrekitException;
  41. import org.orekit.errors.OrekitMessages;
  42. import org.orekit.frames.FieldStaticTransform;
  43. import org.orekit.frames.Frame;
  44. import org.orekit.frames.StaticTransform;
  45. import org.orekit.frames.TopocentricFrame;
  46. import org.orekit.propagation.FieldSpacecraftState;
  47. import org.orekit.propagation.SpacecraftState;
  48. import org.orekit.time.AbsoluteDate;
  49. import org.orekit.time.FieldAbsoluteDate;
  50. import org.orekit.time.TimeScale;
  51. import org.orekit.utils.ParameterDriver;
  52. import org.orekit.utils.TimeSpanMap;

  53. /**
  54.  * Global Ionosphere Map (GIM) model.
  55.  * The ionospheric delay is computed according to the formulas:
  56.  * <pre>
  57.  *           40.3
  58.  *    δ =  --------  *  STEC      with, STEC = VTEC * F(elevation)
  59.  *            f²
  60.  * </pre>
  61.  * With:
  62.  * <ul>
  63.  * <li>f: The frequency of the signal in Hz.</li>
  64.  * <li>STEC: The Slant Total Electron Content in TECUnits.</li>
  65.  * <li>VTEC: The Vertical Total Electron Content in TECUnits.</li>
  66.  * <li>F(elevation): A mapping function which depends on satellite elevation.</li>
  67.  * </ul>
  68.  * The VTEC is read from a IONEX file. A stream contains, for a given day, the values of the TEC for each hour of the day.
  69.  * Values are given on a global 2.5° x 5.0° (latitude x longitude) grid.
  70.  * <p>
  71.  * A bilinear interpolation is performed the case of the user initialize the latitude and the
  72.  * longitude with values that are not contained in the stream.
  73.  * </p><p>
  74.  * A temporal interpolation is also performed to compute the VTEC at the desired date.
  75.  * </p><p>
  76.  * IONEX files are obtained from
  77.  * <a href="ftp://cddis.nasa.gov/gnss/products/ionex/"> The Crustal Dynamics Data Information System</a>.
  78.  * </p><p>
  79.  * The files have to be extracted to UTF-8 text files before being read by this loader.
  80.  * </p><p>
  81.  * Example of file:
  82.  * </p>
  83.  * <pre>
  84.  *      1.0            IONOSPHERE MAPS     GPS                 IONEX VERSION / TYPE
  85.  * BIMINX V5.3         AIUB                16-JAN-19 07:26     PGM / RUN BY / DATE
  86.  * BROADCAST IONOSPHERE MODEL FOR DAY 015, 2019                COMMENT
  87.  *   2019     1    15     0     0     0                        EPOCH OF FIRST MAP
  88.  *   2019     1    16     0     0     0                        EPOCH OF LAST MAP
  89.  *   3600                                                      INTERVAL
  90.  *     25                                                      # OF MAPS IN FILE
  91.  *   NONE                                                      MAPPING FUNCTION
  92.  *      0.0                                                    ELEVATION CUTOFF
  93.  *                                                             OBSERVABLES USED
  94.  *   6371.0                                                    BASE RADIUS
  95.  *      2                                                      MAP DIMENSION
  96.  *    350.0 350.0   0.0                                        HGT1 / HGT2 / DHGT
  97.  *     87.5 -87.5  -2.5                                        LAT1 / LAT2 / DLAT
  98.  *   -180.0 180.0   5.0                                        LON1 / LON2 / DLON
  99.  *     -1                                                      EXPONENT
  100.  * TEC/RMS values in 0.1 TECU; 9999, if no value available     COMMENT
  101.  *                                                             END OF HEADER
  102.  *      1                                                      START OF TEC MAP
  103.  *   2019     1    15     0     0     0                        EPOCH OF CURRENT MAP
  104.  *     87.5-180.0 180.0   5.0 350.0                            LAT/LON1/LON2/DLON/H
  105.  *    92   92   92   92   92   92   92   92   92   92   92   92   92   92   92   92
  106.  *    92   92   92   92   92   92   92   92   92   92   92   92   92   92   92   92
  107.  *    92   92   92   92   92   92   92   92   92   92   92   92   92   92   92   92
  108.  *    92   92   92   92   92   92   92   92   92   92   92   92   92   92   92   92
  109.  *    92   92   92   92   92   92   92   92   92
  110.  *    ...
  111.  * </pre>
  112.  *
  113.  * @see "Schaer, S., W. Gurtner, and J. Feltens, 1998, IONEX: The IONosphere Map EXchange
  114.  *       Format Version 1, February 25, 1998, Proceedings of the IGS AC Workshop
  115.  *       Darmstadt, Germany, February 9–11, 1998"
  116.  *
  117.  * @author Bryan Cazabonne
  118.  *
  119.  */
  120. public class GlobalIonosphereMapModel implements IonosphericModel, IonosphericDelayModel {

  121.     /** Map of interpolable TEC. */
  122.     private final TimeSpanMap<TECMapPair> tecMap;

  123.     /** UTC time scale. */
  124.     private final TimeScale utc;

  125.     /** Loaded IONEX files.
  126.      * @since 12.0
  127.      */
  128.     private String names;

  129.     /**
  130.      * Constructor with supported names given by user. This constructor uses the {@link
  131.      * DataContext#getDefault() default data context}.
  132.      *
  133.      * @param supportedNames regular expression that matches the names of the IONEX files
  134.      *                       to be loaded. See {@link DataProvidersManager#feed(String,
  135.      *                       DataLoader)}.
  136.      * @see #GlobalIonosphereMapModel(String, DataProvidersManager, TimeScale)
  137.      */
  138.     @DefaultDataContext
  139.     public GlobalIonosphereMapModel(final String supportedNames) {
  140.         this(supportedNames,
  141.                 DataContext.getDefault().getDataProvidersManager(),
  142.                 DataContext.getDefault().getTimeScales().getUTC());
  143.     }

  144.     /**
  145.      * Constructor that uses user defined supported names and data context.
  146.      *
  147.      * @param supportedNames       regular expression that matches the names of the IONEX
  148.      *                             files to be loaded. See {@link DataProvidersManager#feed(String,
  149.      *                             DataLoader)}.
  150.      * @param dataProvidersManager provides access to auxiliary data files.
  151.      * @param utc                  UTC time scale.
  152.      * @since 10.1
  153.      */
  154.     public GlobalIonosphereMapModel(final String supportedNames,
  155.                                     final DataProvidersManager dataProvidersManager,
  156.                                     final TimeScale utc) {
  157.         this.utc    = utc;
  158.         this.tecMap = new TimeSpanMap<>(null);
  159.         this.names  = "";

  160.         // Read files
  161.         dataProvidersManager.feed(supportedNames, new Parser());

  162.     }

  163.     /**
  164.      * Constructor that uses user defined data sources.
  165.      *
  166.      * @param utc   UTC time scale.
  167.      * @param ionex sources for the IONEX files
  168.      * @since 12.0
  169.      */
  170.     public GlobalIonosphereMapModel(final TimeScale utc,
  171.                                     final DataSource... ionex) {
  172.         try {
  173.             this.utc    = utc;
  174.             this.tecMap = new TimeSpanMap<>(null);
  175.             this.names  = "";
  176.             final Parser parser = new Parser();
  177.             for (final DataSource source : ionex) {
  178.                 try (InputStream is  = source.getOpener().openStreamOnce();
  179.                     BufferedInputStream bis = new BufferedInputStream(is)) {
  180.                     parser.loadData(bis, source.getName());
  181.                 }
  182.             }
  183.         } catch (IOException ioe) {
  184.             throw new OrekitException(ioe, new DummyLocalizable(ioe.getMessage()));
  185.         }
  186.     }

  187.     /**
  188.      * Calculates the ionospheric path delay for the signal path from a ground
  189.      * station to a satellite traversing ionosphere single layer at some pierce point.
  190.      * <p>
  191.      * The path delay can be computed for any elevation angle.
  192.      * </p>
  193.      * @param date current date
  194.      * @param piercePoint ionospheric pierce point
  195.      * @param elevation elevation of the satellite from receiver point in radians
  196.      * @param frequency frequency of the signal in Hz
  197.      * @return the path delay due to the ionosphere in m
  198.      */
  199.     private double pathDelayAtIPP(final AbsoluteDate date, final GeodeticPoint piercePoint,
  200.                                   final double elevation, final double frequency) {
  201.         // TEC in TECUnits
  202.         final TECMapPair pair = getPairAtDate(date);
  203.         final double tec = pair.getTEC(date, piercePoint);
  204.         // Square of the frequency
  205.         final double freq2 = frequency * frequency;
  206.         // "Slant" Total Electron Content
  207.         final double stec;
  208.         // Check if a mapping factor is needed
  209.         if (pair.mapping) {
  210.             stec = tec;
  211.         } else {
  212.             // Mapping factor
  213.             final double fz = mappingFunction(elevation, pair.r0, pair.h);
  214.             stec = tec * fz;
  215.         }
  216.         // Delay computation
  217.         final double alpha  = 40.3e16 / freq2;
  218.         return alpha * stec;
  219.     }

  220.     @Override
  221.     public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
  222.                             final double frequency, final double[] parameters) {
  223.         return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
  224.     }

  225.     @Override
  226.     public double pathDelay(final SpacecraftState state,
  227.                             final TopocentricFrame baseFrame, final AbsoluteDate receptionDate,
  228.                             final double frequency, final double[] parameters) {

  229.         // we use transform from body frame to inert frame and invert it
  230.         // because base frame could be peered with inertial frame (hence improving performances with caching)
  231.         // but the reverse is almost never used
  232.         final Frame           bodyFrame  = baseFrame.getParentShape().getBodyFrame();
  233.         final StaticTransform body2Inert = bodyFrame.getStaticTransformTo(state.getFrame(), receptionDate);
  234.         final Vector3D        satPoint   = body2Inert.getInverse().transformPosition(state.getPosition());

  235.         // Elevation in radians
  236.         final double   elevation = bodyFrame.
  237.                                    getStaticTransformTo(baseFrame, receptionDate).
  238.                                    transformPosition(satPoint).
  239.                                    getDelta();

  240.         // Only consider measures above the horizon
  241.         if (elevation > 0.0) {
  242.             // Normalized Line Of Sight in body frame
  243.             final Vector3D los = satPoint.subtract(baseFrame.getCartesianPoint()).normalize();
  244.             // Ionosphere Pierce Point
  245.             final GeodeticPoint ipp = piercePoint(state.getDate(), baseFrame.getCartesianPoint(), los, baseFrame.getParentShape());
  246.             if (ipp != null) {
  247.                 // Delay
  248.                 return pathDelayAtIPP(state.getDate(), ipp, elevation, frequency);
  249.             }
  250.         }

  251.         return 0.0;

  252.     }

  253.     /**
  254.      * Calculates the ionospheric path delay for the signal path from a ground
  255.      * station to a satellite traversing ionosphere single layer at some pierce point.
  256.      * <p>
  257.      * The path delay can be computed for any elevation angle.
  258.      * </p>
  259.      * @param <T> type of the elements
  260.      * @param date current date
  261.      * @param piercePoint ionospheric pierce point
  262.      * @param elevation elevation of the satellite from receiver point in radians
  263.      * @param frequency frequency of the signal in Hz
  264.      * @return the path delay due to the ionosphere in m
  265.      */
  266.     private <T extends CalculusFieldElement<T>> T pathDelayAtIPP(final FieldAbsoluteDate<T> date,
  267.                                                                  final GeodeticPoint piercePoint,
  268.                                                                  final T elevation, final double frequency) {
  269.         // TEC in TECUnits
  270.         final TECMapPair pair = getPairAtDate(date.toAbsoluteDate());
  271.         final T tec = pair.getTEC(date, piercePoint);
  272.         // Square of the frequency
  273.         final double freq2 = frequency * frequency;
  274.         // "Slant" Total Electron Content
  275.         final T stec;
  276.         // Check if a mapping factor is needed
  277.         if (pair.mapping) {
  278.             stec = tec;
  279.         } else {
  280.             // Mapping factor
  281.             final T fz = mappingFunction(elevation, pair.r0, pair.h);
  282.             stec = tec.multiply(fz);
  283.         }
  284.         // Delay computation
  285.         final double alpha  = 40.3e16 / freq2;
  286.         return stec.multiply(alpha);
  287.     }

  288.     @Override
  289.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
  290.                                                            final double frequency, final T[] parameters) {
  291.         return pathDelay(state, baseFrame, state.getDate(), frequency, parameters);
  292.     }

  293.     @Override
  294.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state,
  295.                                                            final TopocentricFrame baseFrame, final FieldAbsoluteDate<T> receptionDate,
  296.                                                            final double frequency, final T[] parameters) {

  297.         // we use transform from body frame to inert frame and invert it
  298.         // because base frame could be peered with inertial frame (hence improving performances with caching)
  299.         // but the reverse is almost never used
  300.         final Frame                   bodyFrame  = baseFrame.getParentShape().getBodyFrame();
  301.         final FieldStaticTransform<T> body2Inert = bodyFrame.getStaticTransformTo(state.getFrame(), receptionDate);
  302.         final FieldVector3D<T>        satPoint   = body2Inert.getInverse().transformPosition(state.getPosition());

  303.         // Elevation in radians
  304.         final T                elevation = bodyFrame.
  305.                                            getStaticTransformTo(baseFrame, receptionDate).
  306.                                            transformPosition(satPoint).
  307.                                            getDelta();

  308.         // Only consider measures above the horizon
  309.         if (elevation.getReal() > 0.0) {
  310.             // Normalized Line Of Sight in body frame
  311.             final Vector3D los = satPoint.toVector3D().subtract(baseFrame.getCartesianPoint()).normalize();
  312.             // Ionosphere Pierce Point
  313.             final GeodeticPoint ipp = piercePoint(state.getDate().toAbsoluteDate(), baseFrame.getCartesianPoint(), los, baseFrame.getParentShape());
  314.             if (ipp != null) {
  315.                 // Delay
  316.                 return pathDelayAtIPP(state.getDate(), ipp, elevation, frequency);
  317.             }
  318.         }

  319.         return elevation.getField().getZero();

  320.     }

  321.     /** Get the pair valid at date.
  322.      * @param date computation date
  323.      * @return pair valid at date
  324.      * @since 12.0
  325.      */
  326.     private TECMapPair getPairAtDate(final AbsoluteDate date) {
  327.         final TECMapPair pair = tecMap.get(date);
  328.         if (pair == null) {
  329.             final TimeSpanMap.Transition<TECMapPair> lastTransition = tecMap.getLastTransition();
  330.             if (lastTransition != null && lastTransition.getDate().equals(date)) {
  331.                 // we consider the transition date is in the validity range of the last span
  332.                 return lastTransition.getBefore();
  333.             }
  334.             throw new OrekitException(OrekitMessages.NO_TEC_DATA_IN_FILES_FOR_DATE,
  335.                                       names, date);
  336.         }
  337.         return pair;
  338.     }

  339.     @Override
  340.     public List<ParameterDriver> getParametersDrivers() {
  341.         return Collections.emptyList();
  342.     }

  343.     /**
  344.      * Computes the ionospheric mapping function.
  345.      * @param elevation the elevation of the satellite in radians
  346.      * @param r0 mean Earth radius
  347.      * @param h height of the ionospheric layer
  348.      * @return the mapping function
  349.      */
  350.     private double mappingFunction(final double elevation,
  351.                                    final double r0, final double h) {
  352.         // Calculate the zenith angle from the elevation
  353.         final double z = FastMath.abs(0.5 * FastMath.PI - elevation);
  354.         // Distance ratio
  355.         final double ratio = r0 / (r0 + h);
  356.         // Mapping function
  357.         final double coef = FastMath.sin(z) * ratio;
  358.         final double fz = 1.0 / FastMath.sqrt(1.0 - coef * coef);
  359.         return fz;
  360.     }

  361.     /**
  362.      * Computes the ionospheric mapping function.
  363.      * @param <T> type of the elements
  364.      * @param elevation the elevation of the satellite in radians
  365.      * @param r0 mean Earth radius
  366.      * @param h height of the ionospheric layer
  367.      * @return the mapping function
  368.      */
  369.     private <T extends CalculusFieldElement<T>> T mappingFunction(final T elevation,
  370.                                                                   final double r0, final double h) {
  371.         // Calculate the zenith angle from the elevation
  372.         final T z = FastMath.abs(elevation.getPi().multiply(0.5).subtract(elevation));
  373.         // Distance ratio
  374.         final double ratio = r0 / (r0 + h);
  375.         // Mapping function
  376.         final T coef = FastMath.sin(z).multiply(ratio);
  377.         final T fz = FastMath.sqrt(coef.multiply(coef).negate().add(1.0)).reciprocal();
  378.         return fz;
  379.     }

  380.     /** Compute Ionospheric Pierce Point.
  381.      * <p>
  382.      * The pierce point is computed assuming a spherical ionospheric shell above mean Earth radius.
  383.      * </p>
  384.      * @param date computation date
  385.      * @param recPoint point at receiver station in body frame
  386.      * @param los normalized line of sight in body frame
  387.      * @param bodyShape shape of the body
  388.      * @return pierce point, or null if recPoint is above ionosphere single layer
  389.      * @since 12.0
  390.      */
  391.     private GeodeticPoint piercePoint(final AbsoluteDate date,
  392.                                       final Vector3D recPoint, final Vector3D los,
  393.                                       final BodyShape bodyShape) {

  394.         final TECMapPair pair = getPairAtDate(date);
  395.         final double     r    = pair.r0 + pair.h;
  396.         final double     r2   = r * r;
  397.         final double     p2   = recPoint.getNormSq();
  398.         if (p2 >= r2) {
  399.             // we are above ionosphere single layer
  400.             return null;
  401.         }

  402.         // compute positive k such that recPoint + k los is on the spherical shell at radius r
  403.         final double dot = Vector3D.dotProduct(recPoint, los);
  404.         final double k   = FastMath.sqrt(dot * dot + r2 - p2) - dot;

  405.         // Ionosphere Pierce Point in body frame
  406.         final Vector3D ipp = new Vector3D(1, recPoint, k, los);

  407.         return bodyShape.transform(ipp, bodyShape.getBodyFrame(), null);

  408.     }

  409.     /** Parser for IONEX files. */
  410.     private class Parser implements DataLoader {

  411.         /** String for the end of a TEC map. */
  412.         private static final String END = "END OF TEC MAP";

  413.         /** String for the epoch of a TEC map. */
  414.         private static final String EPOCH = "EPOCH OF CURRENT MAP";

  415.         /** Index of label in data lines. */
  416.         private static final int LABEL_START = 60;

  417.         /** Kilometers to meters conversion factor. */
  418.         private static final double KM_TO_M = 1000.0;

  419.         /** Header of the IONEX file. */
  420.         private IONEXHeader header;

  421.         /** List of TEC Maps. */
  422.         private List<TECMap> maps;

  423.         @Override
  424.         public boolean stillAcceptsData() {
  425.             return true;
  426.         }

  427.         @Override
  428.         public void loadData(final InputStream input, final String name)
  429.             throws IOException {

  430.             maps = new ArrayList<>();

  431.             // Open stream and parse data
  432.             int   lineNumber = 0;
  433.             String line      = null;
  434.             try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
  435.                  BufferedReader    br  = new BufferedReader(isr)) {

  436.                 // Placeholders for parsed data
  437.                 int               nbOfMaps    = 1;
  438.                 int               exponent    = -1;
  439.                 double            baseRadius  = Double.NaN;
  440.                 double            hIon        = Double.NaN;
  441.                 boolean           mappingF    = false;
  442.                 boolean           inTEC       = false;
  443.                 double[]          latitudes   = null;
  444.                 double[]          longitudes  = null;
  445.                 AbsoluteDate      firstEpoch  = null;
  446.                 AbsoluteDate      lastEpoch   = null;
  447.                 AbsoluteDate      epoch       = firstEpoch;
  448.                 ArrayList<Double> values      = new ArrayList<>();

  449.                 for (line = br.readLine(); line != null; line = br.readLine()) {
  450.                     ++lineNumber;
  451.                     if (line.length() > LABEL_START) {
  452.                         switch (line.substring(LABEL_START).trim()) {
  453.                             case "EPOCH OF FIRST MAP" :
  454.                                 firstEpoch = parseDate(line);
  455.                                 break;
  456.                             case "EPOCH OF LAST MAP" :
  457.                                 lastEpoch = parseDate(line);
  458.                                 break;
  459.                             case "INTERVAL" :
  460.                                 // ignored;
  461.                                 break;
  462.                             case "# OF MAPS IN FILE" :
  463.                                 nbOfMaps = parseInt(line, 2, 4);
  464.                                 break;
  465.                             case "BASE RADIUS" :
  466.                                 // Value is in kilometers
  467.                                 baseRadius = parseDouble(line, 2, 6) * KM_TO_M;
  468.                                 break;
  469.                             case "MAPPING FUNCTION" :
  470.                                 mappingF = !parseString(line, 2, 4).equals("NONE");
  471.                                 break;
  472.                             case "EXPONENT" :
  473.                                 exponent = parseInt(line, 4, 2);
  474.                                 break;
  475.                             case "HGT1 / HGT2 / DHGT" :
  476.                                 if (parseDouble(line, 17, 3) == 0.0) {
  477.                                     // Value is in kilometers
  478.                                     hIon = parseDouble(line, 3, 5) * KM_TO_M;
  479.                                 }
  480.                                 break;
  481.                             case "LAT1 / LAT2 / DLAT" :
  482.                                 latitudes = parseCoordinate(line);
  483.                                 break;
  484.                             case "LON1 / LON2 / DLON" :
  485.                                 longitudes = parseCoordinate(line);
  486.                                 break;
  487.                             case "END OF HEADER" :
  488.                                 // Check that latitude and longitude boundaries were found
  489.                                 if (latitudes == null || longitudes == null) {
  490.                                     throw new OrekitException(OrekitMessages.NO_LATITUDE_LONGITUDE_BOUNDARIES_IN_IONEX_HEADER, name);
  491.                                 }
  492.                                 // Check that first and last epochs were found
  493.                                 if (firstEpoch == null || lastEpoch == null) {
  494.                                     throw new OrekitException(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, name);
  495.                                 }
  496.                                 // At the end of the header, we build the IONEXHeader object
  497.                                 header = new IONEXHeader(nbOfMaps, baseRadius, hIon, mappingF);
  498.                                 break;
  499.                             case "START OF TEC MAP" :
  500.                                 inTEC = true;
  501.                                 break;
  502.                             case END :
  503.                                 final BilinearInterpolatingFunction tec = interpolateTEC(values, exponent, latitudes, longitudes);
  504.                                 final TECMap map = new TECMap(epoch, tec);
  505.                                 maps.add(map);
  506.                                 // Reset parameters
  507.                                 inTEC  = false;
  508.                                 values = new ArrayList<>();
  509.                                 epoch  = null;
  510.                                 break;
  511.                             default :
  512.                                 if (inTEC) {
  513.                                     // Date
  514.                                     if (line.endsWith(EPOCH)) {
  515.                                         epoch = parseDate(line);
  516.                                     }
  517.                                     // Fill TEC values list
  518.                                     if (!line.endsWith("LAT/LON1/LON2/DLON/H") &&
  519.                                         !line.endsWith(END) &&
  520.                                         !line.endsWith(EPOCH)) {
  521.                                         for (int fieldStart = 0; fieldStart < line.length(); fieldStart += 5) {
  522.                                             values.add((double) Integer.parseInt(line.substring(fieldStart, fieldStart + 5).trim()));
  523.                                         }
  524.                                     }
  525.                                 }
  526.                                 break;
  527.                         }
  528.                     } else {
  529.                         if (inTEC) {
  530.                             // Here, we are parsing the last line of TEC data for a given latitude
  531.                             // The size of this line is lower than 60.
  532.                             for (int fieldStart = 0; fieldStart < line.length(); fieldStart += 5) {
  533.                                 values.add((double) Integer.parseInt(line.substring(fieldStart, fieldStart + 5).trim()));
  534.                             }
  535.                         }
  536.                     }

  537.                 }

  538.                 // Close the stream after reading
  539.                 input.close();

  540.             } catch (NumberFormatException nfe) {
  541.                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  542.                                           lineNumber, name, line);
  543.             }

  544.             // TEC map
  545.             if (maps.size() != header.getTECMapsNumer()) {
  546.                 throw new OrekitException(OrekitMessages.INCONSISTENT_NUMBER_OF_TEC_MAPS_IN_FILE,
  547.                                           maps.size(), header.getTECMapsNumer());
  548.             }
  549.             TECMap previous = null;
  550.             for (TECMap current : maps) {
  551.                 if (previous != null) {
  552.                     tecMap.addValidBetween(new TECMapPair(previous, current,
  553.                                                           header.getEarthRadius(), header.getHIon(), header.isMappingFunction()),
  554.                                            previous.date, current.date);
  555.                 }
  556.                 previous = current;
  557.             }

  558.             names = names.isEmpty() ? name : (names + ", " + name);

  559.         }

  560.         /** Extract a string from a line.
  561.          * @param line to parse
  562.          * @param start start index of the string
  563.          * @param length length of the string
  564.          * @return parsed string
  565.          */
  566.         private String parseString(final String line, final int start, final int length) {
  567.             return line.substring(start, FastMath.min(line.length(), start + length)).trim();
  568.         }

  569.         /** Extract an integer from a line.
  570.          * @param line to parse
  571.          * @param start start index of the integer
  572.          * @param length length of the integer
  573.          * @return parsed integer
  574.          */
  575.         private int parseInt(final String line, final int start, final int length) {
  576.             return Integer.parseInt(parseString(line, start, length));
  577.         }

  578.         /** Extract a double from a line.
  579.          * @param line to parse
  580.          * @param start start index of the real
  581.          * @param length length of the real
  582.          * @return parsed real
  583.          */
  584.         private double parseDouble(final String line, final int start, final int length) {
  585.             return Double.parseDouble(parseString(line, start, length));
  586.         }

  587.         /** Extract a date from a parsed line.
  588.          * @param line to parse
  589.          * @return an absolute date
  590.          */
  591.         private AbsoluteDate parseDate(final String line) {
  592.             return new AbsoluteDate(parseInt(line, 0, 6),
  593.                                     parseInt(line, 6, 6),
  594.                                     parseInt(line, 12, 6),
  595.                                     parseInt(line, 18, 6),
  596.                                     parseInt(line, 24, 6),
  597.                                     parseDouble(line, 30, 13),
  598.                                     utc);
  599.         }

  600.         /** Build the coordinate array from a parsed line.
  601.          * @param line to parse
  602.          * @return an array of coordinates in radians
  603.          */
  604.         private double[] parseCoordinate(final String line) {
  605.             final double a = parseDouble(line, 2, 6);
  606.             final double b = parseDouble(line, 8, 6);
  607.             final double c = parseDouble(line, 14, 6);
  608.             final double[] coordinate = new double[((int) FastMath.abs((a - b) / c)) + 1];
  609.             int i = 0;
  610.             for (double cor = FastMath.min(a, b); cor <= FastMath.max(a, b); cor += FastMath.abs(c)) {
  611.                 coordinate[i] = FastMath.toRadians(cor);
  612.                 i++;
  613.             }
  614.             return coordinate;
  615.         }

  616.         /** Interpolate the TEC in latitude and longitude.
  617.          * @param exponent exponent defining the unit of the values listed in the data blocks
  618.          * @param values TEC values
  619.          * @param latitudes array containing the different latitudes in radians
  620.          * @param longitudes array containing the different latitudes in radians
  621.          * @return the interpolating TEC functiopn in TECUnits
  622.          */
  623.         private BilinearInterpolatingFunction interpolateTEC(final ArrayList<Double> values, final double exponent,
  624.                                                              final double[] latitudes, final double[] longitudes) {
  625.             // Array dimensions
  626.             final int dimLat = latitudes.length;
  627.             final int dimLon = longitudes.length;

  628.             // Build the array of TEC data
  629.             final double factor = FastMath.pow(10.0, exponent);
  630.             final double[][] fvalTEC = new double[dimLat][dimLon];
  631.             int index = dimLon * dimLat;
  632.             for (int x = 0; x < dimLat; x++) {
  633.                 for (int y = dimLon - 1; y >= 0; y--) {
  634.                     index = index - 1;
  635.                     fvalTEC[x][y] = values.get(index) * factor;
  636.                 }
  637.             }

  638.             // Build Bilinear Interpolation function
  639.             return new BilinearInterpolatingFunction(latitudes, longitudes, fvalTEC);

  640.         }

  641.     }

  642.     /**
  643.      * Container for IONEX data.
  644.      * <p>
  645.      * The TEC contained in the map is previously interpolated
  646.      * according to the latitude and the longitude given by the user.
  647.      * </p>
  648.      */
  649.     private static class TECMap {

  650.         /** Date of the TEC Map. */
  651.         private AbsoluteDate date;

  652.         /** Interpolated TEC [TECUnits]. */
  653.         private BilinearInterpolatingFunction tec;

  654.         /**
  655.          * Constructor.
  656.          * @param date date of the TEC map
  657.          * @param tec interpolated tec
  658.          */
  659.         TECMap(final AbsoluteDate date, final BilinearInterpolatingFunction tec) {
  660.             this.date = date;
  661.             this.tec  = tec;
  662.         }

  663.     }

  664.     /** Container for a consecutive pair of TEC maps.
  665.      * @since 12.0
  666.      */
  667.     private static class TECMapPair {

  668.         /** First snapshot. */
  669.         private final TECMap first;

  670.         /** Second snapshot. */
  671.         private final TECMap second;

  672.         /** Mean earth radius [m]. */
  673.         private double r0;

  674.         /** Height of the ionospheric single layer [m]. */
  675.         private double h;

  676.         /** Flag for mapping function computation. */
  677.         private boolean mapping;

  678.         /** Simple constructor.
  679.          * @param first first snapshot
  680.          * @param second second snapshot
  681.          * @param r0 mean Earth radius
  682.          * @param h height of the ionospheric layer
  683.          * @param mapping flag for mapping computation
  684.          */
  685.         TECMapPair(final TECMap first, final TECMap second,
  686.                    final double r0, final double h, final boolean mapping) {
  687.             this.first   = first;
  688.             this.second  = second;
  689.             this.r0      = r0;
  690.             this.h       = h;
  691.             this.mapping = mapping;
  692.         }

  693.         /** Get TEC at pierce point.
  694.          * @param date date
  695.          * @param ipp Ionospheric Pierce Point
  696.          * @return TEC
  697.          */
  698.         public double getTEC(final AbsoluteDate date, final GeodeticPoint ipp) {
  699.             // Get the TEC values at the two closest dates
  700.             final AbsoluteDate t1   = first.date;
  701.             final double       tec1 = first.tec.value(ipp.getLatitude(), ipp.getLongitude());
  702.             final AbsoluteDate t2   = second.date;
  703.             final double       tec2 = second.tec.value(ipp.getLatitude(), ipp.getLongitude());
  704.             final double       dt   = t2.durationFrom(t1);

  705.             // Perform temporal interpolation (Ref, Eq. 2)
  706.             return (t2.durationFrom(date) / dt) * tec1 + (date.durationFrom(t1) / dt) * tec2;

  707.         }

  708.         /** Get TEC at pierce point.
  709.          * @param date date
  710.          * @param ipp Ionospheric Pierce Point
  711.          * @param <T> type of the field elements
  712.          * @return TEC
  713.          */
  714.         public <T extends CalculusFieldElement<T>> T getTEC(final FieldAbsoluteDate<T> date, final GeodeticPoint ipp) {

  715.             // Get the TEC values at the two closest dates
  716.             final AbsoluteDate t1   = first.date;
  717.             final double       tec1 = first.tec.value(ipp.getLatitude(), ipp.getLongitude());
  718.             final AbsoluteDate t2   = second.date;
  719.             final double       tec2 = second.tec.value(ipp.getLatitude(), ipp.getLongitude());
  720.             final double       dt   = t2.durationFrom(t1);

  721.             // Perform temporal interpolation (Ref, Eq. 2)
  722.             return date.durationFrom(t2).negate().divide(dt).multiply(tec1).add(date.durationFrom(t1).divide(dt).multiply(tec2));

  723.         }

  724.     }

  725.     /** Container for IONEX header. */
  726.     private static class IONEXHeader {

  727.         /** Number of maps contained in the IONEX file. */
  728.         private int nbOfMaps;

  729.         /** Mean earth radius [m]. */
  730.         private double baseRadius;

  731.         /** Height of the ionospheric single layer [m]. */
  732.         private double hIon;

  733.         /** Flag for mapping function adopted for TEC determination. */
  734.         private boolean isMappingFunction;

  735.         /**
  736.          * Constructor.
  737.          * @param nbOfMaps number of TEC maps contained in the file
  738.          * @param baseRadius mean earth radius in meters
  739.          * @param hIon height of the ionospheric single layer in meters
  740.          * @param mappingFunction flag for mapping function adopted for TEC determination
  741.          */
  742.         IONEXHeader(final int nbOfMaps,
  743.                     final double baseRadius, final double hIon,
  744.                     final boolean mappingFunction) {
  745.             this.nbOfMaps          = nbOfMaps;
  746.             this.baseRadius        = baseRadius;
  747.             this.hIon              = hIon;
  748.             this.isMappingFunction = mappingFunction;
  749.         }

  750.         /**
  751.          * Get the number of TEC maps contained in the file.
  752.          * @return the number of TEC maps
  753.          */
  754.         public int getTECMapsNumer() {
  755.             return nbOfMaps;
  756.         }

  757.         /**
  758.          * Get the mean earth radius in meters.
  759.          * @return the mean earth radius
  760.          */
  761.         public double getEarthRadius() {
  762.             return baseRadius;
  763.         }

  764.         /**
  765.          * Get the height of the ionospheric single layer in meters.
  766.          * @return the height of the ionospheric single layer
  767.          */
  768.         public double getHIon() {
  769.             return hIon;
  770.         }

  771.         /**
  772.          * Get the mapping function flag.
  773.          * @return false if mapping function computation is needed
  774.          */
  775.         public boolean isMappingFunction() {
  776.             return isMappingFunction;
  777.         }

  778.     }

  779. }