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;
18
19 import java.util.Collections;
20 import java.util.List;
21
22 import org.hipparchus.CalculusFieldElement;
23 import org.hipparchus.Field;
24 import org.hipparchus.analysis.interpolation.LinearInterpolator;
25 import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
26 import org.hipparchus.util.FastMath;
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.HeightDependentPressureTemperatureHumidityConverter;
31 import org.orekit.models.earth.weather.PressureTemperatureHumidity;
32 import org.orekit.models.earth.weather.PressureTemperatureHumidityProvider;
33 import org.orekit.time.AbsoluteDate;
34 import org.orekit.time.FieldAbsoluteDate;
35 import org.orekit.utils.FieldTrackingCoordinates;
36 import org.orekit.utils.ParameterDriver;
37 import org.orekit.utils.TrackingCoordinates;
38
39 /** The canonical Saastamoinen model.
40 * <p>
41 * Estimates the path delay imposed to
42 * electro-magnetic signals by the troposphere according to the formula:
43 * \[
44 * \delta = \frac{0.002277}{\cos z} \left[P+(\frac{1255}{T}+0.005)e - B(h) \tan^2 z\right]
45 * \]
46 * with the following input data provided to the model:
47 * <ul>
48 * <li>z: zenith angle</li>
49 * <li>P: atmospheric pressure</li>
50 * <li>T: temperature</li>
51 * <li>e: partial pressure of water vapor</li>
52 * </ul>
53 * @author Luc Maisonobe
54 * @see "J Saastamoinen, Atmospheric Correction for the Troposphere and Stratosphere in Radio
55 * Ranging of Satellites"
56 * @since 12.1
57 */
58 public class CanonicalSaastamoinenModel implements TroposphericModel {
59
60 /** Default lowest acceptable elevation angle [rad]. */
61 public static final double DEFAULT_LOW_ELEVATION_THRESHOLD = 0.05;
62
63 /** Base delay coefficient. */
64 private static final double L0 = 2.2768e-5;
65
66 /** Temperature numerator. */
67 private static final double T_NUM = 1255;
68
69 /** Wet offset. */
70 private static final double WET_OFFSET = 0.05;
71
72 /** X values for the B function (table 1 in reference paper). */
73 private static final double[] X_VALUES_FOR_B = {
74 0.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1500.0, 2000.0, 2500.0, 3000.0, 4000.0, 5000.0, 6000.0
75 };
76
77 /** Y values for the B function (table 1 in reference paper).
78 * <p>
79 * The values have been scaled up by a factor 100.0 due to conversion from hPa to Pa.
80 * </p>
81 */
82 private static final double[] Y_VALUES_FOR_B = {
83 116.0, 113.0, 110.0, 107.0, 104.0, 101.0, 94.0, 88.0, 82.0, 76.0, 66.0, 57.0, 49.0
84 };
85
86 /** Interpolation function for the B correction term. */
87 private static final PolynomialSplineFunction B_FUNCTION;
88
89 static {
90 B_FUNCTION = new LinearInterpolator().interpolate(X_VALUES_FOR_B, Y_VALUES_FOR_B);
91 }
92
93 /** Lowest acceptable elevation angle [rad]. */
94 private double lowElevationThreshold;
95
96 /** Provider for pressure, temperature and humidity. */
97 private final PressureTemperatureHumidityProvider pthProvider;
98
99 /**
100 * Create a new Saastamoinen model for the troposphere using the given environmental
101 * conditions and table from the reference book.
102 *
103 * @see HeightDependentPressureTemperatureHumidityConverter
104 */
105 public CanonicalSaastamoinenModel() {
106 this(TroposphericModelUtils.STANDARD_ATMOSPHERE_PROVIDER);
107 }
108
109 /**
110 * Create a new Saastamoinen model for the troposphere using the given environmental
111 * conditions and table from the reference book.
112 * @param pthProvider provider for pressure, temperature and humidity
113 * @since 13.0
114 */
115 public CanonicalSaastamoinenModel(final PressureTemperatureHumidityProvider pthProvider) {
116 this.lowElevationThreshold = DEFAULT_LOW_ELEVATION_THRESHOLD;
117 this.pthProvider = pthProvider;
118 }
119
120 /** {@inheritDoc}
121 * <p>
122 * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
123 * reasons, we use the value for h = 0 when altitude is negative.
124 * </p>
125 * <p>
126 * There are also numerical issues for elevation angles close to zero. For continuity reasons,
127 * elevations lower than a threshold will use the value obtained
128 * for the threshold itself.
129 * </p>
130 * @see #getLowElevationThreshold()
131 * @see #setLowElevationThreshold(double)
132 */
133 @Override
134 public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
135 final double[] parameters, final AbsoluteDate date) {
136
137 // there are no data in the model for negative altitudes and altitude bigger than 6000 m
138 // limit the height to a range of [0, 5000] m
139 final double fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
140 X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);
141
142 // interpolate the b correction term
143 final double B = B_FUNCTION.value(fixedHeight);
144
145 // calculate the zenith angle from the elevation
146 final double z = FastMath.abs(0.5 * FastMath.PI - FastMath.max(trackingCoordinates.getElevation(),
147 lowElevationThreshold));
148
149 // compute weather parameters
150 final PressureTemperatureHumidity weather = pthProvider.getWeatherParameters(point, date);
151
152 // calculate the path delay
153 final double invCos = 1.0 / FastMath.cos(z);
154 final double tan = FastMath.tan(z);
155 final double zh = L0 * weather.getPressure();
156 final double zw = L0 * (T_NUM / weather.getTemperature() + WET_OFFSET) * weather.getWaterVaporPressure();
157 final double sh = zh * invCos;
158 final double sw = (zw - L0 * B * tan * tan) * invCos;
159 return new TroposphericDelay(zh, zw, sh, sw);
160
161 }
162
163 /** {@inheritDoc}
164 * <p>
165 * The Saastamoinen model is not defined for altitudes below 0.0. for continuity
166 * reasons, we use the value for h = 0 when altitude is negative.
167 * </p>
168 * <p>
169 * There are also numerical issues for elevation angles close to zero. For continuity reasons,
170 * elevations lower than a threshold will use the value obtained
171 * for the threshold itself.
172 * </p>
173 * @see #getLowElevationThreshold()
174 * @see #setLowElevationThreshold(double)
175 */
176 @Override
177 public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
178 final FieldGeodeticPoint<T> point,
179 final T[] parameters, final FieldAbsoluteDate<T> date) {
180
181 final Field<T> field = date.getField();
182 final T zero = field.getZero();
183
184 // there are no data in the model for negative altitudes and altitude bigger than 5000 m
185 // limit the height to a range of [0, 5000] m
186 final T fixedHeight = FastMath.min(FastMath.max(point.getAltitude(), X_VALUES_FOR_B[0]),
187 X_VALUES_FOR_B[X_VALUES_FOR_B.length - 1]);
188
189 // interpolate the b correction term
190 final T B = B_FUNCTION.value(fixedHeight);
191
192 // calculate the zenith angle from the elevation
193 final T z = FastMath.abs(zero.getPi().multiply(0.5).
194 subtract(FastMath.max(trackingCoordinates.getElevation(), lowElevationThreshold)));
195
196 // compute weather parameters
197 final FieldPressureTemperatureHumidity<T> weather = pthProvider.getWeatherParameters(point, date);
198
199 // calculate the path delay in m
200 final T invCos = FastMath.cos(z).reciprocal();
201 final T tan = FastMath.tan(z);
202 final T zh = weather.getPressure().multiply(L0);
203 final T zw = weather.getTemperature().reciprocal().multiply(T_NUM).add(WET_OFFSET).
204 multiply(weather.getWaterVaporPressure()).multiply(L0);
205 final T sh = zh.multiply(invCos);
206 final T sw = zw.subtract(B.multiply(tan).multiply(tan).multiply(L0)).multiply(invCos);
207 return new FieldTroposphericDelay<>(zh, zw, sh, sw);
208
209 }
210
211 /** {@inheritDoc} */
212 @Override
213 public List<ParameterDriver> getParametersDrivers() {
214 return Collections.emptyList();
215 }
216
217 /** Get the low elevation threshold value for path delay computation.
218 * @return low elevation threshold, in rad.
219 * @see #pathDelay(TrackingCoordinates, GeodeticPoint, double[], AbsoluteDate)
220 * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
221 */
222 public double getLowElevationThreshold() {
223 return lowElevationThreshold;
224 }
225
226 /** Set the low elevation threshold value for path delay computation.
227 * @param lowElevationThreshold The new value for the threshold [rad]
228 * @see #pathDelay(TrackingCoordinates, GeodeticPoint, double[], AbsoluteDate)
229 * @see #pathDelay(FieldTrackingCoordinates, FieldGeodeticPoint, CalculusFieldElement[], FieldAbsoluteDate)
230 */
231 public void setLowElevationThreshold(final double lowElevationThreshold) {
232 this.lowElevationThreshold = lowElevationThreshold;
233 }
234
235 }
236