GlobalIonosphereMapModel.java

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

  31. import org.hipparchus.Field;
  32. import org.hipparchus.CalculusFieldElement;
  33. import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
  34. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  35. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  36. import org.hipparchus.util.FastMath;
  37. import org.hipparchus.util.MathUtils;
  38. import org.orekit.annotation.DefaultDataContext;
  39. import org.orekit.bodies.GeodeticPoint;
  40. import org.orekit.data.AbstractSelfFeedingLoader;
  41. import org.orekit.data.DataContext;
  42. import org.orekit.data.DataLoader;
  43. import org.orekit.data.DataProvidersManager;
  44. import org.orekit.errors.OrekitException;
  45. import org.orekit.errors.OrekitInternalError;
  46. import org.orekit.errors.OrekitMessages;
  47. import org.orekit.frames.TopocentricFrame;
  48. import org.orekit.propagation.FieldSpacecraftState;
  49. import org.orekit.propagation.SpacecraftState;
  50. import org.orekit.time.AbsoluteDate;
  51. import org.orekit.time.DateTimeComponents;
  52. import org.orekit.time.FieldAbsoluteDate;
  53. import org.orekit.time.TimeComponents;
  54. import org.orekit.time.TimeScale;
  55. import org.orekit.utils.ParameterDriver;

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

  125.     /** Serializable UID. */
  126.     private static final long serialVersionUID = 201928052L;

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

  129.     /** Threshold for latitude and longitude difference. */
  130.     private static final double THRESHOLD = 0.001;

  131.     /** Geodetic site latitude, radians.*/
  132.     private double latitude;

  133.     /** Geodetic site longitude, radians.*/
  134.     private double longitude;

  135.     /** Mean earth radius [m]. */
  136.     private double r0;

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

  139.     /** Time interval between two TEC maps [s]. */
  140.     private double dt;

  141.     /** Number of TEC maps as read on the header of the file. */
  142.     private int nbMaps;

  143.     /** Flag for mapping function computation. */
  144.     private boolean mapping;

  145.     /** Epoch of the first TEC map as read in the header of the IONEX file. */
  146.     private AbsoluteDate startDate;

  147.     /** Epoch of the last TEC map as read in the header of the IONEX file. */
  148.     private AbsoluteDate endDate;

  149.     /** Map of interpolated TEC at a specific date. */
  150.     private Map<AbsoluteDate, Double> tecMap;

  151.     /** UTC time scale. */
  152.     private final TimeScale utc;

  153.     /**
  154.      * Constructor with supported names given by user. This constructor uses the {@link
  155.      * DataContext#getDefault() default data context}.
  156.      *
  157.      * @param supportedNames regular expression that matches the names of the IONEX files
  158.      *                       to be loaded. See {@link DataProvidersManager#feed(String,
  159.      *                       DataLoader)}.
  160.      * @see #GlobalIonosphereMapModel(String, DataProvidersManager, TimeScale)
  161.      */
  162.     @DefaultDataContext
  163.     public GlobalIonosphereMapModel(final String supportedNames) {
  164.         this(supportedNames,
  165.                 DataContext.getDefault().getDataProvidersManager(),
  166.                 DataContext.getDefault().getTimeScales().getUTC());
  167.     }

  168.     /**
  169.      * Constructor that uses user defined supported names and data context.
  170.      *
  171.      * @param supportedNames       regular expression that matches the names of the IONEX
  172.      *                             files to be loaded. See {@link DataProvidersManager#feed(String,
  173.      *                             DataLoader)}.
  174.      * @param dataProvidersManager provides access to auxiliary data files.
  175.      * @param utc                  UTC time scale.
  176.      * @since 10.1
  177.      */
  178.     public GlobalIonosphereMapModel(final String supportedNames,
  179.                                     final DataProvidersManager dataProvidersManager,
  180.                                     final TimeScale utc) {
  181.         super(supportedNames, dataProvidersManager);
  182.         this.latitude       = Double.NaN;
  183.         this.longitude      = Double.NaN;
  184.         this.tecMap         = new HashMap<>();
  185.         this.utc = utc;
  186.     }

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

  219.     @Override
  220.     public double pathDelay(final SpacecraftState state, final TopocentricFrame baseFrame,
  221.                             final double frequency, final double[] parameters) {

  222.         // Elevation in radians
  223.         final Vector3D position  = state.getPVCoordinates(baseFrame).getPosition();
  224.         final double   elevation = position.getDelta();

  225.         // Only consider measures above the horizon
  226.         if (elevation > 0.0) {
  227.             // Date
  228.             final AbsoluteDate date = state.getDate();
  229.             // Geodetic point
  230.             final GeodeticPoint geo = baseFrame.getPoint();
  231.             // Delay
  232.             return pathDelay(date, geo, elevation, frequency);
  233.         }

  234.         return 0.0;

  235.     }

  236.     /**
  237.      * Calculates the ionospheric path delay for the signal path from a ground
  238.      * station to a satellite.
  239.      * <p>
  240.      * The path delay can be computed for any elevation angle.
  241.      * </p>
  242.      * @param <T> type of the elements
  243.      * @param date current date
  244.      * @param geo geodetic point of receiver/station
  245.      * @param elevation elevation of the satellite in radians
  246.      * @param frequency frequency of the signal in Hz
  247.      * @return the path delay due to the ionosphere in m
  248.      */
  249.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldAbsoluteDate<T> date, final GeodeticPoint geo,
  250.                                                        final T elevation, final double frequency) {
  251.         // TEC in TECUnits
  252.         final T tec = getTEC(date, geo);
  253.         // Square of the frequency
  254.         final double freq2 = frequency * frequency;
  255.         // "Slant" Total Electron Content
  256.         final T stec;
  257.         // Check if a mapping factor is needed
  258.         if (mapping) {
  259.             stec = tec;
  260.         } else {
  261.             // Mapping factor
  262.             final T fz = mappingFunction(elevation);
  263.             stec = tec.multiply(fz);
  264.         }
  265.         // Delay computation
  266.         final double alpha  = 40.3e16 / freq2;
  267.         return stec.multiply(alpha);
  268.     }

  269.     @Override
  270.     public <T extends CalculusFieldElement<T>> T pathDelay(final FieldSpacecraftState<T> state, final TopocentricFrame baseFrame,
  271.                                                        final double frequency, final T[] parameters) {

  272.         // Elevation in radians
  273.         final FieldVector3D<T> position = state.getPVCoordinates(baseFrame).getPosition();
  274.         final T elevation = position.getDelta();

  275.         // Only consider measures above the horizon
  276.         if (elevation.getReal() > 0.0) {
  277.             // Date
  278.             final FieldAbsoluteDate<T> date = state.getDate();
  279.             // Geodetic point
  280.             final GeodeticPoint geo = baseFrame.getPoint();
  281.             // Delay
  282.             return pathDelay(date, geo, elevation, frequency);
  283.         }

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

  285.     }

  286.     /**
  287.      * Computes the Total Electron Content (TEC) at a given date by performing a
  288.      * temporal interpolation with the two closest date in the IONEX file.
  289.      * @param date current date
  290.      * @param recPoint geodetic point of receiver/station
  291.      * @return the TEC after a temporal interpolation, in TECUnits
  292.      */
  293.     public double getTEC(final AbsoluteDate date, final GeodeticPoint recPoint) {

  294.         // Load TEC data only if needed
  295.         loadsIfNeeded(recPoint);

  296.         // Check if the date is out of range
  297.         checkDate(date);

  298.         // Date and Time components
  299.         final DateTimeComponents dateTime = date.getComponents(utc);
  300.         // Find the two closest dates of the current date
  301.         final double secInDay   = dateTime.getTime().getSecondsInLocalDay();
  302.         final double ratio      = FastMath.floor(secInDay / dt) * dt;
  303.         final AbsoluteDate tI   = new AbsoluteDate(dateTime.getDate(),
  304.                                                    new TimeComponents(ratio),
  305.                                                    utc);
  306.         final AbsoluteDate tIp1 = tI.shiftedBy(dt);

  307.         // Get the TEC values at the two closest dates
  308.         final double tecI   = tecMap.get(tI);
  309.         final double tecIp1 = tecMap.get(tIp1);

  310.         // Perform temporal interpolation (Ref, Eq. 2)
  311.         final double tec = (tIp1.durationFrom(date) / dt) * tecI + (date.durationFrom(tI) / dt) * tecIp1;
  312.         return tec;
  313.     }

  314.     /**
  315.      * Computes the Total Electron Content (TEC) at a given date by performing a
  316.      * temporal interpolation with the two closest date in the IONEX file.
  317.      * @param <T> type of the elements
  318.      * @param date current date
  319.      * @param recPoint geodetic point of receiver/station
  320.      * @return the TEC after a temporal interpolation, in TECUnits
  321.      */
  322.     public <T extends CalculusFieldElement<T>> T getTEC(final FieldAbsoluteDate<T> date, final GeodeticPoint recPoint) {

  323.         // Load TEC data only if needed
  324.         loadsIfNeeded(recPoint);

  325.         // Check if the date is out of range
  326.         checkDate(date.toAbsoluteDate());

  327.         // Field
  328.         final Field<T> field = date.getField();

  329.         // Date and Time components
  330.         final DateTimeComponents dateTime = date.getComponents(utc);
  331.         // Find the two closest dates of the current date
  332.         final double secInDay           = dateTime.getTime().getSecondsInLocalDay();
  333.         final double ratio              = FastMath.floor(secInDay / dt) * dt;
  334.         final FieldAbsoluteDate<T> tI   = new FieldAbsoluteDate<>(field, dateTime.getDate(),
  335.                                                                   new TimeComponents(ratio),
  336.                                                                   utc);
  337.         final FieldAbsoluteDate<T> tIp1 = tI.shiftedBy(dt);

  338.         // Get the TEC values at the two closest dates
  339.         final double tecI   = tecMap.get(tI.toAbsoluteDate());
  340.         final double tecIp1 = tecMap.get(tIp1.toAbsoluteDate());

  341.         // Perform temporal interpolation (Ref, Eq. 2)
  342.         final T tec = tIp1.durationFrom(date).divide(dt).multiply(tecI).add(date.durationFrom(tI).divide(dt).multiply(tecIp1));
  343.         return tec;
  344.     }

  345.     @Override
  346.     public List<ParameterDriver> getParametersDrivers() {
  347.         return Collections.emptyList();
  348.     }

  349.     /**
  350.      * Computes the ionospheric mapping function.
  351.      * @param elevation the elevation of the satellite in radians
  352.      * @return the mapping function
  353.      */
  354.     private double mappingFunction(final double elevation) {
  355.         // Calculate the zenith angle from the elevation
  356.         final double z = FastMath.abs(0.5 * FastMath.PI - elevation);
  357.         // Distance ratio
  358.         final double ratio = r0 / (r0 + h);
  359.         // Mapping function
  360.         final double coef = FastMath.sin(z) * ratio;
  361.         final double fz = 1.0 / FastMath.sqrt(1.0 - coef * coef);
  362.         return fz;
  363.     }

  364.     /**
  365.      * Computes the ionospheric mapping function.
  366.      * @param <T> type of the elements
  367.      * @param elevation the elevation of the satellite in radians
  368.      * @return the mapping function
  369.      */
  370.     private <T extends CalculusFieldElement<T>> T mappingFunction(final T elevation) {
  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.     /**
  381.      * Lazy loading of TEC data.
  382.      * @param recPoint geodetic point of receiver/station
  383.      */
  384.     private void loadsIfNeeded(final GeodeticPoint recPoint) {

  385.         // Current latitude and longitude of the geodetic point
  386.         final double lat = recPoint.getLatitude();
  387.         final double lon = MathUtils.normalizeAngle(recPoint.getLongitude(), 0.0);

  388.         // Read the file only if the TEC map is empty or if the geodetic point displacement is
  389.         // greater than 0.001 radians (in latitude or longitude)
  390.         if (tecMap.isEmpty() || FastMath.abs(lat - latitude) > THRESHOLD ||  FastMath.abs(lon - longitude) > THRESHOLD) {
  391.             this.latitude  = lat;
  392.             this.longitude = lon;

  393.             // Read file
  394.             final Parser parser = new Parser();
  395.             feed(parser);

  396.             // File header
  397.             final IONEXHeader top = parser.getIONEXHeader();
  398.             this.startDate  = top.getFirstDate();
  399.             this.endDate    = top.getLastDate();
  400.             this.dt         = top.getInterval();
  401.             this.nbMaps     = top.getTECMapsNumer();
  402.             this.r0         = top.getEarthRadius();
  403.             this.h          = top.getHIon();
  404.             this.mapping    = top.isMappingFunction();

  405.             // TEC map
  406.             for (TECMap map : parser.getTECMaps()) {
  407.                 tecMap.put(map.getDate(), map.getTEC());
  408.             }
  409.         }
  410.         checkSize();
  411.     }

  412.     /**
  413.      * Check if the current date is between the startDate and
  414.      * the endDate of the IONEX file.
  415.      * @param date current date
  416.      */
  417.     private void checkDate(final AbsoluteDate date) {
  418.         if (startDate.durationFrom(date) > 0 || date.durationFrom(endDate) > 0) {
  419.             throw new OrekitException(OrekitMessages.NO_TEC_DATA_IN_FILE_FOR_DATE,
  420.                     getSupportedNames(), date);
  421.         }
  422.     }

  423.     /**
  424.      * Check if the number of parsed TEC maps is consistent with the header specification.
  425.      */
  426.     private void checkSize() {
  427.         if (tecMap.size() != nbMaps) {
  428.             throw new OrekitException(OrekitMessages.INCONSISTENT_NUMBER_OF_TEC_MAPS_IN_FILE, tecMap.size(), nbMaps);
  429.         }
  430.     }

  431.     /** Replace the instance with a data transfer object for serialization.
  432.      * @return data transfer object that will be serialized
  433.      */
  434.     @DefaultDataContext
  435.     private Object writeReplace() {
  436.         return new DataTransferObject(getSupportedNames());
  437.     }

  438.     /** Parser for IONEX files. */
  439.     private class Parser implements DataLoader {

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

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

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

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

  448.         /** Header of the IONEX file. */
  449.         private IONEXHeader header;

  450.         /** List of TEC Maps. */
  451.         private List<TECMap> maps;

  452.         @Override
  453.         public boolean stillAcceptsData() {
  454.             return true;
  455.         }

  456.         @Override
  457.         public void loadData(final InputStream input, final String name)
  458.             throws IOException,  ParseException {

  459.             maps = new ArrayList<>();

  460.             // Open stream and parse data
  461.             int   lineNumber = 0;
  462.             String line      = null;
  463.             try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
  464.                  BufferedReader    br = new BufferedReader(isr)) {

  465.                 // Placeholders for parsed data
  466.                 int               interval    = 3600;
  467.                 int               nbOfMaps    = 1;
  468.                 int               exponent    = -1;
  469.                 double            baseRadius  = 6371.0e3;
  470.                 double            hIon        = 350e3;
  471.                 boolean           mappingF    = false;
  472.                 boolean           inTEC       = false;
  473.                 double[]          latitudes   = null;
  474.                 double[]          longitudes  = null;
  475.                 AbsoluteDate      firstEpoch  = null;
  476.                 AbsoluteDate      lastEpoch   = null;
  477.                 AbsoluteDate      epoch       = firstEpoch;
  478.                 ArrayList<Double> values      = new ArrayList<>();

  479.                 for (line = br.readLine(); line != null; line = br.readLine()) {
  480.                     ++lineNumber;
  481.                     if (line.length() > LABEL_START) {
  482.                         switch (line.substring(LABEL_START).trim()) {
  483.                             case "EPOCH OF FIRST MAP" :
  484.                                 firstEpoch = parseDate(line);
  485.                                 break;
  486.                             case "EPOCH OF LAST MAP" :
  487.                                 lastEpoch = parseDate(line);
  488.                                 break;
  489.                             case "INTERVAL" :
  490.                                 interval = parseInt(line, 2, 4);
  491.                                 break;
  492.                             case "# OF MAPS IN FILE" :
  493.                                 nbOfMaps = parseInt(line, 2, 4);
  494.                                 break;
  495.                             case "BASE RADIUS" :
  496.                                 // Value is in kilometers
  497.                                 baseRadius = parseDouble(line, 2, 6) * KM_TO_M;
  498.                                 break;
  499.                             case "MAPPING FUNCTION" :
  500.                                 mappingF = !parseString(line, 2, 4).equals("NONE");
  501.                                 break;
  502.                             case "EXPONENT" :
  503.                                 exponent = parseInt(line, 4, 2);
  504.                                 break;
  505.                             case "HGT1 / HGT2 / DHGT" :
  506.                                 if (parseDouble(line, 17, 3) == 0.0) {
  507.                                     // Value is in kilometers
  508.                                     hIon = parseDouble(line, 3, 5) * KM_TO_M;
  509.                                 }
  510.                                 break;
  511.                             case "LAT1 / LAT2 / DLAT" :
  512.                                 latitudes = parseCoordinate(line);
  513.                                 break;
  514.                             case "LON1 / LON2 / DLON" :
  515.                                 longitudes = parseCoordinate(line);
  516.                                 break;
  517.                             case "END OF HEADER" :
  518.                                 // Check that latitude and longitude bondaries were found
  519.                                 if (latitudes == null || longitudes == null) {
  520.                                     throw new OrekitException(OrekitMessages.NO_LATITUDE_LONGITUDE_BONDARIES_IN_IONEX_HEADER, getSupportedNames());
  521.                                 }
  522.                                 // Check that first and last epochs were found
  523.                                 if (firstEpoch == null || lastEpoch == null) {
  524.                                     throw new OrekitException(OrekitMessages.NO_EPOCH_IN_IONEX_HEADER, getSupportedNames());
  525.                                 }
  526.                                 // At the end of the header, we build the IONEXHeader object
  527.                                 header = new IONEXHeader(firstEpoch, lastEpoch, interval, nbOfMaps,
  528.                                                          baseRadius, hIon, mappingF);
  529.                                 break;
  530.                             case "START OF TEC MAP" :
  531.                                 inTEC = true;
  532.                                 break;
  533.                             case END :
  534.                                 final double tec = interpolateTEC(values, exponent, latitudes, longitudes);
  535.                                 final TECMap map = new TECMap(epoch, tec);
  536.                                 maps.add(map);
  537.                                 // Reset parameters
  538.                                 inTEC  = false;
  539.                                 values = new ArrayList<>();
  540.                                 epoch  = null;
  541.                                 break;
  542.                             default :
  543.                                 if (inTEC) {
  544.                                     // Date
  545.                                     if (line.endsWith(EPOCH)) {
  546.                                         epoch = parseDate(line);
  547.                                     }
  548.                                     // Fill TEC values list
  549.                                     if (!line.endsWith("LAT/LON1/LON2/DLON/H") &&
  550.                                         !line.endsWith(END) &&
  551.                                         !line.endsWith(EPOCH)) {
  552.                                         line = line.trim();
  553.                                         final String[] readLine = SEPARATOR.split(line);
  554.                                         for (final String s : readLine) {
  555.                                             values.add(Double.valueOf(s));
  556.                                         }
  557.                                     }
  558.                                 }
  559.                                 break;
  560.                         }
  561.                     } else {
  562.                         if (inTEC) {
  563.                             // Here, we are parsing the last line of TEC data for a given latitude
  564.                             // The size of this line is lower than 60.
  565.                             line = line.trim();
  566.                             final String[] readLine = SEPARATOR.split(line);
  567.                             for (final String s : readLine) {
  568.                                 values.add(Double.valueOf(s));
  569.                             }
  570.                         }
  571.                     }

  572.                 }

  573.                 // Close the stream after reading
  574.                 input.close();

  575.             } catch (NumberFormatException nfe) {
  576.                 throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  577.                                           lineNumber, name, line);
  578.             }

  579.         }

  580.         /**
  581.          * Get the header of the IONEX file.
  582.          * @return the header of the IONEX file
  583.          */
  584.         public IONEXHeader getIONEXHeader() {
  585.             return header;
  586.         }

  587.         /**
  588.          * Get the list of the TEC maps.
  589.          * @return the list of TEC maps.
  590.          */
  591.         public List<TECMap> getTECMaps() {
  592.             return maps;
  593.         }

  594.         /** Extract a string from a line.
  595.          * @param line to parse
  596.          * @param start start index of the string
  597.          * @param length length of the string
  598.          * @return parsed string
  599.          */
  600.         private String parseString(final String line, final int start, final int length) {
  601.             return line.substring(start, FastMath.min(line.length(), start + length)).trim();
  602.         }

  603.         /** Extract an integer from a line.
  604.          * @param line to parse
  605.          * @param start start index of the integer
  606.          * @param length length of the integer
  607.          * @return parsed integer
  608.          */
  609.         private int parseInt(final String line, final int start, final int length) {
  610.             return Integer.parseInt(parseString(line, start, length));
  611.         }

  612.         /** Extract a double from a line.
  613.          * @param line to parse
  614.          * @param start start index of the real
  615.          * @param length length of the real
  616.          * @return parsed real
  617.          */
  618.         private double parseDouble(final String line, final int start, final int length) {
  619.             return Double.parseDouble(parseString(line, start, length));
  620.         }

  621.         /** Extract a date from a parsed line.
  622.          * @param line to parse
  623.          * @return an absolute date
  624.          */
  625.         private AbsoluteDate parseDate(final String line) {
  626.             return new AbsoluteDate(parseInt(line, 0, 6),
  627.                                     parseInt(line, 6, 6),
  628.                                     parseInt(line, 12, 6),
  629.                                     parseInt(line, 18, 6),
  630.                                     parseInt(line, 24, 6),
  631.                                     parseDouble(line, 30, 13),
  632.                                     utc);
  633.         }

  634.         /** Build the coordinate array from a parsed line.
  635.          * @param line to parse
  636.          * @return an array of coordinates in radians
  637.          */
  638.         private double[] parseCoordinate(final String line) {
  639.             final double a = parseDouble(line, 2, 6);
  640.             final double b = parseDouble(line, 8, 6);
  641.             final double c = parseDouble(line, 14, 6);
  642.             final double[] coordinate = new double[((int) FastMath.abs((a - b) / c)) + 1];
  643.             int i = 0;
  644.             for (double cor = FastMath.min(a, b); cor <= FastMath.max(a, b); cor += FastMath.abs(c)) {
  645.                 coordinate[i] = FastMath.toRadians(cor);
  646.                 i++;
  647.             }
  648.             return coordinate;
  649.         }

  650.         /** Interpolate the TEC in latitude and longitude.
  651.          * @param exponent exponent defining the unit of the values listed in the data blocks
  652.          * @param values TEC values
  653.          * @param latitudes array containing the different latitudes in radians
  654.          * @param longitudes array containing the different latitudes in radians
  655.          * @return the interpolated TEC in TECUnits
  656.          */
  657.         private double interpolateTEC(final ArrayList<Double> values, final double exponent,
  658.                                       final double[] latitudes, final double[] longitudes) {
  659.             // Array dimensions
  660.             final int dimLat = latitudes.length;
  661.             final int dimLon = longitudes.length;

  662.             // Build the array of TEC data
  663.             final double[][] fvalTEC = new double[dimLat][dimLon];
  664.             int index = dimLon * dimLat;
  665.             for (int x = 0; x < dimLat; x++) {
  666.                 for (int y = dimLon - 1; y >= 0; y--) {
  667.                     index = index - 1;
  668.                     fvalTEC[x][y] = values.get(index);
  669.                 }
  670.             }

  671.             // Build Bilinear Interpolation function
  672.             final BilinearInterpolatingFunction functionTEC = new BilinearInterpolatingFunction(latitudes, longitudes, fvalTEC);
  673.             final double tec = functionTEC.value(latitude, longitude) * FastMath.pow(10.0, exponent);
  674.             return tec;
  675.         }
  676.     }

  677.     /**
  678.      * Container for IONEX data.
  679.      * <p>
  680.      * The TEC contained in the map is previously interpolated
  681.      * according to the latitude and the longitude given by the user.
  682.      * </p>
  683.      */
  684.     private static class TECMap {

  685.         /** Date of the TEC Map. */
  686.         private AbsoluteDate date;

  687.         /** Interpolated TEC [TECUnits]. */
  688.         private double tec;

  689.         /**
  690.          * Constructor.
  691.          * @param date date of the TEC map
  692.          * @param tec interpolated tec
  693.          */
  694.         TECMap(final AbsoluteDate date, final double tec) {
  695.             this.date = date;
  696.             this.tec  = tec;
  697.         }

  698.         /**
  699.          * Get the date of the TEC map.
  700.          * @return the date
  701.          */
  702.         public AbsoluteDate getDate() {
  703.             return date;
  704.         }

  705.         /**
  706.          * Get the value of the interpolated TEC.
  707.          * @return the TEC in TECUnits
  708.          */
  709.         public double getTEC() {
  710.             return tec;
  711.         }

  712.     }

  713.     /** Container for IONEX header. */
  714.     private static class IONEXHeader {

  715.         /** Epoch of the first TEC map. */
  716.         private AbsoluteDate firstDate;

  717.         /** Epoch of the last TEC map. */
  718.         private AbsoluteDate lastDate;

  719.         /** Interval between two maps [s]. */
  720.         private int interval;

  721.         /** Number of maps contained in the IONEX file. */
  722.         private int nbOfMaps;

  723.         /** Mean earth radius [m]. */
  724.         private double baseRadius;

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

  727.         /** Flag for mapping function adopted for TEC determination. */
  728.         private boolean isMappingFunction;

  729.         /**
  730.          * Constructor.
  731.          * @param firstDate epoch of the first TEC map.
  732.          * @param lastDate epoch of the last TEC map.
  733.          * @param nbOfMaps number of TEC maps contained in the file
  734.          * @param interval number of seconds between two tec maps.
  735.          * @param baseRadius mean earth radius in meters
  736.          * @param hIon height of the ionospheric single layer in meters
  737.          * @param mappingFunction flag for mapping function adopted for TEC determination
  738.          */
  739.         IONEXHeader(final AbsoluteDate firstDate, final AbsoluteDate lastDate,
  740.                     final int interval, final int nbOfMaps,
  741.                     final double baseRadius, final double hIon,
  742.                     final boolean mappingFunction) {
  743.             this.firstDate         = firstDate;
  744.             this.lastDate          = lastDate;
  745.             this.interval          = interval;
  746.             this.nbOfMaps          = nbOfMaps;
  747.             this.baseRadius        = baseRadius;
  748.             this.hIon              = hIon;
  749.             this.isMappingFunction = mappingFunction;
  750.         }

  751.         /**
  752.          * Get the first date of the IONEX file.
  753.          * @return the first date of the IONEX file
  754.          */
  755.         public AbsoluteDate getFirstDate() {
  756.             return firstDate;
  757.         }

  758.         /**
  759.          * Get the last date of the IONEX file.
  760.          * @return the last date of the IONEX file
  761.          */
  762.         public AbsoluteDate getLastDate() {
  763.             return lastDate;
  764.         }

  765.         /**
  766.          * Get the time interval between two TEC maps.
  767.          * @return the interval between two TEC maps
  768.          */
  769.         public int getInterval() {
  770.             return interval;
  771.         }

  772.         /**
  773.          * Get the number of TEC maps contained in the file.
  774.          * @return the number of TEC maps
  775.          */
  776.         public int getTECMapsNumer() {
  777.             return nbOfMaps;
  778.         }

  779.         /**
  780.          * Get the mean earth radius in meters.
  781.          * @return the mean earth radius
  782.          */
  783.         public double getEarthRadius() {
  784.             return baseRadius;
  785.         }

  786.         /**
  787.          * Get the height of the ionospheric single layer in meters.
  788.          * @return the height of the ionospheric single layer
  789.          */
  790.         public double getHIon() {
  791.             return hIon;
  792.         }

  793.         /**
  794.          * Get the mapping function flag.
  795.          * @return false if mapping function computation is needed
  796.          */
  797.         public boolean isMappingFunction() {
  798.             return isMappingFunction;
  799.         }

  800.     }

  801.     /** Internal class used only for serialization. */
  802.     @DefaultDataContext
  803.     private static class DataTransferObject implements Serializable {

  804.         /** Serializable UID. */
  805.         private static final long serialVersionUID = 201928052L;

  806.         /** Regular expression that matches the names of the IONEX files. */
  807.         private final String supportedNames;

  808.         /** Simple constructor.
  809.          * @param supportedNames regular expression that matches the names of the IONEX files
  810.          */
  811.         DataTransferObject(final String supportedNames) {
  812.             this.supportedNames = supportedNames;
  813.         }

  814.         /** Replace the deserialized data transfer object with a {@link GlobalIonosphereMapModel}.
  815.          * @return replacement {@link GlobalIonosphereMapModel}
  816.          */
  817.         private Object readResolve() {
  818.             try {
  819.                 return new GlobalIonosphereMapModel(supportedNames);
  820.             } catch (OrekitException oe) {
  821.                 throw new OrekitInternalError(oe);
  822.             }
  823.         }

  824.     }

  825. }