/* DOD World Magnetic Model 2005-2010 

See bottom of file for info from original Fortran source code 

Reminder - this is just a model! See the USGS and BGS sites for predicted accuracy.

Usage is 

    var wmm = new WorldMagneticModel();

    then
    
    var dec = wmm.declination(0.0, 59.0, -2.0, 2008.5); 
    
    parameters are 
        altitude in kilometres
        decimal degree latitude, -ve for south
        decimal degree longitude, -ve for west
        date as year + fraction of year
        
    return is declination angle in decimal degrees, 
    +ve for Magnetic North East of True North
    (-999 for < 2005.0 and >= 2010.0)
    
A new set of coefficients for 2010-2015 should be available from USGS/BGS in late 2009.
        
The method knownAnswerTest yields a maximum declination error of 0.12% on the USGS test data set. 
This is assumed to be down to the use of double rather than single precision floating point arithmetic. 
The maximum error is at 0 latitude, 120W longitude. The value produced by this code at this point is 9.191 degrees.
This is the same answer as produced by the BGS calculator at http://www.geomag.bgs.ac.uk/gifs/wmm_calc.html

This javascript port by Bill Chadwick, 27-Oct-2008 
email: w.chadwick<at>sky.com

*/

function WorldMagneticModel (){

/* 2005 - 2010 coefficients from WMM.COF */

this.coff = [
"  1,  0,  -29556.8,       0.0,        8.0,        0.0",
"  1,  1,   -1671.7,    5079.8,       10.6,      -20.9",
"  2,  0,   -2340.6,       0.0,      -15.1,        0.0",
"  2,  1,    3046.9,   -2594.7,       -7.8,      -23.2",
"  2,  2,    1657.0,    -516.7,       -0.8,      -14.6",
"  3,  0,    1335.4,       0.0,        0.4,        0.0",
"  3,  1,   -2305.1,    -199.9,       -2.6,        5.0",
"  3,  2,    1246.7,     269.3,       -1.2,       -7.0",
"  3,  3,     674.0,    -524.2,       -6.5,       -0.6",
"  4,  0,     919.8,       0.0,       -2.5,        0.0",
"  4,  1,     798.1,     281.5,        2.8,        2.2",
"  4,  2,     211.3,    -226.0,       -7.0,        1.6",
"  4,  3,    -379.4,     145.8,        6.2,        5.8",
"  4,  4,     100.0,    -304.7,       -3.8,        0.1",
"  5,  0,    -227.4,       0.0,       -2.8,        0.0",
"  5,  1,     354.6,      42.4,        0.7,        0.0",
"  5,  2,     208.7,     179.8,       -3.2,        1.7",
"  5,  3,    -136.5,    -123.0,       -1.1,        2.1",
"  5,  4,    -168.3,     -19.5,        0.1,        4.8",
"  5,  5,     -14.1,     103.6,       -0.8,       -1.1",
"  6,  0,      73.2,       0.0,       -0.7,        0.0",
"  6,  1,      69.7,     -20.3,        0.4,       -0.6",
"  6,  2,      76.7,      54.7,       -0.3,       -1.9",
"  6,  3,    -151.2,      63.6,        2.3,       -0.4",
"  6,  4,     -14.9,     -63.4,       -2.1,       -0.5",
"  6,  5,      14.6,      -0.1,       -0.6,       -0.3",
"  6,  6,     -86.3,      50.4,        1.4,        0.7",
"  7,  0,      80.1,       0.0,        0.2,        0.0",
"  7,  1,     -74.5,     -61.5,       -0.1,        0.6",
"  7,  2,      -1.4,     -22.4,       -0.3,        0.4",
"  7,  3,      38.5,       7.2,        1.1,        0.2",
"  7,  4,      12.4,      25.4,        0.6,        0.3",
"  7,  5,       9.5,      11.0,        0.5,       -0.8",
"  7,  6,       5.7,     -26.4,       -0.4,       -0.2",
"  7,  7,       1.8,      -5.1,        0.6,        0.1",
"  8,  0,      24.9,       0.0,        0.1,        0.0",
"  8,  1,       7.7,      11.2,        0.3,       -0.2",
"  8,  2,     -11.6,     -21.0,       -0.4,        0.1",
"  8,  3,      -6.9,       9.6,        0.3,        0.3",
"  8,  4,     -18.2,     -19.8,       -0.3,        0.4",
"  8,  5,      10.0,      16.1,        0.2,        0.1",
"  8,  6,       9.2,       7.7,        0.4,       -0.2",
"  8,  7,     -11.6,     -12.9,       -0.7,        0.4",
"  8,  8,      -5.2,      -0.2,        0.4,        0.4",
"  9,  0,       5.6,       0.0,        0.0,        0.0",
"  9,  1,       9.9,     -20.1,        0.0,        0.0",
"  9,  2,       3.5,      12.9,        0.0,        0.0",
"  9,  3,      -7.0,      12.6,        0.0,        0.0",
"  9,  4,       5.1,      -6.7,        0.0,        0.0",
"  9,  5,     -10.8,      -8.1,        0.0,        0.0",
"  9,  6,      -1.3,       8.0,        0.0,        0.0",
"  9,  7,       8.8,       2.9,        0.0,        0.0",
"  9,  8,      -6.7,      -7.9,        0.0,        0.0",
"  9,  9,      -9.1,       6.0,        0.0,        0.0",
" 10,  0,      -2.3,       0.0,        0.0,        0.0",
" 10,  1,      -6.3,       2.4,        0.0,        0.0",
" 10,  2,       1.6,       0.2,        0.0,        0.0",
" 10,  3,      -2.6,       4.4,        0.0,        0.0",
" 10,  4,       0.0,       4.8,        0.0,        0.0",
" 10,  5,       3.1,      -6.5,        0.0,        0.0",
" 10,  6,       0.4,      -1.1,        0.0,        0.0",
" 10,  7,       2.1,      -3.4,        0.0,        0.0",
" 10,  8,       3.9,      -0.8,        0.0,        0.0",
" 10,  9,      -0.1,      -2.3,        0.0,        0.0",
" 10, 10,      -2.3,      -7.9,        0.0,        0.0",
" 11,  0,       2.8,       0.0,        0.0,        0.0",
" 11,  1,      -1.6,       0.3,        0.0,        0.0",
" 11,  2,      -1.7,       1.2,        0.0,        0.0",
" 11,  3,       1.7,      -0.8,        0.0,        0.0",
" 11,  4,      -0.1,      -2.5,        0.0,        0.0",
" 11,  5,       0.1,       0.9,        0.0,        0.0",
" 11,  6,      -0.7,      -0.6,        0.0,        0.0",
" 11,  7,       0.7,      -2.7,        0.0,        0.0",
" 11,  8,       1.8,      -0.9,        0.0,        0.0",
" 11,  9,       0.0,      -1.3,        0.0,        0.0",
" 11, 10,       1.1,      -2.0,        0.0,        0.0",
" 11, 11,       4.1,      -1.2,        0.0,        0.0",
" 12,  0,      -2.4,       0.0,        0.0,        0.0",
" 12,  1,      -0.4,      -0.4,        0.0,        0.0",
" 12,  2,       0.2,       0.3,        0.0,        0.0",
" 12,  3,       0.8,       2.4,        0.0,        0.0",
" 12,  4,      -0.3,      -2.6,        0.0,        0.0",
" 12,  5,       1.1,       0.6,        0.0,        0.0",
" 12,  6,      -0.5,       0.3,        0.0,        0.0",
" 12,  7,       0.4,       0.0,        0.0,        0.0",
" 12,  8,      -0.3,       0.0,        0.0,        0.0",
" 12,  9,      -0.3,       0.3,        0.0,        0.0",
" 12, 10,      -0.1,      -0.9,        0.0,        0.0",
" 12, 11,      -0.3,      -0.4,        0.0,        0.0",
" 12, 12,      -0.1,       0.8,        0.0,        0.0"
];

/* static variables */

/* some 13x13 2D arrays */
this.c = new Array(13); 
this.cd = new Array(13); 
this.tc = new Array(13); 
this.dp = new Array(13); 
this.k = new Array(13); 

for(var i=0; i< 13; i++)
{
    this.c[i] = new Array(13); 
    this.cd[i] = new Array(13); 
    this.tc[i] = new Array(13); 
    this.dp[i] = new Array(13); 
    this.k[i] = new Array(13); 
}

/* some 1D arrays */
this.snorm = new Array(169); 
this.sp = new Array(13); 
this.cp = new Array(13); 
this.fn = new Array(13); 
this.fm = new Array(13); 
this.pp = new Array(13); 


/* locals */

var maxdeg = 12;
var maxord;
var i,j,D1,D2,n,m;
var a,b,a2,b2,c2,a4,b4,c4,re;
var gnm,hnm,dgnm,dhnm,flnmj;
var c_str;
var c_flds;

/* INITIALIZE CONSTANTS */

maxord = maxdeg;
this.sp[0] = 0.0;
this.cp[0] = this.snorm[0] = this.pp[0] = 1.0;
this.dp[0][0] = 0.0;
a = 6378.137;
b = 6356.7523142;
re = 6371.2;
a2 = a*a;
b2 = b*b;
c2 = a2-b2;
a4 = a2*a2;
b4 = b2*b2;
c4 = a4 - b4;

/* READ WORLD MAGNETIC MODEL SPHERICAL HARMONIC COEFFICIENTS */
this.c[0][0] = 0.0;
this.cd[0][0] = 0.0;

for(i=0; i<this.coff.length; i++)
{
  c_str = this.coff[i];
  c_flds = c_str.split(",");

  n = parseInt(c_flds[0]);   
  m = parseInt(c_flds[1]);   
  gnm = parseFloat(c_flds[2]);
  hnm = parseFloat(c_flds[3]);
  dgnm = parseFloat(c_flds[4]);
  dhnm = parseFloat(c_flds[5]);
     
  if (m <= n)
  {
	this.c[m][n] = gnm;
	this.cd[m][n] = dgnm;
	if (m != 0) 
	{
	  this.c[n][m-1] = hnm;
	  this.cd[n][m-1] = dhnm;
	}
  }
}

/* CONVERT SCHMIDT NORMALIZED GAUSS COEFFICIENTS TO UNNORMALIZED */

this.snorm[0] = 1.0;
for (n=1; n<=maxord; n++) 
{
    this.snorm[n] = this.snorm[n-1]*(2*n-1)/n;
    j = 2;
    for (m=0,D1=1,D2=(n-m+D1)/D1; D2>0; D2--,m+=D1) 
    {
        this.k[m][n] = (((n-1)*(n-1))-(m*m))/((2*n-1)*(2*n-3));
        if (m > 0) 
        {
        flnmj = ((n-m+1)*j)/(n+m);
        this.snorm[n+m*13] = this.snorm[n+(m-1)*13]*Math.sqrt(flnmj);
        j = 1;
        this.c[n][m-1] = this.snorm[n+m*13]*this.c[n][m-1];
        this.cd[n][m-1] = this.snorm[n+m*13]*this.cd[n][m-1];
        }
        this.c[m][n] = this.snorm[n+m*13]*this.c[m][n];
        this.cd[m][n] = this.snorm[n+m*13]*this.cd[m][n];
    }
    this.fn[n] = (n+1);
    this.fm[n] = n;
}
this.k[1][1] = 0.0;
this.fm[0] = 0.0;// !!!!!! WMM C and Fortran both have a bug in that fm[0] is not initialised 

}

WorldMagneticModel.prototype.declination = function(altitudeKm, latitudeDegrees, longitudeDegrees, yearFloat){
		


/* locals */

var a = 6378.137;
var b = 6356.7523142;
var re = 6371.2;
var a2 = a*a;
var b2 = b*b;
var c2 = a2-b2;
var a4 = a2*a2;
var b4 = b2*b2;
var c4 = a4 - b4;
var D3, D4;
var dip, ti, gv, dec;
var n,m;
	  
var pi, dt, rlon, rlat, srlon, srlat, crlon, crlat,srlat2,
crlat2, q, q1, q2, ct, d,aor,ar, br, r2, bpp, par,
temp1,parp,temp2,bx,by,bz,bh,dtr,bp,bt, st,ca,sa;

var maxord = 12;
var alt = altitudeKm;
var glon = longitudeDegrees;
var glat = latitudeDegrees;

/*************************************************************************/

dt = yearFloat - 2005.0;
//if more then 6 years has passed since last epoch update then return invalid
if ((dt < 0.0) || (dt > 5.0)) 
    return -999;


pi = 3.14159265359;
dtr = pi/180.0;
rlon = glon*dtr;
rlat = glat*dtr;
srlon = Math.sin(rlon);
srlat = Math.sin(rlat);
crlon = Math.cos(rlon);
crlat = Math.cos(rlat);
srlat2 = srlat*srlat;
crlat2 = crlat*crlat;
this.sp[1] = srlon;
this.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++)
{
    this.sp[m] = this.sp[1]*this.cp[m-1]+this.cp[1]*this.sp[m-1];
    this.cp[m] = this.cp[1]*this.cp[m-1]-this.sp[1]*this.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,D3=1,D4=(n+m+D3)/D3; D4>0; D4--,m+=D3) 
    {
        /*
           COMPUTE UNNORMALIZED ASSOCIATED LEGENDRE POLYNOMIALS
           AND DERIVATIVES VIA RECURSION RELATIONS
        */
            
          
        if (n == m) 
        {
          this.snorm[n+m*13] = st*this.snorm[n-1+(m-1)*13];
          this.dp[m][n] = st*this.dp[m-1][n-1]+ct*this.snorm[n-1+(m-1)*13];
        }
        else if (n == 1 && m == 0) 
        {
          this.snorm[n+m*13] = ct*this.snorm[n-1+m*13];
          this.dp[m][n] = ct*this.dp[m][n-1]-st*this.snorm[n-1+m*13];
        }
        else if (n > 1 && n != m) 
        {
          if (m > n-2) this.snorm[n-2+m*13] = 0.0;
          if (m > n-2) this.dp[m][n-2] = 0.0;
          this.snorm[n+m*13] = ct*this.snorm[n-1+m*13]-this.k[m][n]*this.snorm[n-2+m*13];
          this.dp[m][n] = ct*this.dp[m][n-1] - st*this.snorm[n-1+m*13]-this.k[m][n]*this.dp[m][n-2];
         }
         
        /*
        TIME ADJUST THE GAUSS COEFFICIENTS
        */
        this.tc[m][n] = this.c[m][n]+dt*this.cd[m][n];
        if (m != 0) this.tc[n][m-1] = this.c[n][m-1]+dt*this.cd[n][m-1];
        
        /*
        ACCUMULATE TERMS OF THE SPHERICAL HARMONIC EXPANSIONS
        */
        par = ar*this.snorm[n+m*13];
        if (m == 0) 
        {
            temp1 = this.tc[m][n]*this.cp[m];
            temp2 = this.tc[m][n]*this.sp[m];
        }
        else 
        {
            temp1 = this.tc[m][n]*this.cp[m]+this.tc[n][m-1]*this.sp[m];
            temp2 = this.tc[m][n]*this.sp[m]-this.tc[n][m-1]*this.cp[m];
        }
        bt = bt-ar*temp1*this.dp[m][n];
        bp += (this.fm[m]*temp2*par);
        br += (this.fn[n]*temp1*par);
        /*
        SPECIAL CASE:  NORTH/SOUTH GEOGRAPHIC POLES
        */
        if (st == 0.0 && m == 1) 
        {
            if (n == 1) this.pp[n] = this.pp[n-1];
            else this.pp[n] = this.ct*this.pp[n-1]-this.k[m][n]*this.pp[n-2];
            parp = ar*this.pp[n];
            bpp += (this.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));
ti = Math.sqrt((bh*bh)+(bz*bz));
dec = Math.atan2(by,bx)/dtr;
dip = Math.atan2(bz,bh)/dtr;
/*
    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
*/
gv = -999.0;
if (Math.abs(glat) >= 55.0) 
{
    if (glat > 0.0 && glon >= 0.0) gv = dec-glon;
    if (glat > 0.0 && glon < 0.0) gv = dec+Math.abs(glon);
    if (glat < 0.0 && glon >= 0.0) gv = dec+glon;
    if (glat < 0.0 && glon < 0.0) gv = dec-Math.abs(glon);
    if (gv > +180.0) gv -= 360.0;
    if (gv < -180.0) gv += 360.0;
}
      
return dec;
}

WorldMagneticModel.prototype.knownAnswerTest = function(){

/* http://www.ngdc.noaa.gov/geomag/WMM/data/2006TestValues_WMM2005.pdf */
/* Lat	Lon Dec	    */
/* Lon 240 = 120W, Lon 300 = 60W */
var kat = [
"80.00	,0.00	 ,-7.55	    ",	   
"80.00	,60.00	 ,35.53	    ",	   
"80.00	,120.00	 ,0.85	    ",	   
"80.00	,180.00	 ,9.32	    ",	   
"80.00	,240.00	 ,36.35	    ",	   
"80.00	,300.00	 ,-55.92	",	   
"40.00	,0.00	 ,-1.13	    ",	   
"40.00	,60.00	 ,4.98	    ",	   
"40.00	,120.00	 ,-7.27	    ",	   
"40.00	,180.00	 ,5.85	    ",	   
"40.00	,240.00	 ,14.93	    ",	   
"40.00	,300.00	 ,-17.98	",	   
"0.00	,0.00	 ,-6.60	    ",	   
"0.00	,60.00	 ,-4.20	    ",	   
"0.00	,120.00	 ,1.12	    ",	   
"0.00	,180.00	 ,9.63	    ",	   
"0.00	,240.00	 ,9.18	    ",	   
"0.00	,300.00	 ,-14.48	",	   
"-40.00	,0.00	 ,-23.47	",	   
"-40.00	,60.00	 ,-45.20	",	   
"-40.00	,120.00	 ,-3.42	    ",	   
"-40.00	,180.00	 ,21.73	    ",	   
"-40.00	,240.00	 ,22.43	    ",	   
"-40.00	,300.00	 ,-2.45	    ",	   
"-80.00	,0.00	 ,-21.65	",	   
"-80.00	,60.00	 ,-74.13	",	   
"-80.00	,120.00	 ,-140.50	",	   
"-80.00	,180.00	 ,131.73	",	   
"-80.00	,240.00	 ,70.50	    ",	   
"-80.00	,300.00	 ,23.95	    "
];

var maxErr = 0.0;

for(var i=0; i<kat.length; i++){

  var c_str = kat[i];
  var c_flds = c_str.split(",");
 
  var lat = parseFloat(c_flds[0]);
  var lon = parseFloat(c_flds[1]);
  var exp = parseFloat(c_flds[2]);
  var maxExp;

  var dec = this.declination(0,lat,lon,2006.0); 
  if(Math.abs(dec-exp) > maxErr){
      maxErr = Math.abs(dec-exp);
      maxExp = exp;
  } 

}

return maxErr * 100 / maxExp;//max % error

}

/*

C***********************************************************************
C
C
C     SUBROUTINE GEOMAG (GEOMAGNETIC FIELD COMPUTATION)
C
C
C***********************************************************************
C
C     GEOMAG IS A NATIONAL GEOSPATIAL INTELLIGENCE AGENCY (NGA) STANDARD
C     PRODUCT.  IT IS COVERED UNDER NGA MILITARY SPECIFICATION:
C     MIL-W-89500 (1993).
C
C***********************************************************************
C     Contact Information
C
C     Software and Model Support
C     	National Geophysical Data Center
C     	NOAA EGC/2
C     	325 Broadway
C     	Boulder, CO 80303 USA
C     	Attn: Susan McLean or Stefan Maus
C     	Phone:  (303) 497-6478 or -6522
C     	Email:  Susan.McLean@noaa.gov or Stefan.Maus@noaa.gov
C		Web: http://www.ngdc.noaa.gov/seg/WMM/
C
C     Sponsoring Government Agency
C	   National Geospatial-Intelligence Agency
C    	   PRG / CSAT, M.S. L-41
C    	   3838 Vogel Road
C    	   Arnold, MO 63010
C    	   Attn: Craig Rollins
C    	   Phone:  (314) 263-4186
C    	   Email:  Craig.M.Rollins@Nga.Mil
C
C      Original Program By:
C        Dr. John Quinn
C        FLEET PRODUCTS DIVISION, CODE N342
C        NAVAL OCEANOGRAPHIC OFFICE (NAVOCEANO)
C        STENNIS SPACE CENTER (SSC), MS 39522-5001
C
C***********************************************************************
C
C     PURPOSE:  THIS ROUTINE COMPUTES THE DECLINATION (DEC),
C               INCLINATION (DIP), TOTAL INTENSITY (TI) AND
C               GRID VARIATION (GV - POLAR REGIONS ONLY, REFERENCED
C               TO GRID NORTH OF A STEREOGRAPHIC PROJECTION) OF THE
C               EARTH'S MAGNETIC FIELD IN GEODETIC COORDINATES
C               FROM THE COEFFICIENTS OF THE CURRENT OFFICIAL
C               DEPARTMENT OF DEFENSE (DOD) SPHERICAL HARMONIC WORLD
C               MAGNETIC MODEL (WMM.COF).  THE WMM SERIES OF MODELS IS
C               UPDATED EVERY 5 YEARS ON JANUARY 1ST OF THOSE YEARS
C               WHICH ARE DIVISIBLE BY 5 (I.E. 2000, 2005, 2010 ETC.)
C               BY NOAA'S NATIONAL GEOPHYSICAL DATA CENTER IN
C               COOPERATION WITH THE BRITISH GEOLOGICAL SURVEY (BGS).
C               THE MODEL IS BASED ON GEOMAGNETIC FIELD MEASUREMENTS
C               FROM SATELLITE AND GROUND OBSERVATORIES.
C
C***********************************************************************
C
C     MODEL:  THE WMM SERIES GEOMAGNETIC MODELS ARE COMPOSED
C             OF TWO PARTS:  THE MAIN FIELD MODEL, WHICH IS
C             VALID AT THE BASE EPOCH OF THE CURRENT MODEL AND
C             A SECULAR VARIATION MODEL, WHICH ACCOUNTS FOR SLOW
C             TEMPORAL VARIATIONS IN THE MAIN GEOMAGNETIC FIELD
C             FROM THE BASE EPOCH TO A MAXIMUM OF 5 YEARS BEYOND
C             THE BASE EPOCH.  FOR EXAMPLE, THE BASE EPOCH OF
C             THE WMM-2005 MODEL IS 2005.0.  THIS MODEL IS THEREFORE
C             CONSIDERED VALID BETWEEN 2005.0 AND 2010.0. THE
C             COMPUTED MAGNETIC PARAMETERS ARE REFERENCED TO THE
C             WGS-84 ELLIPSOID.
C
C***********************************************************************
C
C     ACCURACY:  IN OCEAN AREAS AT THE EARTH'S SURFACE OVER THE
C                ENTIRE 5 YEAR LIFE OF THE DEGREE AND ORDER 12
C                SPHERICAL HARMONIC MODEL WMM-2005, THE ESTIMATED
C                MAXIMUM RMS ERRORS FOR THE VARIOUS MAGNETIC COMPONENTS
C                ARE:
C
C                DEC  -   0.5 Degrees
C                DIP  -   0.5 Degrees
C                TI   - 280.0 nanoTeslas (nT)
C                GV   -   0.5 Degrees
C
C                OTHER MAGNETIC COMPONENTS THAT CAN BE DERIVED FROM
C                THESE FOUR BY SIMPLE TRIGONOMETRIC RELATIONS WILL
C                HAVE THE FOLLOWING APPROXIMATE ERRORS OVER OCEAN AREAS:
C
C                X    - 140 nT (North)
C                Y    - 140 nT (East)
C                Z    - 200 nT (Vertical) Positive is down
C                H    - 200 nT (Horizontal)
C
C                OVER LAND THE MAXIMUM RMS ERRORS ARE EXPECTED TO BE
C                HIGHER, ALTHOUGH THE RMS ERRORS FOR DEC, DIP, AND GV
C                ARE STILL ESTIMATED TO BE LESS THAN 1.0 DEGREE, FOR
C                THE ENTIRE 5-YEAR LIFE OF THE MODEL AT THE EARTH's
C                SURFACE.  THE OTHER COMPONENT ERRORS OVER LAND ARE
C                MORE DIFFICULT TO ESTIMATE AND SO ARE NOT GIVEN.
C
C                THE ACCURACY AT ANY GIVEN TIME FOR ALL OF THESE
C                GEOMAGNETIC PARAMETERS DEPENDS ON THE GEOMAGNETIC
C                LATITUDE.  THE ERRORS ARE LEAST FROM THE EQUATOR TO
C                MID-LATITUDES AND GREATEST NEAR THE MAGNETIC POLES.
C
C                IT IS VERY IMPORTANT TO NOTE THAT A DEGREE AND
C                ORDER 12 MODEL, SUCH AS WMM-2005, DESCRIBES ONLY
C                THE LONG WAVELENGTH SPATIAL MAGNETIC FLUCTUATIONS
C                DUE TO EARTH'S CORE.  NOT INCLUDED IN THE WMM SERIES
C                MODELS ARE INTERMEDIATE AND SHORT WAVELENGTH
C                SPATIAL FLUCTUATIONS OF THE GEOMAGNETIC FIELD
C                WHICH ORIGINATE IN THE EARTH'S MANTLE AND CRUST.
C                CONSEQUENTLY, ISOLATED ANGULAR ERRORS AT VARIOUS
C                POSITIONS ON THE SURFACE (PRIMARILY OVER LAND, IN
C                CONTINENTAL MARGINS AND OVER OCEANIC SEAMOUNTS,
C                RIDGES AND TRENCHES) OF SEVERAL DEGREES MAY BE
C                EXPECTED. ALSO NOT INCLUDED IN THE MODEL ARE
C                NONSECULAR TEMPORAL FLUCTUATIONS OF THE GEOMAGNETIC
C                FIELD OF MAGNETOSPHERIC AND IONOSPHERIC ORIGIN.
C                DURING MAGNETIC STORMS, TEMPORAL FLUCTUATIONS CAN
C                CAUSE SUBSTANTIAL DEVIATIONS OF THE GEOMAGNETIC
C                FIELD FROM MODEL VALUES.  IN ARCTIC AND ANTARCTIC
C                REGIONS, AS WELL AS IN EQUATORIAL REGIONS, DEVIATIONS
C                FROM MODEL VALUES ARE BOTH FREQUENT AND PERSISTENT.
C
C                IF THE REQUIRED DECLINATION ACCURACY IS MORE
C                STRINGENT THAN THE WMM SERIES OF MODELS PROVIDE, THEN
C                THE USER IS ADVISED TO REQUEST SPECIAL (REGIONAL OR
C                LOCAL) SURVEYS BE PERFORMED AND MODELS PREPARED.
C                REQUESTS OF THIS NATURE SHOULD BE MADE TO NIMA
C                AT THE ADDRESS ABOVE.
C
C***********************************************************************
C
C     USAGE:  THIS ROUTINE IS BROKEN UP INTO TWO PARTS:
C
C             A) AN INITIALIZATION MODULE, WHICH IS CALLED ONLY
C                ONCE AT THE BEGINNING OF THE MAIN (CALLING)
C                PROGRAM
C             B) A PROCESSING MODULE, WHICH COMPUTES THE MAGNETIC
C                FIELD PARAMETERS FOR EACH SPECIFIED GEODETIC
C                POSITION (ALTITUDE, LATITUDE, LONGITUDE) AND TIME
C
C             INITIALIZATION IS MADE VIA A SINGLE CALL TO THE MAIN
C             ENTRY POINT (GEOMAG), WHILE SUBSEQUENT PROCESSING
C             CALLS ARE MADE THROUGH THE SECOND ENTRY POINT (GEOMG1).
C             ONE CALL TO THE PROCESSING MODULE IS REQUIRED FOR EACH
C             POSITION AND TIME.
C
C             THE VARIABLE MAXDEG IN THE INITIALIZATION CALL IS THE
C             MAXIMUM DEGREE TO WHICH THE SPHERICAL HARMONIC MODEL
C             IS TO BE COMPUTED.  IT MUST BE SPECIFIED BY THE USER
C             IN THE CALLING ROUTINE.  NORMALLY IT IS 12 BUT IT MAY
C             BE SET LESS THAN 12 TO INCREASE COMPUTATIONAL SPEED AT
C             THE EXPENSE OF REDUCED ACCURACY.
C
C             THE PC VERSION OF THIS SUBROUTINE MUST BE COMPILED
C             WITH A FORTRAN 77 COMPATIBLE COMPILER SUCH AS THE
C             MICROSOFT OPTIMIZING FORTRAN COMPILER VERSION 4.1
C             OR LATER.
C
C**********************************************************************
C
C     REFERENCES:
C
C       JOHN M. QUINN, DAVID J. KERRIDGE AND DAVID R. BARRACLOUGH,
C            WORLD MAGNETIC CHARTS FOR 1985 - SPHERICAL HARMONIC
C            MODELS OF THE GEOMAGNETIC FIELD AND ITS SECULAR
C            VARIATION, GEOPHYS. J. R. ASTR. SOC. (1986) 87,
C            PP 1143-1157
C
C       DEFENSE MAPPING AGENCY TECHNICAL REPORT, TR 8350.2:
C            DEPARTMENT OF DEFENSE WORLD GEODETIC SYSTEM 1984,
C            SEPT. 30 (1987)
C
C       JOHN M. QUINN, RACHEL J. COLEMAN, MICHAEL R. PECK, AND
C            STEPHEN E. LAUBER; THE JOINT US/UK 1990 EPOCH
C            WORLD MAGNETIC MODEL, TECHNICAL REPORT NO. 304,
C            NAVAL OCEANOGRAPHIC OFFICE (1991)
C
C       JOHN M. QUINN, RACHEL J. COLEMAN, DONALD L. SHIEL, AND
C            JOHN M. NIGRO; THE JOINT US/UK 1995 EPOCH WORLD
C            MAGNETIC MODEL, TECHNICAL REPORT NO. 314, NAVAL
C            OCEANOGRAPHIC OFFICE (1995)
C
C            SUSAN AMCMILLAN, DAVID R. BARRACLOUGH, JOHN M. QUINN, AND
C            RACHEL J. COLEMAN;  THE 1995 REVISION OF THE JOINT US/UK
C            GEOMAGNETIC FIELD MODELS - I. SECULAR VARIATION, JOURNAL OF
C            GEOMAGNETISM AND GEOELECTRICITY, VOL. 49, PP. 229-243
C            (1997)
C
C            JOHN M. QUINN, RACHEL J. COELMAN, SUSAM MACMILLAN, AND
C            DAVID R. BARRACLOUGH;  THE 1995 REVISION OF THE JOINT
C            US/UK GEOMAGNETIC FIELD MODELS: II. MAIN FIELD,JOURNAL OF
C            GEOMAGNETISM AND GEOELECTRICITY, VOL. 49, PP. 245 - 261
C            (1997)
C
C***********************************************************************
C
C     PARAMETER DESCRIPTIONS:
C
C       A      - SEMIMAJOR AXIS OF WGS-84 ELLIPSOID (KM)
C       B      - SEMIMINOR AXIS OF WGS-84 ELLIPSOID (KM)
C       RE     - MEAN RADIUS OF IAU-66 ELLIPSOID (KM)
C       SNORM  - SCHMIDT NORMALIZATION FACTORS
C       C      - GAUSS COEFFICIENTS OF MAIN GEOMAGNETIC MODEL (NT)
C       CD     - GAUSS COEFFICIENTS OF SECULAR GEOMAGNETIC MODEL (NT/YR)
C       TC     - TIME ADJUSTED GEOMAGNETIC GAUSS COEFFICIENTS (NT)
C       OTIME  - TIME ON PREVIOUS CALL TO GEOMAG (YRS)
C       OALT   - GEODETIC ALTITUDE ON PREVIOUS CALL TO GEOMAG (YRS)
C       OLAT   - GEODETIC LATITUDE ON PREVIOUS CALL TO GEOMAG (DEG.)
C       TIME   - COMPUTATION TIME (YRS)                        (INPUT)
C                (EG. 1 JULY 1995 = 1995.500)
C       ALT    - GEODETIC ALTITUDE (KM)                        (INPUT)
C       GLAT   - GEODETIC LATITUDE (DEG.)                      (INPUT)
C       GLON   - GEODETIC LONGITUDE (DEG.)                     (INPUT)
C       EPOCH  - BASE TIME OF GEOMAGNETIC MODEL (YRS)
C       DTR    - DEGREE TO RADIAN CONVERSION
C       SP(M)  - SINE OF (M*SPHERICAL COORD. LONGITUDE)
C       CP(M)  - COSINE OF (M*SPHERICAL COORD. LONGITUDE)
C       ST     - SINE OF (SPHERICAL COORD. LATITUDE)
C       CT     - COSINE OF (SPHERICAL COORD. LATITUDE)
C       R      - SPHERICAL COORDINATE RADIAL POSITION (KM)
C       CA     - COSINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE
C       SA     - SINE OF SPHERICAL TO GEODETIC VECTOR ROTATION ANGLE
C       BR     - RADIAL COMPONENT OF GEOMAGNETIC FIELD (NT)
C       BT     - THETA COMPONENT OF GEOMAGNETIC FIELD (NT)
C       BP     - PHI COMPONENT OF GEOMAGNETIC FIELD (NT)
C       P(N,M) - ASSOCIATED LEGENDRE POLYNOMIALS (UNNORMALIZED)
C       PP(N)  - ASSOCIATED LEGENDRE POLYNOMIALS FOR M=1 (UNNORMALIZED)
C       DP(N,M)- THETA DERIVATIVE OF P(N,M) (UNNORMALIZED)
C       BX     - NORTH GEOMAGNETIC COMPONENT (NT)
C       BY     - EAST GEOMAGNETIC COMPONENT (NT)
C       BZ     - VERTICALLY DOWN GEOMAGNETIC COMPONENT (NT)
C       BH     - HORIZONTAL GEOMAGNETIC COMPONENT (NT)
C       DEC    - GEOMAGNETIC DECLINATION (DEG.)                (OUTPUT)
C                  EAST=POSITIVE ANGLES
C                  WEST=NEGATIVE ANGLES
C       DIP    - GEOMAGNETIC INCLINATION (DEG.)                (OUTPUT)
C                  DOWN=POSITIVE ANGLES
C                    UP=NEGATIVE ANGLES
C       TI     - GEOMAGNETIC TOTAL INTENSITY (NT)              (OUTPUT)
C       GV     - GEOMAGNETIC GRID VARIATION (DEG.)             (OUTPUT)
C                REFERENCED TO GRID NORTH
C                GRID NORTH REFERENCED TO 0 MERIDIAN
C                OF A POLAR STEREOGRAPHIC PROJECTION
C                (ARCTIC/ANTARCTIC ONLY)
C       MAXDEG - MAXIMUM DEGREE OF SPHERICAL HARMONIC MODEL    (INPUT)
C       MOXORD - MAXIMUM ORDER OF SPHERICAL HARMONIC MODEL
C
C***********************************************************************
C
C     NOTE:  THIS VERSION OF GEOMAG USES A WMM SERIES GEOMAGNETIC
C            FIELS MODEL REFERENCED TO THE WGS-84 GRAVITY MODEL
C            ELLIPSOID
C


*/