;+

; FUNCTION:	dipole_offset

;

; PURPOSE:	Calculate the location of the dipole using the IGRF

;		coefficients and their secular variation.  IGRF coefficients

;		are the coefficients used in the spherical harmonic

;		expansion of the earth's magnetic field and there is a

;		standard formula for converting the first few of these

;		coefficients to dipole offset.

;

;		This version is accurate for years 1995-2000

;

;		Returns a 3-element vector in kilometers.

;

; INPUT:	Optional argument is 4-digit year

;

; Algorithm by:	Mike Temerin

; Written by:	J.Rauchleiba  4/8/97

;-

function dipole_offset, refyear



if NOT keyword_set(refyear) then refyear=strmid(systime(0), 20, 4)

refyear = fix(refyear)

revision = 1995

year = refyear - revision



g10 = -29682.0 + 17.6*year

g11 = -1789.0 + 13.0*year

g20 = -2197.0 - 31.2*year

g21 = 3074.0 + 3.7*year

g22 = 1685.0 - .8*year



h11 = 5318.0 - 18.3*year

h21 = -2356.0 - 15.0*year

h22 = -425.0 - 8.8*year

h2 = g10*g10 + g11*g11 + h11*h11



sq3 = sqrt(3.0)



L0 = 2.*g10*g20 + (g11*g21 + h11*h21)*sq3

L2 = -h11*g20 + (g10*h21 - h11*g22 + g11*h22)*sq3

L1 = -g11*g20 + (g10*g21 + g11*g22 + h11*h22)*sq3



e = (L0*g10 + L1*g11 + L2*h11)/(4.*h2)



MEAN_RAD = 6372.1 ;km



magdisp = dblarr(3)



magdisp(0) = MEAN_RAD*(L1 - g11*e)/(3.*h2)

magdisp(1) = MEAN_RAD*(L2 - h11*e)/(3.*h2)

magdisp(2) = MEAN_RAD*(L0 - g10*e)/(3.*h2)



return, magdisp



end