ModifiedHopfieldModel.java

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

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

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

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

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

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

  49.     /** Temperature reference. */
  50.     private static final double T0 = 273.16;

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

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

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

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

  59.     /** BAse radius. */
  60.     private static final double RE = 6378137.0;

  61.     /** Create a new Hopfield model.
  62.      */
  63.     public ModifiedHopfieldModel() {
  64.         // nothing to do
  65.     }

  66.     /** {@inheritDoc} */
  67.     @Override
  68.     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
  69.                                        final GeodeticPoint point,
  70.                                        final PressureTemperatureHumidity weather,
  71.                                        final double[] parameters, final AbsoluteDate date) {

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

  74.         // dry component
  75.         final double hd  = HD0 + HD1 * (weather.getTemperature() - T0);
  76.         final double nd  = ND * TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()) /
  77.                            weather.getTemperature();

  78.         // wet component
  79.         final double hw  = HW0;
  80.         final double nw  = (NW1 + NW2 / weather.getTemperature()) / weather.getTemperature();

  81.         return  new TroposphericDelay(delay(0.0,         hd, nd),
  82.                                       delay(0.0,         hw, nw),
  83.                                       delay(zenithAngle, hd, nd),
  84.                                       delay(zenithAngle, hw, nw));

  85.     }

  86.     /** {@inheritDoc}
  87.      * <p>
  88.      * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
  89.      * reasons, we use the value for h = 0 when altitude is negative.
  90.      * </p>
  91.      * <p>
  92.      * There are also numerical issues for elevation angles close to zero. For continuity reasons,
  93.      * elevations lower than a threshold will use the value obtained
  94.      * for the threshold itself.
  95.      * </p>
  96.      */
  97.     @Override
  98.     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
  99.                                                                                    final FieldGeodeticPoint<T> point,
  100.                                                                                    final FieldPressureTemperatureHumidity<T> weather,
  101.                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {

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

  104.         // dry component
  105.         final T hd = weather.getTemperature().subtract(T0).multiply(HD1).add(HD0);
  106.         final T nd = TroposphericModelUtils.HECTO_PASCAL.fromSI(weather.getPressure()).
  107.                      multiply(ND).
  108.                      divide(weather.getTemperature());

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

  112.         return  new FieldTroposphericDelay<>(delay(date.getField().getZero(), hd, nd),
  113.                                              delay(date.getField().getZero(), hw, nw),
  114.                                              delay(zenithAngle,               hd, nd),
  115.                                              delay(zenithAngle,               hw, nw));

  116.     }

  117.     /** {@inheritDoc} */
  118.     @Override
  119.     public List<ParameterDriver> getParametersDrivers() {
  120.         return Collections.emptyList();
  121.     }

  122.     /** Compute the 9 terms sum delay.
  123.      * @param zenithAngle zenith angle
  124.      * @param hi altitude effect
  125.      * @param ni delay factor
  126.      * @return 9 terms sum delay
  127.      */
  128.     private double delay(final double zenithAngle, final double hi, final double ni) {

  129.         // equation 5.107
  130.         final SinCos scZ   = FastMath.sinCos(zenithAngle);
  131.         final double rePhi = RE + hi;
  132.         final double reS   = RE * scZ.sin();
  133.         final double reC   = RE * scZ.cos();
  134.         final double ri    = FastMath.sqrt(rePhi * rePhi - reS * reS) - reC;

  135.         final double ai    = -scZ.cos() / hi;
  136.         final double bi    = -scZ.sin() * scZ.sin() / (2 * hi * RE);
  137.         final double ai2   = ai * ai;
  138.         final double bi2   = bi * bi;

  139.         final double f1i   = 1;
  140.         final double f2i   = 4 * ai;
  141.         final double f3i   = 6 * ai2 + 4 * bi;
  142.         final double f4i   = 4 * ai * (ai2 + 3 * bi);
  143.         final double f5i   = ai2 * ai2 + 12 * ai2 * bi + 6 * bi2;
  144.         final double f6i   = 4 * ai * bi * (ai2 + 3 * bi);
  145.         final double f7i   = bi2 * (6 * ai2 + 4 * bi);
  146.         final double f8i   = 4 * ai * bi * bi2;
  147.         final double f9i   = bi2 * bi2;

  148.         return ni * (ri * (f1i +
  149.                            ri * (f2i / 2 +
  150.                                  ri * (f3i / 3 +
  151.                                        ri * (f4i / 4 +
  152.                                              ri * (f5i / 5 +
  153.                                                    ri * (f6i / 6 +
  154.                                                           ri * (f7i / 7 +
  155.                                                                 ri * (f8i / 8 +
  156.                                                                       ri * f9i / 9)))))))));

  157.     }

  158.     /** Compute the 9 terms sum delay.
  159.      * @param <T> type of the elements
  160.      * @param zenithAngle zenith angle
  161.      * @param hi altitude effect
  162.      * @param ni delay factor
  163.      * @return 9 terms sum delay
  164.      */
  165.     private <T extends CalculusFieldElement<T>> T delay(final T zenithAngle, final T hi, final T ni) {

  166.         // equation 5.107
  167.         final FieldSinCos<T> scZ   = FastMath.sinCos(zenithAngle);
  168.         final T rePhi = hi.add(RE);
  169.         final T reS   = scZ.sin().multiply(RE);
  170.         final T reC   = scZ.cos().multiply(RE);
  171.         final T ri    = FastMath.sqrt(rePhi.multiply(rePhi).subtract(reS.multiply(reS))).subtract(reC);

  172.         final T ai    = scZ.cos().negate().divide(hi);
  173.         final T bi    = scZ.sin().multiply(scZ.sin()).negate().divide(hi.add(hi).multiply(RE));
  174.         final T ai2   = ai.multiply(ai);
  175.         final T bi2   = bi.multiply(bi);

  176.         final T f1i   = ai.getField().getOne();
  177.         final T f2i   = ai.multiply(4);
  178.         final T f3i   = ai2.multiply(6).add(bi.multiply(4));
  179.         final T f4i   = ai.multiply(4).multiply(ai2.add(bi.multiply(3)));
  180.         final T f5i   = ai2.multiply(ai2).add(ai2.multiply(12).multiply(bi)).add(bi2.multiply(6));
  181.         final T f6i   = ai.multiply(4).multiply(bi).multiply(ai2.add(bi.multiply(3)));
  182.         final T f7i   = bi2.multiply(ai2.multiply(6).add(bi.multiply(4)));
  183.         final T f8i   = ai.multiply(4).multiply(bi).multiply(bi2);
  184.         final T f9i   = bi2.multiply(bi2);

  185.         return ni.
  186.                multiply(ri.
  187.                         multiply(f1i.
  188.                                  add(ri.
  189.                                      multiply(f2i.divide(2).
  190.                                               add(ri.
  191.                                                   multiply(f3i.divide(3).
  192.                                                            add(ri.
  193.                                                                multiply(f4i.divide(4).
  194.                                                                         add(ri.
  195.                                                                             multiply(f5i.divide(5).
  196.                                                                                      add(ri.
  197.                                                                                          multiply(f6i.divide(6).
  198.                                                                                                   add(ri.
  199.                                                                                                       multiply(f7i.divide(7).
  200.                                                                                                                add(ri.
  201.                                                                                                                    multiply(f8i.divide(8).
  202.                                                                                                                             add(ri.multiply(f9i.divide(9)))))))))))))))))));

  203.     }

  204. }