1 /* 2 * Copyright (C) 2009 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.hardware; 18 19 import java.util.Calendar; 20 import java.util.TimeZone; 21 22 /** 23 * Estimates magnetic field at a given point on 24 * Earth, and in particular, to compute the magnetic declination from true 25 * north. 26 * 27 * <p>This uses the World Magnetic Model produced by the United States National 28 * Geospatial-Intelligence Agency. More details about the model can be found at 29 * <a href="http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml">http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml</a>. 30 * This class currently uses WMM-2020 which is valid until 2025, but should 31 * produce acceptable results for several years after that. Future versions of 32 * Android may use a newer version of the model. 33 */ 34 public class GeomagneticField { 35 // The magnetic field at a given point, in nanoteslas in geodetic 36 // coordinates. 37 private float mX; 38 private float mY; 39 private float mZ; 40 41 // Geocentric coordinates -- set by computeGeocentricCoordinates. 42 private float mGcLatitudeRad; 43 private float mGcLongitudeRad; 44 private float mGcRadiusKm; 45 46 // Constants from WGS84 (the coordinate system used by GPS) 47 static private final float EARTH_SEMI_MAJOR_AXIS_KM = 6378.137f; 48 static private final float EARTH_SEMI_MINOR_AXIS_KM = 6356.7523142f; 49 static private final float EARTH_REFERENCE_RADIUS_KM = 6371.2f; 50 51 // These coefficients and the formulae used below are from: 52 // NOAA Technical Report: The US/UK World Magnetic Model for 2020-2025 53 static private final float[][] G_COEFF = new float[][]{ 54 {0.0f}, 55 {-29404.5f, -1450.7f}, 56 {-2500.0f, 2982.0f, 1676.8f}, 57 {1363.9f, -2381.0f, 1236.2f, 525.7f}, 58 {903.1f, 809.4f, 86.2f, -309.4f, 47.9f}, 59 {-234.4f, 363.1f, 187.8f, -140.7f, -151.2f, 13.7f}, 60 {65.9f, 65.6f, 73.0f, -121.5f, -36.2f, 13.5f, -64.7f}, 61 {80.6f, -76.8f, -8.3f, 56.5f, 15.8f, 6.4f, -7.2f, 9.8f}, 62 {23.6f, 9.8f, -17.5f, -0.4f, -21.1f, 15.3f, 13.7f, -16.5f, -0.3f}, 63 {5.0f, 8.2f, 2.9f, -1.4f, -1.1f, -13.3f, 1.1f, 8.9f, -9.3f, -11.9f}, 64 {-1.9f, -6.2f, -0.1f, 1.7f, -0.9f, 0.6f, -0.9f, 1.9f, 1.4f, -2.4f, -3.9f}, 65 {3.0f, -1.4f, -2.5f, 2.4f, -0.9f, 0.3f, -0.7f, -0.1f, 1.4f, -0.6f, 0.2f, 3.1f}, 66 {-2.0f, -0.1f, 0.5f, 1.3f, -1.2f, 0.7f, 0.3f, 0.5f, -0.2f, -0.5f, 0.1f, -1.1f, -0.3f}}; 67 68 static private final float[][] H_COEFF = new float[][]{ 69 {0.0f}, 70 {0.0f, 4652.9f}, 71 {0.0f, -2991.6f, -734.8f}, 72 {0.0f, -82.2f, 241.8f, -542.9f}, 73 {0.0f, 282.0f, -158.4f, 199.8f, -350.1f}, 74 {0.0f, 47.7f, 208.4f, -121.3f, 32.2f, 99.1f}, 75 {0.0f, -19.1f, 25.0f, 52.7f, -64.4f, 9.0f, 68.1f}, 76 {0.0f, -51.4f, -16.8f, 2.3f, 23.5f, -2.2f, -27.2f, -1.9f}, 77 {0.0f, 8.4f, -15.3f, 12.8f, -11.8f, 14.9f, 3.6f, -6.9f, 2.8f}, 78 {0.0f, -23.3f, 11.1f, 9.8f, -5.1f, -6.2f, 7.8f, 0.4f, -1.5f, 9.7f}, 79 {0.0f, 3.4f, -0.2f, 3.5f, 4.8f, -8.6f, -0.1f, -4.2f, -3.4f, -0.1f, -8.8f}, 80 {0.0f, 0.0f, 2.6f, -0.5f, -0.4f, 0.6f, -0.2f, -1.7f, -1.6f, -3.0f, -2.0f, -2.6f}, 81 {0.0f, -1.2f, 0.5f, 1.3f, -1.8f, 0.1f, 0.7f, -0.1f, 0.6f, 0.2f, -0.9f, 0.0f, 0.5f}}; 82 83 static private final float[][] DELTA_G = new float[][]{ 84 {0.0f}, 85 {6.7f, 7.7f}, 86 {-11.5f, -7.1f, -2.2f}, 87 {2.8f, -6.2f, 3.4f, -12.2f}, 88 {-1.1f, -1.6f, -6.0f, 5.4f, -5.5f}, 89 {-0.3f, 0.6f, -0.7f, 0.1f, 1.2f, 1.0f}, 90 {-0.6f, -0.4f, 0.5f, 1.4f, -1.4f, 0.0f, 0.8f}, 91 {-0.1f, -0.3f, -0.1f, 0.7f, 0.2f, -0.5f, -0.8f, 1.0f}, 92 {-0.1f, 0.1f, -0.1f, 0.5f, -0.1f, 0.4f, 0.5f, 0.0f, 0.4f}, 93 {-0.1f, -0.2f, 0.0f, 0.4f, -0.3f, 0.0f, 0.3f, 0.0f, 0.0f, -0.4f}, 94 {0.0f, 0.0f, 0.0f, 0.2f, -0.1f, -0.2f, 0.0f, -0.1f, -0.2f, -0.1f, 0.0f}, 95 {0.0f, -0.1f, 0.0f, 0.0f, 0.0f, -0.1f, 0.0f, 0.0f, -0.1f, -0.1f, -0.1f, -0.1f}, 96 {0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, -0.1f}}; 97 98 static private final float[][] DELTA_H = new float[][]{ 99 {0.0f}, 100 {0.0f, -25.1f}, 101 {0.0f, -30.2f, -23.9f}, 102 {0.0f, 5.7f, -1.0f, 1.1f}, 103 {0.0f, 0.2f, 6.9f, 3.7f, -5.6f}, 104 {0.0f, 0.1f, 2.5f, -0.9f, 3.0f, 0.5f}, 105 {0.0f, 0.1f, -1.8f, -1.4f, 0.9f, 0.1f, 1.0f}, 106 {0.0f, 0.5f, 0.6f, -0.7f, -0.2f, -1.2f, 0.2f, 0.3f}, 107 {0.0f, -0.3f, 0.7f, -0.2f, 0.5f, -0.3f, -0.5f, 0.4f, 0.1f}, 108 {0.0f, -0.3f, 0.2f, -0.4f, 0.4f, 0.1f, 0.0f, -0.2f, 0.5f, 0.2f}, 109 {0.0f, 0.0f, 0.1f, -0.3f, 0.1f, -0.2f, 0.1f, 0.0f, -0.1f, 0.2f, 0.0f}, 110 {0.0f, 0.0f, 0.1f, 0.0f, 0.2f, 0.0f, 0.0f, 0.1f, 0.0f, -0.1f, 0.0f, 0.0f}, 111 {0.0f, 0.0f, 0.0f, -0.1f, 0.1f, 0.0f, 0.0f, 0.0f, 0.1f, 0.0f, 0.0f, 0.0f, -0.1f}}; 112 113 static private final long BASE_TIME = new Calendar.Builder() 114 .setTimeZone(TimeZone.getTimeZone("UTC")) 115 .setDate(2020, Calendar.JANUARY, 1) 116 .build() 117 .getTimeInMillis(); 118 119 // The ratio between the Gauss-normalized associated Legendre functions and 120 // the Schmid quasi-normalized ones. Compute these once staticly since they 121 // don't depend on input variables at all. 122 static private final float[][] SCHMIDT_QUASI_NORM_FACTORS = 123 computeSchmidtQuasiNormFactors(G_COEFF.length); 124 125 /** 126 * Estimate the magnetic field at a given point and time. 127 * 128 * @param gdLatitudeDeg 129 * Latitude in WGS84 geodetic coordinates -- positive is east. 130 * @param gdLongitudeDeg 131 * Longitude in WGS84 geodetic coordinates -- positive is north. 132 * @param altitudeMeters 133 * Altitude in WGS84 geodetic coordinates, in meters. 134 * @param timeMillis 135 * Time at which to evaluate the declination, in milliseconds 136 * since January 1, 1970. (approximate is fine -- the declination 137 * changes very slowly). 138 */ GeomagneticField(float gdLatitudeDeg, float gdLongitudeDeg, float altitudeMeters, long timeMillis)139 public GeomagneticField(float gdLatitudeDeg, 140 float gdLongitudeDeg, 141 float altitudeMeters, 142 long timeMillis) { 143 final int MAX_N = G_COEFF.length; // Maximum degree of the coefficients. 144 145 // We don't handle the north and south poles correctly -- pretend that 146 // we're not quite at them to avoid crashing. 147 gdLatitudeDeg = Math.min(90.0f - 1e-5f, 148 Math.max(-90.0f + 1e-5f, gdLatitudeDeg)); 149 computeGeocentricCoordinates(gdLatitudeDeg, 150 gdLongitudeDeg, 151 altitudeMeters); 152 153 assert G_COEFF.length == H_COEFF.length; 154 155 // Note: LegendreTable computes associated Legendre functions for 156 // cos(theta). We want the associated Legendre functions for 157 // sin(latitude), which is the same as cos(PI/2 - latitude), except the 158 // derivate will be negated. 159 LegendreTable legendre = 160 new LegendreTable(MAX_N - 1, 161 (float) (Math.PI / 2.0 - mGcLatitudeRad)); 162 163 // Compute a table of (EARTH_REFERENCE_RADIUS_KM / radius)^n for i in 164 // 0..MAX_N-2 (this is much faster than calling Math.pow MAX_N+1 times). 165 float[] relativeRadiusPower = new float[MAX_N + 2]; 166 relativeRadiusPower[0] = 1.0f; 167 relativeRadiusPower[1] = EARTH_REFERENCE_RADIUS_KM / mGcRadiusKm; 168 for (int i = 2; i < relativeRadiusPower.length; ++i) { 169 relativeRadiusPower[i] = relativeRadiusPower[i - 1] * 170 relativeRadiusPower[1]; 171 } 172 173 // Compute tables of sin(lon * m) and cos(lon * m) for m = 0..MAX_N -- 174 // this is much faster than calling Math.sin and Math.com MAX_N+1 times. 175 float[] sinMLon = new float[MAX_N]; 176 float[] cosMLon = new float[MAX_N]; 177 sinMLon[0] = 0.0f; 178 cosMLon[0] = 1.0f; 179 sinMLon[1] = (float) Math.sin(mGcLongitudeRad); 180 cosMLon[1] = (float) Math.cos(mGcLongitudeRad); 181 182 for (int m = 2; m < MAX_N; ++m) { 183 // Standard expansions for sin((m-x)*theta + x*theta) and 184 // cos((m-x)*theta + x*theta). 185 int x = m >> 1; 186 sinMLon[m] = sinMLon[m-x] * cosMLon[x] + cosMLon[m-x] * sinMLon[x]; 187 cosMLon[m] = cosMLon[m-x] * cosMLon[x] - sinMLon[m-x] * sinMLon[x]; 188 } 189 190 float inverseCosLatitude = 1.0f / (float) Math.cos(mGcLatitudeRad); 191 float yearsSinceBase = 192 (timeMillis - BASE_TIME) / (365f * 24f * 60f * 60f * 1000f); 193 194 // We now compute the magnetic field strength given the geocentric 195 // location. The magnetic field is the derivative of the potential 196 // function defined by the model. See NOAA Technical Report: The US/UK 197 // World Magnetic Model for 2020-2025 for the derivation. 198 float gcX = 0.0f; // Geocentric northwards component. 199 float gcY = 0.0f; // Geocentric eastwards component. 200 float gcZ = 0.0f; // Geocentric downwards component. 201 202 for (int n = 1; n < MAX_N; n++) { 203 for (int m = 0; m <= n; m++) { 204 // Adjust the coefficients for the current date. 205 float g = G_COEFF[n][m] + yearsSinceBase * DELTA_G[n][m]; 206 float h = H_COEFF[n][m] + yearsSinceBase * DELTA_H[n][m]; 207 208 // Negative derivative with respect to latitude, divided by 209 // radius. This looks like the negation of the version in the 210 // NOAA Technical report because that report used 211 // P_n^m(sin(theta)) and we use P_n^m(cos(90 - theta)), so the 212 // derivative with respect to theta is negated. 213 gcX += relativeRadiusPower[n+2] 214 * (g * cosMLon[m] + h * sinMLon[m]) 215 * legendre.mPDeriv[n][m] 216 * SCHMIDT_QUASI_NORM_FACTORS[n][m]; 217 218 // Negative derivative with respect to longitude, divided by 219 // radius. 220 gcY += relativeRadiusPower[n+2] * m 221 * (g * sinMLon[m] - h * cosMLon[m]) 222 * legendre.mP[n][m] 223 * SCHMIDT_QUASI_NORM_FACTORS[n][m] 224 * inverseCosLatitude; 225 226 // Negative derivative with respect to radius. 227 gcZ -= (n + 1) * relativeRadiusPower[n+2] 228 * (g * cosMLon[m] + h * sinMLon[m]) 229 * legendre.mP[n][m] 230 * SCHMIDT_QUASI_NORM_FACTORS[n][m]; 231 } 232 } 233 234 // Convert back to geodetic coordinates. This is basically just a 235 // rotation around the Y-axis by the difference in latitudes between the 236 // geocentric frame and the geodetic frame. 237 double latDiffRad = Math.toRadians(gdLatitudeDeg) - mGcLatitudeRad; 238 mX = (float) (gcX * Math.cos(latDiffRad) 239 + gcZ * Math.sin(latDiffRad)); 240 mY = gcY; 241 mZ = (float) (- gcX * Math.sin(latDiffRad) 242 + gcZ * Math.cos(latDiffRad)); 243 } 244 245 /** 246 * @return The X (northward) component of the magnetic field in nanoteslas. 247 */ getX()248 public float getX() { 249 return mX; 250 } 251 252 /** 253 * @return The Y (eastward) component of the magnetic field in nanoteslas. 254 */ getY()255 public float getY() { 256 return mY; 257 } 258 259 /** 260 * @return The Z (downward) component of the magnetic field in nanoteslas. 261 */ getZ()262 public float getZ() { 263 return mZ; 264 } 265 266 /** 267 * @return The declination of the horizontal component of the magnetic 268 * field from true north, in degrees (i.e. positive means the 269 * magnetic field is rotated east that much from true north). 270 */ getDeclination()271 public float getDeclination() { 272 return (float) Math.toDegrees(Math.atan2(mY, mX)); 273 } 274 275 /** 276 * @return The inclination of the magnetic field in degrees -- positive 277 * means the magnetic field is rotated downwards. 278 */ getInclination()279 public float getInclination() { 280 return (float) Math.toDegrees(Math.atan2(mZ, 281 getHorizontalStrength())); 282 } 283 284 /** 285 * @return Horizontal component of the field strength in nanoteslas. 286 */ getHorizontalStrength()287 public float getHorizontalStrength() { 288 return (float) Math.hypot(mX, mY); 289 } 290 291 /** 292 * @return Total field strength in nanoteslas. 293 */ getFieldStrength()294 public float getFieldStrength() { 295 return (float) Math.sqrt(mX * mX + mY * mY + mZ * mZ); 296 } 297 298 /** 299 * @param gdLatitudeDeg 300 * Latitude in WGS84 geodetic coordinates. 301 * @param gdLongitudeDeg 302 * Longitude in WGS84 geodetic coordinates. 303 * @param altitudeMeters 304 * Altitude above sea level in WGS84 geodetic coordinates. 305 * @return Geocentric latitude (i.e. angle between closest point on the 306 * equator and this point, at the center of the earth. 307 */ computeGeocentricCoordinates(float gdLatitudeDeg, float gdLongitudeDeg, float altitudeMeters)308 private void computeGeocentricCoordinates(float gdLatitudeDeg, 309 float gdLongitudeDeg, 310 float altitudeMeters) { 311 float altitudeKm = altitudeMeters / 1000.0f; 312 float a2 = EARTH_SEMI_MAJOR_AXIS_KM * EARTH_SEMI_MAJOR_AXIS_KM; 313 float b2 = EARTH_SEMI_MINOR_AXIS_KM * EARTH_SEMI_MINOR_AXIS_KM; 314 double gdLatRad = Math.toRadians(gdLatitudeDeg); 315 float clat = (float) Math.cos(gdLatRad); 316 float slat = (float) Math.sin(gdLatRad); 317 float tlat = slat / clat; 318 float latRad = 319 (float) Math.sqrt(a2 * clat * clat + b2 * slat * slat); 320 321 mGcLatitudeRad = (float) Math.atan(tlat * (latRad * altitudeKm + b2) 322 / (latRad * altitudeKm + a2)); 323 324 mGcLongitudeRad = (float) Math.toRadians(gdLongitudeDeg); 325 326 float radSq = altitudeKm * altitudeKm 327 + 2 * altitudeKm * (float) Math.sqrt(a2 * clat * clat + 328 b2 * slat * slat) 329 + (a2 * a2 * clat * clat + b2 * b2 * slat * slat) 330 / (a2 * clat * clat + b2 * slat * slat); 331 mGcRadiusKm = (float) Math.sqrt(radSq); 332 } 333 334 335 /** 336 * Utility class to compute a table of Gauss-normalized associated Legendre 337 * functions P_n^m(cos(theta)) 338 */ 339 static private class LegendreTable { 340 // These are the Gauss-normalized associated Legendre functions -- that 341 // is, they are normal Legendre functions multiplied by 342 // (n-m)!/(2n-1)!! (where (2n-1)!! = 1*3*5*...*2n-1) 343 public final float[][] mP; 344 345 // Derivative of mP, with respect to theta. 346 public final float[][] mPDeriv; 347 348 /** 349 * @param maxN 350 * The maximum n- and m-values to support 351 * @param thetaRad 352 * Returned functions will be Gauss-normalized 353 * P_n^m(cos(thetaRad)), with thetaRad in radians. 354 */ LegendreTable(int maxN, float thetaRad)355 public LegendreTable(int maxN, float thetaRad) { 356 // Compute the table of Gauss-normalized associated Legendre 357 // functions using standard recursion relations. Also compute the 358 // table of derivatives using the derivative of the recursion 359 // relations. 360 float cos = (float) Math.cos(thetaRad); 361 float sin = (float) Math.sin(thetaRad); 362 363 mP = new float[maxN + 1][]; 364 mPDeriv = new float[maxN + 1][]; 365 mP[0] = new float[] { 1.0f }; 366 mPDeriv[0] = new float[] { 0.0f }; 367 for (int n = 1; n <= maxN; n++) { 368 mP[n] = new float[n + 1]; 369 mPDeriv[n] = new float[n + 1]; 370 for (int m = 0; m <= n; m++) { 371 if (n == m) { 372 mP[n][m] = sin * mP[n - 1][m - 1]; 373 mPDeriv[n][m] = cos * mP[n - 1][m - 1] 374 + sin * mPDeriv[n - 1][m - 1]; 375 } else if (n == 1 || m == n - 1) { 376 mP[n][m] = cos * mP[n - 1][m]; 377 mPDeriv[n][m] = -sin * mP[n - 1][m] 378 + cos * mPDeriv[n - 1][m]; 379 } else { 380 assert n > 1 && m < n - 1; 381 float k = ((n - 1) * (n - 1) - m * m) 382 / (float) ((2 * n - 1) * (2 * n - 3)); 383 mP[n][m] = cos * mP[n - 1][m] - k * mP[n - 2][m]; 384 mPDeriv[n][m] = -sin * mP[n - 1][m] 385 + cos * mPDeriv[n - 1][m] - k * mPDeriv[n - 2][m]; 386 } 387 } 388 } 389 } 390 } 391 392 /** 393 * Compute the ration between the Gauss-normalized associated Legendre 394 * functions and the Schmidt quasi-normalized version. This is equivalent to 395 * sqrt((m==0?1:2)*(n-m)!/(n+m!))*(2n-1)!!/(n-m)! 396 */ 397 private static float[][] computeSchmidtQuasiNormFactors(int maxN) { 398 float[][] schmidtQuasiNorm = new float[maxN + 1][]; 399 schmidtQuasiNorm[0] = new float[] { 1.0f }; 400 for (int n = 1; n <= maxN; n++) { 401 schmidtQuasiNorm[n] = new float[n + 1]; 402 schmidtQuasiNorm[n][0] = 403 schmidtQuasiNorm[n - 1][0] * (2 * n - 1) / (float) n; 404 for (int m = 1; m <= n; m++) { 405 schmidtQuasiNorm[n][m] = schmidtQuasiNorm[n][m - 1] 406 * (float) Math.sqrt((n - m + 1) * (m == 1 ? 2 : 1) 407 / (float) (n + m)); 408 } 409 } 410 return schmidtQuasiNorm; 411 } 412 } 413