GlobalPressureTemperature2Model.java

/* Copyright 2002-2019 CS Systèmes d'Information
 * Licensed to CS Systèmes d'Information (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.models.earth;

import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.nio.charset.StandardCharsets;
import java.text.ParseException;
import java.util.ArrayList;
import java.util.List;
import java.util.SortedSet;
import java.util.TreeSet;
import java.util.concurrent.atomic.AtomicReference;
import java.util.function.ToDoubleFunction;

import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.SinCos;
import org.orekit.data.DataLoader;
import org.orekit.data.DataProvidersManager;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;

/** The Global Pressure and Temperature 2 (GPT2) model.
 * This model is an empirical model that provides the temperature, the pressure and the water vapor pressure
 * of a site depending its latitude and  longitude. This model also provides the a<sub>h</sub>
 * and a<sub>w</sub> coefficients used for the {@link ViennaOneModel Vienna 1} model.
 * <p>
 * The requisite coefficients for the computation of the weather parameters are provided by the
 * Department of Geodesy and Geoinformation of the Vienna University. They are based on an
 * external grid file like "gpt2_1.grd" (1° x 1°) or "gpt2_5.grd" (5° x 5°) available at:
 * <a href="http://vmf.geo.tuwien.ac.at/codes/"> link</a>
 * </p>
 * <p>
 * A bilinear interpolation is performed in order to obtained the correct values of the weather parameters.
 * </p>
 * <p>
 * The format is always the same, with and example shown below for the pressure and the temperature.
 * <p>
 * Example:
 * </p>
 * <pre>
 * %  lat    lon   p:a0    A1   B1   A2   B2  T:a0    A1   B1   A2   B2
 *   87.5    2.5 101421    21  409 -217 -122 259.2 -13.2 -6.1  2.6  0.3
 *   87.5    7.5 101416    21  411 -213 -120 259.3 -13.1 -6.1  2.6  0.3
 *   87.5   12.5 101411    22  413 -209 -118 259.3 -13.1 -6.1  2.6  0.3
 *   87.5   17.5 101407    23  415 -205 -116 259.4 -13.0 -6.1  2.6  0.3
 *   ...
 * </pre>
 *
 * @see K. Lagler, M. Schindelegger, J. Böhm, H. Krasna, T. Nilsson (2013),
 * GPT2: empirical slant delay model for radio space geodetic techniques. Geophys
 * Res Lett 40(6):1069–1073. doi:10.1002/grl.50288
 *
 * @author Bryan Cazabonne
 *
 */
public class GlobalPressureTemperature2Model implements WeatherModel {

    /** Default supported files name pattern. */
    public static final String DEFAULT_SUPPORTED_NAMES = "gpt2_\\d+.grd";

    /** Standard gravity constant [m/s²]. */
    private static final double G = Constants.G0_STANDARD_GRAVITY;

    /** Ideal gas constant for dry air [J/kg/K]. */
    private static final double R = 287.0;

    /** Conversion factor from degrees to mill arcseconds. */
    private static final int DEG_TO_MAS = 3600000;

    /** Shared lazily loaded grid. */
    private static final AtomicReference<Grid> SHARED_GRID = new AtomicReference<>(null);

    /** South-West grid entry. */
    private final GridEntry southWest;

    /** South-East grid entry. */
    private final GridEntry southEast;

    /** North-West grid entry. */
    private final GridEntry northWest;

    /** North-East grid entry. */
    private final GridEntry northEast;

    /** The hydrostatic and wet a coefficients loaded. */
    private double[] coefficientsA;

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

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

    /** Temperature site, in kelvins. */
    private double temperature;

    /** Pressure site, in hPa. */
    private double pressure;

    /** water vapour pressure, in hPa. */
    private double e0;

    /** The height of the station in m. */
    private double height;

    /** Geoid used to compute the undulations. */
    private final Geoid geoid;

    /** Current date. */
    private AbsoluteDate date;

    /** Constructor with supported names given by user.
     * @param supportedNames supported names
     * @param latitude geodetic latitude of the station, in radians
     * @param longitude longitude geodetic longitude of the station, in radians
     * @param geoid level surface of the gravity potential of a body
     */
    public GlobalPressureTemperature2Model(final String supportedNames, final double latitude,
                                           final double longitude, final Geoid geoid) {
        this.coefficientsA = null;
        this.temperature   = Double.NaN;
        this.pressure      = Double.NaN;
        this.e0            = Double.NaN;
        this.geoid         = geoid;
        this.latitude      = latitude;

        // get the lazily loaded shared grid
        Grid grid = SHARED_GRID.get();
        if (grid == null) {
            // this is the first instance we create, we need to load the grid data
            final Parser parser = new Parser();
            DataProvidersManager.getInstance().feed(supportedNames, parser);
            SHARED_GRID.compareAndSet(null, parser.grid);
            grid = parser.grid;
        }

        // Normalize longitude according to the grid
        this.longitude = MathUtils.normalizeAngle(longitude, grid.entries[0][0].longitude + FastMath.PI);

        final int southIndex = grid.getSouthIndex(this.latitude);
        final int westIndex  = grid.getWestIndex(this.longitude);
        this.southWest = grid.entries[southIndex    ][westIndex    ];
        this.southEast = grid.entries[southIndex    ][westIndex + 1];
        this.northWest = grid.entries[southIndex + 1][westIndex    ];
        this.northEast = grid.entries[southIndex + 1][westIndex + 1];

    }

    /** Constructor with default supported names.
     * @param latitude geodetic latitude of the station, in radians
     * @param longitude geodetic latitude of the station, in radians
     * @param geoid level surface of the gravity potential of a body
     */
    public GlobalPressureTemperature2Model(final double latitude, final double longitude, final Geoid geoid) {
        this(DEFAULT_SUPPORTED_NAMES, latitude, longitude, geoid);
    }

    /** Returns the a coefficients array.
     * <ul>
     * <li>double[0] = a<sub>h</sub>
     * <li>double[1] = a<sub>w</sub>
     * </ul>
     * @return the a coefficients array
     */
    public double[] getA() {
        return coefficientsA.clone();
    }

    /** Returns the temperature at the station [K].
     * @return the temperature at the station [K]
     */
    public double getTemperature() {
        return temperature;
    }

    /** Returns the pressure at the station [hPa].
     * @return the pressure at the station [hPa]
     */
    public double getPressure() {
        return pressure;
    }

    /** Returns the water vapor pressure at the station [hPa].
     * @return the water vapor pressure at the station [hPa]
     */
    public double getWaterVaporPressure() {
        return e0;
    }

    @Override
    public void weatherParameters(final double stationHeight, final AbsoluteDate currentDate) {

        this.date      = currentDate;
        this.height    = stationHeight;

        final int dayOfYear = currentDate.getComponents(TimeScalesFactory.getUTC()).getDate().getDayOfYear();

        // ah and aw coefficients
        coefficientsA = new double[] {
            interpolate(e -> evaluate(dayOfYear, e.ah)) * 0.001,
            interpolate(e -> evaluate(dayOfYear, e.aw)) * 0.001
        };

        // Corrected height (can be negative)
        final double undu            = geoid.getUndulation(latitude, longitude, date);
        final double correctedheight = height - undu - interpolate(e -> e.hS);

        // Temperature gradient [K/m]
        final double dTdH = interpolate(e -> evaluate(dayOfYear, e.dT)) * 0.001;

        // Specific humidity
        final double qv = interpolate(e -> evaluate(dayOfYear, e.qv0)) * 0.001;

        // For the computation of the temperature and the pressure, we use
        // the standard ICAO atmosphere formulas.

        // Temperature [K]
        final double t0 = interpolate(e -> evaluate(dayOfYear, e.temperature0));
        this.temperature = t0 + dTdH * correctedheight;

        // Pressure [hPa]
        final double p0 = interpolate(e -> evaluate(dayOfYear, e.pressure0));
        final double exponent = G / (dTdH * R);
        this.pressure = p0 * FastMath.pow(1 - (dTdH / t0) * correctedheight, exponent) * 0.01;

        // Water vapor pressure [hPa]
        this.e0 = qv * pressure / (0.622 + 0.378 * qv);

    }

    /** Interpolate a grid function.
     * @param gridGetter getter for the grid function
     * @return interpolated function"
     */
    private double interpolate(final ToDoubleFunction<GridEntry> gridGetter) {

        // cell surrounding the point
        final double[] xVal = new double[] {
            southWest.longitude, southEast.longitude
        };
        final double[] yVal = new double[] {
            southWest.latitude, northWest.latitude
        };

        // evaluate grid points at specified day
        final double[][] fval = new double[][] {
            {
                gridGetter.applyAsDouble(southWest),
                gridGetter.applyAsDouble(northWest)
            }, {
                gridGetter.applyAsDouble(southEast),
                gridGetter.applyAsDouble(northEast)
            }
        };

        // perform interpolation in the grid
        return new BilinearInterpolatingFunction(xVal, yVal, fval).value(longitude, latitude);

    }

    /** Evaluate a model for some day.
     * @param dayOfYear day to evaluate
     * @param model model array
     * @return model value at specified day
     */
    private double evaluate(final int dayOfYear, final double[] model) {

        final double coef = (dayOfYear / 365.25) * 2 * FastMath.PI;
        final SinCos sc1  = FastMath.sinCos(coef);
        final SinCos sc2  = FastMath.sinCos(2.0 * coef);

        return model[0] +
               model[1] * sc1.cos() + model[2] * sc1.sin() +
               model[3] * sc2.cos() + model[4] * sc2.sin();

    }

    /** Parser for GPT2 grid files. */
    private static class Parser implements DataLoader {

        /** Grid entries. */
        private Grid grid;

        @Override
        public boolean stillAcceptsData() {
            return grid == null;
        }

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

            final SortedSet<Integer> latSample = new TreeSet<>();
            final SortedSet<Integer> lonSample = new TreeSet<>();
            final List<GridEntry>    entries   = new ArrayList<>();

            // Open stream and parse data
            int   lineNumber = 0;
            String line      = null;
            try (InputStreamReader isr = new InputStreamReader(input, StandardCharsets.UTF_8);
                 BufferedReader    br = new BufferedReader(isr)) {
                final String splitter = "\\s+";

                for (line = br.readLine(); line != null; line = br.readLine()) {
                    ++lineNumber;
                    line = line.trim();

                    // read grid data
                    if (line.length() > 0 && !line.startsWith("%")) {
                        final GridEntry entry = new GridEntry(line.split(splitter));
                        latSample.add(entry.latKey);
                        lonSample.add(entry.lonKey);
                        entries.add(entry);
                    }

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

            // organize entries in a grid that wraps arouns Earth in longitude
            grid = new Grid(latSample, lonSample, entries, name);

        }

    }

    /** Container for complete grid. */
    private static class Grid {

        /** Latitude sample. */
        private final SortedSet<Integer> latitudeSample;

        /** Longitude sample. */
        private final SortedSet<Integer> longitudeSample;

        /** Grid entries. */
        private final GridEntry[][] entries;

        /** Simple constructor.
         * @param latitudeSample latitude sample
         * @param longitudeSample longitude sample
         * @param loadedEntries loaded entries, organized as a simple list
         * @param name file name
         */
        Grid(final SortedSet<Integer> latitudeSample, final SortedSet<Integer> longitudeSample,
             final List<GridEntry> loadedEntries, final String name) {

            final int nA         = latitudeSample.size();
            final int nO         = longitudeSample.size() + 1; // we add one here for wrapping the grid
            this.entries         = new GridEntry[nA][nO];
            this.latitudeSample  = latitudeSample;
            this.longitudeSample = longitudeSample;

            // organize entries in the regular grid
            for (final GridEntry entry : loadedEntries) {
                final int latitudeIndex  = latitudeSample.headSet(entry.latKey + 1).size() - 1;
                final int longitudeIndex = longitudeSample.headSet(entry.lonKey + 1).size() - 1;
                entries[latitudeIndex][longitudeIndex] = entry;
            }

            // finalize the grid
            for (final GridEntry[] row : entries) {

                // check for missing entries
                for (int longitudeIndex = 0; longitudeIndex < nO - 1; ++longitudeIndex) {
                    if (row[longitudeIndex] == null) {
                        throw new OrekitException(OrekitMessages.IRREGULAR_OR_INCOMPLETE_GRID, name);
                    }
                }

                // wrap the grid around the Earth in longitude
                row[nO - 1] = new GridEntry(row[0].latitude, row[0].latKey,
                                            row[0].longitude + 2 * FastMath.PI,
                                            row[0].lonKey + DEG_TO_MAS * 360,
                                            row[0].hS, row[0].pressure0, row[0].temperature0,
                                            row[0].qv0, row[0].dT, row[0].ah, row[0].aw);

            }

        }

        /** Get index of South entries in the grid.
         * @param latitude latitude to locate (radians)
         * @return index of South entries in the grid
         */
        public int getSouthIndex(final double latitude) {

            final int latKey = (int) FastMath.rint(FastMath.toDegrees(latitude) * DEG_TO_MAS);
            final int index  = latitudeSample.headSet(latKey + 1).size() - 1;

            // make sure we have at least one point remaining on North by clipping to size - 2
            return FastMath.min(index, latitudeSample.size() - 2);

        }

        /** Get index of West entries in the grid.
         * @param longitude longitude to locate (radians)
         * @return index of West entries in the grid
         */
        public int getWestIndex(final double longitude) {

            final int lonKey = (int) FastMath.rint(FastMath.toDegrees(longitude) * DEG_TO_MAS);
            final int index  = longitudeSample.headSet(lonKey + 1).size() - 1;

            // we don't do clipping in longitude because we have added a row to wrap around the Earth
            return index;

        }

    }

    /** Container for grid entries. */
    private static class GridEntry {

        /** Latitude (radian). */
        private final double latitude;

        /** Latitude key (mas). */
        private final int latKey;

        /** Longitude (radian). */
        private final double longitude;

        /** Longitude key (mas). */
        private final int lonKey;

        /** Height correction. */
        private final double hS;

        /** Pressure model. */
        private final double[] pressure0;

        /** Temperature model. */
        private final double[] temperature0;

        /** Specific humidity model. */
        private final double[] qv0;

        /** Temperature gradient model. */
        private final double[] dT;

        /** ah coefficient model. */
        private final double[] ah;

        /** aw coefficient model. */
        private final double[] aw;

        /** Build an entry from a parsed line.
         * @param fields line fields
         */
        GridEntry(final String[] fields) {

            final double latDegree = Double.parseDouble(fields[0]);
            final double lonDegree = Double.parseDouble(fields[1]);
            latitude     = FastMath.toRadians(latDegree);
            longitude    = FastMath.toRadians(lonDegree);
            latKey       = (int) FastMath.rint(latDegree * DEG_TO_MAS);
            lonKey       = (int) FastMath.rint(lonDegree * DEG_TO_MAS);

            hS           = Double.valueOf(fields[23]);

            pressure0    = createModel(fields, 2);
            temperature0 = createModel(fields, 7);
            qv0          = createModel(fields, 12);
            dT           = createModel(fields, 17);
            ah           = createModel(fields, 24);
            aw           = createModel(fields, 29);

        }

        /** Build an entry from its components.
         * @param latitude latitude (radian)
         * @param latKey latitude key (mas)
         * @param longitude longitude (radian)
         * @param lonKey longitude key (mas)
         * @param hS height correction
         * @param pressure0 pressure model
         * @param temperature0 temperature model
         * @param qv0 specific humidity model
         * @param dT temperature gradient model
         * @param ah ah coefficient model
         * @param aw aw coefficient model
         */
        GridEntry(final double latitude, final int latKey, final double longitude, final int lonKey,
                  final double hS, final double[] pressure0, final double[] temperature0,
                  final double[] qv0, final double[] dT, final double[] ah, final double[] aw) {

            this.latitude     = latitude;
            this.latKey       = latKey;
            this.longitude    = longitude;
            this.lonKey       = lonKey;
            this.hS           = hS;
            this.pressure0    = pressure0.clone();
            this.temperature0 = temperature0.clone();
            this.qv0          = qv0.clone();
            this.dT           = dT.clone();
            this.ah           = ah.clone();
            this.aw           = aw.clone();
        }

        /** Create a time model array.
         * @param fields line fields
         * @param first index of the first component of the model
         * @return time model array
         */
        private double[] createModel(final String[] fields, final int first) {
            return new double[] {
                Double.parseDouble(fields[first    ]),
                Double.parseDouble(fields[first + 1]),
                Double.parseDouble(fields[first + 2]),
                Double.parseDouble(fields[first + 3]),
                Double.parseDouble(fields[first + 4])
            };
        }

    }

}