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.           * @since 7.1
  107.      */
  108.     public SaastamoinenModel(final double t0, final double p0, final double r0,
  109.                              final String deltaRFileName) {
  110.         this(t0, p0, r0,
  111.              deltaRFileName == null ? defaultDeltaR() : loadDeltaR(deltaRFileName));
  112.     }

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

  131.     /** Create a new Saastamoinen model using a standard atmosphere model.
  132.      *
  133.      * <ul>
  134.      * <li>temperature: 18 degree Celsius
  135.      * <li>pressure: 1013.25 mbar
  136.      * <li>humidity: 50%
  137.      * </ul>
  138.      *
  139.      * @return a Saastamoinen model with standard environmental values
  140.      */
  141.     public static SaastamoinenModel getStandardModel() {
  142.         return new SaastamoinenModel(273.16 + 18, 1013.25, 0.5, (String) null);
  143.     }

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

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

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

  164.         // interpolate the b correction term
  165.         final double B = bFunction.value(fixedHeight / 1e3);
  166.         // calculate e
  167.         final double e = R * FastMath.exp(eFunction.value(T));

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

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

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

  180.         return delta;
  181.     }

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

  197.     /** Load δR function.
  198.      * @param deltaRFileName regular expression for filename containing δR
  199.      * correction term table
  200.      * @return δR function
  201.      */
  202.     private static BilinearInterpolatingFunction loadDeltaR(final String deltaRFileName) {

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

  219.     /** Create the default δR function.
  220.      * @return δR function
  221.      */
  222.     private static BilinearInterpolatingFunction defaultDeltaR() {

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

  228.         // the height in m
  229.         final double xValForR[] = {
  230.             0, 500, 1000, 1500, 2000, 3000, 4000, 5000
  231.         };

  233.         // the zenith angle
  234.         final double yValForR[] = {
  235.             FastMath.toRadians( 0.00), FastMath.toRadians(60.00), FastMath.toRadians(66.00), FastMath.toRadians(70.00),
  236.             FastMath.toRadians(73.00), FastMath.toRadians(75.00), FastMath.toRadians(76.00), FastMath.toRadians(77.00),
  237.             FastMath.toRadians(78.00), FastMath.toRadians(78.50), FastMath.toRadians(79.00), FastMath.toRadians(79.50),
  238.             FastMath.toRadians(79.75), FastMath.toRadians(80.00), FastMath.toRadians(90.00)
  239.         };

  241.         final double[][] fval = new double[][] {
  242.             {
  243.                 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
  244.             }, {
  245.                 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
  246.             }, {
  247.                 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
  248.             }, {
  249.                 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
  250.             }, {
  251.                 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
  252.             }, {
  253.                 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
  254.             }, {
  255.                 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
  256.             }, {
  257.                 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
  258.             }
  259.         };

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

  264.     }

  266.     /** Replace the instance with a data transfer object for serialization.
  267.      * @return data transfer object that will be serialized
  268.      */
  269.     private Object writeReplace() {
  270.         return new DataTransferObject(this);
  271.     }

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

  276.         /** Serializable UID. */
  277.         private static final long serialVersionUID = 20160126L;

  279.         /** The temperature at the station [K]. */
  280.         private double t0;

  282.         /** The atmospheric pressure [mbar]. */
  283.         private double p0;

  285.         /** The humidity [percent]. */
  286.         private double r0;

  288.         /** Samples x-coordinates. */
  289.         private final double[] xval;

  291.         /** Samples y-coordinates. */
  292.         private final double[] yval;

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

  297.         /** Simple constructor.
  298.          * @param model model to serialize
  299.          */
  300.         DataTransferObject(final SaastamoinenModel model) {
  301.             this.t0 = model.t0;
  302.             this.p0 = model.p0;
  303.             this.r0 = model.r0;
  304.             this.xval = model.deltaRFunction.xval.clone();
  305.             this.yval = model.deltaRFunction.yval.clone();
  306.             this.fval = model.deltaRFunction.fval.clone();
  307.         }

  309.         /** Replace the deserialized data transfer object with a {@link SaastamoinenModel}.
  310.          * @return replacement {@link SaastamoinenModel}
  311.          */
  312.         private Object readResolve() {
  313.             return new SaastamoinenModel(t0, p0, r0,
  314.                                          new BilinearInterpolatingFunction(xval, yval, fval));
  315.         }

  317.     }

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

  329.         /**
  330.          * The minimum number of points that are needed to compute the
  331.          * function.
  332.          */
  333.         private static final int MIN_NUM_POINTS = 2;

  335.         /** Samples x-coordinates. */
  336.         private final double[] xval;

  338.         /** Samples y-coordinates. */
  339.         private final double[] yval;

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

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

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

  361.             final int xLen = x.length;
  362.             final int yLen = y.length;

  364.             if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
  365.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
  366.             }

  368.             if (xLen < MIN_NUM_POINTS || yLen < MIN_NUM_POINTS || f.length < MIN_NUM_POINTS ||
  369.                             f[0].length < MIN_NUM_POINTS) {
  370.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DATA);
  371.             }

  373.             if (xLen != f.length) {
  374.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  375.                                                        xLen, f.length);
  376.             }

  378.             if (yLen != f[0].length) {
  379.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  380.                                                        yLen, f[0].length);
  381.             }

  383.             MathArrays.checkOrder(x);
  384.             MathArrays.checkOrder(y);

  386.             xval = x.clone();
  387.             yval = y.clone();
  388.             fval = f.clone();
  389.         }

  391.         @Override
  392.         public double value(final double x, final double y) {
  393.             final int offset = 1;
  394.             final int count = offset + 1;
  395.             final int i = searchIndex(x, xval, offset, count);
  396.             final int j = searchIndex(y, yval, offset, count);

  398.             final double x1 = xval[i];
  399.             final double x2 = xval[i + 1];
  400.             final double y1 = yval[j];
  401.             final double y2 = yval[j + 1];
  402.             final double fQ11 = fval[i][j];
  403.             final double fQ21 = fval[i + 1][j];
  404.             final double fQ12 = fval[i][j + 1];
  405.             final double fQ22 = fval[i + 1][j + 1];

  407.             final double f = (fQ11 * (x2 - x)  * (y2 - y) +
  408.                             fQ21 * (x  - x1) * (y2 - y) +
  409.                             fQ12 * (x2 - x)  * (y  - y1) +
  410.                             fQ22 * (x  - x1) * (y  - y1)) /
  411.                             ((x2 - x1) * (y2 - y1));

  413.             return f;
  414.         }

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

  432.             if (r == -1 || r == -val.length - 1) {
  433.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
  434.                                                        c, val[0], val[val.length - 1]);
  435.             }

  437.             if (r < 0) {
  438.                 // "c" in within an interpolation sub-interval, which
  439.                 // returns
  440.                 // negative
  441.                 // need to remove the negative sign for consistency
  442.                 r = -r - offset - 1;
  443.             } else {
  444.                 r -= offset;
  445.             }

  447.             if (r < 0) {
  448.                 r = 0;
  449.             }

  451.             if ((r + count) >= val.length) {
  452.                 // "c" is the last sample of the range: Return the index
  453.                 // of the sample at the lower end of the last sub-interval.
  454.                 r = val.length - count;
  455.             }

  457.             return r;
  458.         }

  460.     }

  462. }
Created with JaCoCo 0.8.2.201808211720