NiellMappingFunctionModel.java

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

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

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

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

  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.     };

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

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

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

  72.     /** X values for the ch amplitude function. */
  73.     private static final double[] VALUES_FOR_CH_AMPLITUDE = {
  74.         0.0, 9.0128400e-5, 4.3497037e-5, 84.795348e-5, 170.37206e-5
  75.     };

  76.     /** Values for the aw function. */
  77.     private static final double[] VALUES_FOR_AW = {
  78.         5.8021897e-4, 5.6794847e-4, 5.8118019e-4, 5.9727542e-4, 6.1641693e-4
  79.     };

  80.     /** Values for the bw function. */
  81.     private static final double[] VALUES_FOR_BW = {
  82.         1.4275268e-3, 1.5138625e-3, 1.4572752e-3, 1.5007428e-3, 1.7599082e-3
  83.     };

  84.     /** Values for the cw function. */
  85.     private static final double[] VALUES_FOR_CW = {
  86.         4.3472961e-2, 4.6729510e-2, 4.3908931e-2, 4.4626982e-2, 5.4736038e-2
  87.     };

  88.     /** Values for the cw function. */
  89.     private static final double[] LATITUDE_VALUES = {
  90.         FastMath.toRadians(15.0), FastMath.toRadians(30.0), FastMath.toRadians(45.0), FastMath.toRadians(60.0), FastMath.toRadians(75.0),
  91.     };

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

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

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

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

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

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

  104.     /** Interpolation function for the aw term. */
  105.     private final UnivariateFunction awFunction;

  106.     /** Interpolation function for the bw term. */
  107.     private final UnivariateFunction bwFunction;

  108.     /** Interpolation function for the cw term. */
  109.     private final UnivariateFunction cwFunction;

  110.     /** UTC time scale. */
  111.     private final TimeScale utc;

  112.     /** Builds a new instance.
  113.      *
  114.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  115.      *
  116.      * @see #NiellMappingFunctionModel(TimeScale)
  117.      */
  118.     @DefaultDataContext
  119.     public NiellMappingFunctionModel() {
  120.         this(DataContext.getDefault().getTimeScales().getUTC());
  121.     }

  122.     /** Builds a new instance.
  123.      * @param utc UTC time scale.
  124.      * @since 10.1
  125.      */
  126.     public NiellMappingFunctionModel(final TimeScale utc) {
  127.         this.utc = utc;
  128.         // Interpolation functions for hydrostatic coefficients
  129.         this.ahAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AH_AVERAGE);
  130.         this.bhAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BH_AVERAGE);
  131.         this.chAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CH_AVERAGE);
  132.         this.ahAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AH_AMPLITUDE);
  133.         this.bhAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BH_AMPLITUDE);
  134.         this.chAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CH_AMPLITUDE);

  135.         // Interpolation functions for wet coefficients
  136.         this.awFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AW);
  137.         this.bwFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BW);
  138.         this.cwFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CW);
  139.     }

  140.     /** {@inheritDoc} */
  141.     @Override
  142.     @Deprecated
  143.     public double[] mappingFactors(final double elevation, final GeodeticPoint point,
  144.                                    final AbsoluteDate date) {
  145.         return mappingFactors(new TrackingCoordinates(0.0, elevation, 0.0),
  146.                               point,
  147.                               TroposphericModelUtils.STANDARD_ATMOSPHERE,
  148.                               date);
  149.     }

  150.     /** {@inheritDoc} */
  151.     @Override
  152.     public double[] mappingFactors(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
  153.                                    final PressureTemperatureHumidity weather,
  154.                                    final AbsoluteDate date) {
  155.         // Day of year computation
  156.         final DateTimeComponents dtc = date.getComponents(utc);
  157.         final int dofyear = dtc.getDate().getDayOfYear();

  158.         // Temporal factor
  159.         double t0 = 28;
  160.         if (point.getLatitude() < 0) {
  161.             // southern hemisphere: t0 = 28 + an integer half of year
  162.             t0 += 183;
  163.         }
  164.         final double coef    = 2 * FastMath.PI * ((dofyear - t0) / 365.25);
  165.         final double cosCoef = FastMath.cos(coef);

  166.         // Compute ah, bh and ch Eq. 5
  167.         double absLatidude = FastMath.abs(point.getLatitude());
  168.         // there are no data in the model for latitudes lower than 15°
  169.         absLatidude = FastMath.max(FastMath.toRadians(15.0), absLatidude);
  170.         // there are no data in the model for latitudes greater than 75°
  171.         absLatidude = FastMath.min(FastMath.toRadians(75.0), absLatidude);
  172.         final double ah = ahAverageFunction.value(absLatidude) - ahAmplitudeFunction.value(absLatidude) * cosCoef;
  173.         final double bh = bhAverageFunction.value(absLatidude) - bhAmplitudeFunction.value(absLatidude) * cosCoef;
  174.         final double ch = chAverageFunction.value(absLatidude) - chAmplitudeFunction.value(absLatidude) * cosCoef;

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

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

  178.         // Wet mapping factor
  179.         function[1] = TroposphericModelUtils.mappingFunction(awFunction.value(absLatidude),
  180.                                                              bwFunction.value(absLatidude),
  181.                                                              cwFunction.value(absLatidude),
  182.                                                              trackingCoordinates.getElevation());

  183.         // Apply height correction
  184.         final double correction = TroposphericModelUtils.computeHeightCorrection(trackingCoordinates.getElevation(),
  185.                                                                                  point.getAltitude());
  186.         function[0] = function[0] + correction;

  187.         return function;
  188.     }

  189.     /** {@inheritDoc} */
  190.     @Override
  191.     @Deprecated
  192.     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final T elevation, final FieldGeodeticPoint<T> point,
  193.                                                                   final FieldAbsoluteDate<T> date) {
  194.         return mappingFactors(new FieldTrackingCoordinates<>(date.getField().getZero(), elevation, date.getField().getZero()),
  195.                               point,
  196.                               new FieldPressureTemperatureHumidity<>(date.getField(),
  197.                                                                      TroposphericModelUtils.STANDARD_ATMOSPHERE),
  198.                               date);
  199.     }

  200.     /** {@inheritDoc} */
  201.     @Override
  202.     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
  203.                                                                   final FieldGeodeticPoint<T> point,
  204.                                                                   final FieldPressureTemperatureHumidity<T> weather,
  205.                                                                   final FieldAbsoluteDate<T> date) {
  206.         final Field<T> field = date.getField();
  207.         final T zero = field.getZero();

  208.         // Day of year computation
  209.         final DateTimeComponents dtc = date.getComponents(utc);
  210.         final int dofyear = dtc.getDate().getDayOfYear();

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

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

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

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

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

  235.         // Apply height correction
  236.         final T correction = TroposphericModelUtils.computeHeightCorrection(trackingCoordinates.getElevation(),
  237.                                                                             point.getAltitude(),
  238.                                                                             field);
  239.         function[0] = function[0].add(correction);

  240.         return function;
  241.     }

  242. }