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

  19. import java.util.Collections;
  20. import java.util.List;

  22. import org.hipparchus.CalculusFieldElement;
  23. import org.hipparchus.util.FastMath;
  24. import org.hipparchus.util.FieldSinCos;
  25. import org.hipparchus.util.MathUtils;
  26. import org.hipparchus.util.SinCos;
  27. import org.orekit.bodies.FieldGeodeticPoint;
  28. import org.orekit.bodies.GeodeticPoint;
  29. import org.orekit.models.earth.weather.FieldPressureTemperatureHumidity;
  30. import org.orekit.models.earth.weather.PressureTemperatureHumidity;
  31. import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
  32. import org.orekit.time.AbsoluteDate;
  33. import org.orekit.time.FieldAbsoluteDate;
  34. import org.orekit.utils.FieldTrackingCoordinates;
  35. import org.orekit.utils.ParameterDriver;
  36. import org.orekit.utils.TrackingCoordinates;

  38. /** The modified Hopfield model.
  39.  * <p>
  40.  * This model from Hopfield 1969, 1970, 1972 is described in equations
  41.  * 5.105, 5.106, 5.107 and 5.108 in Guochang Xu, GPS - Theory, Algorithms
  42.  * and Applications, Springer, 2007.
  43.  * </p>
  44.  * @author Luc Maisonobe
  45.  * @see "Guochang Xu, GPS - Theory, Algorithms and Applications, Springer, 2007"
  46.  * @since 12.1
  47.  */
  48. public class ModifiedHopfieldModel implements TroposphericModel {

  50.     /** Constant for dry altitude effect. */
  51.     private static final double HD0 = 40136.0;

  53.     /** Slope for dry altitude effect. */
  54.     private static final double HD1 = 148.72;

  56.     /** Temperature reference. */
  57.     private static final double T0 = 273.16;

  59.     /** Constant for wet altitude effect. */
  60.     private static final double HW0 = 11000.0;

  62.     /** Dry delay factor. */
  63.     private static final double ND = 77.64e-6;

  65.     /** Wet delay factor, degree 1. */
  66.     private static final double NW1 = -12.96e-6;

  68.     /** Wet delay factor, degree 2. */
  69.     private static final double NW2 = 0.371800;

  71.     /** BAse radius. */
  72.     private static final double RE = 6378137.0;

  74.     /** Provider for pressure, temperature and humidity.
  75.      * @since 13.0
  76.      */
  77.     private final PressureTemperatureHumidityProvider pthProvider;

  79.     /** Create a new Hopfield model.
  80.      * @param pthProvider provider for pressure, temperature and humidity
  81.      * @since 13.0
  82.      */
  83.     public ModifiedHopfieldModel(final PressureTemperatureHumidityProvider pthProvider) {
  84.         this.pthProvider = pthProvider;
  85.     }

  87.     /** {@inheritDoc} */
  88.     @Override
  89.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
  90.                                        final GeodeticPoint point,
  91.                                        final double[] parameters, final AbsoluteDate date) {

  93.         // zenith angle
  94.         final double zenithAngle = MathUtils.SEMI_PI - trackingCoordinates.getElevation();

  96.         // compute weather parameters
  97.         final PressureTemperatureHumidity weather = pthProvider.getWeatherParameters(point, date);

  99.         // dry component
  100.         final double hd  = HD0 + HD1 * (weather.getTemperature() - T0);
  101.         final double nd  = ND * TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()) /
  102.                            weather.getTemperature();

  104.         // wet component
  105.         final double hw  = HW0;
  106.         final double nw  = (NW1 + NW2 / weather.getTemperature()) / weather.getTemperature();

  108.         return  new TroposphericDelay(delay(0.0,         hd, nd),
  109.                                       delay(0.0,         hw, nw),
  110.                                       delay(zenithAngle, hd, nd),
  111.                                       delay(zenithAngle, hw, nw));

  113.     }

  115.     /** {@inheritDoc}
  116.      * <p>
  117.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  118.      * reasons, we use the value for h = 0 when altitude is negative.
  119.      * </p>
  120.      * <p>
  121.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  122.      * elevations lower than a threshold will use the value obtained
  123.      * for the threshold itself.
  124.      * </p>
  125.      */
  126.     @Override
  127.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  128.                                                                                    final FieldGeodeticPoint<T> point,
  129.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

  131.         // zenith angle
  132.         final T zenithAngle = trackingCoordinates.getElevation().negate().add(MathUtils.SEMI_PI);

  134.         // compute weather parameters
  135.         final FieldPressureTemperatureHumidity<T> weather = pthProvider.getWeatherParameters(point, date);

  137.         // dry component
  138.         final T hd = weather.getTemperature().subtract(T0).multiply(HD1).add(HD0);
  139.         final T nd = TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()).
  140.                      multiply(ND).
  141.                      divide(weather.getTemperature());

  143.         // wet component
  144.         final T hw = date.getField().getZero().newInstance(HW0);
  145.         final T nw = weather.getTemperature().reciprocal().multiply(NW2).add(NW1).divide(weather.getTemperature());

  147.         return  new FieldTroposphericDelay<>(delay(date.getField().getZero(), hd, nd),
  148.                                              delay(date.getField().getZero(), hw, nw),
  149.                                              delay(zenithAngle,               hd, nd),
  150.                                              delay(zenithAngle,               hw, nw));

  152.     }

  154.     /** {@inheritDoc} */
  155.     @Override
  156.     public List<ParameterDriver> getParametersDrivers() {
  157.         return Collections.emptyList();
  158.     }

  160.     /** Compute the 9 terms sum delay.
  161.      * @param zenithAngle zenith angle
  162.      * @param hi altitude effect
  163.      * @param ni delay factor
  164.      * @return 9 terms sum delay
  165.      */
  166.     private double delay(final double zenithAngle, final double hi, final double ni) {

  168.         // equation 5.107
  169.         final SinCos scZ   = FastMath.sinCos(zenithAngle);
  170.         final double rePhi = RE + hi;
  171.         final double reS   = RE * scZ.sin();
  172.         final double reC   = RE * scZ.cos();
  173.         final double ri    = FastMath.sqrt(rePhi * rePhi - reS * reS) - reC;

  175.         final double ai    = -scZ.cos() / hi;
  176.         final double bi    = -scZ.sin() * scZ.sin() / (2 * hi * RE);
  177.         final double ai2   = ai * ai;
  178.         final double bi2   = bi * bi;

  180.         final double f1i   = 1;
  181.         final double f2i   = 4 * ai;
  182.         final double f3i   = 6 * ai2 + 4 * bi;
  183.         final double f4i   = 4 * ai * (ai2 + 3 * bi);
  184.         final double f5i   = ai2 * ai2 + 12 * ai2 * bi + 6 * bi2;
  185.         final double f6i   = 4 * ai * bi * (ai2 + 3 * bi);
  186.         final double f7i   = bi2 * (6 * ai2 + 4 * bi);
  187.         final double f8i   = 4 * ai * bi * bi2;
  188.         final double f9i   = bi2 * bi2;

  190.         return ni * (ri * (f1i +
  191.                            ri * (f2i / 2 +
  192.                                  ri * (f3i / 3 +
  193.                                        ri * (f4i / 4 +
  194.                                              ri * (f5i / 5 +
  195.                                                    ri * (f6i / 6 +
  196.                                                           ri * (f7i / 7 +
  197.                                                                 ri * (f8i / 8 +
  198.                                                                       ri * f9i / 9)))))))));

  200.     }

  202.     /** Compute the 9 terms sum delay.
  203.      * @param <T> type of the elements
  204.      * @param zenithAngle zenith angle
  205.      * @param hi altitude effect
  206.      * @param ni delay factor
  207.      * @return 9 terms sum delay
  208.      */
  209.     private <T extends CalculusFieldElement<T>> T delay(final T zenithAngle, final T hi, final T ni) {

  211.         // equation 5.107
  212.         final FieldSinCos<T> scZ   = FastMath.sinCos(zenithAngle);
  213.         final T rePhi = hi.add(RE);
  214.         final T reS   = scZ.sin().multiply(RE);
  215.         final T reC   = scZ.cos().multiply(RE);
  216.         final T ri    = FastMath.sqrt(rePhi.multiply(rePhi).subtract(reS.multiply(reS))).subtract(reC);

  218.         final T ai    = scZ.cos().negate().divide(hi);
  219.         final T bi    = scZ.sin().multiply(scZ.sin()).negate().divide(hi.add(hi).multiply(RE));
  220.         final T ai2   = ai.multiply(ai);
  221.         final T bi2   = bi.multiply(bi);

  223.         final T f1i   = ai.getField().getOne();
  224.         final T f2i   = ai.multiply(4);
  225.         final T f3i   = ai2.multiply(6).add(bi.multiply(4));
  226.         final T f4i   = ai.multiply(4).multiply(ai2.add(bi.multiply(3)));
  227.         final T f5i   = ai2.multiply(ai2).add(ai2.multiply(12).multiply(bi)).add(bi2.multiply(6));
  228.         final T f6i   = ai.multiply(4).multiply(bi).multiply(ai2.add(bi.multiply(3)));
  229.         final T f7i   = bi2.multiply(ai2.multiply(6).add(bi.multiply(4)));
  230.         final T f8i   = ai.multiply(4).multiply(bi).multiply(bi2);
  231.         final T f9i   = bi2.multiply(bi2);

  233.         return ni.
  234.                multiply(ri.
  235.                         multiply(f1i.
  236.                                  add(ri.
  237.                                      multiply(f2i.divide(2).
  238.                                               add(ri.
  239.                                                   multiply(f3i.divide(3).
  240.                                                            add(ri.
  241.                                                                multiply(f4i.divide(4).
  242.                                                                         add(ri.
  243.                                                                             multiply(f5i.divide(5).
  244.                                                                                      add(ri.
  245.                                                                                          multiply(f6i.divide(6).
  246.                                                                                                   add(ri.
  247.                                                                                                       multiply(f7i.divide(7).
  248.                                                                                                                add(ri.
  249.                                                                                                                    multiply(f8i.divide(8).
  250.                                                                                                                             add(ri.multiply(f9i.divide(9)))))))))))))))))));

  252.     }

  254. }
Created with JaCoCo 0.8.12.202403310830