SaastamoinenModel.java

  1. /* Copyright 2011-2012 Space Applications Services
  2.  * Licensed to CS Communication & Systèmes (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;

  19. import java.io.Serializable;
  20. import java.util.Arrays;

  22. import org.hipparchus.analysis.BivariateFunction;
  23. import org.hipparchus.analysis.UnivariateFunction;
  24. import org.hipparchus.analysis.interpolation.LinearInterpolator;
  25. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  26. import org.hipparchus.exception.LocalizedCoreFormats;
  27. import org.hipparchus.exception.MathIllegalArgumentException;
  28. import org.hipparchus.util.FastMath;
  29. import org.hipparchus.util.MathArrays;
  30. import org.orekit.data.DataProvidersManager;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitMessages;
  33. import org.orekit.utils.InterpolationTableLoader;

  35. /** The modified Saastamoinen model. Estimates the path delay imposed to
  36.  * electro-magnetic signals by the troposphere according to the formula:
  37.  * <pre>
  38.  * δ = 2.277e-3 / cos z * (P + (1255 / T + 0.05) * e - B * tan²
  39.  * z) + δR
  40.  * </pre>
  41.  * with the following input data provided to the model:
  42.  * <ul>
  43.  * <li>z: zenith angle</li>
  44.  * <li>P: atmospheric pressure</li>
  45.  * <li>T: temperature</li>
  46.  * <li>e: partial pressure of water vapour</li>
  47.  * <li>B, δR: correction terms</li>
  48.  * </ul>
  49.  * <p>
  50.  * The model supports custom δR correction terms to be read from a
  51.  * configuration file (saastamoinen-correction.txt) via the
  52.  * {@link DataProvidersManager}.
  53.  * </p>
  54.  * @author Thomas Neidhart
  55.  * @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
  56.  */
  57. public class SaastamoinenModel implements TroposphericModel {

  59.     /** Default file name for δR correction term table. */
  60.     public static final String DELTA_R_FILE_NAME = "^saastamoinen-correction\\.txt$";

  62.     /** Serializable UID. */
  63.     private static final long serialVersionUID = 20160126L;

  65.     /** X values for the B function. */
  66.     private static final double[] X_VALUES_FOR_B = {
  67.         0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0
  68.     };

  70.     /** E values for the B function. */
  71.     private static final double[] Y_VALUES_FOR_B = {
  72.         1.156, 1.079, 1.006, 0.938, 0.874, 0.813, 0.757, 0.654, 0.563
  73.     };

  75.     /** Coefficients for the partial pressure of water vapor polynomial. */
  76.     private static final double[] E_COEFFICIENTS = {
  77.         -37.2465, 0.213166, -0.000256908
  78.     };

  80.     /** Interpolation function for the B correction term. */
  81.     private final transient UnivariateFunction bFunction;

  83.     /** Polynomial function for the e term. */
  84.     private final transient PolynomialFunction eFunction;

  86.     /** Interpolation function for the delta R correction term. */
  87.     private final transient BilinearInterpolatingFunction deltaRFunction;

  89.     /** The temperature at the station [K]. */
  90.     private double t0;

  92.     /** The atmospheric pressure [mbar]. */
  93.     private double p0;

  95.     /** The humidity [percent]. */
  96.     private double r0;

  98.     /** Create a new Saastamoinen model for the troposphere using the given
  99.      * environmental conditions.
  100.      * @param t0 the temperature at the station [K]
  101.      * @param p0 the atmospheric pressure at the station [mbar]
  102.      * @param r0 the humidity at the station [fraction] (50% -&gt; 0.5)
  103.      * @param deltaRFileName regular expression for filename containing δR
  104.      * correction term table (typically {@link #DELTA_R_FILE_NAME}), if null
  105.      * default values from the reference book are used
  106.      * @exception OrekitException if δR correction term table cannot be loaded
  107.      * @since 7.1
  108.      */
  109.     public SaastamoinenModel(final double t0, final double p0, final double r0,
  110.                              final String deltaRFileName)
  111.         throws OrekitException {
  112.         this(t0, p0, r0,
  113.              deltaRFileName == null ? defaultDeltaR() : loadDeltaR(deltaRFileName));
  114.     }

  116.     /** Create a new Saastamoinen model.
  117.      * @param t0 the temperature at the station [K]
  118.      * @param p0 the atmospheric pressure at the station [mbar]
  119.      * @param r0 the humidity at the station [fraction] (50% -> 0.5)
  120.      * @param deltaR δR correction term function
  121.      * @since 7.1
  122.      */
  123.     private SaastamoinenModel(final double t0, final double p0, final double r0,
  124.                               final BilinearInterpolatingFunction deltaR) {
  125.         this.t0             = t0;
  126.         this.p0             = p0;
  127.         this.r0             = r0;
  128.         this.bFunction      = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
  129.         this.eFunction      = new PolynomialFunction(E_COEFFICIENTS);
  130.         this.deltaRFunction = deltaR;
  131.     }

  133.     /** Create a new Saastamoinen model using a standard atmosphere model.
  134.      *
  135.      * <ul>
  136.      * <li>temperature: 18 degree Celsius
  137.      * <li>pressure: 1013.25 mbar
  138.      * <li>humidity: 50%
  139.      * </ul>
  140.      *
  141.      * @return a Saastamoinen model with standard environmental values
  142.      * @exception OrekitException if δR correction term table cannot be loaded
  143.      */
  144.     public static SaastamoinenModel getStandardModel()
  145.         throws OrekitException {
  146.         return new SaastamoinenModel(273.16 + 18, 1013.25, 0.5, (String) null);
  147.     }

  149.     /** {@inheritDoc}
  150.      * <p>
  151.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  152.      * reasons, we use the value for h = 0 when altitude is negative.
  153.      * </p>
  154.      */
  155.     public double pathDelay(final double elevation, final double height) {

  157.         // there are no data in the model for negative altitudes
  158.         // we use the data for the lowest available altitude: 0.0
  159.         final double fixedHeight = FastMath.max(0.0, height);

  161.         // the corrected temperature using a temperature gradient of -6.5 K/km
  162.         final double T = t0 - 6.5e-3 * fixedHeight;
  163.         // the corrected pressure
  164.         final double P = p0 * FastMath.pow(1.0 - 2.26e-5 * fixedHeight, 5.225);
  165.         // the corrected humidity
  166.         final double R = r0 * FastMath.exp(-6.396e-4 * fixedHeight);

  168.         // interpolate the b correction term
  169.         final double B = bFunction.value(fixedHeight / 1e3);
  170.         // calculate e
  171.         final double e = R * FastMath.exp(eFunction.value(T));

  173.         // calculate the zenith angle from the elevation
  174.         final double z = FastMath.abs(0.5 * FastMath.PI - elevation);

  176.         // get correction factor
  177.         final double deltaR = getDeltaR(fixedHeight, z);

  179.         // calculate the path delay in m
  180.         final double tan = FastMath.tan(z);
  181.         final double delta = 2.277e-3 / FastMath.cos(z) *
  182.                              (P + (1255d / T + 5e-2) * e - B * tan * tan) + deltaR;

  184.         return delta;
  185.     }

  187.     /** Calculates the delta R correction term using linear interpolation.
  188.      * @param height the height of the station in m
  189.      * @param zenith the zenith angle of the satellite
  190.      * @return the delta R correction term in m
  191.      */
  192.     private double getDeltaR(final double height, final double zenith) {
  193.         // limit the height to a range of [0, 5000] m
  194.         final double h = FastMath.min(FastMath.max(0, height), 5000);
  195.         // limit the zenith angle to 90 degree
  196.         // Note: the function is symmetric for negative zenith angles
  197.         final double z = FastMath.min(Math.abs(zenith), 0.5 * FastMath.PI);
  198.         return deltaRFunction.value(h, z);
  199.     }

  201.     /** Load δR function.
  202.      * @param deltaRFileName regular expression for filename containing δR
  203.      * correction term table
  204.      * @return δR function
  205.      * @exception OrekitException if table cannot be loaded
  206.      */
  207.     private static BilinearInterpolatingFunction loadDeltaR(final String deltaRFileName)
  208.         throws OrekitException {

  210.         // read the δR interpolation function from the config file
  211.         final InterpolationTableLoader loader = new InterpolationTableLoader();
  212.         DataProvidersManager.getInstance().feed(deltaRFileName, loader);
  213.         if (!loader.stillAcceptsData()) {
  214.             final double[] elevations = loader.getOrdinateGrid();
  215.             for (int i = 0; i < elevations.length; ++i) {
  216.                 elevations[i] = FastMath.toRadians(elevations[i]);
  217.             }
  218.             return new BilinearInterpolatingFunction(loader.getAbscissaGrid(), elevations,
  219.                                                      loader.getValuesSamples());
  220.         }
  221.         throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE,
  222.                                   deltaRFileName.replaceAll("^\\^", "").replaceAll("\\$$", ""));
  223.     }

  225.     /** Create the default δR function.
  226.      * @return δR function
  227.      */
  228.     private static BilinearInterpolatingFunction defaultDeltaR() {

  230.         // the correction table in the referenced book only contains values for an angle of 60 - 80
  231.         // degree, thus for 0 degree, the correction term is assumed to be 0, for degrees > 80 it
  232.         // is assumed to be the same value as for 80.

  234.         // the height in m
  235.         final double xValForR[] = {
  236.             0, 500, 1000, 1500, 2000, 3000, 4000, 5000
  237.         };

  239.         // the zenith angle
  240.         final double yValForR[] = {
  241.             FastMath.toRadians( 0.00), FastMath.toRadians(60.00), FastMath.toRadians(66.00), FastMath.toRadians(70.00),
  242.             FastMath.toRadians(73.00), FastMath.toRadians(75.00), FastMath.toRadians(76.00), FastMath.toRadians(77.00),
  243.             FastMath.toRadians(78.00), FastMath.toRadians(78.50), FastMath.toRadians(79.00), FastMath.toRadians(79.50),
  244.             FastMath.toRadians(79.75), FastMath.toRadians(80.00), FastMath.toRadians(90.00)
  245.         };

  247.         final double[][] fval = new double[][] {
  248.             {
  249.                 0.000, 0.003, 0.006, 0.012, 0.020, 0.031, 0.039, 0.050, 0.065, 0.075, 0.087, 0.102, 0.111, 0.121, 0.121
  250.             }, {
  251.                 0.000, 0.003, 0.006, 0.011, 0.018, 0.028, 0.035, 0.045, 0.059, 0.068, 0.079, 0.093, 0.101, 0.110, 0.110
  252.             }, {
  253.                 0.000, 0.002, 0.005, 0.010, 0.017, 0.025, 0.032, 0.041, 0.054, 0.062, 0.072, 0.085, 0.092, 0.100, 0.100
  254.             }, {
  255.                 0.000, 0.002, 0.005, 0.009, 0.015, 0.023, 0.029, 0.037, 0.049, 0.056, 0.065, 0.077, 0.083, 0.091, 0.091
  256.             }, {
  257.                 0.000, 0.002, 0.004, 0.008, 0.013, 0.021, 0.026, 0.033, 0.044, 0.051, 0.059, 0.070, 0.076, 0.083, 0.083
  258.             }, {
  259.                 0.000, 0.002, 0.003, 0.006, 0.011, 0.017, 0.021, 0.027, 0.036, 0.042, 0.049, 0.058, 0.063, 0.068, 0.068
  260.             }, {
  261.                 0.000, 0.001, 0.003, 0.005, 0.009, 0.014, 0.017, 0.022, 0.030, 0.034, 0.040, 0.047, 0.052, 0.056, 0.056
  262.             }, {
  263.                 0.000, 0.001, 0.002, 0.004, 0.007, 0.011, 0.014, 0.018, 0.024, 0.028, 0.033, 0.039, 0.043, 0.047, 0.047
  264.             }
  265.         };

  267.         // the actual delta R is interpolated using a a bilinear interpolator
  268.         return new BilinearInterpolatingFunction(xValForR, yValForR, fval);

  270.     }

  272.     /** Replace the instance with a data transfer object for serialization.
  273.      * @return data transfer object that will be serialized
  274.      */
  275.     private Object writeReplace() {
  276.         return new DataTransferObject(this);
  277.     }

  279.     /** Specialization of the data transfer object for serialization. */
  280.     private static class DataTransferObject implements Serializable {

  282.         /** Serializable UID. */
  283.         private static final long serialVersionUID = 20160126L;

  285.         /** The temperature at the station [K]. */
  286.         private double t0;

  288.         /** The atmospheric pressure [mbar]. */
  289.         private double p0;

  291.         /** The humidity [percent]. */
  292.         private double r0;

  294.         /** Samples x-coordinates. */
  295.         private final double[] xval;

  297.         /** Samples y-coordinates. */
  298.         private final double[] yval;

  300.         /** Set of cubic splines patching the whole data grid. */
  301.         private final double[][] fval;

  303.         /** Simple constructor.
  304.          * @param model model to serialize
  305.          */
  306.         DataTransferObject(final SaastamoinenModel model) {
  307.             this.t0 = model.t0;
  308.             this.p0 = model.p0;
  309.             this.r0 = model.r0;
  310.             this.xval = model.deltaRFunction.xval.clone();
  311.             this.yval = model.deltaRFunction.yval.clone();
  312.             this.fval = model.deltaRFunction.fval.clone();
  313.         }

  315.         /** Replace the deserialized data transfer object with a {@link SaastamoinenModel}.
  316.          * @return replacement {@link SaastamoinenModel}
  317.          */
  318.         private Object readResolve() {
  319.             return new SaastamoinenModel(t0, p0, r0,
  320.                                          new BilinearInterpolatingFunction(xval, yval, fval));
  321.         }

  323.     }

  325.     /**
  326.      * Function that implements a standard bilinear interpolation.
  327.      * The interpolation as found
  328.      * in the Wikipedia reference <a href =
  329.      * "http://en.wikipedia.org/wiki/Bilinear_interpolation">BiLinear
  330.      * Interpolation</a>. This is a stand-in until Apache Math has a
  331.      * bilinear interpolator
  332.      */
  333.     private static class BilinearInterpolatingFunction implements BivariateFunction {

  335.         /**
  336.          * The minimum number of points that are needed to compute the
  337.          * function.
  338.          */
  339.         private static final int MIN_NUM_POINTS = 2;

  341.         /** Samples x-coordinates. */
  342.         private final double[] xval;

  344.         /** Samples y-coordinates. */
  345.         private final double[] yval;

  347.         /** Set of cubic splines patching the whole data grid. */
  348.         private final double[][] fval;

  350.         /**
  351.          * @param x Sample values of the x-coordinate, in increasing order.
  352.          * @param y Sample values of the y-coordinate, in increasing order.
  353.          * @param f Values of the function on every grid point. the expected
  354.          *        number of elements.
  355.          * @throws MathIllegalArgumentException if the length of x and y don't
  356.          *         match the row, column height of f, or if any of the arguments
  357.          *         are null, or if any of the arrays has zero length, or if
  358.          *         {@code x} or {@code y} are not strictly increasing.
  359.          */
  360.         BilinearInterpolatingFunction(final double[] x, final double[] y, final double[][] f)
  361.                         throws MathIllegalArgumentException {

  363.             if (x == null || y == null || f == null || f[0] == null) {
  364.                 throw new IllegalArgumentException("All arguments must be non-null");
  365.             }

  367.             final int xLen = x.length;
  368.             final int yLen = y.length;

  370.             if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
  371.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
  372.             }

  374.             if (xLen < MIN_NUM_POINTS || yLen < MIN_NUM_POINTS || f.length < MIN_NUM_POINTS ||
  375.                             f[0].length < MIN_NUM_POINTS) {
  376.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DATA);
  377.             }

  379.             if (xLen != f.length) {
  380.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  381.                                                        xLen, f.length);
  382.             }

  384.             if (yLen != f[0].length) {
  385.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  386.                                                        yLen, f[0].length);
  387.             }

  389.             MathArrays.checkOrder(x);
  390.             MathArrays.checkOrder(y);

  392.             xval = x.clone();
  393.             yval = y.clone();
  394.             fval = f.clone();
  395.         }

  397.         @Override
  398.         public double value(final double x, final double y) {
  399.             final int offset = 1;
  400.             final int count = offset + 1;
  401.             final int i = searchIndex(x, xval, offset, count);
  402.             final int j = searchIndex(y, yval, offset, count);

  404.             final double x1 = xval[i];
  405.             final double x2 = xval[i + 1];
  406.             final double y1 = yval[j];
  407.             final double y2 = yval[j + 1];
  408.             final double fQ11 = fval[i][j];
  409.             final double fQ21 = fval[i + 1][j];
  410.             final double fQ12 = fval[i][j + 1];
  411.             final double fQ22 = fval[i + 1][j + 1];

  413.             final double f = (fQ11 * (x2 - x)  * (y2 - y) +
  414.                             fQ21 * (x  - x1) * (y2 - y) +
  415.                             fQ12 * (x2 - x)  * (y  - y1) +
  416.                             fQ22 * (x  - x1) * (y  - y1)) /
  417.                             ((x2 - x1) * (y2 - y1));

  419.             return f;
  420.         }

  422.         /**
  423.          * @param c Coordinate.
  424.          * @param val Coordinate samples.
  425.          * @param offset how far back from found value to offset for
  426.          *        querying
  427.          * @param count total number of elements forward from beginning that
  428.          *        will be queried
  429.          * @return the index in {@code val} corresponding to the interval
  430.          *         containing {@code c}.
  431.          * @throws MathIllegalArgumentException if {@code c} is out of the range
  432.          *         defined by the boundary values of {@code val}.
  433.          */
  434.         private int searchIndex(final double c, final double[] val, final int offset, final int count)
  435.             throws MathIllegalArgumentException {
  436.             int r = Arrays.binarySearch(val, c);

  438.             if (r == -1 || r == -val.length - 1) {
  439.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  440.                                                        c, val[0], val[val.length - 1]);
  441.             }

  443.             if (r < 0) {
  444.                 // "c" in within an interpolation sub-interval, which
  445.                 // returns
  446.                 // negative
  447.                 // need to remove the negative sign for consistency
  448.                 r = -r - offset - 1;
  449.             } else {
  450.                 r -= offset;
  451.             }

  453.             if (r < 0) {
  454.                 r = 0;
  455.             }

  457.             if ((r + count) >= val.length) {
  458.                 // "c" is the last sample of the range: Return the index
  459.                 // of the sample at the lower end of the last sub-interval.
  460.                 r = val.length - count;
  461.             }

  463.             return r;
  464.         }

  466.     }

  468. }
Created with JaCoCo 0.7.9.201702052155