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 java.util.Collections;
20  import java.util.List;
21  
22  import org.hipparchus.CalculusFieldElement;
23  import org.hipparchus.util.FastMath;
24  import org.orekit.annotation.DefaultDataContext;
25  import org.orekit.bodies.FieldGeodeticPoint;
26  import org.orekit.bodies.GeodeticPoint;
27  import org.orekit.models.earth.weather.ConstantPressureTemperatureHumidityProvider;
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  
35  /** An estimated tropospheric model. The tropospheric delay is computed according to the formula:
36   * <p>
37   * δ = δ<sub>h</sub> * m<sub>h</sub> + (δ<sub>t</sub> - δ<sub>h</sub>) * m<sub>w</sub>
38   * <p>
39   * With:
40   * <ul>
41   * <li>δ<sub>h</sub>: Tropospheric zenith hydro-static delay.</li>
42   * <li>δ<sub>t</sub>: Tropospheric total zenith delay.</li>
43   * <li>m<sub>h</sub>: Hydro-static mapping function.</li>
44   * <li>m<sub>w</sub>: Wet mapping function.</li>
45   * </ul>
46   * <p>
47   * The mapping functions m<sub>h</sub>(e) and m<sub>w</sub>(e) are
48   * computed thanks to a {@link #model} initialized by the user.
49   * The user has the possibility to use several mapping function models for the computations:
50   * the {@link GlobalMappingFunctionModel Global Mapping Function}, or
51   * the {@link NiellMappingFunctionModel Niell Mapping Function}
52   * </p> <p>
53   * The tropospheric zenith delay δ<sub>h</sub> is computed empirically with a
54   * {@link TroposphericModel tropospheric model}
55   * while the tropospheric total zenith delay δ<sub>t</sub> is estimated as a {@link ParameterDriver},
56   * hence the wet part is the difference between the two.
57   * @since 12.1
58   */
59  public class EstimatedModel implements TroposphericModel {
60  
61      /** Name of the parameter of this model: the total zenith delay. */
62      public static final String TOTAL_ZENITH_DELAY = "total zenith delay";
63  
64      /** Mapping Function model. */
65      private final TroposphereMappingFunction model;
66  
67      /** Driver for the tropospheric zenith total delay.*/
68      private final ParameterDriver totalZenithDelay;
69  
70      /** Model for hydrostatic component. */
71      private final TroposphericModel hydrostatic;
72  
73      /** Build a new instance using the given environmental conditions.
74       * <p>
75       * This constructor uses a {@link ModifiedSaastamoinenModel} for the hydrostatic contribution.
76       * </p>
77       * @param h0 altitude of the station [m]
78       * @param t0 the temperature at the station [K]
79       * @param p0 the atmospheric pressure at the station [mbar]
80       * @param model mapping function model.
81       * @param totalDelay initial value for the tropospheric zenith total delay [m]
82       */
83      @DefaultDataContext
84      public EstimatedModel(final double h0, final double t0, final double p0,
85                            final TroposphereMappingFunction model, final double totalDelay) {
86          this(new ModifiedSaastamoinenModel(new ConstantPressureTemperatureHumidityProvider(new PressureTemperatureHumidity(h0,
87                                                                                                                             TroposphericModelUtils.HECTO_PASCAL.toSI(p0),
88                                                                                                                             t0,
89                                                                                                                             0.0,
90                                                                                                                             Double.NaN,
91                                                                                                                             Double.NaN))),
92               model, totalDelay);
93      }
94  
95      /** Build a new instance using the given environmental conditions.
96       * @param hydrostatic model for hydrostatic component
97       * @param model mapping function model.
98       * @param totalDelay initial value for the tropospheric zenith total delay [m]
99       * @since 12.1
100      */
101     public EstimatedModel(final TroposphericModel hydrostatic,
102                           final TroposphereMappingFunction model,
103                           final double totalDelay) {
104 
105         totalZenithDelay = new ParameterDriver(EstimatedModel.TOTAL_ZENITH_DELAY,
106                                                totalDelay, FastMath.scalb(1.0, 0), 0.0, Double.POSITIVE_INFINITY);
107 
108         this.hydrostatic = hydrostatic;
109         this.model = model;
110     }
111 
112     /** Build a new instance using a standard atmosphere model.
113      * <ul>
114      * <li>altitude: 0m
115      * <li>temperature: 18 degree Celsius
116      * <li>pressure: 1013.25 mbar
117      * </ul>
118      * @param model mapping function model.
119      * @param totalDelay initial value for the tropospheric zenith total delay [m]
120      */
121     @DefaultDataContext
122     public EstimatedModel(final TroposphereMappingFunction model, final double totalDelay) {
123         this(0.0, 273.15 + 18.0, 1013.25, model, totalDelay);
124     }
125 
126     /** {@inheritDoc} */
127     @Override
128     public List<ParameterDriver> getParametersDrivers() {
129         return Collections.singletonList(totalZenithDelay);
130     }
131 
132     /** {@inheritDoc} */
133     @Override
134     public TroposphericDelay pathDelay(final TrackingCoordinates trackingCoordinates,
135                                        final GeodeticPoint point,
136                                        final double[] parameters, final AbsoluteDate date) {
137 
138         // zenith hydrostatic delay
139         final double zd = hydrostatic.pathDelay(trackingCoordinates, point, parameters, date).getZh();
140 
141         // zenith wet delay
142         final double wd = parameters[0] - zd;
143 
144         // mapping functions
145         final double[] mf = model.mappingFactors(trackingCoordinates, point, date);
146 
147         // composite delay
148         return new TroposphericDelay(zd, wd, mf[0] * zd, mf[1] * wd);
149 
150     }
151 
152     /** {@inheritDoc} */
153     @Override
154     public <T extends CalculusFieldElement<T>> FieldTroposphericDelay<T> pathDelay(final FieldTrackingCoordinates<T> trackingCoordinates,
155                                                                                    final FieldGeodeticPoint<T> point,
156                                                                                    final T[] parameters, final FieldAbsoluteDate<T> date) {
157 
158         // zenith hydrostatic delay
159         final T zd = hydrostatic.pathDelay(trackingCoordinates, point, parameters, date).getZh();
160 
161         // zenith wet delay
162         final T wd = parameters[0].subtract(zd);
163 
164         // mapping functions
165         final T[] mf = model.mappingFactors(trackingCoordinates, point, date);
166 
167         // composite delay
168         return new FieldTroposphericDelay<>(zd, wd, mf[0].multiply(zd), mf[1].multiply(wd));
169 
170     }
171 
172 }