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;
18  
19  import org.hipparchus.Field;
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.util.FastMath;
22  import org.orekit.models.earth.weather.ConstantPressureTemperatureHumidityProvider;
23  import org.orekit.models.earth.weather.PressureTemperatureHumidity;
24  import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
25  import org.orekit.models.earth.weather.water.CIPM2007;
26  import org.orekit.utils.units.Unit;
27  
28  /**
29   * Utility class for tropospheric models.
30   * @author Bryan Cazabonne
31   * @since 11.0
32   */
33  public class TroposphericModelUtils {
34  
35      /** Nanometers unit.
36       * @since 12.1
37       */
38      public static final Unit NANO_M = Unit.parse("nm");
39  
40      /** Micrometers unit.
41       * @since 12.1
42       */
43      public static final Unit MICRO_M = Unit.parse("µm");
44  
45      /** HectoPascal unit.
46       * @since 12.1
47       */
48      public static final Unit HECTO_PASCAL = Unit.parse("hPa");
49  
50      /** Standard atmosphere.
51       * <ul>
52       * <li>altitude: 0m</li>
53       * <li>temperature: 20 degree Celsius</li>
54       * <li>pressure: 1013.25 mbar</li>
55       * <li>humidity: 50%</li>
56       * </ul>
57       * @see #STANDARD_ATMOSPHERE_PROVIDER
58       * @since 12.1
59       */
60      public static final PressureTemperatureHumidity STANDARD_ATMOSPHERE;
61  
62      /** Provider for {@link #STANDARD_ATMOSPHERE standard atmosphere}.
63       * @since 12.1
64       */
65      public static final PressureTemperatureHumidityProvider STANDARD_ATMOSPHERE_PROVIDER;
66  
67      static {
68          final double h  = 0.0;
69          final double p  = HECTO_PASCAL.toSI(1013.25);
70          final double t  = 273.15 + 20;
71          final double rh = 0.5;
72          STANDARD_ATMOSPHERE = new PressureTemperatureHumidity(h, p, t,
73                                                                new CIPM2007().waterVaporPressure(p, t, rh),
74                                                                Double.NaN, Double.NaN);
75          STANDARD_ATMOSPHERE_PROVIDER =
76                          new ConstantPressureTemperatureHumidityProvider(STANDARD_ATMOSPHERE);
77      }
78  
79      /**
80       * Private constructor as class is a utility.
81       */
82      private TroposphericModelUtils() {
83          // Nothing to do
84      }
85  
86      /** Compute the mapping function related to the coefficient values and the elevation.
87       * @param a a coefficient
88       * @param b b coefficient
89       * @param c c coefficient
90       * @param elevation the elevation of the satellite, in radians.
91       * @return the value of the function at a given elevation
92       */
93      public static double mappingFunction(final double a, final double b, final double c, final double elevation) {
94          final double sinE = FastMath.sin(elevation);
95          // Numerator
96          final double numMP = 1 + a / (1 + b / (1 + c));
97          // Denominator
98          final double denMP = sinE + a / (sinE + b / (sinE + c));
99  
100         final double fElevation = numMP / denMP;
101 
102         return fElevation;
103     }
104 
105     /** Compute the mapping function related to the coefficient values and the elevation.
106      * @param <T> type of the elements
107      * @param a a coefficient
108      * @param b b coefficient
109      * @param c c coefficient
110      * @param elevation the elevation of the satellite, in radians.
111      * @return the value of the function at a given elevation
112      */
113     public static <T extends CalculusFieldElement<T>> T mappingFunction(final T a, final T b, final T c, final T elevation) {
114         final T sinE = FastMath.sin(elevation);
115         // Numerator
116         final T numMP = a.divide(b.divide(c.add(1.0)).add(1.0)).add(1.0);
117         // Denominator
118         final T denMP = a.divide(b.divide(c.add(sinE)).add(sinE)).add(sinE);
119 
120         final T fElevation = numMP.divide(denMP);
121 
122         return fElevation;
123     }
124 
125     /** This method computes the height correction for the hydrostatic
126      *  component of the mapping function.
127      *  The formulas are given by Neill's paper, 1996:
128      *<p>
129      *      Niell A. E. (1996)
130      *      "Global mapping functions for the atmosphere delay of radio wavelengths,”
131      *      J. Geophys. Res., 101(B2), pp.  3227–3246, doi:  10.1029/95JB03048.
132      *</p>
133      * @param elevation the elevation of the satellite, in radians.
134      * @param height the height of the station in m above sea level.
135      * @return the height correction, in m
136      */
137     public static double computeHeightCorrection(final double elevation, final double height) {
138         final double fixedHeight = FastMath.max(0.0, height);
139         final double sinE = FastMath.sin(elevation);
140         // Ref: Eq. 4
141         final double function = TroposphericModelUtils.mappingFunction(2.53e-5, 5.49e-3, 1.14e-3, elevation);
142         // Ref: Eq. 6
143         final double dmdh = (1 / sinE) - function;
144         // Ref: Eq. 7
145         final double correction = dmdh * (fixedHeight / 1000.0);
146         return correction;
147     }
148 
149     /** This method computes the height correction for the hydrostatic
150      *  component of the mapping function.
151      *  The formulas are given by Neill's paper, 1996:
152      *<p>
153      *      Niell A. E. (1996)
154      *      "Global mapping functions for the atmosphere delay of radio wavelengths,”
155      *      J. Geophys. Res., 101(B2), pp.  3227–3246, doi:  10.1029/95JB03048.
156      *</p>
157      * @param <T> type of the elements
158      * @param elevation the elevation of the satellite, in radians.
159      * @param height the height of the station in m above sea level.
160      * @param field field to which the elements belong
161      * @return the height correction, in m
162      */
163     public static <T extends CalculusFieldElement<T>> T computeHeightCorrection(final T elevation, final T height, final Field<T> field) {
164         final T zero = field.getZero();
165         final T fixedHeight = FastMath.max(zero, height);
166         final T sinE = FastMath.sin(elevation);
167         // Ref: Eq. 4
168         final T function = TroposphericModelUtils.mappingFunction(zero.newInstance(2.53e-5), zero.newInstance(5.49e-3), zero.newInstance(1.14e-3), elevation);
169         // Ref: Eq. 6
170         final T dmdh = sinE.reciprocal().subtract(function);
171         // Ref: Eq. 7
172         final T correction = dmdh.multiply(fixedHeight.divide(1000.0));
173         return correction;
174     }
175 
176 }