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.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