;+
; NAME: 
;       LEG
; 
; PURPOSE:
;       This procedure calculates the associated Legendre functions P(l,n)
;       at a given point x, up to a maximum order P(nmax,nmax). It is based
;       on the varying order recurrence relation given by Abramovitz & Stegun,
;       Handbook of Mathematical Functions, 1970 (eq. 8.5.3). The
;       Schmidt normalization is used here , in accordance with usual practice
;       in geomagnetic studies. The procedure also gives d P(n,m) / d (theta), 
;       where theta = arccos(x), by taking the derivative of the same
;       recurrence relation. Finally, the procedure gives the value of
;	P(n,m)/sin(theta).
;
; CATEGORY:
;       Mathematical procedures.
;
; CALLING SEQUENCE:
;       LEG, x, nmax, p, q, s
;
; INPUTS:
;          X: the point at which P(n,m) is to be evaluated.
;
;       NMAX: maximum order of P(n,m).
;
; KEYWORDS:
;       None.
;
; OUTPUTS:
;       P: an array of size (nmax + 1) x (nmax + 1), containing the function values
;          (P(n,m) = 0 for n < m).
;       Q: an array of size (nmax + 1) x (nmax + 1), containing the values of 
;          d P(n,m) / d theta.
;	S: an array of size (nmax + 1) x (nmax + 1), containing the values of
;	   P(n,m)/sin(theta).
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None.
;
; RESTRICTIONS:
;       None.
;
; MODIFICATION HISTORY:
;       September 6, 1993. Written by A. Louro.
;-

pro leg, x, nmax, p, q, s

p = fltarr(nmax + 1, nmax + 1)
q = fltarr(nmax + 1, nmax + 1)
s = fltarr(nmax + 1, nmax + 1)
aux = sqrt(1. - x * x)
sr2 = sqrt(2.)

p(0,0) = sr2 				; Set initial values for
p(1,0) = sr2 * x			; recurrence procedure
q(0,0) = 0.
q(1,0) = - sr2 * aux

for n = 1, nmax do begin 		; Calculate diagonal values P(n,n)
  c1 = sqrt(1. - 1. / (2. * n))
  p(n,n) = c1 * aux * p(n - 1, n - 1)
  q(n,n) = c1 * (x * p(n - 1, n - 1) + aux * q(n -1, n - 1))
  s(n,n) = c1 * p(n - 1, n - 1)
endfor

for m = 0, nmax - 1 do begin		; Fill in rest of array
  mm = max([m, 1])
  for n = mm, nmax - 1 do begin
    c2 = 2. * n + 1.
    c3 = sqrt((n + m) * (n - m))
    c4 = sqrt((n - m + 1) * (n + m + 1))
    p(n + 1, m) = (c2 * x * p(n,m) - c3 * p(n - 1, m)) / c4
    q(n + 1, m) = (c2 * (x * q(n,m) - aux * p(n,m)) - $
		   c3 * q(n - 1, m)) / c4
    s(n + 1, m) = (c2 * x * s(n,m) - c3 * s(n - 1, m)) / c4
  endfor
endfor

for n = 0, nmax do begin		; "Correct" (n,0) terms by
  p(n, 0) = p(n, 0) / sr2		; dividing by sqrt(2)
  q(n, 0) = q(n, 0) / sr2
endfor
return
end

;+
; NAME:
;	IGRF
; PURPOSE:
;       Computation of magnetic field vectors according to the
;       series of International Geomagnetic Reference Fields (1945
;	- 1990).
; CATEGORY:
;	Magnetic field modelling.
; CALLING SEQUENCE:
;	IGRF, LAT, LON, HEIGHT, YEAR, X, Y, Z, HOR, DIP, DEC
; INPUTS:
;       LAT:    The geographic latitude (degrees) of the point for
;               which the magnetic field vector is to be computed.
;	LON:	The geographic longitude of the point.
;	HEIGHT:	The altitude (km) of the point above mean sea level.
;	YEAR:	The year of observations (either 2 or 4 digits).
; OPTIONAL INPUTS:
;	None.
; KEYWORD PARAMETERS:
;	None.
; OUTPUTS:
;       X:      The geographic northward component of the magnetic
;		field vector (gauss; 1 gauss = 1E-4 T).
;	Y:	The geographic eastward component of the vector.
;	Z:	The geographic downward component of the vector.
;	HOR:	The geographic horizontal component of the vector.
;       DIP:    The angle of the field vector below the horizontal
;		(degrees).
;       DEC:    The direction of the horizontal component (degrees
;		east of north).
; OPTIONAL OUTPUTS:
;	None.
; COMMON BLOCKS:
;	None.
; SIDE EFFECTS:
;       Files containing the spherical harmonic coefficients of
;	the field expansion are read.
; RESTRICTIONS:
;       Only scalar arguments are permitted.  The YEAR must be in
;	the range 1945 - 1995 (for now).
; PROCEDURE:
;	See Akasofu & Chapman, Solar Terrestrial Physics, pp. 68ff.
; EXAMPLE:
;	
; SEE ALSO:
;	
; MODIFICATION HISTORY:
; 	Written by:	Dr. Alfredo Louro, ISR, 1991-1995.
;-

pro igrf,lat,lon,height,year,x,y,z,hor,dip,dec

IF N_PARAMS() LT 7 THEN BEGIN
    doc_library,'igrf'
    RETURN
    ENDIF

; Set up OS-specific parameters.
@isitdos

nlat=N_ELEMENTS(lat)
nlon=N_ELEMENTS(lon)
nht=N_ELEMENTS(height)
nyr=N_ELEMENTS(year)
IF nlat>nlon>nht>nyr GT 1 THEN BEGIN
    MESSAGE,'Only scalar arguments permitted',/INFORMATIONAL
    RETURN
    ENDIF

lat=lat(0)
lon=lon(0)
height=height(0)
year=year(0)

r0 = 6371.2
ra = 6378.16
rb= 6356.775

; Ensure the year is in the correct format (last two digits only).
yr=year MOD 100

; convert lat, lan from degrees to radians.
conv=!pi/180.
;lat_s=90.-lat
lat_s=lat*conv
lon_s=lon*conv

; convert (oblate-spheroidal) geographic coordinates to spherical.

ra2=ra*ra
ra4=ra2*ra2
rb2=rb*rb
rb4=rb2*rb2
sin2lat=sin(lat_s)*sin(lat_s)
aux2=ra2-(ra2-rb2)*sin2lat
aux4=ra4-(ra4-rb4)*sin2lat
r=sqrt(height*height+2.*height*sqrt(aux2)+aux4/aux2)
f=(height*sqrt(aux2)+ra2)^2/(height*sqrt(aux2)+rb2)^2
tlat=acos(sin(lat_s)/sqrt(f+sin2lat*(1.-f)))

; read spherical harmonic coefficients

g=fltarr(11,11)
h=fltarr(11,11)

yearlo=floor(yr/5.)*5
yearhi=yearlo+5
if yr ge 90 then begin
  gdot=fltarr(11,11)
  hdot=fltarr(11,11)
  openr,ilun1,igrfdir+dd+'igrf90.dat',/GET_LUN
  openr,ilun2,igrfdir+dd+'igrf90s.dat',/GET_LUN
  WHILE NOT EOF(ilun1) DO BEGIN
    READF,ilun1,n,m,g_aux,h_aux
    g(n,m)=g_aux
    h(n,m)=h_aux
  ENDWHILE
  WHILE NOT EOF(ilun2) DO BEGIN
    readf,ilun2,n,m,gdot_aux,hdot_aux
    gdot(n,m)=gdot_aux
    hdot(n,m)=hdot_aux
  endwhile
  FREE_LUN,ilun1
  FREE_LUN,ilun2
  g=g+(yr-90)*gdot
  h=h+(yr-90)*hdot
endif else begin
  glo=fltarr(11,11)
  ghi=fltarr(11,11)
  hlo=fltarr(11,11)
  hhi=fltarr(11,11)
  filelo=igrfdir+dd+'dgrf'+strcompress(yearlo,/remove_all)+'.dat'
  if yearhi eq 90 then filehi=igrfdir+dd+'igrf90.dat' else $
    filehi=igrfdir+dd+'dgrf'+strcompress(yearhi,/remove_all)+'.dat'
  openr,ilun1,filelo,/GET_LUN
  openr,ilun2,filehi,/GET_LUN
  WHILE NOT EOF(ilun1) DO BEGIN
    READF,ilun1,n,m,glo_aux,hlo_aux
    glo(n,m)=glo_aux
    hlo(n,m)=hlo_aux
  ENDWHILE
  WHILE NOT EOF(ilun2) DO BEGIN
    READF,ilun2,n,m,ghi_aux,hhi_aux
    ghi(n,m)=ghi_aux
    hhi(n,m)=hhi_aux
  ENDWHILE
  FREE_LUN,ilun1
  FREE_LUN,ilun2

  g=(yr*(ghi-glo)+(yearhi*glo-yearlo*ghi))/5.
  h=(yr*(hhi-hlo)+(yearhi*hlo-yearlo*hhi))/5.
endelse
; create phi-dependent factor of Tn and T'n

phix=fltarr(11,11)
phiy=fltarr(11,11)
phiz=fltarr(11,11)
for n=1,10 do begin
  fact=(r0/r)^(n+2)
  for m=0,n do begin
    phix(n,m)=fact*g(n,m)*cos(m*lon_s)+h(n,m)*sin(m*lon_s)
    phiy(n,m)=fact*m*(g(n,m)*sin(m*lon_s)-h(n,m)*cos(m*lon_s))
    phiz(n,m)=-phix(n,m)*(n+1)
  endfor
endfor

; calculate the associated Legendre functions P(n,m)

leg,cos(tlat),10,p,der_p,p_sin

; calculate Z and H components of geomagnetic field

x=0.
y=0.
z=0.
for n=1,10 do begin
  for m=0,n do begin
    z=z+phiz(n,m)*p(n,m)
    x=x+phix(n,m)*der_p(n,m)
    y=y+phiy(n,m)*p_sin(n,m)
  endfor
endfor

delta=tlat+lat_s-!pi/2.
xc=x
zc=z
x=xc*cos(delta)+zc*sin(delta)
z=-xc*sin(delta)+zc*cos(delta)

z=z(0)*1.e-5
x=x(0)*1.e-5
y=y(0)*1.e-5

; calculate dip angle, and declination
; and convert back from radians to degrees 

hor=sqrt(x^2+y^2)

dip=atan(z,hor)/conv
dec=atan(y,x)/conv
finish:
return
end