// $Revision: 1.1 $ // Derived from geomag.c, by Dean Nelson // http://geostarslib.sourceforge.net/ // Java conversion contributed by Edward Falk package com.android.flightdeck; import java.lang.Math; import java.util.Calendar; import java.util.GregorianCalendar; public class GeoMag { static final int wmm_n[] = { 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12}; static final int wmm_m[] = { 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12}; static final double wmm_gnm[] = { -29556.8, -1671.7, -2340.6, 3046.9, 1657.0, 1335.4, -2305.1, 1246.7, 674.0, 919.8, 798.1, 211.3, -379.4, 100.0, -227.4, 354.6, 208.7, -136.5, -168.3, -14.1, 73.2, 69.7, 76.7, -151.2, -14.9, 14.6, -86.3, 80.1, -74.5, -1.4, 38.5, 12.4, 9.5, 5.7, 1.8, 24.9, 7.7, -11.6, -6.9, -18.2, 10.0, 9.2, -11.6, -5.2, 5.6, 9.9, 3.5, -7.0, 5.1, -10.8, -1.3, 8.8, -6.7, -9.1, -2.3, -6.3, 1.6, -2.6, 0.0, 3.1, 0.4, 2.1, 3.9, -0.1, -2.3, 2.8, -1.6, -1.7, 1.7, -0.1, 0.1, -0.7, 0.7, 1.8, 0.0, 1.1, 4.1, -2.4, -0.4, 0.2, 0.8, -0.3, 1.1, -0.5, 0.4, -0.3, -0.3, -0.1, -0.3, -0.1}; static final double wmm_hnm[] = { 0.0, 5079.8, 0.0, -2594.7, -516.7, 0.0, -199.9, 269.3, -524.2, 0.0, 281.5, -226.0, 145.8, -304.7, 0.0, 42.4, 179.8, -123.0, -19.5, 103.6, 0.0, -20.3, 54.7, 63.6, -63.4, -0.1, 50.4, 0.0, -61.5, -22.4, 7.2, 25.4, 11.0, -26.4, -5.1, 0.0, 11.2, -21.0, 9.6, -19.8, 16.1, 7.7, -12.9, -0.2, 0.0, -20.1, 12.9, 12.6, -6.7, -8.1, 8.0, 2.9, -7.9, 6.0, 0.0, 2.4, 0.2, 4.4, 4.8, -6.5, -1.1, -3.4, -0.8, -2.3, -7.9, 0.0, 0.3, 1.2, -0.8, -2.5, 0.9, -0.6, -2.7, -0.9, -1.3, -2.0, -1.2, 0.0, -0.4, 0.3, 2.4, -2.6, 0.6, 0.3, 0.0, 0.0, 0.3, -0.9, -0.4, 0.8}; static final double wmm_dgnm[] = { 8.0, 10.6, -15.1, -7.8, -0.8, 0.4, -2.6, -1.2, -6.5, -2.5, 2.8, -7.0, 6.2, -3.8, -2.8, 0.7, -3.2, -1.1, 0.1, -0.8, -0.7, 0.4, -0.3, 2.3, -2.1, -0.6, 1.4, 0.2, -0.1, -0.3, 1.1, 0.6, 0.5, -0.4, 0.6, 0.1, 0.3, -0.4, 0.3, -0.3, 0.2, 0.4, -0.7, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; static final double wmm_dhnm[] = { 0.0, -20.9, 0.0, -23.2, -14.6, 0.0, 5.0, -7.0, -0.6, 0.0, 2.2, 1.6, 5.8, 0.1, 0.0, 0.0, 1.7, 2.1, 4.8, -1.1, 0.0, -0.6, -1.9, -0.4, -0.5, -0.3, 0.7, 0.0, 0.6, 0.4, 0.2, 0.3, -0.8, -0.2, 0.1, 0.0, -0.2, 0.1, 0.3, 0.4, 0.1, -0.2, 0.4, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; static private final int MAXDEG = 12; static private final double DEG_TO_RAD = Math.PI/180; static private final double RAD_TO_DEG = 180./Math.PI; static private final double DEG_TO_MIN = 60; static private final double CIRCLE = 360.; static private final double HALF_CIRCLE = CIRCLE / 2.0; static private final double ellips_a = 6378137.0; // From WGS 1984 static private final double ellips_fl = 298.257223563; static private boolean started = false; static private double c[][] = new double[13][13]; // Gaus Coeffs of main model (NT) static private double cd[][] = new double[13][13]; // Gaus Coeffs of secular model (NT/YR) static private double tc[][] = new double[13][13]; // Time adjusted gaus coeffs (NT) static private double dp[][] = new double[13][13]; // theta derivative of P(n,m) static private double p[][] = new double[13][13]; // Schmidt normalization factors static private double sp[] = new double[13]; // Sine(M*spherical coord longitude) static private double cp[] = new double[13]; // Cosine(M*spherical coord longitude) static private double fn[] = new double[13]; static private double fm[] = new double[13]; static private double pp[] = new double[13]; // Legendre polynomials for m = 1 static private double k[][] = new double[13][13]; static private double epoch; static private int maxord; // Max order of spherical harmonic model static class GeoMagValues { public double dec, dip, ti, gv, adec, adip, ati, x, y, z, h, ax, ay, az, ah; } static private void geoMagInit() { int n, m, i, j; double flnmj; // Only do this once if (started) return; started = true; // INITIALIZE CONSTANTS maxord = MAXDEG; maxord = MAXDEG; sp[0] = 0.0; cp[0] = p[0][0] = pp[0] = 1.0; dp[0][0] = 0.0; // READ WORLD MAGNETIC MODEL SPHERICAL HARMONIC COEFFICIENTS c[0][0] = 0.0; cd[0][0] = 0.0; // Load the WMM Coefficients epoch=2005.0; for( i=0;i<90;i++) { n=wmm_n[i]; m=wmm_m[i]; if (m <= n) { c[m][n] = wmm_gnm[i]; cd[m][n] = wmm_dgnm[i]; if (m != 0) { c[n][m-1] = wmm_hnm[i]; cd[n][m-1] = wmm_dhnm[i]; } } } //end for /* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */ p[0][0] = 1.0; for (n=1; n<=maxord; n++) { p[0][n] = p[0][n-1]*(double)(2*n-1)/(double)n; j = 2; for (m=0;m<=n;m++) { k[m][n] = (double)(((n-1)*(n-1))-(m*m))/(double)((2*n-1)*(2*n-3)); if (m > 0) { flnmj = (double)((n-m+1)*j)/(double)(n+m); p[m][n] = p[m-1][n]*Math.sqrt(flnmj); j = 1; c[n][m-1] = p[m][n]*c[n][m-1]; cd[n][m-1] = p[m][n]*cd[n][m-1]; } c[m][n] = p[m][n]*c[m][n]; cd[m][n] = p[m][n]*cd[m][n]; } //end for m fn[n] = (double)(n+1); fm[n] = (double)n; } //end for n k[1][1] = 0.0; } // Main GeoMag routine /*! \brief This routine computes all of the relevant geomagnetic data \param double alt : in: altitude in meters \param double glat : in: latitude in decimal degrees \param double glon : in: longitude in decimal degrees \param double time : in: time in decimal years \param double *dec :out: declination in degrees \param double *dip :out: dip in degrees \param double *ti :out: total intensity in nT \param double *gv :out: geomagnetic grid variation \param double *adec :out: Annual declination in degrees \param double *adip :out: Annual dip in degrees \param double *ati :out: Annual total intensity in nT \param double *x :out: X Component of the magnetic field \param double *y :out: Y Component of the magnetic field \param double *z :out: Z Component of the magnetic field \param double *h :out: h Component of the magnetic field \param double *ax :out: Annual X Component of the magnetic field \param double *ay :out: Annual Y Component of the magnetic field \param double *az :out: Annual Z Component of the magnetic field \param double *ah :out: Annual H Component of the magnetic field \retval values on success \retval null on error */ static public GeoMagValues geoMag( double alt, double glat, double glon, double time) { double time1; double rdec, rdip, rdec2, rdip2; GeoMagValues values = new GeoMagValues(); GeoMagValues values2 = new GeoMagValues(); geoMagInit(); geomg1(alt,glat,glon,time, values); time1 = time + 1.0; // add a year to time geomg1(alt,glat,glon,time1, values2); /*COMPUTE X, Y, Z, AND H COMPONENTS OF THE MAGNETIC FIELD*/ rdec = values.dec * DEG_TO_RAD; rdip = values.dip * DEG_TO_RAD; rdec2 = values2.dec * DEG_TO_RAD; rdip2 = values2.dip * DEG_TO_RAD; values.x = values.ti * (Math.cos(rdec) * Math.cos(rdip)); values.y = values.ti * (Math.cos(rdip) * Math.sin(rdec)); values.z = values.ti * (Math.sin(rdip)); values.h = values.ti * (Math.cos(rdip)); values2.x = values2.ti * (Math.cos(rdec2) * Math.cos(rdip2)); values2.y = values2.ti * (Math.cos(rdip2) * Math.sin(rdec2)); values2.z = values2.ti * (Math.sin(rdip2)); values2.h = values2.ti * (Math.cos(rdip2)); /* COMPUTE ANNUAL CHANGE FOR TOTAL INTENSITY */ values.ati = values2.ti - values.ti; /* COMPUTE ANNUAL CHANGE FOR DIP & DEC (in minutes) */ values.adip = (rdip - rdip) * DEG_TO_MIN; values.adec = (rdec - rdec) * DEG_TO_MIN; /* COMPUTE ANNUAL CHANGE FOR X, Y, Z, AND H */ values.ax = values2.x - values.x; values.ay = values2.y - values.y; values.az = values2.z - values.z; values.ah = values2.h - values.h; return values; } static void geomg1( double alt, double glat, double glon, double time, GeoMagValues values) { /*! This is documentation that came with the WMM. It is included here for completeness. \verbatim // MAXDEG - MAXIMUM DEGREE OF SPHERICAL HARMONIC MODEL (INPUT) // TIME - COMPUTATION TIME (YRS) (INPUT) // (EG. 1 JULY 1995 = 1995.500) // ALT - GEODETIC ALTITUDE (M) (INPUT) // GLAT - GEODETIC LATITUDE (DEG.) (INPUT) // GLON - GEODETIC LONGITUDE (DEG.) (INPUT) // EPOCH - BASE TIME OF GEOMAGNETIC MODEL (YRS) // P(N,M) - ASSOCIATED LEGENDRE POLYNOMIALS (UNNORMALIZED) // DEC - GEOMAGNETIC DECLINATION (DEG.) (OUTPUT) // EAST=POSITIVE ANGLES // WEST=NEGATIVE ANGLES // DIP - GEOMAGNETIC INCLINATION (DEG.) (OUTPUT) // DOWN=POSITIVE ANGLES // UP=NEGATIVE ANGLES // TI - GEOMAGNETIC TOTAL INTENSITY (NT) (OUTPUT) // GV - GEOMAGNETIC GRID VARIATION (DEG.) (OUTPUT) // REFERENCED TO GRID NORTH // GRID NORTH REFERENCED TO 0 MERIDIAN // OF A POLAR STEREOGRAPHIC PROJECTION // (ARCTIC/ANTARCTIC ONLY) // A - SEMIMAJOR AXIS OF WGS-84 ELLIPSOID (KM) // B - SEMIMINOR AXIS OF WGS-84 ELLIPSOID (KM) // RE - MEAN RADIUS OF IAU-66 ELLIPSOID (KM) // a2, // a * a // b2, // b * b // c2, // a2 - b2 // a4, // a2 * a2 // b4, // b2 * b2 // c4, // a4 - b4 // ct, //C CT - COSINE OF (SPHERICAL COORD. LATITUDE) // st, //C ST - SINE OF (SPHERICAL COORD. LATITUDE) // r2, // r, //C R - SPHERICAL COORDINATE RADIAL POSITION (KM) // ca, //C CA - COSINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE // sa, //C SA - SINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE // br, //C BR - RADIAL COMPONENT OF GEOMAGNETIC FIELD (NT) // bt, //C BT - THETA COMPONENT OF GEOMAGNETIC FIELD (NT) // bp, //C BP - PHI COMPONENT OF GEOMAGNETIC FIELD (NT) // bx, //C BX - NORTH GEOMAGNETIC COMPONENT (NT) // by, //C BY - EAST GEOMAGNETIC COMPONENT (NT) // bz, //C BZ - VERTICALLY DOWN GEOMAGNETIC COMPONENT (NT) // bh; //C BH - HORIZONTAL GEOMAGNETIC COMPONENT (NT) \endverbatim */ int m,n; double temp1, temp2,d,ar,aor,dt,bpp, par, parp,r,r2,sa,ca,st,ct,br,bt,bp,bx,by,bz,bh; double rlon,rlat, srlon, srlat, crlon, crlat, srlat2, crlat2, a, // A - SEMIMAJOR AXIS OF WGS-84 ELLIPSOID (KM) b, // B - SEMIMINOR AXIS OF WGS-84 ELLIPSOID (KM) re, // RE - MEAN RADIUS OF IAU-66 ELLIPSOID (KM) a2, // a * a b2, // b * b c2, // a2 - b2 a4, // a2 * a2 b4, // b2 * b2 c4, // a4 - b4 q, q1, q2; a = ellips_a; //6378137.0; b = GEO_B(ellips_a, (ellips_fl)); //6356752.3142; //wgs-84 re = 6371200.0; a2 = a*a; b2 = b*b; c2 = a2-b2; a4 = a2*a2; b4 = b2*b2; c4 = a4 - b4; dt = time - epoch; // if (dt < 0.0 || dt > 5.0) // printf("\n WARNING geoMag - TIME EXTENDS BEYOND MODEL 5-YEAR LIFE SPAN (dt=%f, time=%f)\n",dt,time); rlon = glon*DEG_TO_RAD; rlat = glat*DEG_TO_RAD; srlon = Math.sin(rlon); srlat = Math.sin(rlat); crlon = Math.cos(rlon); crlat = Math.cos(rlat); srlat2 = srlat*srlat; crlat2 = crlat*crlat; sp[1] = srlon; cp[1] = crlon; /* CONVERT FROM GEODETIC COORDS. TO SPHERICAL COORDS. */ q = Math.sqrt(a2-c2*srlat2); q1 = alt*q; q2 = ((q1+a2)/(q1+b2))*((q1+a2)/(q1+b2)); ct = srlat/Math.sqrt(q2*crlat2+srlat2); st = Math.sqrt(1.0-(ct*ct)); r2 = (alt*alt)+2.0*q1+(a4-c4*srlat2)/(q*q); r = Math.sqrt(r2); d = Math.sqrt(a2*crlat2+b2*srlat2); ca = (alt+d)/r; sa = c2*crlat*srlat/(r*d); for (m=2; m<=maxord; m++) { sp[m] = sp[1]*cp[m-1]+cp[1]*sp[m-1]; cp[m] = cp[1]*cp[m-1]-sp[1]*sp[m-1]; } aor = re/r; ar = aor*aor; br = bt = bp = bpp = 0.0; for (n=1; n<=maxord; n++) { ar = ar*aor; for (m=0;m<=n;m++) { /* COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS AND DERIVATIVES VIA RECURSION RELATIONS */ if (n == m) { p[m][n] = st * p[m-1][n-1]; dp[m][n] = st * dp[m-1][n-1] + ct * p[m-1][n-1]; } else if (n == 1 && m == 0) { p[m][n] = ct * p[m][n-1]; dp[m][n] = ct * dp[m][n-1] - st * p[m][n-1]; } else if (n > 1 && n != m) { if (m > n-2) p[m][n-2] = 0.0; if (m > n-2) dp[m][n-2] = 0.0; p[m][n] = ct * p[m][n-1] - k[m][n] * p[m][n-2]; dp[m][n] = ct * dp[m][n-1] - st * p[m][n-1] - k[m][n] * dp[m][n-2]; } // TIME ADJUST THE GAUSS COEFFICIENTS tc[m][n] = c[m][n]+dt*cd[m][n]; if (m != 0) tc[n][m-1] = c[n][m-1]+dt*cd[n][m-1]; // ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS par = ar * p[m][n]; if (m == 0) { temp1 = tc[m][n]*cp[m]; temp2 = tc[m][n]*sp[m]; } else { temp1 = tc[m][n]*cp[m]+tc[n][m-1]*sp[m]; temp2 = tc[m][n]*sp[m]-tc[n][m-1]*cp[m]; } bt = bt-ar*temp1*dp[m][n]; bp += (fm[m]*temp2*par); br += (fn[n]*temp1*par); // SPECIAL CASE: NORTH/SOUTH GEOGRAPHIC POLES if (st == 0.0 && m == 1) { if (n == 1) pp[n] = pp[n-1]; else pp[n] = ct*pp[n-1]-k[m][n]*pp[n-2]; parp = ar*pp[n]; bpp += (fm[m]*temp2*parp); } } } if (st == 0.0) bp = bpp; else bp /= st; /* ROTATE MAGNETIC VECTOR COMPONENTS FROM SPHERICAL TO GEODETIC COORDINATES */ bx = -bt*ca-br*sa; by = bp; bz = bt*sa-br*ca; /* ** COMPUTE DECLINATION (DEC), INCLINATION (DIP) AND ** TOTAL INTENSITY (TI) */ bh = Math.sqrt((bx*bx)+(by*by)); values.ti = Math.sqrt((bh*bh)+(bz*bz)); values.dec = Math.atan2(by,bx)*RAD_TO_DEG; values.dip = Math.atan2(bz,bh)*RAD_TO_DEG; /* ** COMPUTE MAGNETIC GRID VARIATION IF THE CURRENT ** GEODETIC POSITION IS IN THE ARCTIC OR ANTARCTIC ** (I.E. GLAT > +55 DEGREES OR GLAT < -55 DEGREES) ** ** OTHERWISE, SET MAGNETIC GRID VARIATION TO -999.0 */ values.gv = -999.0; if (Math.abs(glat) >= 55.) { if (glat > 0.0 && glon >= 0.0) values.gv = values.dec-glon; if (glat > 0.0 && glon < 0.0) values.gv = values.dec+Math.abs(glon); if (glat < 0.0 && glon >= 0.0) values.gv = values.dec+glon; if (glat < 0.0 && glon < 0.0) values.gv = values.dec-Math.abs(glon); if (values.gv > +HALF_CIRCLE) values.gv -= CIRCLE; if (values.gv < -HALF_CIRCLE) values.gv += CIRCLE; } } static private double GEO_B(double a, double fl) { return a*(1.0-(1.0/fl)); } /*************************************************************************/ // Short GeoMag routine /*! \brief This routine computes magnetic declination \param double lat : in: latitude in decimal degrees \param double lon : in: longitude in decimal degrees \param double hgt : in: altitude in meters \param double time : in: time in decimal years */ static double geoMagGetDecRet(double lat, double lon, double hgt, double time) { GeoMagValues values = new GeoMagValues(); /* Determine declination */ geoMagInit(); geomg1(hgt,lat,lon,time, values); return values.dec; } /*************************************************************************/ // Short GeoMag routine /*! \brief This routine computes magnetic declination \param double lat : in: latitude in decimal degrees \param double lon : in: longitude in decimal degrees \param double hgt : in: altitude in meters \param int month : in: Month of the year \param int day : in: Day of the month \param int year : in: Year */ static double geoMagGetDecRet(double lat, double lon, double hgt, int month, int day, int year) { /* Days in each month 1-12 */ int jday[] = {0,0,31,59,90,120,151,180,211,242,272,303,333 }; double time; /* Determine decimal years from date */ time = ((jday[month]+day) / 365.0) + year; return geoMagGetDecRet(lat, lon, hgt, time); } /*************************************************************************/ // Short GeoMag routine /*! \brief This routine computes magnetic declination \param double lat : in: latitude in decimal degrees \param double lon : in: longitude in decimal degrees \param double hgt : in: altitude in meters */ static double geoMagGetDecRet(double lat, double lon, double hgt) { GregorianCalendar date = new GregorianCalendar(); return geoMagGetDecRet(lat, lon, hgt, date.get(Calendar.YEAR)); } }