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.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;
37  
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 {
49  
50      /** Constant for dry altitude effect. */
51      private static final double HD0 = 40136.0;
52  
53      /** Slope for dry altitude effect. */
54      private static final double HD1 = 148.72;
55  
56      /** Temperature reference. */
57      private static final double T0 = 273.16;
58  
59      /** Constant for wet altitude effect. */
60      private static final double HW0 = 11000.0;
61  
62      /** Dry delay factor. */
63      private static final double ND = 77.64e-6;
64  
65      /** Wet delay factor, degree 1. */
66      private static final double NW1 = -12.96e-6;
67  
68      /** Wet delay factor, degree 2. */
69      private static final double NW2 = 0.371800;
70  
71      /** BAse radius. */
72      private static final double RE = 6378137.0;
73  
74      /** Provider for pressure, temperature and humidity.
75       * @since 13.0
76       */
77      private final PressureTemperatureHumidityProvider pthProvider;
78  
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      }
86  
87      /** {@inheritDoc} */
88      @Override
89      public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
90                                         final GeodeticPoint point,
91                                         final double[] parameters, final AbsoluteDate date) {
92  
93          // zenith angle
94          final double zenithAngle = MathUtils.SEMI_PI - trackingCoordinates.getElevation();
95  
96          // compute weather parameters
97          final PressureTemperatureHumidity weather = pthProvider.getWeatherParameters(point, date);
98  
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();
103 
104         // wet component
105         final double hw  = HW0;
106         final double nw  = (NW1 + NW2 / weather.getTemperature()) / weather.getTemperature();
107 
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));
112 
113     }
114 
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) {
130 
131         // zenith angle
132         final T zenithAngle = trackingCoordinates.getElevation().negate().add(MathUtils.SEMI_PI);
133 
134         // compute weather parameters
135         final FieldPressureTemperatureHumidity<T> weather = pthProvider.getWeatherParameters(point, date);
136 
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());
142 
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());
146 
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));
151 
152     }
153 
154     /** {@inheritDoc} */
155     @Override
156     public List<ParameterDriver> getParametersDrivers() {
157         return Collections.emptyList();
158     }
159 
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) {
167 
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;
174 
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;
179 
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;
189 
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)))))))));
199 
200     }
201 
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) {
210 
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);
217 
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);
222 
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);
232 
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)))))))))))))))))));
251 
252     }
253 
254 }