1 /*
2  * Copyright (C) 2017 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *      http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 package android.location.cts.gnss.pseudorange;
18 
19 /**
20  * Calculate the troposheric delay based on the ENGOS Tropospheric model.
21  *
22  * <p>The tropospheric delay is modeled as a combined effect of the delay experienced due to
23  * hyrostatic (dry) and wet components of the troposphere. Both delays experienced at zenith are
24  * scaled with a mapping function to get the delay at any specific elevation.
25  *
26  * <p>The tropospheric model algorithm of EGNOS model by Penna, N., A. Dodson and W. Chen (2001)
27  * (http://espace.library.curtin.edu.au/cgi-bin/espace.pdf?file=/2008/11/13/file_1/18917) is used
28  * for calculating the zenith delays. In this model, the weather parameters are extracted using
29  * interpolation from lookup table derived from the US Standard Atmospheric Supplements, 1966.
30  *
31  * <p>A close form mapping function is built using Guo and Langley, 2003
32  * (http://gauss2.gge.unb.ca/papers.pdf/iongpsgnss2003.guo.pdf) which is able to calculate accurate
33  * mapping down to 2 degree elevations.
34  *
35  * <p>Sources:
36  * <p>http://espace.library.curtin.edu.au/cgi-bin/espace.pdf?file=/2008/11/13/file_1/18917
37  * <p>- http://www.academia.edu/3512180/Assessment_of_UNB3M_neutral
38  * _atmosphere_model_and_EGNOS_model_for_near-equatorial-tropospheric_delay_correction
39  * <p>- http://gauss.gge.unb.ca/papers.pdf/ion52am.collins.pdf
40  * <p>- http://www.navipedia.net/index.php/Tropospheric_Delay#cite_ref-3
41  * <p>Hydrostatic and non-hydrostatic mapping functions are obtained from:
42  * http://gauss2.gge.unb.ca/papers.pdf/iongpsgnss2003.guo.pdf
43  *
44  */
45 public class TroposphericModelEgnos {
46   // parameters of the EGNOS models
47   private static final int INDEX_15_DEGREES = 0;
48   private static final int INDEX_75_DEGREES = 4;
49   private static final int LATITUDE_15_DEGREES = 15;
50   private static final int LATITUDE_75_DEGREES = 75;
51   // Lookup Average parameters
52   // Troposphere average presssure mbar
53   private static final double[] latDegreeToPressureMbarAvgMap =
54     {1013.25,  1017.25, 1015.75, 1011.75, 1013.0};
55   // Troposphere average temperature Kelvin
56   private static final double[] latDegreeToTempKelvinAvgMap =
57     {299.65, 294.15, 283.15, 272.15, 263.65};
58   // Troposphere average wator vapor pressure
59   private static final double[] latDegreeToWVPressureMbarAvgMap = {26.31, 21.79, 11.66, 6.78, 4.11};
60   // Troposphere average temperature lapse rate K/m
61   private static final double[] latDegreeToBetaAvgMapKPM =
62     {6.30e-3, 6.05e-3, 5.58e-3, 5.39e-3, 4.53e-3};
63   // Troposphere average water vapor lapse rate (dimensionless)
64   private static final double[] latDegreeToLampdaAvgMap = {2.77, 3.15, 2.57, 1.81, 1.55};
65 
66   // Lookup Amplitude parameters
67   // Troposphere amplitude presssure mbar
68   private static final double[] latDegreeToPressureMbarAmpMap = {0.0, -3.75, -2.25, -1.75, -0.5};
69   // Troposphere amplitude temperature Kelvin
70   private static final double[] latDegreeToTempKelvinAmpMap = {0.0, 7.0, 11.0, 15.0, 14.5};
71   // Troposphere amplitude wator vapor pressure
72   private static final double[] latDegreeToWVPressureMbarAmpMap = {0.0, 8.85, 7.24, 5.36, 3.39};
73   // Troposphere amplitude temperature lapse rate K/m
74   private static final double[] latDegreeToBetaAmpMapKPM =
75     {0.0, 0.25e-3, 0.32e-3, 0.81e-3, 0.62e-3};
76   // Troposphere amplitude water vapor lapse rate (dimensionless)
77   private static final double[] latDegreeToLampdaAmpMap = {0.0, 0.33, 0.46, 0.74, 0.30};
78   // Zenith delay dry constant K/mbar
79   private static final double K1 = 77.604;
80   // Zenith delay wet constant K^2/mbar
81   private static final double K2 = 382000.0;
82   // gas constant for dry air J/kg/K
83   private static final double RD = 287.054;
84   // Acceleration of gravity at the atmospheric column centroid m/s^-2
85   private static final double GM = 9.784;
86   // Gravity m/s^2
87   private static final double GRAVITY_MPS2 = 9.80665;
88 
89   private static final double MINIMUM_INTERPOLATION_THRESHOLD = 1e-25;
90   private static final double B_HYDROSTATIC = 0.0035716;
91   private static final double C_HYDROSTATIC = 0.082456;
92   private static final double B_NON_HYDROSTATIC = 0.0018576;
93   private static final double C_NON_HYDROSTATIC = 0.062741;
94   private static final double SOUTHERN_HEMISPHERE_DMIN = 211.0;
95   private static final double NORTHERN_HEMISPHERE_DMIN = 28.0;
96   // Days recalling that every fourth year is a leap year and has an extra day - February 29th
97   private static final double DAYS_PER_YEAR = 365.25;
98 
99   /**
100    * Compute the tropospheric correction in meters given the satellite elevation in radians, the
101    * user latitude in radians, the user Orthometric height above sea level in meters and the day of
102    * the year.
103    *
104    * <p>Dry and wet delay zenith delay components are calculated and then scaled with the mapping
105    * function at the given satellite elevation.
106    *
107    */
calculateTropoCorrectionMeters(double satElevationRadians, double userLatitudeRadian, double heightMetersAboveSeaLevel, int dayOfYear1To366)108   public static double calculateTropoCorrectionMeters(double satElevationRadians,
109       double userLatitudeRadian, double heightMetersAboveSeaLevel, int dayOfYear1To366) {
110     DryAndWetMappingValues dryAndWetMappingValues =
111         computeDryAndWetMappingValuesUsingUNBabcMappingFunction(satElevationRadians,
112             userLatitudeRadian, heightMetersAboveSeaLevel);
113     DryAndWetZenithDelays dryAndWetZenithDelays = calculateZenithDryAndWetDelaysSec
114         (userLatitudeRadian, heightMetersAboveSeaLevel, dayOfYear1To366);
115 
116     double drydelaySeconds =
117         dryAndWetZenithDelays.dryZenithDelaySec * dryAndWetMappingValues.dryMappingValue;
118     double wetdelaySeconds =
119         dryAndWetZenithDelays.wetZenithDelaySec * dryAndWetMappingValues.wetMappingValue;
120     return drydelaySeconds + wetdelaySeconds;
121   }
122 
123   /**
124    * Compute the dry and wet mapping values based on the University of Brunswick UNBabc model. The
125    * mapping function inputs are satellite elevation in radians, user latitude in radians and user
126    * orthometric height above sea level in meters. The function returns
127    * {@code DryAndWetMappingValues} containing dry and wet mapping values.
128    *
129    * <p>From the many dry and wet mapping functions of components of the troposphere, the method
130    * from the University of Brunswick in Canada was selected due to its reasonable computation time
131    * and accuracy with satellites as low as 2 degrees elevation.
132    * <p>Source: http://gauss2.gge.unb.ca/papers.pdf/iongpsgnss2003.guo.pdf
133    */
computeDryAndWetMappingValuesUsingUNBabcMappingFunction( double satElevationRadians, double userLatitudeRadians, double heightMetersAboveSeaLevel)134   private static DryAndWetMappingValues computeDryAndWetMappingValuesUsingUNBabcMappingFunction(
135       double satElevationRadians, double userLatitudeRadians, double heightMetersAboveSeaLevel) {
136 
137     if (satElevationRadians > Math.PI / 2.0) {
138       satElevationRadians = Math.PI / 2.0;
139     } else if (satElevationRadians < 2.0 * Math.PI / 180.0) {
140       satElevationRadians = Math.toRadians(2.0);
141     }
142 
143     // dry components mapping parameters
144     double aHidrostatic = (1.18972 - 0.026855 * heightMetersAboveSeaLevel / 1000.0 + 0.10664
145         * Math.cos(userLatitudeRadians)) / 1000.0;
146 
147 
148     double numeratorDry = 1.0 + (aHidrostatic / (1.0 + (B_HYDROSTATIC / (1.0 + C_HYDROSTATIC))));
149     double denominatorDry = Math.sin(satElevationRadians) + (aHidrostatic / (
150         Math.sin(satElevationRadians)
151         + (B_HYDROSTATIC / (Math.sin(satElevationRadians) + C_HYDROSTATIC))));
152 
153     double drymap = numeratorDry / denominatorDry;
154 
155     // wet components mapping parameters
156     double aNonHydrostatic = (0.61120 - 0.035348 * heightMetersAboveSeaLevel / 1000.0 - 0.01526
157         * Math.cos(userLatitudeRadians)) / 1000.0;
158 
159 
160     double numeratorWet =
161         1.0 + (aNonHydrostatic / (1.0 + (B_NON_HYDROSTATIC / (1.0 + C_NON_HYDROSTATIC))));
162     double denominatorWet = Math.sin(satElevationRadians) + (aNonHydrostatic / (
163         Math.sin(satElevationRadians)
164         + (B_NON_HYDROSTATIC / (Math.sin(satElevationRadians) + C_NON_HYDROSTATIC))));
165 
166     double wetmap = numeratorWet / denominatorWet;
167     return new DryAndWetMappingValues(drymap, wetmap);
168   }
169 
170   /**
171    * Compute the combined effect of the delay at zenith experienced due to hyrostatic (dry) and wet
172    * components of the troposphere. The function inputs are the user latitude in radians, user
173    * orthometric height above sea level in meters and the day of the year (1-366). The function
174    * returns a {@code DryAndWetZenithDelays} containing dry and wet delays at zenith.
175    *
176    * <p>EGNOS Tropospheric model by Penna et al. (2001) is used in this case.
177    * (http://espace.library.curtin.edu.au/cgi-bin/espace.pdf?file=/2008/11/13/file_1/18917)
178    *
179    */
calculateZenithDryAndWetDelaysSec(double userLatitudeRadians, double heightMetersAboveSeaLevel, int dayOfyear1To366)180   private static DryAndWetZenithDelays calculateZenithDryAndWetDelaysSec(double userLatitudeRadians,
181       double heightMetersAboveSeaLevel, int dayOfyear1To366) {
182     // interpolated meteorological values
183     double pressureMbar;
184     double tempKelvin;
185     double waterVaporPressureMbar;
186     // temperature lapse rate, [K/m]
187     double beta;
188     // water vapor lapse rate, dimensionless
189     double lambda;
190 
191     double absLatitudeDeg = Math.toDegrees(Math.abs(userLatitudeRadians));
192     // day of year min constant
193     double dmin;
194     if (userLatitudeRadians < 0) {
195       dmin = SOUTHERN_HEMISPHERE_DMIN;
196     } else {
197       dmin = NORTHERN_HEMISPHERE_DMIN;
198 
199     }
200     double amplitudeScalefactor = Math.cos((2 * Math.PI * (dayOfyear1To366 - dmin))
201         / DAYS_PER_YEAR);
202 
203     if (absLatitudeDeg <= LATITUDE_15_DEGREES) {
204       pressureMbar = latDegreeToPressureMbarAvgMap[INDEX_15_DEGREES]
205           - latDegreeToPressureMbarAmpMap[INDEX_15_DEGREES] * amplitudeScalefactor;
206       tempKelvin = latDegreeToTempKelvinAvgMap[INDEX_15_DEGREES]
207           - latDegreeToTempKelvinAmpMap[INDEX_15_DEGREES] * amplitudeScalefactor;
208       waterVaporPressureMbar = latDegreeToWVPressureMbarAvgMap[INDEX_15_DEGREES]
209           - latDegreeToWVPressureMbarAmpMap[INDEX_15_DEGREES] * amplitudeScalefactor;
210       beta = latDegreeToBetaAvgMapKPM[INDEX_15_DEGREES] - latDegreeToBetaAmpMapKPM[INDEX_15_DEGREES]
211           * amplitudeScalefactor;
212       lambda = latDegreeToLampdaAmpMap[INDEX_15_DEGREES] - latDegreeToLampdaAmpMap[INDEX_15_DEGREES]
213           * amplitudeScalefactor;
214     } else if (absLatitudeDeg > LATITUDE_15_DEGREES && absLatitudeDeg < LATITUDE_75_DEGREES) {
215       int key = (int) (absLatitudeDeg / LATITUDE_15_DEGREES);
216 
217       double averagePressureMbar = interpolate(key * LATITUDE_15_DEGREES,
218           latDegreeToPressureMbarAvgMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
219           latDegreeToPressureMbarAvgMap[key], absLatitudeDeg);
220       double amplitudePressureMbar = interpolate(key * LATITUDE_15_DEGREES,
221           latDegreeToPressureMbarAmpMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
222           latDegreeToPressureMbarAmpMap[key], absLatitudeDeg);
223       pressureMbar = averagePressureMbar - amplitudePressureMbar * amplitudeScalefactor;
224 
225       double averageTempKelvin = interpolate(key * LATITUDE_15_DEGREES,
226           latDegreeToTempKelvinAvgMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
227           latDegreeToTempKelvinAvgMap[key], absLatitudeDeg);
228       double amplitudeTempKelvin = interpolate(key * LATITUDE_15_DEGREES,
229           latDegreeToTempKelvinAmpMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
230           latDegreeToTempKelvinAmpMap[key], absLatitudeDeg);
231       tempKelvin = averageTempKelvin - amplitudeTempKelvin * amplitudeScalefactor;
232 
233       double averageWaterVaporPressureMbar = interpolate(key * LATITUDE_15_DEGREES,
234           latDegreeToWVPressureMbarAvgMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
235           latDegreeToWVPressureMbarAvgMap[key], absLatitudeDeg);
236       double amplitudeWaterVaporPressureMbar = interpolate(key * LATITUDE_15_DEGREES,
237           latDegreeToWVPressureMbarAmpMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
238           latDegreeToWVPressureMbarAmpMap[key], absLatitudeDeg);
239       waterVaporPressureMbar =
240           averageWaterVaporPressureMbar - amplitudeWaterVaporPressureMbar * amplitudeScalefactor;
241 
242       double averageBeta = interpolate(key * LATITUDE_15_DEGREES, latDegreeToBetaAvgMapKPM[key - 1],
243           (key + 1) * LATITUDE_15_DEGREES, latDegreeToBetaAvgMapKPM[key], absLatitudeDeg);
244       double amplitudeBeta = interpolate(key * LATITUDE_15_DEGREES,
245           latDegreeToBetaAmpMapKPM[key - 1], (key + 1) * LATITUDE_15_DEGREES,
246           latDegreeToBetaAmpMapKPM[key], absLatitudeDeg);
247       beta = averageBeta - amplitudeBeta * amplitudeScalefactor;
248 
249       double averageLambda = interpolate(key * LATITUDE_15_DEGREES,
250           latDegreeToLampdaAvgMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
251           latDegreeToLampdaAvgMap[key], absLatitudeDeg);
252       double amplitudeLambda = interpolate(key * LATITUDE_15_DEGREES,
253           latDegreeToLampdaAmpMap[key - 1], (key + 1) * LATITUDE_15_DEGREES,
254           latDegreeToLampdaAmpMap[key], absLatitudeDeg);
255       lambda = averageLambda - amplitudeLambda * amplitudeScalefactor;
256     } else {
257       pressureMbar = latDegreeToPressureMbarAvgMap[INDEX_75_DEGREES]
258           - latDegreeToPressureMbarAmpMap[INDEX_75_DEGREES] * amplitudeScalefactor;
259       tempKelvin = latDegreeToTempKelvinAvgMap[INDEX_75_DEGREES]
260           - latDegreeToTempKelvinAmpMap[INDEX_75_DEGREES] * amplitudeScalefactor;
261       waterVaporPressureMbar = latDegreeToWVPressureMbarAvgMap[INDEX_75_DEGREES]
262           - latDegreeToWVPressureMbarAmpMap[INDEX_75_DEGREES] * amplitudeScalefactor;
263       beta = latDegreeToBetaAvgMapKPM[INDEX_75_DEGREES] - latDegreeToBetaAmpMapKPM[INDEX_75_DEGREES]
264           * amplitudeScalefactor;
265       lambda = latDegreeToLampdaAmpMap[INDEX_75_DEGREES] - latDegreeToLampdaAmpMap[INDEX_75_DEGREES]
266           * amplitudeScalefactor;
267     }
268 
269     double zenithDryDelayAtSeaLevelSeconds = (1.0e-6 * K1 * RD * pressureMbar) / GM;
270     double zenithWetDelayAtSeaLevelSeconds = (((1.0e-6 * K2 * RD)
271         / (GM * (lambda + 1.0) - beta * RD)) * (waterVaporPressureMbar / tempKelvin));
272     double commonBase = 1.0 - ((beta * heightMetersAboveSeaLevel) / tempKelvin);
273 
274     double powerDry = (GRAVITY_MPS2 / (RD * beta));
275     double powerWet = (((lambda + 1.0) * GRAVITY_MPS2) / (RD * beta)) - 1.0;
276     double zenithDryDelaySeconds = zenithDryDelayAtSeaLevelSeconds * Math.pow(commonBase, powerDry);
277     double zenithWetDelaySeconds = zenithWetDelayAtSeaLevelSeconds * Math.pow(commonBase, powerWet);
278     return new DryAndWetZenithDelays(zenithDryDelaySeconds, zenithWetDelaySeconds);
279   }
280 
281   /**
282    * Interpolate linearly given two points (point1X, point1Y) and (point2X, point2Y). Given the
283    * desired value of x (xInterpolated), an interpolated value of y shall be computed and returned.
284    */
interpolate(double point1X, double point1Y, double point2X, double point2Y, double xOutput)285   private static double interpolate(double point1X, double point1Y, double point2X, double point2Y,
286       double xOutput) {
287     // Check that xOutput is between the two interpolation points.
288     if ((point1X < point2X && (xOutput < point1X || xOutput > point2X))
289         || (point2X < point1X && (xOutput < point2X || xOutput > point1X))) {
290       throw new IllegalArgumentException("Interpolated value is outside the interpolated region");
291     }
292     double deltaX = point2X - point1X;
293     double yOutput;
294 
295     if (Math.abs(deltaX) > MINIMUM_INTERPOLATION_THRESHOLD) {
296       yOutput = point1Y + (xOutput - point1X) / deltaX * (point2Y - point1Y);
297     } else {
298       yOutput = point1Y;
299     }
300     return yOutput;
301   }
302 
303   /* A class containing dry and wet mapping values */
304   private static class DryAndWetMappingValues {
305     public double dryMappingValue;
306     public double wetMappingValue;
307 
DryAndWetMappingValues(double dryMappingValue, double wetMappingValue)308     public DryAndWetMappingValues(double dryMappingValue, double wetMappingValue) {
309       this.dryMappingValue = dryMappingValue;
310       this.wetMappingValue = wetMappingValue;
311     }
312   }
313 
314   /* A class containing dry and wet delays in seconds experienced at zenith */
315   private static class DryAndWetZenithDelays {
316     public double dryZenithDelaySec;
317     public double wetZenithDelaySec;
318 
DryAndWetZenithDelays(double dryZenithDelay, double wetZenithDelay)319     public DryAndWetZenithDelays(double dryZenithDelay, double wetZenithDelay) {
320       this.dryZenithDelaySec = dryZenithDelay;
321       this.wetZenithDelaySec = wetZenithDelay;
322     }
323   }
324 }
325