NiellMappingFunctionModel.java

  1. /* Copyright 2002-2025 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.troposphere;

  19. import org.hipparchus.Field;
  20. import org.hipparchus.CalculusFieldElement;
  21. import org.hipparchus.analysis.UnivariateFunction;
  22. import org.hipparchus.analysis.interpolation.LinearInterpolator;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.MathArrays;
  25. import org.orekit.annotation.DefaultDataContext;
  26. import org.orekit.bodies.FieldGeodeticPoint;
  27. import org.orekit.bodies.GeodeticPoint;
  28. import org.orekit.data.DataContext;
  29. import org.orekit.time.AbsoluteDate;
  30. import org.orekit.time.FieldAbsoluteDate;
  31. import org.orekit.time.TimeScale;
  32. import org.orekit.utils.FieldTrackingCoordinates;
  33. import org.orekit.utils.TrackingCoordinates;

  35. /** The Niell Mapping Function  model for radio wavelengths.
  36.  *  This model is an empirical mapping function. It only needs the
  37.  *  values of the station latitude, height and the date for the computations.
  38.  *  <p>
  39.  *  With this model, the hydrostatic mapping function is time and latitude dependent
  40.  *  whereas the wet mapping function is only latitude dependent.
  41.  *  </p>
  42.  *
  43.  * @see "A. E. Niell(1996), Global mapping functions for the atmosphere delay of radio wavelengths,
  44.  *      J. Geophys. Res., 101(B2), pp.  3227–3246, doi:  10.1029/95JB03048."
  45.  *
  46.  * @author Bryan Cazabonne
  47.  *
  48.  */
  49. public class NiellMappingFunctionModel implements TroposphereMappingFunction {

  51.     /** Values for the ah average function. */
  52.     private static final double[] VALUES_FOR_AH_AVERAGE = {
  53.         1.2769934e-3, 1.2683230e-3, 1.2465397e-3, 1.2196049e-3, 1.2045996e-3
  54.     };

  56.     /** Values for the bh average function. */
  57.     private static final double[] VALUES_FOR_BH_AVERAGE = {
  58.         2.9153695e-3, 2.9152299e-3, 2.9288445e-3, 2.9022565e-3, 2.9024912e-3
  59.     };

  61.     /** Values for the ch average function. */
  62.     private static final double[] VALUES_FOR_CH_AVERAGE = {
  63.         62.610505e-3, 62.837393e-3, 63.721774e-3, 63.824265e-3, 64.258455e-3
  64.     };

  66.     /** Values for the ah amplitude function. */
  67.     private static final double[] VALUES_FOR_AH_AMPLITUDE = {
  68.         0.0, 1.2709626e-5, 2.6523662e-5, 3.4000452e-5, 4.1202191e-5
  69.     };

  71.     /** Values for the bh amplitude function. */
  72.     private static final double[] VALUES_FOR_BH_AMPLITUDE = {
  73.         0.0, 2.1414979e-5, 3.0160779e-5, 7.2562722e-5, 11.723375e-5
  74.     };

  76.     /** X values for the ch amplitude function. */
  77.     private static final double[] VALUES_FOR_CH_AMPLITUDE = {
  78.         0.0, 9.0128400e-5, 4.3497037e-5, 84.795348e-5, 170.37206e-5
  79.     };

  81.     /** Values for the aw function. */
  82.     private static final double[] VALUES_FOR_AW = {
  83.         5.8021897e-4, 5.6794847e-4, 5.8118019e-4, 5.9727542e-4, 6.1641693e-4
  84.     };

  86.     /** Values for the bw function. */
  87.     private static final double[] VALUES_FOR_BW = {
  88.         1.4275268e-3, 1.5138625e-3, 1.4572752e-3, 1.5007428e-3, 1.7599082e-3
  89.     };

  91.     /** Values for the cw function. */
  92.     private static final double[] VALUES_FOR_CW = {
  93.         4.3472961e-2, 4.6729510e-2, 4.3908931e-2, 4.4626982e-2, 5.4736038e-2
  94.     };

  96.     /** Values for the cw function. */
  97.     private static final double[] LATITUDE_VALUES = {
  98.         FastMath.toRadians(15.0), FastMath.toRadians(30.0), FastMath.toRadians(45.0), FastMath.toRadians(60.0), FastMath.toRadians(75.0),
  99.     };

  101.     /** Interpolation function for the ah (average) term. */
  102.     private final UnivariateFunction ahAverageFunction;

  104.     /** Interpolation function for the bh (average) term. */
  105.     private final UnivariateFunction bhAverageFunction;

  107.     /** Interpolation function for the ch (average) term. */
  108.     private final UnivariateFunction chAverageFunction;

  110.     /** Interpolation function for the ah (amplitude) term. */
  111.     private final UnivariateFunction ahAmplitudeFunction;

  113.     /** Interpolation function for the bh (amplitude) term. */
  114.     private final UnivariateFunction bhAmplitudeFunction;

  116.     /** Interpolation function for the ch (amplitude) term. */
  117.     private final UnivariateFunction chAmplitudeFunction;

  119.     /** Interpolation function for the aw term. */
  120.     private final UnivariateFunction awFunction;

  122.     /** Interpolation function for the bw term. */
  123.     private final UnivariateFunction bwFunction;

  125.     /** Interpolation function for the cw term. */
  126.     private final UnivariateFunction cwFunction;

  128.     /** UTC time scale. */
  129.     private final TimeScale utc;

  131.     /** Builds a new instance.
  132.      *
  133.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  134.      *
  135.      * @see #NiellMappingFunctionModel(TimeScale)
  136.      */
  137.     @DefaultDataContext
  138.     public NiellMappingFunctionModel() {
  139.         this(DataContext.getDefault().getTimeScales().getUTC());
  140.     }

  142.     /** Builds a new instance.
  143.      * @param utc UTC time scale.
  144.      * @since 10.1
  145.      */
  146.     public NiellMappingFunctionModel(final TimeScale utc) {
  147.         this.utc = utc;
  148.         // Interpolation functions for hydrostatic coefficients
  149.         this.ahAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AH_AVERAGE);
  150.         this.bhAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BH_AVERAGE);
  151.         this.chAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CH_AVERAGE);
  152.         this.ahAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AH_AMPLITUDE);
  153.         this.bhAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BH_AMPLITUDE);
  154.         this.chAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CH_AMPLITUDE);

  156.         // Interpolation functions for wet coefficients
  157.         this.awFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AW);
  158.         this.bwFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BW);
  159.         this.cwFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CW);
  160.     }

  162.     /** {@inheritDoc} */
  163.     @Override
  164.     public double[] mappingFactors(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
  165.                                    final AbsoluteDate date) {

  167.         // Temporal factor
  168.         double t0 = 28;
  169.         if (point.getLatitude() < 0) {
  170.             // southern hemisphere: t0 = 28 + an integer half of year
  171.             t0 += 183;
  172.         }
  173.         final double coef    = 2 * FastMath.PI * ((date.getDayOfYear(utc) - t0) / 365.25);
  174.         final double cosCoef = FastMath.cos(coef);

  176.         // Compute ah, bh and ch Eq. 5
  177.         double absLatidude = FastMath.abs(point.getLatitude());
  178.         // there are no data in the model for latitudes lower than 15°
  179.         absLatidude = FastMath.max(FastMath.toRadians(15.0), absLatidude);
  180.         // there are no data in the model for latitudes greater than 75°
  181.         absLatidude = FastMath.min(FastMath.toRadians(75.0), absLatidude);
  182.         final double ah = ahAverageFunction.value(absLatidude) - ahAmplitudeFunction.value(absLatidude) * cosCoef;
  183.         final double bh = bhAverageFunction.value(absLatidude) - bhAmplitudeFunction.value(absLatidude) * cosCoef;
  184.         final double ch = chAverageFunction.value(absLatidude) - chAmplitudeFunction.value(absLatidude) * cosCoef;

  186.         final double[] function = new double[2];

  188.         // Hydrostatic mapping factor
  189.         function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch, trackingCoordinates.getElevation());

  191.         // Wet mapping factor
  192.         function[1] = TroposphericModelUtils.mappingFunction(awFunction.value(absLatidude),
  193.                                                              bwFunction.value(absLatidude),
  194.                                                              cwFunction.value(absLatidude),
  195.                                                              trackingCoordinates.getElevation());

  197.         // Apply height correction
  198.         final double correction = TroposphericModelUtils.computeHeightCorrection(trackingCoordinates.getElevation(),
  199.                                                                                  point.getAltitude());
  200.         function[0] = function[0] + correction;

  202.         return function;
  203.     }

  205.     /** {@inheritDoc} */
  206.     @Override
  207.     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
  208.                                                                   final FieldGeodeticPoint<T> point,
  209.                                                                   final FieldAbsoluteDate<T> date) {
  210.         final Field<T> field = date.getField();
  211.         final T zero = field.getZero();

  213.         // Temporal factor
  214.         double t0 = 28;
  215.         if (point.getLatitude().getReal() < 0) {
  216.             // southern hemisphere: t0 = 28 + an integer half of year
  217.             t0 += 183;
  218.         }
  219.         final T coef    = zero.getPi().multiply(2.0).multiply((date.getDayOfYear(utc).subtract(t0)).divide(365.25));
  220.         final T cosCoef = FastMath.cos(coef);

  222.         // Compute ah, bh and ch Eq. 5
  223.         double absLatidude = FastMath.abs(point.getLatitude().getReal());
  224.         // there are no data in the model for latitudes lower than 15°
  225.         absLatidude = FastMath.max(FastMath.toRadians(15.0), absLatidude);
  226.         // there are no data in the model for latitudes greater than 75°
  227.         absLatidude = FastMath.min(FastMath.toRadians(75.0), absLatidude);
  228.         final T ah = cosCoef.multiply(ahAmplitudeFunction.value(absLatidude)).negate().add(ahAverageFunction.value(absLatidude));
  229.         final T bh = cosCoef.multiply(bhAmplitudeFunction.value(absLatidude)).negate().add(bhAverageFunction.value(absLatidude));
  230.         final T ch = cosCoef.multiply(chAmplitudeFunction.value(absLatidude)).negate().add(chAverageFunction.value(absLatidude));

  232.         final T[] function = MathArrays.buildArray(field, 2);

  234.         // Hydrostatic mapping factor
  235.         function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch,
  236.                                                              trackingCoordinates.getElevation());

  238.         // Wet mapping factor
  239.         function[1] = TroposphericModelUtils.mappingFunction(zero.newInstance(awFunction.value(absLatidude)), zero.newInstance(bwFunction.value(absLatidude)),
  240.                                                              zero.newInstance(cwFunction.value(absLatidude)), trackingCoordinates.getElevation());

  242.         // Apply height correction
  243.         final T correction = TroposphericModelUtils.computeHeightCorrection(trackingCoordinates.getElevation(),
  244.                                                                             point.getAltitude(),
  245.                                                                             field);
  246.         function[0] = function[0].add(correction);

  248.         return function;
  249.     }

  251. }
Created with JaCoCo 0.8.12.202403310830