/*  module magfield.c */

/* Module to calculate magnetic variation and field given position,
**               altitude, and date
** Implements the NIMA (formerly DMA) WMM and IGRF models
**
**    http://www.nima.mil/GandG/ngdc-wmm2000.html
**    For WMM2000 coefficients:
**    ftp://ftp.ngdc.noaa.gov/Solid_Earth/Mainfld_Mag/DoD_Model/wmm.cof
**    For IGRF/DGRF coefficients:
**    http://swdcdb.kugi.kyoto-u.ac.jp/igrf/coef/igrfall.d
**
** Copyright (C) 2000  Edward A Williams <Ed_Williams@compuserve.com>
**
**  The routine uses a spherical harmonic expansion of the magnetic
** potential up to twelfth order, together with its time variation, as
** described in Chapter 4 of "Geomagnetism, Vol 1, Ed. J.A.Jacobs,
** Academic Press (London 1987)". The program first converts geodetic
** coordinates (lat/long on elliptic earth and altitude) to spherical
** geocentric (spherical lat/long and radius) coordinates. Using this,
** the spherical (B_r, B_theta, B_phi) magnetic field components are
** computed from the model. These are finally referred to surface (X, Y,
** Z) coordinates.
**
**   Fields are accurate to better than 200nT, variation and dip to
** better than 0.5 degrees, with the exception of the declination near
** the magnetic poles (where it is ill-defined) where the error may reach
** 4 degrees or more.
**
**   Variation is undefined at both the geographic and  
** magnetic poles, even though the field itself is well-behaved. To
** avoid the routine blowing up, latitude entries corresponding to
** the geographic poles are slightly offset. At the magnetic poles,
** the routine returns zero variation.
**
** HISTORY
** Adapted from EAW Excel 3.0 version 3/27/94 EAW
** Recoded in C++ by Starry Chan
** WMM95 added and rearranged in ANSI-C EAW 7/9/95
** Put shell around program and made Borland & GCC compatible EAW 11/22/95
** IGRF95 added 2/96 EAW
** WMM2000 IGR2000 added 2/00 EAW
** Released under GPL  3/26/00 EAW
** Adaptions and modifications for the SimGear project  3/27/2000 CLO
** Removed all pow() calls and made static roots[][] arrays to
** save many sqrt() calls on subsequent invocations
** 3/28/2000  Norman Vine -- nhv@yahoo.com
** Put in some bullet-proofing to handle magnetic and geographic poles.
** 3/28/2000 EAW
** Added missing comment close, the lack of which caused the altitude 
** correction to be omitted.
** 01/31/01 Jim Seymour (jseymour@LinxNet.com)
**
** This program is free software; you can redistribute it and/or
** modify it under the terms of the GNU General Public License as
** published by the Free Software Foundation; either version 2 of the
** License, or (at your option) any later version.
**
** This program is distributed in the hope that it will be useful, but
** WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
** General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**
*/ 

#include	<stdio.h>
#include	<stdlib.h>
#include	<math.h>
#include        "magfield.h"

#define	max(a,b)	(((a) > (b)) ? (a) : (b))

static const double pi = 3.14159265358979;
static const double a = 6378.16;	/* major radius (km) IAU66 ellipsoid */
static const double f = 1.0 / 298.25;	/* inverse flattening IAU66 ellipsoid */
static const double b = 6378.16 * (1.0 -1.0 / 298.25 );
	/* minor radius b=a*(1-f) */
static const double r_0 = 6371.2;	/* "mean radius" for spherical harmonic expansion */

/*IGRF90 constants */
static const double gnm_igrf90[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29775.4, -1851, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-2135.8, 3058.2, 1693.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{1314.6, -2240.2, 1245.6, 806.5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{938.9, 782.3, 323.9, -422.7, 141.7, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-211, 352.5, 243.8, -110.8, -165.6, -37, 0, 0, 0, 0, 0, 0, 0}, 
	{60.7, 63.9, 60.4, -177.5, 2, 16.7, -96.3, 0, 0, 0, 0, 0, 0}, 
	{76.6, -64.2, 3.7, 27.5, 0.9, 5.7, 9.8, -0.5, 0, 0, 0, 0, 0}, 
	{22.4, 5.1, -0.9, -10.8, -12.4, 3.8, 3.8, 2.6, -6, 0, 0, 0, 0}, 
	{4.4, 9.9, 0.8, -12, 9.3, -3.9, -1.4, 7.3, 1.5, -5.5, 0, 0, 0}, 
	{-3.6, -3.9, 2.4, -5.3, -2.4, 4.4, 3, 1.2, 2.2, 2.9, 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 const double hnm_igrf90[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 5410.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -2277.7, -380, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -286.5, 293.3, -348.5, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 248.1, -239.5, 87, -299.4, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 47.2, 153.5, -154.4, -69.2, 97.7, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -15.8, 82.7, 68.3, -52.5, 1.8, 26.9, 0, 0, 0, 0, 0, 0},
	{0, -81.1, -27.3, 0.6, 20.4, 16.4, -22.6, -5, 0, 0, 0, 0, 0}, 
	{0, 9.7, -19.9, 7.1, -22.1, 11.9, 11, -16, -10.7, 0, 0, 0, 0}, 
	{0, -20.8, 15.4, 9.5, -5.7, -6.4, 8.6, 9.1, -6.6, 1.9, 0, 0, 0}, 
	{0, 1.3, 0.4, 3.1, 5.6, -4.2, -0.5, -1.5, 3.8, -0.5, -6.2, 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 const double gtnm_igrf90[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{18, 10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-12.9, 2.4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{3.3, -6.7, 0, -5.9, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0.5, 0.6, -7, 0.5, -5.5, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0.6, -0.1, -1.6, -3.1, 0, 2.3, 0, 0, 0, 0, 0, 0, 0}, 
	{1.3, -0.2, 1.8, 1.3, -0.2, 0.1, 1.2, 0, 0, 0, 0, 0, 0}, 
	{0.6, -0.5, -0.3, 0.6, 1.6, 0.2, 0.2, 0.3, 0, 0, 0, 0, 0}, 
	{0.2, -0.7, -0.2, 0.1, -1.1, 0, 0, -0.5, -0.6, 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 const double htnm_igrf90[13][13] =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -16.1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -15.8, -13.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 4.4, 1.6, -10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 2.6, 1.8, 3.1, -1.4, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -0.1, 0.5, 0.4, 1.7, 0.4, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 0.2, -1.3, 0, -0.9, 0.5, 1.2, 0, 0, 0, 0, 0, 0}, 
	{0, 0.6, 0.2, 0.8, -0.5, -0.2, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 0.5, -0.2, 0.3, 0.3, 0.4, -0.5, -0.3, 0.6, 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}};
/* IGRF95 */
static const double gnm_igrf95[13][13] =
{
  {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},
  {-29682.0, -1789.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {-2197.0, 3074.0, 1685.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {1329.0, -2268.0, 1249.0, 769.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {941.0, 782.0, 291.0, -421.0, 116.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
  {-210.0, 352.0, 237.0, -122.0, -167.0, -26.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0,0.0},
  {66.0, 64.0, 65.0, -172.0, 2.0, 17.0, -94.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {78.0, -67.0, 1.0, 29.0, 4.0, 8.0, 10.0, -2.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {24.0, 4.0, -1.0, -9.0, -14.0, 4.0, 5.0, 0.0, -7.0, 0.0,0.0, 0.0, 0.0},
  {4.0, 9.0, 1.0, -12.0, 9.0, -4.0, -2.0, 7.0, 0.0, -6.0, 0.0, 0.0, 0.0},
  {-3.0, -4.0, 2.0, -5.0, -2.0, 4.0, 3.0, 1.0, 3.0, 3.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 const double hnm_igrf95[13][13]=
{
  {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, 5318.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, -2356.0, -425.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -263.0, 302.0, -406.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 262.0, -232.0, 98.0, -301.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 44.0, 157.0, -152.0, -64.0, 99.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -16.0, 77.0, 67.0, -57.0, 4.0, 28.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -77.0, -25.0, 3.0, 22.0, 16.0, -23.0, -3.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 12.0, -20.0, 7.0, -21.0, 12.0, 10.0, -17.0, -10.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -19.0, 15.0, 11.0, -7.0, -7.0, 9.0, 7.0, -8.0, 1.0, 0.0, 0.0, 0.0},
  {0.0, 2.0, 1.0, 3.0, 6.0, -4.0, 0.0, -2.0, 3.0, -1.0, -6.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 const double gtnm_igrf95[13][13]=
{
  {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},
  {17.6, 13.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {-13.2, 3.7, -0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {1.5, -6.4, -0.2, -8.1, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.8, 0.9, -6.9, 0.5, -4.6, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.8, 0.1, -1.5, -2.0, -0.1, 2.3, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.5, -0.4, 0.6, 1.9, -0.2, -0.2, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {-0.2, -0.8, -0.6, 0.6, 1.2, 0.1, 0.2, -0.6, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.3, -0.2, 0.1, 0.4, -1.1, 0.3, 0.2, -0.9, -0.3, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0}
 };

static const double htnm_igrf95[13][13]=
{
  {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, -18.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, -15.0, -8.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 4.1, 2.2, -12.1, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 1.8, 1.2, 2.7, -1.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.2, 1.2, 0.3, 1.8, 0.9, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.3, -1.6, -0.2, -0.9, 1.0, 2.2, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.8, 0.2, 0.6, -0.4, 0.0, -0.3, 0.0, 0.0,0.0, 0.0, 0.0, 0.0},
  {0.0, 0.4, -0.2, 0.2, 0.7, 0.0, -1.2, -0.7, -0.6, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0}
};

/*WMM85 constants */
static const double gnm_wmm85[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29879.8, -1903.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-2070.6, 3040.8, 1696.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{1303.9, -2203, 1241.7, 839.4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{933.8, 781.8, 359, -424.5, 164.5, 0, 0, 0, 0, 0, 0, 0, 0},
	{-216.4, 353, 254.3, -93.7, -157.5, -45.2, 0, 0, 0, 0, 0, 0, 0}, 
	{53.2, 63.8, 51.3, -188.4, 3.3, 20.3, -101.7, 0, 0, 0, 0, 0, 0}, 
	{76.9, -60.7, 0.7, 25.4, -8.1, 6.9, 7, -4.4, 0, 0, 0, 0, 0}, 
	{18.4, 5.1, 1.2, -12, -9.1, 0.1, 4.7, 6.5, -9.5, 0, 0, 0, 0}, 
	{5.7, 10.9, 0.9, -12.2, 9.5, -3.3, -1, 6.5, 1.5, -4.8, 0, 0, 0}, 
	{-3.4, -4.7, 2.5, -5.5, -2.1, 4.6, 3.2, 0.6, 1.9, 2.8, -0.2, 0, 0}, 
	{2.3, -0.8, -2, 2.1, 0.2, -0.4, -0.4, 1.6, 1.5, -0.7, 2.3, 3.5, 0}, 
	{-1.8, 0, 0.1, -0.3, 0.5, 0.5, -0.6, -0.4, 0, -0.5, 0, 0.7, -0.2}};
static const double hnm_wmm85[13][13] = 
	{
	{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 5490.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -2189.1, -312, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -310.3, 282.6, -299.2, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 227.2, -246.7, 72.5, -299.1, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 43.4, 148.2, -154.8, -71.8, 91.5, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -12.3, 87.9, 67.8, -51.1, -4, 20.8, 0, 0, 0, 0, 0, 0}, 
	{0, -80.1, -25.9, -0.9, 21.6, 18.5, -20, -7.7, 0, 0, 0, 0, 0}, 
	{0, 3.8, -20.2, 5, -24.2, 12.2, 7.6, -16.3, -10.9, 0, 0, 0, 0}, 
	{0, -20.8, 15.8, 9, -5, -6.4, 9.1, 9.9, -5.8, 2.3, 0, 0, 0}, 
	{0, 1.2, 0.4, 2.5, 5.6, -4.4, -0.5, -1.6, 3.7, -0.5, -6.1, 0, 0}, 
	{0, 1.3, 2, -1.1, -2.8, 0.7, -0.1, -2.4, -0.4, -1.5, -1.5, 0.7, 0},
	{0, 0.3, 0.6, 2.5, -1.7, 0.3, 0.2, -0.1, 0.1, 0.1, -1.4, 0.4, 0.7}};
static const double gtnm_wmm85[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{21.9, 10.6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-11.2, 1.8, 9.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{8.3, -2, -0.6, 2.4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-1.2, 0.1, -9.7, -1.7, -9.3, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{1.4, -0.5, -1.2, -2.2, 0.9, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{3.1, 0, 1.8, -0.2, -0.4, 2.4, 1.8, 0, 0, 0, 0, 0, 0}, 
	{-0.1, -0.8, -1.2, 1.1, 0, 0.6, -1.8, -1.2, 0, 0, 0, 0, 0},
	{0.2, 0, 0.7, 0.1, 0.2, -0.3, -0.1, 0.2, -2.2, 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 const double htnm_wmm85[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -31.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -9.7, -19.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 6.1, 1.3, -13, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 1.3, 3.6, 2.5, 0.6, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -0.9, 0.6, 0.3, 2.4, -1.4, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 0.7, -2.1, -1.4, -4.3, -0.7, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 0, 1.2, 2, 2.6, 0.9, 0.8, 0.4, 0, 0, 0, 0, 0}, 
	{0, -0.6, -1.5, 0.1, -1.1, 0.4, -2, 0.9, 1.5, 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}};

/* wmm90 constants */
static const double gnm_wmm90[13][13] =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29780.5, -1851.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-2134.3, 3062.2, 1691.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{1312.9, -2244.7, 1246.8, 808.6, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{933.5, 784.9, 323.5, -421.7, 139.2, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-208.3, 352.2, 246.5, -110.8, -162.3, -37.2, 0, 0, 0, 0, 0, 0, 0}, 
	{59, 63.7, 60, -181.3, 0.4, 15.4, -96, 0, 0, 0, 0, 0, 0}, 
	{76.1, -62.1, 1.3, 30.2, 4.7, 7.9, 10.1, 1.9, 0, 0, 0, 0, 0}, 
	{22.9, 2.3, -1.2, -11.7, -17.5, 2.2, 5.7, 3, -7, 0, 0, 0, 0}, 
	{3.6, 9.5, -0.9, -10.7, 10.7, -3.2, -1.4, 6.3, 0.8, -5.5, 0, 0, 0},
	{-3.3, -2.6, 4.5, -5.6, -3.6, 3.9, 3.2, 1.7, 3, 3.7, 0.7, 0, 0}, 
	{1.3, -1.4, -2.5, 3.2, 0.2, -1.1, 0.3, -0.3, 0.9, -1.1, 2.4, 3, 0}, 
	{-1.3, 0.1, 0.5, 0.7, 0.4, -0.2, -1.1, 0.9, -0.6, 0.8, 0.2, 0.4, 0.2}};
static const double hnm_wmm90[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, 5407.2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -2278.3, -384.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -284.9, 291.7, -352.4, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 249.4, -232.7, 91.3, -296.5, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 40.8, 148.7, -154.6, -67.6, 97.4, 0, 0, 0, 0, 0, 0, 0},
	{0, -14.7, 82.2, 70, -56.2, -1.4, 24.6, 0, 0, 0, 0, 0, 0}, 
	{0, -78.6, -26.7, 0.1, 19.9, 17.9, -21.5, -6.8, 0, 0, 0, 0, 0}, 
	{0, 9.7, -19.3, 6.6, -20.1, 13.4, 9.8, -19, -9.1, 0, 0, 0, 0}, 
	{0, -21.9, 14.3, 9.5, -6.7, -6.4, 9.1, 8.9, -8, 2.1, 0, 0, 0}, 
	{0, 2.6, 1.2, 2.6, 5.7, -4, -0.4, -1.7, 3.8, -0.8, -6.5, 0, 0}, 
	{0, 0, 1, -1.6, -2.2, 1.1, -0.7, -1.7, -1.5, -1.3, -1.1, 0.6, 0}, 
	{0, 0.7, 0.7, 1.3, -1.5, 0.3, 0.2, -1.1, 1.2, -0.2, -1.3, 0.6, 0.6}};
static const double gtnm_wmm90[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{16, 9.3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-11.7, 3.7, 1.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{2.1, -7.6, 0, -5.8, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{-0.8, 1, -7.4, 0.8, -6.4, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{1.7, 0, 0, -2.7, 0, 3, 0, 0, 0, 0, 0, 0, 0}, 
	{0.8, 0, 1.5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0.5, 0, -0.9, 1.5, 2.7, -1, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -1.1, 0, 0, -2.1, 0, 1, 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 const double htnm_wmm90[13][13] = 
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0, -13.8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -12.8, -14.9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 3.1, 0.8, -11.3, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 3.3, 3.7, 2.8, 0, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 0, -2.1, 1.2, 1.2, 0.6, 0, 0, 0, 0, 0, 0, 0}, 
	{0, -0.6, -0.6, 0, -2.3, 0, 0, 0, 0, 0, 0, 0, 0}, 
	{0, 0.6, 0.8, 0, 0, 0, 0.4, 0, 0, 0, 0, 0, 0}, 
	{0, 0.4, -0.8, 0.5, 0.3, 0.5, 0, -0.7, 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}};
/* wmm95 constants */
static const double gnm_wmm95[13][13] =
        {
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{-29682.1,-1782.2,0,0,0,0,0,0,0,0,0,0,0},
        {-2194.7,3078.6,1685.7,0,0,0,0,0,0,0,0,0,0},
	{1318.8,-2273.6,1246.9,766.3,0,0,0,0,0,0,0,0,0},
        {940,782.9,290.9,-418.9,113.8,0,0,0,0,0,0,0,0},
        {-209.5,354,238.2,-122.1,-162.8,-23.3,0,0,0,0,0,0,0},
        {68.5,65.6,64.1,-169.1,-0.5,16.5,-91,0,0,0,0,0,0},
        {78,-68.1,0.1,29.6,6,8.7,9.2,-2.4,0,0,0,0,0},
        {24.7,3.4,-1.5,-9.6,-16.5,2.6,3.6,-4.9,-8.5,0,0,0,0},
        {2.9,7.5,0.4,-10.3,9.7,-2.3,-2.4,6.8,-0.5,-6.5,0,0,0},
        {-2.9,-3.3,2.8,-4.3,-3.1,2.4,2.8,0.7,4.1,3.6,0.6,0,0},
        {1.7,-1.6,-3.6,1.2,-0.6,0.1,-0.7,-0.8,1.3,-0.3,2.2,4.2,0},
        {-1.8,0.9,-0.1,-0.5,0.8,0.2,0.5,0.4,-0.4,0.3,0.2,0.4,0.6}};
static const double hnm_wmm95[13][13] =
	{
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
        {0,5315.6,0,0,0,0,0,0,0,0,0,0,0},
        {0,-2359.1,-418.6,0,0,0,0,0,0,0,0,0,0},
        {0,-261.1,301,-416.5,0,0,0,0,0,0,0,0,0},
        {0,259.4,-230.9,99.8,-306.1,0,0,0,0,0,0,0,0},
        {0,43.7,157.6,-150.1,-59.2,104.4,0,0,0,0,0,0,0},
        {0,-15.2,74.3,69.4,-55.3,3,33.3,0,0,0,0,0,0},
        {0,-76.1,-24.5,1.6,20,16.5,-23.6,-6.8,0,0,0,0,0},
        {0,14.9,-19.5,6.3,-20.4,12.2,7,-19,-8.8,0,0,0,0},
        {0,-19.8,14.6,10.9,-7.5,-6.8,9.3,7.7,-8.1,2.6,0,0,0},
	{0,3.2,1.7,2.9,5.6,-3.4,-0.7,-2.9,2.3,-1.6,-6.6,0,0},
        {0,0.3,1,-3.6,-1.4,1.9,0.2,-1.3,-2.4,-0.6,-2.2,1.3,0},
        {0,0.3,1.4,0.8,-3,0.7,0.5,-0.8,0.6,0.1,-1.3,-0.4,0.9}};
static const double gtnm_wmm95[13][13] =
        {
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
        {17.6,13.2,0,0,0,0,0,0,0,0,0,0,0},
        {-13.7,4,-0.3,0,0,0,0,0,0,0,0,0,0},
        {0.8,-6.6,-0.5,-8.5,0,0,0,0,0,0,0,0,0},
        {1.2,1.1,-6.8,0.3,-4.5,0,0,0,0,0,0,0,0},
        {0.9,0.5,-1.4,-1.7,0,2.1,0,0,0,0,0,0,0},
	{0.4,-0.3,0.3,2.1,0,-0.4,-0.4,0,0,0,0,0,0},
        {-0.3,-1.1,-0.5,0.5,1.3,0.1,0,-0.9,0,0,0,0,0},
        {0.1,0,0.4,0.3,-1.3,0.5,0.4,-0.9,0.1,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 const double htnm_wmm95[13][13] =
        {
        {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
	{0,-18,0,0,0,0,0,0,0,0,0,0,0},
	{0,-14.6,-7.2,0,0,0,0,0,0,0,0,0,0},
        {0,4,2.2,-12.6,0,0,0,0,0,0,0,0,0},
        {0,1.3,1,2.5,-1.2,0,0,0,0,0,0,0,0},
        {0,0.5,1.5,0.6,1.7,0.6,0,0,0,0,0,0,0},
        {0,0.7,-1.5,-0.5,-0.7,1.1,2.6,0,0,0,0,0,0},
        {0,0.3,0,0.7,-0.6,0.1,-0.6,-0.4,0,0,0,0,0},
        {0,0.4,-0.3,0.1,0.8,-0.1,-1.3,-0.9,-1.1,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 const double gnm_wmm2000[13][13] =
        {
	{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},
	{-29616.0, -1722.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-2266.7, 3070.2, 1677.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1322.4, -2291.5, 1255.9, 724.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{932.1, 786.3, 250.6, -401.5, 106.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-211.9, 351.6, 220.8, -134.5, -168.8, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{73.8, 68.2, 74.1, -163.5, -3.8, 17.1, -85.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{77.4, -73.9, 2.2, 35.7, 7.3, 5.2, 8.4, -1.5, 0.0, 0.0, 0.0, 0.0, 0.0},
	{23.3, 7.3, -8.5, -6.6, -16.9, 8.6, 4.9, -7.8, -7.6, 0.0, 0.0, 0.0, 0.0},
	{5.7, 8.5, 2.0, -9.8, 7.6, -7.0, -2.0, 9.2, -2.2, -6.6, 0.0, 0.0, 0.0},
	{-2.2, -5.7, 1.6, -3.7, -0.6, 4.1, 2.2, 2.2, 4.6, 2.3, 0.1, 0.0, 0.0},
	{3.3, -1.1, -2.4, 2.6, -1.3, -1.7, -0.6, 0.4, 0.7, -0.3, 2.3, 4.2, 0.0},
	{-1.5, -0.2, -0.3, 0.5, 0.2, 0.9, -1.4, 0.6, -0.6, -1.0, -0.3, 0.3, 0.4},
	};

static const double hnm_wmm2000[13][13]=
        {
	{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, 5194.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -2484.8, -467.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -224.7, 293.0, -486.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 273.3, -227.9, 120.9, -302.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 42.0, 173.8, -135.0, -38.6, 105.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -17.4, 61.2, 63.2, -62.9, 0.2, 43.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -62.3, -24.5, 8.9, 23.4, 15.0, -27.6, -7.8, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 12.4, -20.8, 8.4, -21.2, 15.5, 9.1, -15.5, -5.4, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.4, 13.9, 12.0, -6.2, -8.6, 9.4, 5.0, -8.4, 3.2, 0.0, 0.0, 0.0},
	{0.0, 0.9, -0.7, 3.9, 4.8, -5.3, -1.0, -2.4, 1.3, -2.3, -6.4, 0.0, 0.0},
	{0.0, -1.5, 0.7, -1.1, -2.3, 1.3, -0.6, -2.8, -1.6, -0.1, -1.9, 1.4, 0.0},
	{0.0, -1.0, 0.7, 2.2, -2.5, -0.2, 0.0, -0.2, 0.0, 0.2, -0.9, -0.2, 1.0},
	};

static const double gtnm_wmm2000[13][13]=
        {
	{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},
	{14.7, 11.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-13.6, -0.7, -1.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.3, -4.3, 0.9, -8.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-1.6, 0.9, -7.6, 2.2, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.9, -0.2, -2.5, -2.7, -0.9, 1.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1.2, 0.2, 1.7, 1.6, -0.1, -0.3, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.4, -0.8, -0.2, 1.1, 0.4, 0.0, -0.2, -0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.3, 0.6, -0.8, 0.3, -0.2, 0.5, 0.0, -0.6, 0.1, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

static const double htnm_wmm2000[13][13]=
        {
        {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, -20.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, -21.5, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 6.4, -1.3, -13.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 2.3, 0.7, 3.7, -0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.0, 2.1, 2.3, 3.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.3, -1.7, -0.9, -1.0, -0.1, 1.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 1.4, 0.2, 0.7, 0.4, -0.3, -0.8, -0.1, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.5, 0.1, -0.2, 0.0, 0.1, -0.1, 0.3, 0.2, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
        };

static const double gnm_igrf2000[13][13] =
        {
	{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},
	{-29615.0, -1728.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-2267.0, 3072.0, 1672.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1341.0, -2290.0, 1253.0, 715.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{935.0, 787.0, 251.0, -405.0, 110.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-217.0, 351.0, 222.0, -131.0, -169.0, -12.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{72.0, 68.0, 74.0, -161.0, -5.0, 17.0, -91.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{79.0, -74.0, 0.0, 33.0, 9.0, 7.0, 8.0, -2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{25.0, 6.0, -9.0, -8.0, -17.0, 9.0, 7.0, -8.0, -7.0, 0.0, 0.0, 0.0, 0.0},
	{5.0, 9.0, 3.0, -8.0, 6.0, -9.0, -2.0, 9.0, -4.0, -8.0, 0.0, 0.0, 0.0},
	{-2.0, -6.0, 2.0, -3.0, 0.0, 4.0, 1.0, 2.0, 4.0, 0.0, -1.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 const double hnm_igrf2000[13][13]=
        {
	{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, 5186.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, -2478.0, -458.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -227.0, 296.0, -492.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 272.0, -232.0, 119.0, -304.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 44.0, 172.0, -134.0, -40.0, 107.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -17.0, 64.0, 65.0, -61.0, 1.0, 44.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -65.0, -24.0, 6.0, 24.0, 15.0, -25.0, -6.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 12.0, -22.0, 8.0, -21.0, 15.0, 9.0, -16.0, -3.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.0, 13.0, 12.0, -6.0, -8.0, 9.0, 4.0, -8.0, 5.0, 0.0, 0.0, 0.0},
	{0.0, 1.0, 0.0, 4.0, 5.0, -6.0, -1.0, -3.0, 0.0, -2.0, -8.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 const double gtnm_igrf2000[13][13]=
        {
	{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},
	{14.6, 10.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-12.4, 1.1, -1.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.7, -5.4, 0.9, -7.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-1.3, 1.6, -7.3, 2.9, -3.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.7, -2.1, -2.8, -0.8, 2.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{1.0, -0.4, 0.9, 2.0, -0.6, -0.3, 1.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.4, -0.4, -0.3, 1.1, 1.1, -0.2, 0.6, -0.9, 0.0, 0.0, 0.0, 0.0, 0.0},
	{-0.3, 0.2, -0.3, 0.4, -1.0, 0.3, -0.5, -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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};

static const double htnm_igrf2000[13][13]=
        {
	{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, -22.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -20.6, -9.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 6.0, -0.1, -14.2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 2.1, 1.3, 5.0, 0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.1, 0.6, 1.7, 1.9, 0.1, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, -0.2, -1.4, 0.0, -0.8, 0.0, 0.9, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 1.1, 0.0, 0.3, -0.1, -0.6, -0.7, 0.2, 0.0, 0.0, 0.0, 0.0, 0.0},
	{0.0, 0.1, 0.0, 0.0, 0.3, 0.6, -0.4, 0.3, 0.7, 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, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
	};



static const int nmax = 12;

double P[13][13];
double DP[13][13];
double gnm[13][13];
double hnm[13][13];
double sm[13];
double cm[13];

static double root[13];
static double roots[13][13][2];

/* Convert date to Julian day    1950-2049 */
unsigned long int yymmdd_to_julian_days(int yy,int mm,int dd)
{
  unsigned long jd;
 
  yy = (yy < 50) ? (2000 + yy) : (1900 + yy);
  jd = dd - 32075L + 1461L * (yy + 4800L + (mm - 14) / 12 ) / 4;
  jd = jd + 367L * (mm - 2 - (mm - 14) / 12*12) / 12;
  jd = jd - 3 * ((yy + 4900L + (mm - 14) / 12) / 100) / 4;

  return(jd);
} 

/* Convert degrees to radians */
double deg_to_rad(double deg)
{
return deg*pi/180.;
}

/* Convert radians to degrees */
double rad_to_deg(double rad)
{
return rad*180./pi;
}
	     



/* 
 * return variation (in radians) given geodetic latitude (radians), longitude
 * (radians) ,height (km) and (Julian) date
 * model=1 is IGRF90, 2 is WMM85, 3 is WMM90, 4 is WMM95, 
 * 5 is IGRF95, 6 is WMM2000, 7 is IGRF2000
 * N and E lat and long are positive, S and W negative
*/

double SGMagVar( double lat, double lon, double h, long dat, int model, double* field )
{
    /* output field B_r,B_th,B_phi,B_x,B_y,B_z */
    int n,m;
    double yearfrac,sr,r,theta,c,s,psi,fn,fn_0,B_r,B_theta,B_phi,X,Y,Z;
    double sinpsi, cospsi, inv_s;

    static int been_here = 0;
	
    double sinlat = sin(lat);
    double coslat = cos(lat);

    /* convert to geocentric */ 
    sr = sqrt(a*a*coslat*coslat + b*b*sinlat*sinlat);
    /* sr is effective radius */
    theta = atan2(coslat * (h*sr + a*a), sinlat * (h*sr + b*b));
		  
    /* theta is geocentric co-latitude */

    r = h*h + 2.0*h * sr +
	(a*a*a*a - ( a*a*a*a - b*b*b*b ) * sinlat*sinlat ) / 
	(a*a - (a*a - b*b) * sinlat*sinlat );

    r = sqrt(r);

    /* r is geocentric radial distance */
    c = cos(theta);
    s = sin(theta);
      /* protect against zero divide at geographic poles */
    inv_s =  1.0 / (s + (s == 0.)*1.0e-8); 

    /*zero out arrays */
    for ( n = 0; n <= nmax; n++ ) {
      for ( m = 0; m <= n; m++ ) {
	P[n][m] = 0;
	DP[n][m] = 0;
      }
    }

    /* diagonal elements */
    P[0][0] = 1;
    P[1][1] = s;
    DP[0][0] = 0;
    DP[1][1] = c;
    P[1][0] = c ;
    DP[1][0] = -s;

    /* these values will not change for subsequent function calls */
    if( !been_here ) {
      for ( n = 2; n <= nmax; n++ ) {
	root[n] = sqrt((2.0*n-1) / (2.0*n));
      }
	
      for ( m = 0; m <= nmax; m++ ) {
	double mm = m*m;
	for ( n = max(m + 1, 2); n <= nmax; n++ ) {
	  roots[m][n][0] = sqrt((n-1)*(n-1) - mm);
	  roots[m][n][1] = 1.0 / sqrt( n*n - mm);
	}
      }
      been_here = 1;
    }
	
    for ( n=2; n <= nmax; n++ ) {
      /*  double root = sqrt((2.0*n-1) / (2.0*n)); */
      P[n][n] = P[n-1][n-1] * s * root[n];
      DP[n][n] = (DP[n-1][n-1] * s + P[n-1][n-1] * c) * root[n];
    }

    /* lower triangle */
    for ( m = 0; m <= nmax; m++ ) {
      /*  double mm = m*m;  */
      for ( n = max(m + 1, 2); n <= nmax; n++ ) {
	/* double root1 = sqrt((n-1)*(n-1) - mm); */
	/* double root2 = 1.0 / sqrt( n*n - mm);  */
	P[n][m] = (P[n-1][m] * c * (2.0*n-1) -
		   P[n-2][m] * roots[m][n][0]) * roots[m][n][1];
	DP[n][m] = ((DP[n-1][m] * c - P[n-1][m] * s) *
		    (2.0*n-1) - DP[n-2][m] * roots[m][n][0]) * roots[m][n][1];
      }
    }

    /* compute gnm, hnm at dat */
    switch(model) {
    case 1:	/* IGRF90 */
      yearfrac = (dat - yymmdd_to_julian_days(90,1,1)) / 365.25;
      for (n=1;n<=nmax;n++)
        for (m = 0;m<=nmax;m++) {
	  gnm[n][m] = gnm_igrf90[n][m] + yearfrac * gtnm_igrf90[n][m];
	  hnm[n][m] = hnm_igrf90[n][m] + yearfrac * htnm_igrf90[n][m];
	}
      break;

    case 2:	/* WMM85 */
      yearfrac = (dat - yymmdd_to_julian_days(85,1,1)) / 365.25;
      for (n=1;n<=nmax;n++)
        for (m = 0;m<=nmax;m++) {
	  gnm[n][m] = gnm_wmm85[n][m] + yearfrac * gtnm_wmm85[n][m];
	  hnm[n][m] = hnm_wmm85[n][m] + yearfrac * htnm_wmm85[n][m];
	}
      break;

    case 3:	/* WMM90 */
      yearfrac = (dat - yymmdd_to_julian_days(90,1,1)) / 365.25;
      for (n=1;n<=nmax;n++)
	for (m = 0;m<=nmax;m++) {
	  gnm[n][m] = gnm_wmm90[n][m] + yearfrac * gtnm_wmm90[n][m];
	  hnm[n][m] = hnm_wmm90[n][m] + yearfrac * htnm_wmm90[n][m];
	}
      break;

    case 4:	/* WMM95 */
      yearfrac = (dat - yymmdd_to_julian_days(95,1,1)) / 365.25;
      for (n=1;n<=nmax;n++)
        for (m = 0;m<=nmax;m++) {
           gnm[n][m] = gnm_wmm95[n][m] + yearfrac * gtnm_wmm95[n][m];
	   hnm[n][m] = hnm_wmm95[n][m] + yearfrac * htnm_wmm95[n][m];
	}
      break;
    case 5:	/* IGRF95 */
      yearfrac = (dat - yymmdd_to_julian_days(95,1,1)) / 365.25;
      for (n=1;n<=nmax;n++)
        for (m = 0;m<=nmax;m++) {
           gnm[n][m] = gnm_igrf95[n][m] + yearfrac * gtnm_igrf95[n][m];
	   hnm[n][m] = hnm_igrf95[n][m] + yearfrac * htnm_igrf95[n][m];
	}
      break;
    case 6:      /* WMM2000 */
      yearfrac = (dat - yymmdd_to_julian_days(0,1,1)) / 365.25;
       for (n=1;n<=nmax;n++)
         for (m = 0;m<=nmax;m++) {
           gnm[n][m] = gnm_wmm2000[n][m] + yearfrac * gtnm_wmm2000[n][m];
	   hnm[n][m] = hnm_wmm2000[n][m] + yearfrac * htnm_wmm2000[n][m];
	 }
       break;
    case 7:      /* IGRF2000 */
      yearfrac = (dat - yymmdd_to_julian_days(0,1,1)) / 365.25;
      for (n=1;n<=nmax;n++)
	for (m = 0;m<=nmax;m++) {
	  gnm[n][m] = gnm_igrf2000[n][m] + yearfrac * gtnm_igrf2000[n][m];
	  hnm[n][m] = hnm_igrf2000[n][m] + yearfrac * htnm_igrf2000[n][m];
	}
      break;  
    }

    /* compute sm (sin(m lon) and cm (cos(m lon)) */
    for (m = 0;m<=nmax;m++) {
      sm[m] = sin(m * lon);
      cm[m] = cos(m * lon);
    }

    /* compute B fields */
    B_r = 0.0;
    B_theta = 0.0;
    B_phi = 0.0;
    fn_0 = r_0/r;
    fn = fn_0 * fn_0;
	
    for ( n = 1; n <= nmax; n++ ) {
      double c1_n=0;
      double c2_n=0;
      double c3_n=0;
      for ( m = 0; m <= n; m++ ) {
	double tmp = (gnm[n][m] * cm[m] + hnm[n][m] * sm[m]); 
	c1_n += tmp * P[n][m];
	c2_n += tmp * DP[n][m];
	c3_n +=  m * (gnm[n][m] * sm[m] - hnm[n][m] * cm[m]) * P[n][m];
      }
      /* fn=pow(r_0/r,n+2.0);   */
      fn *= fn_0;
      B_r += (n + 1) * c1_n * fn;
      B_theta -= c2_n * fn;
      B_phi += c3_n * fn * inv_s;
    }



    /* Find geodetic field components: */
    psi = theta - (pi / 2.0 - lat);
    sinpsi = sin(psi);
    cospsi = cos(psi);
    X = -B_theta * cospsi - B_r * sinpsi;
    Y = B_phi;
    Z = B_theta * sinpsi - B_r * cospsi;

    field[0]=B_r;
    field[1]=B_theta;
    field[2]=B_phi;
    field[3]=X;
    field[4]=Y;
    field[5]=Z;   /* output fields */
    /* find variation in radians */
    /* return zero variation at magnetic pole X=Y=0. */
    /* E is positive */
    return (X != 0. || Y != 0.) ? atan2(Y, X) : (double) 0.; 
}