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.iturp834;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
21  import org.hipparchus.util.FastMath;
22  import org.hipparchus.util.FieldSinCos;
23  import org.hipparchus.util.MathArrays;
24  import org.hipparchus.util.MathUtils;
25  import org.hipparchus.util.SinCos;
26  import org.orekit.bodies.FieldGeodeticPoint;
27  import org.orekit.bodies.GeodeticPoint;
28  import org.orekit.errors.OrekitException;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.models.earth.troposphere.TroposphereMappingFunction;
31  import org.orekit.time.AbsoluteDate;
32  import org.orekit.time.FieldAbsoluteDate;
33  import org.orekit.time.TimeScale;
34  import org.orekit.utils.FieldTrackingCoordinates;
35  import org.orekit.utils.TrackingCoordinates;
36  
37  import java.io.BufferedReader;
38  import java.io.IOException;
39  import java.io.InputStream;
40  import java.io.InputStreamReader;
41  import java.nio.charset.StandardCharsets;
42  import java.util.regex.Pattern;
43  
44  /** ITU-R P.834 mapping function.
45   * @see ITURP834PathDelay
46   * @see ITURP834WeatherParametersProvider
47   * @author Luc Maisonobe
48   * @see <a href="https://www.itu.int/rec/R-REC-P.834/en">P.834 : Effects of tropospheric refraction on radiowave propagation</a>
49   * @since 13.0
50   */
51  public class ITURP834MappingFunction implements TroposphereMappingFunction {
52  
53      /** Splitter for fields in lines. */
54      private static final Pattern SPLITTER = Pattern.compile("\\s+");
55  
56      /** Name of data file. */
57      private static final String MAPPING_FUNCTION_NAME = "/assets/org/orekit/ITU-R-P.834/p834_mf_coeff_v1.txt";
58  
59      /** Minimum latitude, including extra margin for dealing with boundaries (degrees). */
60      private static final double MIN_LAT = -92.5;
61  
62      /** Maximum latitude, including extra margin for dealing with boundaries (degrees). */
63      private static final double MAX_LAT =  92.5;
64  
65      /** Grid step in latitude (degrees). */
66      private static final double STEP_LAT =  5.0;
67  
68      /** Minimum longitude, including extra margin for dealing with boundaries (degrees). */
69      private static final double MIN_LON = -182.5;
70  
71      /** Maximum longitude, including extra margin for dealing with boundaries (degrees). */
72      private static final double MAX_LON =  182.5;
73  
74      /** Grid step in longitude (degrees). */
75      private static final double STEP_LON =  5.0;
76  
77      /** Second coefficient for hydrostatic component. */
78      private static final double BH = 0.0029;
79  
80      /** Constants for third coefficient for hydrostatic component, Northern hemisphere. */
81      private static final double[] CH_NORTH = { 0.062, 0.001, 0.005, 0.0 };
82  
83      /** Constants for third coefficient for hydrostatic component, Southern hemisphere. */
84      private static final double[] CH_SOUTH = { 0.062, 0.002, 0.007, FastMath.PI };
85  
86      /** Reference day of year for third coefficient for hydrostatic component. */
87      private static final double REF_DOY = 28;
88  
89      /** Year length (in days). */
90      private static final double YEAR = 365.25;
91  
92      /** Second coefficient for wet component. */
93      private static final double BW = 0.00146;
94  
95      /** Third coefficient for wet component. */
96      private static final double CW = 0.04391;
97  
98      /** Global factor to apply to first coefficient. */
99      private static final double FACTOR = 1.0e-3;
100 
101     /** Interpolator for constant hydrostatic coefficient. */
102     private static final BilinearInterpolatingFunction A0H;
103 
104     /** Interpolator for annual cosine hydrostatic coefficient. */
105     private static final BilinearInterpolatingFunction A1H;
106 
107     /** Interpolator for annual sine hydrostatic coefficient. */
108     private static final BilinearInterpolatingFunction B1H;
109 
110     /** Interpolator for semi-annual cosine hydrostatic coefficient. */
111     private static final BilinearInterpolatingFunction A2H;
112 
113     /** Interpolator for semin-annual sine hydrostatic coefficient. */
114     private static final BilinearInterpolatingFunction B2H;
115 
116     /** Interpolator for constant wet coefficient. */
117     private static final BilinearInterpolatingFunction A0W;
118 
119     /** Interpolator for annual cosine wet coefficient. */
120     private static final BilinearInterpolatingFunction A1W;
121 
122     /** Interpolator for annual sine wet coefficient. */
123     private static final BilinearInterpolatingFunction B1W;
124 
125     /** Interpolator for semi-annual cosine wet coefficient. */
126     private static final BilinearInterpolatingFunction A2W;
127 
128     /** Interpolator for semin-annual sine wet coefficient. */
129     private static final BilinearInterpolatingFunction B2W;
130 
131     // load model data
132     static {
133 
134         // create the various tables
135         // we add extra lines and columns to the official files, for dealing with boundaries
136         final double[]   longitudes = interpolationPoints(MIN_LON, MAX_LON, STEP_LON);
137         final double[]   latitudes  = interpolationPoints(MIN_LAT, MAX_LAT, STEP_LAT);
138         final double[][] a0h        = new double[longitudes.length][latitudes.length];
139         final double[][] a1h        = new double[longitudes.length][latitudes.length];
140         final double[][] b1h        = new double[longitudes.length][latitudes.length];
141         final double[][] a2h        = new double[longitudes.length][latitudes.length];
142         final double[][] b2h        = new double[longitudes.length][latitudes.length];
143         final double[][] a0w        = new double[longitudes.length][latitudes.length];
144         final double[][] a1w        = new double[longitudes.length][latitudes.length];
145         final double[][] b1w        = new double[longitudes.length][latitudes.length];
146         final double[][] a2w        = new double[longitudes.length][latitudes.length];
147         final double[][] b2w        = new double[longitudes.length][latitudes.length];
148 
149         try (InputStream       is     = ITURP834MappingFunction.class.getResourceAsStream(MAPPING_FUNCTION_NAME);
150              InputStreamReader isr    = is  == null ? null : new InputStreamReader(is, StandardCharsets.UTF_8);
151              BufferedReader    reader = isr == null ? null : new BufferedReader(isr)) {
152             if (reader == null) {
153                 // this should never happen with embedded data
154                 throw new OrekitException(OrekitMessages.UNABLE_TO_FIND_FILE, MAPPING_FUNCTION_NAME);
155             }
156             int lineNumber = 0;
157             for (String line = reader.readLine(); line != null; line = reader.readLine()) {
158                 ++lineNumber;
159                 final String[] fields = SPLITTER.split(line.trim());
160                 if (fields.length != 12) {
161                     // this should never happen with the embedded data
162                     throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
163                                               lineNumber, MAPPING_FUNCTION_NAME, line);
164                 }
165 
166                 // parse the fields
167                 final double[] numericFields = new double[fields.length];
168                 for (int i = 0; i < fields.length; ++i) {
169                     numericFields[i] = Double.parseDouble(fields[i]);
170                 }
171 
172                 // find indices in our extended grid
173                 final int longitudeIndex = (int) FastMath.rint((numericFields[1] - MIN_LON) / STEP_LON);
174                 final int latitudeIndex  = (int) FastMath.rint((numericFields[0] - MIN_LAT) / STEP_LAT);
175 
176                 // fill-in tables
177                 a0h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 2];
178                 a1h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 3];
179                 b1h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 4];
180                 a2h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 5];
181                 b2h[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 6];
182                 a0w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 7];
183                 a1w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 8];
184                 b1w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[ 9];
185                 a2w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[10];
186                 b2w[longitudeIndex][latitudeIndex] = FACTOR * numericFields[11];
187 
188             }
189 
190             // extend tables in latitude to cover poles
191             for (int i = 1; i < longitudes.length - 1; ++i) {
192                 a0h[i][0]                    = a0h[i][1];
193                 a0h[i][latitudes.length - 1] = a0h[i][latitudes.length - 2];
194                 a1h[i][0]                    = a1h[i][1];
195                 a1h[i][latitudes.length - 1] = a1h[i][latitudes.length - 2];
196                 b1h[i][0]                    = b1h[i][1];
197                 b1h[i][latitudes.length - 1] = b1h[i][latitudes.length - 2];
198                 a2h[i][0]                    = a2h[i][1];
199                 a2h[i][latitudes.length - 1] = a2h[i][latitudes.length - 2];
200                 b2h[i][0]                    = b2h[i][1];
201                 b2h[i][latitudes.length - 1] = b2h[i][latitudes.length - 2];
202                 a0w[i][0]                    = a0w[i][1];
203                 a0w[i][latitudes.length - 1] = a0w[i][latitudes.length - 2];
204                 a1w[i][0]                    = a1w[i][1];
205                 a1w[i][latitudes.length - 1] = a1w[i][latitudes.length - 2];
206                 b1w[i][0]                    = b1w[i][1];
207                 b1w[i][latitudes.length - 1] = b1w[i][latitudes.length - 2];
208                 a2w[i][0]                    = a2w[i][1];
209                 a2w[i][latitudes.length - 1] = a2w[i][latitudes.length - 2];
210                 b2w[i][0]                    = b2w[i][1];
211                 b2w[i][latitudes.length - 1] = b2w[i][latitudes.length - 2];
212             }
213 
214             // extend tables in longitude to cover anti-meridian
215             for (int j = 0; j < latitudes.length; ++j) {
216                 a0h[0][j]                     = a0h[longitudes.length - 2][j];
217                 a0h[longitudes.length - 1][j] = a0h[1][j];
218                 a1h[0][j]                     = a1h[longitudes.length - 2][j];
219                 a1h[longitudes.length - 1][j] = a1h[1][j];
220                 b1h[0][j]                     = b1h[longitudes.length - 2][j];
221                 b1h[longitudes.length - 1][j] = b1h[1][j];
222                 a2h[0][j]                     = a2h[longitudes.length - 2][j];
223                 a2h[longitudes.length - 1][j] = a2h[1][j];
224                 b2h[0][j]                     = b2h[longitudes.length - 2][j];
225                 b2h[longitudes.length - 1][j] = b2h[1][j];
226                 a0w[0][j]                     = a0w[longitudes.length - 2][j];
227                 a0w[longitudes.length - 1][j] = a0w[1][j];
228                 a1w[0][j]                     = a1w[longitudes.length - 2][j];
229                 a1w[longitudes.length - 1][j] = a1w[1][j];
230                 b1w[0][j]                     = b1w[longitudes.length - 2][j];
231                 b1w[longitudes.length - 1][j] = b1w[1][j];
232                 a2w[0][j]                     = a2w[longitudes.length - 2][j];
233                 a2w[longitudes.length - 1][j] = a2w[1][j];
234                 b2w[0][j]                     = b2w[longitudes.length - 2][j];
235                 b2w[longitudes.length - 1][j] = b2w[1][j];
236             }
237 
238             // build interpolators
239             A0H = new BilinearInterpolatingFunction(longitudes, latitudes, a0h);
240             A1H = new BilinearInterpolatingFunction(longitudes, latitudes, a1h);
241             B1H = new BilinearInterpolatingFunction(longitudes, latitudes, b1h);
242             A2H = new BilinearInterpolatingFunction(longitudes, latitudes, a2h);
243             B2H = new BilinearInterpolatingFunction(longitudes, latitudes, b2h);
244             A0W = new BilinearInterpolatingFunction(longitudes, latitudes, a0w);
245             A1W = new BilinearInterpolatingFunction(longitudes, latitudes, a1w);
246             B1W = new BilinearInterpolatingFunction(longitudes, latitudes, b1w);
247             A2W = new BilinearInterpolatingFunction(longitudes, latitudes, a2w);
248             B2W = new BilinearInterpolatingFunction(longitudes, latitudes, b2w);
249 
250         } catch (IOException ioe) {
251             // this should never happen with the embedded data
252             throw new OrekitException(OrekitMessages.INTERNAL_ERROR, ioe);
253         }
254     }
255 
256     /** UTC time scale. */
257     private final TimeScale utc;
258 
259     /** Simple constructor.
260      * @param utc UTC time scale
261      */
262     public ITURP834MappingFunction(final TimeScale utc) {
263         this.utc = utc;
264     }
265 
266     /** Create interpolation points coordinates.
267      * @param min min angle in degrees (included)
268      * @param max max angle in degrees (included)
269      * @param step step between points
270      * @return interpolation points coordinates (in radians)
271      */
272     private static double[] interpolationPoints(final double min, final double max, final double step) {
273         final double[] points = new double[(int) FastMath.rint((max - min) / step) + 1];
274         for (int i = 0; i < points.length; i++) {
275             points[i] = FastMath.toRadians(min + i * step);
276         }
277         return points;
278     }
279 
280     @Override
281     public double[] mappingFactors(final TrackingCoordinates trackingCoordinates, final GeodeticPoint point,
282                                    final AbsoluteDate date) {
283 
284         final double doy = date.getDayOfYear(utc);
285 
286         // compute third dry coefficient
287         // equation 28c in ITU-R P.834 recommendation
288         final double[] c    = point.getLatitude() >= 0.0 ? CH_NORTH : CH_SOUTH;
289         final double   cosL = FastMath.cos(point.getLatitude());
290         final double   cosD = FastMath.cos((doy - REF_DOY) * MathUtils.TWO_PI / YEAR + c[3]);
291         final double   ch   = c[0] + ((cosD + 1) * c[2] * 0.5 + c[1]) * (1 - cosL);
292 
293         // compute first coefficient
294         final SinCos sc1 = FastMath.sinCos(MathUtils.TWO_PI * doy / YEAR);
295         final SinCos sc2 = SinCos.sum(sc1, sc1);
296 
297         // equation 28d in ITU-R P.834 recommendation
298         final double ah  = A0H.value(point.getLongitude(), point.getLatitude()) +
299                            A1H.value(point.getLongitude(), point.getLatitude()) * sc1.cos() +
300                            B1H.value(point.getLongitude(), point.getLatitude()) * sc1.sin() +
301                            A2H.value(point.getLongitude(), point.getLatitude()) * sc2.cos() +
302                            B2H.value(point.getLongitude(), point.getLatitude()) * sc2.sin();
303 
304         // equation 28e in ITU-R P.834 recommendation
305         final double aw  = A0W.value(point.getLongitude(), point.getLatitude()) +
306                            A1W.value(point.getLongitude(), point.getLatitude()) * sc1.cos() +
307                            B1W.value(point.getLongitude(), point.getLatitude()) * sc1.sin() +
308                            A2W.value(point.getLongitude(), point.getLatitude()) * sc2.cos() +
309                            B2W.value(point.getLongitude(), point.getLatitude()) * sc2.sin();
310 
311         // mapping function, equations 28a and 28b in ITU-R P.834 recommendation
312         final double sinTheta = FastMath.sin(trackingCoordinates.getElevation());
313         final double mh = (1 + ah / (1 + BH / (1 + ch))) /
314                           (sinTheta + ah / (sinTheta + BH / (sinTheta + ch)));
315         final double mw = (1 + aw / (1 + BW / (1 + CW))) /
316                           (sinTheta + aw / (sinTheta + BW / (sinTheta + CW)));
317 
318         return new double[] {
319             mh, mw
320         };
321 
322     }
323 
324     @Override
325     public <T extends CalculusFieldElement<T>> T[] mappingFactors(final FieldTrackingCoordinates<T> trackingCoordinates,
326                                                                   final FieldGeodeticPoint<T> point,
327                                                                   final FieldAbsoluteDate<T> date) {
328 
329         final T doy = date.getDayOfYear(utc);
330 
331         // compute third dry coefficient
332         // equation 28c in ITU-R P.834 recommendation
333         final double[] c    = point.getLatitude().getReal() >= 0.0 ? CH_NORTH : CH_SOUTH;
334         final T        cosL = FastMath.cos(point.getLatitude());
335         final T        cosD = FastMath.cos(doy.subtract(REF_DOY).multiply(MathUtils.TWO_PI / YEAR).add(c[3]));
336         final T        ch   = cosD.add(1).multiply(c[2] * 0.5).add(c[1]).multiply(cosL.subtract(1).negate()).add(c[0]);
337 
338         // compute first coefficient
339         final FieldSinCos<T> sc1 = FastMath.sinCos(doy.multiply(MathUtils.TWO_PI / YEAR));
340         final FieldSinCos<T> sc2 = FieldSinCos.sum(sc1, sc1);
341 
342         // equation 28d in ITU-R P.834 recommendation
343         final T ah  =     A0H.value(point.getLongitude(), point.getLatitude()).
344                       add(A1H.value(point.getLongitude(), point.getLatitude()).multiply(sc1.cos())).
345                       add(B1H.value(point.getLongitude(), point.getLatitude()).multiply(sc1.sin())).
346                       add(A2H.value(point.getLongitude(), point.getLatitude()).multiply(sc2.cos())).
347                       add(B2H.value(point.getLongitude(), point.getLatitude()).multiply(sc2.sin()));
348 
349         // equation 28e in ITU-R P.834 recommendation
350         final T aw  =     A0W.value(point.getLongitude(), point.getLatitude()).
351                       add(A1W.value(point.getLongitude(), point.getLatitude()).multiply(sc1.cos())).
352                       add(B1W.value(point.getLongitude(), point.getLatitude()).multiply(sc1.sin())).
353                       add(A2W.value(point.getLongitude(), point.getLatitude()).multiply(sc2.cos())).
354                       add(B2W.value(point.getLongitude(), point.getLatitude()).multiply(sc2.sin()));
355 
356         // mapping function, equations 28a and 28b in ITU-R P.834 recommendation
357         final T sinTheta = FastMath.sin(trackingCoordinates.getElevation());
358         final T mh = ch.add(1).reciprocal().multiply(BH).add(1).reciprocal().multiply(ah).add(1).
359                      divide(ch.add(sinTheta).reciprocal().multiply(BH).add(sinTheta).reciprocal().multiply(ah).add(sinTheta));
360         final T mw = aw.divide(1 + BW / (1 + CW)).add(1).
361                      divide(sinTheta.add(CW).reciprocal().multiply(BW).add(sinTheta).reciprocal().multiply(aw).add(sinTheta));
362 
363         final T[] m = MathArrays.buildArray(date.getField(), 2);
364         m[0] = mh;
365         m[1] = mw;
366         return m;
367 
368     }
369 
370 }