1   /* Copyright 2002-2021 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.analysis.UnivariateFunction;
22  import org.hipparchus.analysis.interpolation.LinearInterpolator;
23  import org.hipparchus.util.FastMath;
24  import org.hipparchus.util.MathArrays;
25  import org.orekit.annotation.DefaultDataContext;
26  import org.orekit.bodies.FieldGeodeticPoint;
27  import org.orekit.bodies.GeodeticPoint;
28  import org.orekit.data.DataContext;
29  import org.orekit.time.AbsoluteDate;
30  import org.orekit.time.DateTimeComponents;
31  import org.orekit.time.FieldAbsoluteDate;
32  import org.orekit.time.TimeScale;
33  
34  /** The Niell Mapping Function  model for radio wavelengths.
35   *  This model is an empirical mapping function. It only needs the
36   *  values of the station latitude, height and the date for the computations.
37   *  <p>
38   *  With this model, the hydrostatic mapping function is time and latitude dependent
39   *  whereas the wet mapping function is only latitude dependent.
40   *  </p>
41   *
42   * @see "A. E. Niell(1996), Global mapping functions for the atmosphere delay of radio wavelengths,
43   *      J. Geophys. Res., 101(B2), pp.  3227–3246, doi:  10.1029/95JB03048."
44   *
45   * @author Bryan Cazabonne
46   *
47   */
48  public class NiellMappingFunctionModel implements MappingFunction {
49  
50      /** Values for the ah average function. */
51      private static final double[] VALUES_FOR_AH_AVERAGE = {
52          1.2769934e-3, 1.2683230e-3, 1.2465397e-3, 1.2196049e-3, 1.2045996e-3
53      };
54  
55      /** Values for the bh average function. */
56      private static final double[] VALUES_FOR_BH_AVERAGE = {
57          2.9153695e-3, 2.9152299e-3, 2.9288445e-3, 2.9022565e-3, 2.9024912e-3
58      };
59  
60      /** Values for the ch average function. */
61      private static final double[] VALUES_FOR_CH_AVERAGE = {
62          62.610505e-3, 62.837393e-3, 63.721774e-3, 63.824265e-3, 64.258455e-3
63      };
64  
65      /** Values for the ah amplitude function. */
66      private static final double[] VALUES_FOR_AH_AMPLITUDE = {
67          0.0, 1.2709626e-5, 2.6523662e-5, 3.4000452e-5, 4.1202191e-5
68      };
69  
70      /** Values for the bh amplitude function. */
71      private static final double[] VALUES_FOR_BH_AMPLITUDE = {
72          0.0, 2.1414979e-5, 3.0160779e-5, 7.2562722e-5, 11.723375e-5
73      };
74  
75      /** X values for the ch amplitude function. */
76      private static final double[] VALUES_FOR_CH_AMPLITUDE = {
77          0.0, 9.0128400e-5, 4.3497037e-5, 84.795348e-5, 170.37206e-5
78      };
79  
80      /** Values for the aw function. */
81      private static final double[] VALUES_FOR_AW = {
82          5.8021897e-4, 5.6794847e-4, 5.8118019e-4, 5.9727542e-4, 6.1641693e-4
83      };
84  
85      /** Values for the bw function. */
86      private static final double[] VALUES_FOR_BW = {
87          1.4275268e-3, 1.5138625e-3, 1.4572752e-3, 1.5007428e-3, 1.7599082e-3
88      };
89  
90      /** Values for the cw function. */
91      private static final double[] VALUES_FOR_CW = {
92          4.3472961e-2, 4.6729510e-2, 4.3908931e-2, 4.4626982e-2, 5.4736038e-2
93      };
94  
95      /** Values for the cw function. */
96      private static final double[] LATITUDE_VALUES = {
97          FastMath.toRadians(15.0), FastMath.toRadians(30.0), FastMath.toRadians(45.0), FastMath.toRadians(60.0), FastMath.toRadians(75.0),
98      };
99  
100     /** Interpolation function for the ah (average) term. */
101     private final UnivariateFunction ahAverageFunction;
102 
103     /** Interpolation function for the bh (average) term. */
104     private final UnivariateFunction bhAverageFunction;
105 
106     /** Interpolation function for the ch (average) term. */
107     private final UnivariateFunction chAverageFunction;
108 
109     /** Interpolation function for the ah (amplitude) term. */
110     private final UnivariateFunction ahAmplitudeFunction;
111 
112     /** Interpolation function for the bh (amplitude) term. */
113     private final UnivariateFunction bhAmplitudeFunction;
114 
115     /** Interpolation function for the ch (amplitude) term. */
116     private final UnivariateFunction chAmplitudeFunction;
117 
118     /** Interpolation function for the aw term. */
119     private final UnivariateFunction awFunction;
120 
121     /** Interpolation function for the bw term. */
122     private final UnivariateFunction bwFunction;
123 
124     /** Interpolation function for the cw term. */
125     private final UnivariateFunction cwFunction;
126 
127     /** UTC time scale. */
128     private final TimeScale utc;
129 
130     /** Builds a new instance.
131      *
132      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
133      *
134      * @see #NiellMappingFunctionModel(TimeScale)
135      */
136     @DefaultDataContext
137     public NiellMappingFunctionModel() {
138         this(DataContext.getDefault().getTimeScales().getUTC());
139     }
140 
141     /** Builds a new instance.
142      * @param utc UTC time scale.
143      * @since 10.1
144      */
145     public NiellMappingFunctionModel(final TimeScale utc) {
146         this.utc = utc;
147         // Interpolation functions for hydrostatic coefficients
148         this.ahAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AH_AVERAGE);
149         this.bhAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BH_AVERAGE);
150         this.chAverageFunction    = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CH_AVERAGE);
151         this.ahAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AH_AMPLITUDE);
152         this.bhAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BH_AMPLITUDE);
153         this.chAmplitudeFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CH_AMPLITUDE);
154 
155         // Interpolation functions for wet coefficients
156         this.awFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_AW);
157         this.bwFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_BW);
158         this.cwFunction  = new LinearInterpolator().interpolate(LATITUDE_VALUES, VALUES_FOR_CW);
159     }
160 
161     /** {@inheritDoc} */
162     @Override
163     public double[] mappingFactors(final double elevation, final GeodeticPoint point,
164                                    final AbsoluteDate date) {
165         // Day of year computation
166         final DateTimeComponents dtc = date.getComponents(utc);
167         final int dofyear = dtc.getDate().getDayOfYear();
168 
169         // Temporal factor
170         double t0 = 28;
171         if (point.getLatitude() < 0) {
172             // southern hemisphere: t0 = 28 + an integer half of year
173             t0 += 183;
174         }
175         final double coef    = 2 * FastMath.PI * ((dofyear - t0) / 365.25);
176         final double cosCoef = FastMath.cos(coef);
177 
178         // Compute ah, bh and ch Eq. 5
179         double absLatidude = FastMath.abs(point.getLatitude());
180         // there are no data in the model for latitudes lower than 15°
181         absLatidude = FastMath.max(FastMath.toRadians(15.0), absLatidude);
182         // there are no data in the model for latitudes greater than 75°
183         absLatidude = FastMath.min(FastMath.toRadians(75.0), absLatidude);
184         final double ah = ahAverageFunction.value(absLatidude) - ahAmplitudeFunction.value(absLatidude) * cosCoef;
185         final double bh = bhAverageFunction.value(absLatidude) - bhAmplitudeFunction.value(absLatidude) * cosCoef;
186         final double ch = chAverageFunction.value(absLatidude) - chAmplitudeFunction.value(absLatidude) * cosCoef;
187 
188         final double[] function = new double[2];
189 
190         // Hydrostatic mapping factor
191         function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch, elevation);
192 
193         // Wet mapping factor
194         function[1] = TroposphericModelUtils.mappingFunction(awFunction.value(absLatidude), bwFunction.value(absLatidude), cwFunction.value(absLatidude), elevation);
195 
196         // Apply height correction
197         final double correction = TroposphericModelUtils.computeHeightCorrection(elevation, point.getAltitude());
198         function[0] = function[0] + correction;
199 
200         return function;
201     }
202 
203     /** {@inheritDoc} */
204     @Override
205     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final T elevation, final FieldGeodeticPoint<T> point,
206                                                               final FieldAbsoluteDate<T> date) {
207         final Field<T> field = date.getField();
208         final T zero = field.getZero();
209 
210         // Day of year computation
211         final DateTimeComponents dtc = date.getComponents(utc);
212         final int dofyear = dtc.getDate().getDayOfYear();
213 
214         // Temporal factor
215         double t0 = 28;
216         if (point.getLatitude().getReal() < 0) {
217             // southern hemisphere: t0 = 28 + an integer half of year
218             t0 += 183;
219         }
220         final T coef    = zero.getPi().multiply(2.0).multiply((dofyear - t0) / 365.25);
221         final T cosCoef = FastMath.cos(coef);
222 
223         // Compute ah, bh and ch Eq. 5
224         double absLatidude = FastMath.abs(point.getLatitude().getReal());
225         // there are no data in the model for latitudes lower than 15°
226         absLatidude = FastMath.max(FastMath.toRadians(15.0), absLatidude);
227         // there are no data in the model for latitudes greater than 75°
228         absLatidude = FastMath.min(FastMath.toRadians(75.0), absLatidude);
229         final T ah = cosCoef.multiply(ahAmplitudeFunction.value(absLatidude)).negate().add(ahAverageFunction.value(absLatidude));
230         final T bh = cosCoef.multiply(bhAmplitudeFunction.value(absLatidude)).negate().add(bhAverageFunction.value(absLatidude));
231         final T ch = cosCoef.multiply(chAmplitudeFunction.value(absLatidude)).negate().add(chAverageFunction.value(absLatidude));
232 
233         final T[] function = MathArrays.buildArray(field, 2);
234 
235         // Hydrostatic mapping factor
236         function[0] = TroposphericModelUtils.mappingFunction(ah, bh, ch, elevation);
237 
238         // Wet mapping factor
239         function[1] = TroposphericModelUtils.mappingFunction(zero.add(awFunction.value(absLatidude)), zero.add(bwFunction.value(absLatidude)),
240                                                              zero.add(cwFunction.value(absLatidude)), elevation);
241 
242         // Apply height correction
243         final T correction = TroposphericModelUtils.computeHeightCorrection(elevation, point.getAltitude(), field);
244         function[0] = function[0].add(correction);
245 
246         return function;
247     }
248 
249 }