GeomagneticFieldpublic class GeomagneticField extends Object | This class is used to estimated estimate magnetic field at a given point on
Earth, and in particular, to compute the magnetic declination from true
north.
This uses the World Magnetic Model produced by the United States National
Geospatial-Intelligence Agency. More details about the model can be found at
http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml.
This class currently uses WMM-2005 which is valid until 2010, but should
produce acceptable results for several years after that. |
| Fields Summary |
|---|
| private float | mX | | private float | mY | | private float | mZ | | private float | mGcLatitudeRad | | private float | mGcLongitudeRad | | private float | mGcRadiusKm | | private static final float | EARTH_SEMI_MAJOR_AXIS_KM | | private static final float | EARTH_SEMI_MINOR_AXIS_KM | | private static final float | EARTH_REFERENCE_RADIUS_KM | | private static final float[] | G_COEFF | | private static final float[] | H_COEFF | | private static final float[] | DELTA_G | | private static final float[] | DELTA_H | | private static final long | BASE_TIME | | private static final float[] | SCHMIDT_QUASI_NORM_FACTORS |
| Constructors Summary |
|---|
public GeomagneticField(float gdLatitudeDeg, float gdLongitudeDeg, float altitudeMeters, long timeMillis)Estimate the magnetic field at a given point and time.
final int MAX_N = G_COEFF.length; // Maximum degree of the coefficients.
// We don't handle the north and south poles correctly -- pretend that
// we're not quite at them to avoid crashing.
gdLatitudeDeg = Math.min(90.0f - 1e-5f,
Math.max(-90.0f + 1e-5f, gdLatitudeDeg));
computeGeocentricCoordinates(gdLatitudeDeg,
gdLongitudeDeg,
altitudeMeters);
assert G_COEFF.length == H_COEFF.length;
// Note: LegendreTable computes associated Legendre functions for
// cos(theta). We want the associated Legendre functions for
// sin(latitude), which is the same as cos(PI/2 - latitude), except the
// derivate will be negated.
LegendreTable legendre =
new LegendreTable(MAX_N - 1,
(float) (Math.PI / 2.0 - mGcLatitudeRad));
// Compute a table of (EARTH_REFERENCE_RADIUS_KM / radius)^n for i in
// 0..MAX_N-2 (this is much faster than calling Math.pow MAX_N+1 times).
float[] relativeRadiusPower = new float[MAX_N + 2];
relativeRadiusPower[0] = 1.0f;
relativeRadiusPower[1] = EARTH_REFERENCE_RADIUS_KM / mGcRadiusKm;
for (int i = 2; i < relativeRadiusPower.length; ++i) {
relativeRadiusPower[i] = relativeRadiusPower[i - 1] *
relativeRadiusPower[1];
}
// Compute tables of sin(lon * m) and cos(lon * m) for m = 0..MAX_N --
// this is much faster than calling Math.sin and Math.com MAX_N+1 times.
float[] sinMLon = new float[MAX_N];
float[] cosMLon = new float[MAX_N];
sinMLon[0] = 0.0f;
cosMLon[0] = 1.0f;
sinMLon[1] = (float) Math.sin(mGcLongitudeRad);
cosMLon[1] = (float) Math.cos(mGcLongitudeRad);
for (int m = 2; m < MAX_N; ++m) {
// Standard expansions for sin((m-x)*theta + x*theta) and
// cos((m-x)*theta + x*theta).
int x = m >> 1;
sinMLon[m] = sinMLon[m-x] * cosMLon[x] + cosMLon[m-x] * sinMLon[x];
cosMLon[m] = cosMLon[m-x] * cosMLon[x] - sinMLon[m-x] * sinMLon[x];
}
float inverseCosLatitude = 1.0f / (float) Math.cos(mGcLatitudeRad);
float yearsSinceBase =
(timeMillis - BASE_TIME) / (365f * 24f * 60f * 60f * 1000f);
// We now compute the magnetic field strength given the geocentric
// location. The magnetic field is the derivative of the potential
// function defined by the model. See NOAA Technical Report: The US/UK
// World Magnetic Model for 2005-2010 for the derivation.
float gcX = 0.0f; // Geocentric northwards component.
float gcY = 0.0f; // Geocentric eastwards component.
float gcZ = 0.0f; // Geocentric downwards component.
for (int n = 1; n < MAX_N; n++) {
for (int m = 0; m <= n; m++) {
// Adjust the coefficients for the current date.
float g = G_COEFF[n][m] + yearsSinceBase * DELTA_G[n][m];
float h = H_COEFF[n][m] + yearsSinceBase * DELTA_H[n][m];
// Negative derivative with respect to latitude, divided by
// radius. This looks like the negation of the version in the
// NOAA Techincal report because that report used
// P_n^m(sin(theta)) and we use P_n^m(cos(90 - theta)), so the
// derivative with respect to theta is negated.
gcX += relativeRadiusPower[n+2]
* (g * cosMLon[m] + h * sinMLon[m])
* legendre.mPDeriv[n][m]
* SCHMIDT_QUASI_NORM_FACTORS[n][m];
// Negative derivative with respect to longitude, divided by
// radius.
gcY += relativeRadiusPower[n+2] * m
* (g * sinMLon[m] - h * cosMLon[m])
* legendre.mP[n][m]
* SCHMIDT_QUASI_NORM_FACTORS[n][m]
* inverseCosLatitude;
// Negative derivative with respect to radius.
gcZ -= (n + 1) * relativeRadiusPower[n+2]
* (g * cosMLon[m] + h * sinMLon[m])
* legendre.mP[n][m]
* SCHMIDT_QUASI_NORM_FACTORS[n][m];
}
}
// Convert back to geodetic coordinates. This is basically just a
// rotation around the Y-axis by the difference in latitudes between the
// geocentric frame and the geodetic frame.
double latDiffRad = Math.toRadians(gdLatitudeDeg) - mGcLatitudeRad;
mX = (float) (gcX * Math.cos(latDiffRad)
+ gcZ * Math.sin(latDiffRad));
mY = gcY;
mZ = (float) (- gcX * Math.sin(latDiffRad)
+ gcZ * Math.cos(latDiffRad));
|
| Methods Summary |
|---|
| private void | computeGeocentricCoordinates(float gdLatitudeDeg, float gdLongitudeDeg, float altitudeMeters)
float altitudeKm = altitudeMeters / 1000.0f;
float a2 = EARTH_SEMI_MAJOR_AXIS_KM * EARTH_SEMI_MAJOR_AXIS_KM;
float b2 = EARTH_SEMI_MINOR_AXIS_KM * EARTH_SEMI_MINOR_AXIS_KM;
double gdLatRad = Math.toRadians(gdLatitudeDeg);
float clat = (float) Math.cos(gdLatRad);
float slat = (float) Math.sin(gdLatRad);
float tlat = slat / clat;
float latRad =
(float) Math.sqrt(a2 * clat * clat + b2 * slat * slat);
mGcLatitudeRad = (float) Math.atan(tlat * (latRad * altitudeKm + b2)
/ (latRad * altitudeKm + a2));
mGcLongitudeRad = (float) Math.toRadians(gdLongitudeDeg);
float radSq = altitudeKm * altitudeKm
+ 2 * altitudeKm * (float) Math.sqrt(a2 * clat * clat +
b2 * slat * slat)
+ (a2 * a2 * clat * clat + b2 * b2 * slat * slat)
/ (a2 * clat * clat + b2 * slat * slat);
mGcRadiusKm = (float) Math.sqrt(radSq);
| | private static float[][] | computeSchmidtQuasiNormFactors(int maxN)Compute the ration between the Gauss-normalized associated Legendre
functions and the Schmidt quasi-normalized version. This is equivalent to
sqrt((m==0?1:2)*(n-m)!/(n+m!))*(2n-1)!!/(n-m)!
float[][] schmidtQuasiNorm = new float[maxN + 1][];
schmidtQuasiNorm[0] = new float[] { 1.0f };
for (int n = 1; n <= maxN; n++) {
schmidtQuasiNorm[n] = new float[n + 1];
schmidtQuasiNorm[n][0] =
schmidtQuasiNorm[n - 1][0] * (2 * n - 1) / (float) n;
for (int m = 1; m <= n; m++) {
schmidtQuasiNorm[n][m] = schmidtQuasiNorm[n][m - 1]
* (float) Math.sqrt((n - m + 1) * (m == 1 ? 2 : 1)
/ (float) (n + m));
}
}
return schmidtQuasiNorm;
| | public float | getDeclination()
return (float) Math.toDegrees(Math.atan2(mY, mX));
| | public float | getFieldStrength()
return (float) Math.sqrt(mX * mX + mY * mY + mZ * mZ);
| | public float | getHorizontalStrength()
return (float) Math.sqrt(mX * mX + mY * mY);
| | public float | getInclination()
return (float) Math.toDegrees(Math.atan2(mZ,
getHorizontalStrength()));
| | public float | getX()
return mX;
| | public float | getY()
return mY;
| | public float | getZ()
return mZ;
|
|
