ITURP834MappingFunction.java

  1. /* Copyright 2022-2025 Thales Alenia Space
  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.troposphere.iturp834;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
  21. import org.hipparchus.util.FastMath;
  22. import org.hipparchus.util.FieldSinCos;
  23. import org.hipparchus.util.MathArrays;
  24. import org.hipparchus.util.MathUtils;
  25. import org.hipparchus.util.SinCos;
  26. import org.orekit.bodies.FieldGeodeticPoint;
  27. import org.orekit.bodies.GeodeticPoint;
  28. import org.orekit.errors.OrekitException;
  29. import org.orekit.errors.OrekitMessages;
  30. import org.orekit.models.earth.troposphere.TroposphereMappingFunction;
  31. import org.orekit.time.AbsoluteDate;
  32. import org.orekit.time.FieldAbsoluteDate;
  33. import org.orekit.time.TimeScale;
  34. import org.orekit.utils.FieldTrackingCoordinates;
  35. import org.orekit.utils.TrackingCoordinates;

  37. import java.io.BufferedReader;
  38. import java.io.IOException;
  39. import java.io.InputStream;
  40. import java.io.InputStreamReader;
  41. import java.nio.charset.StandardCharsets;
  42. import java.util.regex.Pattern;

  44. /** ITU-R P.834 mapping function.
  45.  * @see ITURP834PathDelay
  46.  * @see ITURP834WeatherParametersProvider
  47.  * @author Luc Maisonobe
  48.  * @see <a href="https://www.itu.int/rec/R-REC-P.834/en">P.834 : Effects of tropospheric refraction on radiowave propagation</a>
  49.  * @since 13.0
  50.  */
  51. public class ITURP834MappingFunction implements TroposphereMappingFunction {

  53.     /** Splitter for fields in lines. */
  54.     private static final Pattern SPLITTER = Pattern.compile("\\s+");

  56.     /** Name of data file. */
  57.     private static final String MAPPING_FUNCTION_NAME = "/assets/org/orekit/ITU-R-P.834/p834_mf_coeff_v1.txt";

  59.     /** Minimum latitude, including extra margin for dealing with boundaries (degrees). */
  60.     private static final double MIN_LAT = -92.5;

  62.     /** Maximum latitude, including extra margin for dealing with boundaries (degrees). */
  63.     private static final double MAX_LAT =  92.5;

  65.     /** Grid step in latitude (degrees). */
  66.     private static final double STEP_LAT =  5.0;

  68.     /** Minimum longitude, including extra margin for dealing with boundaries (degrees). */
  69.     private static final double MIN_LON = -182.5;

  71.     /** Maximum longitude, including extra margin for dealing with boundaries (degrees). */
  72.     private static final double MAX_LON =  182.5;

  74.     /** Grid step in longitude (degrees). */
  75.     private static final double STEP_LON =  5.0;

  77.     /** Second coefficient for hydrostatic component. */
  78.     private static final double BH = 0.0029;

  80.     /** Constants for third coefficient for hydrostatic component, Northern hemisphere. */
  81.     private static final double[] CH_NORTH = { 0.062, 0.001, 0.005, 0.0 };

  83.     /** Constants for third coefficient for hydrostatic component, Southern hemisphere. */
  84.     private static final double[] CH_SOUTH = { 0.062, 0.002, 0.007, FastMath.PI };

  86.     /** Reference day of year for third coefficient for hydrostatic component. */
  87.     private static final double REF_DOY = 28;

  89.     /** Year length (in days). */
  90.     private static final double YEAR = 365.25;

  92.     /** Second coefficient for wet component. */
  93.     private static final double BW = 0.00146;

  95.     /** Third coefficient for wet component. */
  96.     private static final double CW = 0.04391;

  98.     /** Global factor to apply to first coefficient. */
  99.     private static final double FACTOR = 1.0e-3;

  101.     /** Interpolator for constant hydrostatic coefficient. */
  102.     private static final BilinearInterpolatingFunction A0H;

  104.     /** Interpolator for annual cosine hydrostatic coefficient. */
  105.     private static final BilinearInterpolatingFunction A1H;

  107.     /** Interpolator for annual sine hydrostatic coefficient. */
  108.     private static final BilinearInterpolatingFunction B1H;

  110.     /** Interpolator for semi-annual cosine hydrostatic coefficient. */
  111.     private static final BilinearInterpolatingFunction A2H;

  113.     /** Interpolator for semin-annual sine hydrostatic coefficient. */
  114.     private static final BilinearInterpolatingFunction B2H;

  116.     /** Interpolator for constant wet coefficient. */
  117.     private static final BilinearInterpolatingFunction A0W;

  119.     /** Interpolator for annual cosine wet coefficient. */
  120.     private static final BilinearInterpolatingFunction A1W;

  122.     /** Interpolator for annual sine wet coefficient. */
  123.     private static final BilinearInterpolatingFunction B1W;

  125.     /** Interpolator for semi-annual cosine wet coefficient. */
  126.     private static final BilinearInterpolatingFunction A2W;

  128.     /** Interpolator for semin-annual sine wet coefficient. */
  129.     private static final BilinearInterpolatingFunction B2W;

  131.     // load model data
  132.     static {

  134.         // create the various tables
  135.         // we add extra lines and columns to the official files, for dealing with boundaries
  136.         final double[]   longitudes = interpolationPoints(MIN_LON, MAX_LON, STEP_LON);
  137.         final double[]   latitudes  = interpolationPoints(MIN_LAT, MAX_LAT, STEP_LAT);
  138.         final double[][] a0h        = new double[longitudes.length][latitudes.length];
  139.         final double[][] a1h        = new double[longitudes.length][latitudes.length];
  140.         final double[][] b1h        = new double[longitudes.length][latitudes.length];
  141.         final double[][] a2h        = new double[longitudes.length][latitudes.length];
  142.         final double[][] b2h        = new double[longitudes.length][latitudes.length];
  143.         final double[][] a0w        = new double[longitudes.length][latitudes.length];
  144.         final double[][] a1w        = new double[longitudes.length][latitudes.length];
  145.         final double[][] b1w        = new double[longitudes.length][latitudes.length];
  146.         final double[][] a2w        = new double[longitudes.length][latitudes.length];
  147.         final double[][] b2w        = new double[longitudes.length][latitudes.length];

  149.         try (InputStream       is     = ITURP834MappingFunction.class.getResourceAsStream(MAPPING_FUNCTION_NAME);
  150.              InputStreamReader isr    = is  == null ? null : new InputStreamReader(is, StandardCharsets.UTF_8);
  151.              BufferedReader    reader = isr == null ? null : new BufferedReader(isr)) {
  152.             if (reader == null) {
  153.                 // this should never happen with embedded data
  154.                 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, MAPPING_FUNCTION_NAME);
  155.             }
  156.             int lineNumber = 0;
  157.             for (String line = reader.readLine(); line != null; line = reader.readLine()) {
  158.                 ++lineNumber;
  159.                 final String[] fields = SPLITTER.split(line.trim());
  160.                 if (fields.length != 12) {
  161.                     // this should never happen with the embedded data
  162.                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  163.                                               lineNumber, MAPPING_FUNCTION_NAME, line);
  164.                 }

  166.                 // parse the fields
  167.                 final double[] numericFields = new double[fields.length];
  168.                 for (int i = 0; i < fields.length; ++i) {
  169.                     numericFields[i] = Double.parseDouble(fields[i]);
  170.                 }

  172.                 // find indices in our extended grid
  173.                 final int longitudeIndex = (int) FastMath.rint((numericFields[1] - MIN_LON) / STEP_LON);
  174.                 final int latitudeIndex  = (int) FastMath.rint((numericFields[0] - MIN_LAT) / STEP_LAT);

  176.                 // fill-in tables
  177.                 a0h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 2];
  178.                 a1h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 3];
  179.                 b1h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 4];
  180.                 a2h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 5];
  181.                 b2h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 6];
  182.                 a0w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 7];
  183.                 a1w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 8];
  184.                 b1w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 9];
  185.                 a2w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[10];
  186.                 b2w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[11];

  188.             }

  190.             // extend tables in latitude to cover poles
  191.             for (int i = 1; i < longitudes.length - 1; ++i) {
  192.                 a0h[i][0]                    = a0h[i][1];
  193.                 a0h[i][latitudes.length - 1] = a0h[i][latitudes.length - 2];
  194.                 a1h[i][0]                    = a1h[i][1];
  195.                 a1h[i][latitudes.length - 1] = a1h[i][latitudes.length - 2];
  196.                 b1h[i][0]                    = b1h[i][1];
  197.                 b1h[i][latitudes.length - 1] = b1h[i][latitudes.length - 2];
  198.                 a2h[i][0]                    = a2h[i][1];
  199.                 a2h[i][latitudes.length - 1] = a2h[i][latitudes.length - 2];
  200.                 b2h[i][0]                    = b2h[i][1];
  201.                 b2h[i][latitudes.length - 1] = b2h[i][latitudes.length - 2];
  202.                 a0w[i][0]                    = a0w[i][1];
  203.                 a0w[i][latitudes.length - 1] = a0w[i][latitudes.length - 2];
  204.                 a1w[i][0]                    = a1w[i][1];
  205.                 a1w[i][latitudes.length - 1] = a1w[i][latitudes.length - 2];
  206.                 b1w[i][0]                    = b1w[i][1];
  207.                 b1w[i][latitudes.length - 1] = b1w[i][latitudes.length - 2];
  208.                 a2w[i][0]                    = a2w[i][1];
  209.                 a2w[i][latitudes.length - 1] = a2w[i][latitudes.length - 2];
  210.                 b2w[i][0]                    = b2w[i][1];
  211.                 b2w[i][latitudes.length - 1] = b2w[i][latitudes.length - 2];
  212.             }

  214.             // extend tables in longitude to cover anti-meridian
  215.             for (int j = 0; j < latitudes.length; ++j) {
  216.                 a0h[0][j]                     = a0h[longitudes.length - 2][j];
  217.                 a0h[longitudes.length - 1][j] = a0h[1][j];
  218.                 a1h[0][j]                     = a1h[longitudes.length - 2][j];
  219.                 a1h[longitudes.length - 1][j] = a1h[1][j];
  220.                 b1h[0][j]                     = b1h[longitudes.length - 2][j];
  221.                 b1h[longitudes.length - 1][j] = b1h[1][j];
  222.                 a2h[0][j]                     = a2h[longitudes.length - 2][j];
  223.                 a2h[longitudes.length - 1][j] = a2h[1][j];
  224.                 b2h[0][j]                     = b2h[longitudes.length - 2][j];
  225.                 b2h[longitudes.length - 1][j] = b2h[1][j];
  226.                 a0w[0][j]                     = a0w[longitudes.length - 2][j];
  227.                 a0w[longitudes.length - 1][j] = a0w[1][j];
  228.                 a1w[0][j]                     = a1w[longitudes.length - 2][j];
  229.                 a1w[longitudes.length - 1][j] = a1w[1][j];
  230.                 b1w[0][j]                     = b1w[longitudes.length - 2][j];
  231.                 b1w[longitudes.length - 1][j] = b1w[1][j];
  232.                 a2w[0][j]                     = a2w[longitudes.length - 2][j];
  233.                 a2w[longitudes.length - 1][j] = a2w[1][j];
  234.                 b2w[0][j]                     = b2w[longitudes.length - 2][j];
  235.                 b2w[longitudes.length - 1][j] = b2w[1][j];
  236.             }

  238.             // build interpolators
  239.             A0H = new BilinearInterpolatingFunction(longitudes, latitudes, a0h);
  240.             A1H = new BilinearInterpolatingFunction(longitudes, latitudes, a1h);
  241.             B1H = new BilinearInterpolatingFunction(longitudes, latitudes, b1h);
  242.             A2H = new BilinearInterpolatingFunction(longitudes, latitudes, a2h);
  243.             B2H = new BilinearInterpolatingFunction(longitudes, latitudes, b2h);
  244.             A0W = new BilinearInterpolatingFunction(longitudes, latitudes, a0w);
  245.             A1W = new BilinearInterpolatingFunction(longitudes, latitudes, a1w);
  246.             B1W = new BilinearInterpolatingFunction(longitudes, latitudes, b1w);
  247.             A2W = new BilinearInterpolatingFunction(longitudes, latitudes, a2w);
  248.             B2W = new BilinearInterpolatingFunction(longitudes, latitudes, b2w);

  250.         } catch (IOException ioe) {
  251.             // this should never happen with the embedded data
  252.             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, ioe);
  253.         }
  254.     }

  256.     /** UTC time scale. */
  257.     private final TimeScale utc;

  259.     /** Simple constructor.
  260.      * @param utc UTC time scale
  261.      */
  262.     public ITURP834MappingFunction(final TimeScale utc) {
  263.         this.utc = utc;
  264.     }

  266.     /** Create interpolation points coordinates.
  267.      * @param min min angle in degrees (included)
  268.      * @param max max angle in degrees (included)
  269.      * @param step step between points
  270.      * @return interpolation points coordinates (in radians)
  271.      */
  272.     private static double[] interpolationPoints(final double min, final double max, final double step) {
  273.         final double[] points = new double[(int) FastMath.rint((max - min) / step) + 1];
  274.         for (int i = 0; i < points.length; i++) {
  275.             points[i] = FastMath.toRadians(min + i * step);
  276.         }
  277.         return points;
  278.     }

  280.     @Override
  281.     public double[] mappingFactors(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
  282.                                    final AbsoluteDate date) {

  284.         final double doy = date.getDayOfYear(utc);

  286.         // compute third dry coefficient
  287.         // equation 28c in ITU-R P.834 recommendation
  288.         final double[] c    = point.getLatitude() >= 0.0 ? CH_NORTH : CH_SOUTH;
  289.         final double   cosL = FastMath.cos(point.getLatitude());
  290.         final double   cosD = FastMath.cos((doy - REF_DOY) * MathUtils.TWO_PI / YEAR + c[3]);
  291.         final double   ch   = c[0] + ((cosD + 1) * c[2] * 0.5 + c[1]) * (1 - cosL);

  293.         // compute first coefficient
  294.         final SinCos sc1 = FastMath.sinCos(MathUtils.TWO_PI * doy / YEAR);
  295.         final SinCos sc2 = SinCos.sum(sc1, sc1);

  297.         // equation 28d in ITU-R P.834 recommendation
  298.         final double ah  = A0H.value(point.getLongitude(), point.getLatitude()) +
  299.                            A1H.value(point.getLongitude(), point.getLatitude()) * sc1.cos() +
  300.                            B1H.value(point.getLongitude(), point.getLatitude()) * sc1.sin() +
  301.                            A2H.value(point.getLongitude(), point.getLatitude()) * sc2.cos() +
  302.                            B2H.value(point.getLongitude(), point.getLatitude()) * sc2.sin();

  304.         // equation 28e in ITU-R P.834 recommendation
  305.         final double aw  = A0W.value(point.getLongitude(), point.getLatitude()) +
  306.                            A1W.value(point.getLongitude(), point.getLatitude()) * sc1.cos() +
  307.                            B1W.value(point.getLongitude(), point.getLatitude()) * sc1.sin() +
  308.                            A2W.value(point.getLongitude(), point.getLatitude()) * sc2.cos() +
  309.                            B2W.value(point.getLongitude(), point.getLatitude()) * sc2.sin();

  311.         // mapping function, equations 28a and 28b in ITU-R P.834 recommendation
  312.         final double sinTheta = FastMath.sin(trackingCoordinates.getElevation());
  313.         final double mh = (1 + ah / (1 + BH / (1 + ch))) /
  314.                           (sinTheta + ah / (sinTheta + BH / (sinTheta + ch)));
  315.         final double mw = (1 + aw / (1 + BW / (1 + CW))) /
  316.                           (sinTheta + aw / (sinTheta + BW / (sinTheta + CW)));

  318.         return new double[] {
  319.             mh, mw
  320.         };

  322.     }

  324.     @Override
  325.     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
  326.                                                                   final FieldGeodeticPoint<T> point,
  327.                                                                   final FieldAbsoluteDate<T> date) {

  329.         final T doy = date.getDayOfYear(utc);

  331.         // compute third dry coefficient
  332.         // equation 28c in ITU-R P.834 recommendation
  333.         final double[] c    = point.getLatitude().getReal() >= 0.0 ? CH_NORTH : CH_SOUTH;
  334.         final T        cosL = FastMath.cos(point.getLatitude());
  335.         final T        cosD = FastMath.cos(doy.subtract(REF_DOY).multiply(MathUtils.TWO_PI / YEAR).add(c[3]));
  336.         final T        ch   = cosD.add(1).multiply(c[2] * 0.5).add(c[1]).multiply(cosL.subtract(1).negate()).add(c[0]);

  338.         // compute first coefficient
  339.         final FieldSinCos<T> sc1 = FastMath.sinCos(doy.multiply(MathUtils.TWO_PI / YEAR));
  340.         final FieldSinCos<T> sc2 = FieldSinCos.sum(sc1, sc1);

  342.         // equation 28d in ITU-R P.834 recommendation
  343.         final T ah  =     A0H.value(point.getLongitude(), point.getLatitude()).
  344.                       add(A1H.value(point.getLongitude(), point.getLatitude()).multiply(sc1.cos())).
  345.                       add(B1H.value(point.getLongitude(), point.getLatitude()).multiply(sc1.sin())).
  346.                       add(A2H.value(point.getLongitude(), point.getLatitude()).multiply(sc2.cos())).
  347.                       add(B2H.value(point.getLongitude(), point.getLatitude()).multiply(sc2.sin()));

  349.         // equation 28e in ITU-R P.834 recommendation
  350.         final T aw  =     A0W.value(point.getLongitude(), point.getLatitude()).
  351.                       add(A1W.value(point.getLongitude(), point.getLatitude()).multiply(sc1.cos())).
  352.                       add(B1W.value(point.getLongitude(), point.getLatitude()).multiply(sc1.sin())).
  353.                       add(A2W.value(point.getLongitude(), point.getLatitude()).multiply(sc2.cos())).
  354.                       add(B2W.value(point.getLongitude(), point.getLatitude()).multiply(sc2.sin()));

  356.         // mapping function, equations 28a and 28b in ITU-R P.834 recommendation
  357.         final T sinTheta = FastMath.sin(trackingCoordinates.getElevation());
  358.         final T mh = ch.add(1).reciprocal().multiply(BH).add(1).reciprocal().multiply(ah).add(1).
  359.                      divide(ch.add(sinTheta).reciprocal().multiply(BH).add(sinTheta).reciprocal().multiply(ah).add(sinTheta));
  360.         final T mw = aw.divide(1 + BW / (1 + CW)).add(1).
  361.                      divide(sinTheta.add(CW).reciprocal().multiply(BW).add(sinTheta).reciprocal().multiply(aw).add(sinTheta));

  363.         final T[] m = MathArrays.buildArray(date.getField(), 2);
  364.         m[0] = mh;
  365.         m[1] = mw;
  366.         return m;

  368.     }

  370. }
Created with JaCoCo 0.8.12.202403310830