FileDocCategorySizeDatePackage
GeomagneticField.javaAPI DocAndroid 1.5 API18096Wed May 06 22:41:54 BST 2009android.hardware

GeomagneticField

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

param
gdLatitudeDeg Latitude in WGS84 geodetic coordinates -- positive is east.
param
gdLongitudeDeg Longitude in WGS84 geodetic coordinates -- positive is north.
param
altitudeMeters Altitude in WGS84 geodetic coordinates, in meters.
param
timeMillis Time at which to evaluate the declination, in milliseconds since January 1, 1970. (approximate is fine -- the declination changes very slowly).


                                                                                                                                            
      
                             
                             
                              
        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 voidcomputeGeocentricCoordinates(float gdLatitudeDeg, float gdLongitudeDeg, float altitudeMeters)

param
gdLatitudeDeg Latitude in WGS84 geodetic coordinates.
param
gdLongitudeDeg Longitude in WGS84 geodetic coordinates.
param
altitudeMeters Altitude above sea level in WGS84 geodetic coordinates.
return
Geocentric latitude (i.e. angle between closest point on the equator and this point, at the center of the earth.

        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 floatgetDeclination()

return
The declination of the horizontal component of the magnetic field from true north, in degrees (i.e. positive means the magnetic field is rotated east that much from true north).

        return (float) Math.toDegrees(Math.atan2(mY, mX));
    
public floatgetFieldStrength()

return
Total field strength in nanoteslas.

        return (float) Math.sqrt(mX * mX + mY * mY + mZ * mZ);
    
public floatgetHorizontalStrength()

return
Horizontal component of the field strength in nonoteslas.

        return (float) Math.sqrt(mX * mX + mY * mY);
    
public floatgetInclination()

return
The inclination of the magnetic field in degrees -- positive means the magnetic field is rotated downwards.

        return (float) Math.toDegrees(Math.atan2(mZ,
                                                 getHorizontalStrength()));
    
public floatgetX()

return
The X (northward) component of the magnetic field in nanoteslas.

        return mX;
    
public floatgetY()

return
The Y (eastward) component of the magnetic field in nanoteslas.

        return mY;
    
public floatgetZ()

return
The Z (downward) component of the magnetic field in nanoteslas.

        return mZ;