;+
; FUNCTION:	auroral_zone
;
; PURPOSE:	IDL function to calculate auroral zone position as a 
;		function of magnetic local time (mlt: 0.0-24.0 hours), 
;		activity index Q (q: 0-6).
;		Returns corrected geomagnetic colatitude in radians.
;
; OPTIONS:	To return poleward edge, set poleward.
;		To return latitude, set latitude.
;		To return value for southern oval, set south.
;		To return value in degrees, set degrees.
;
; See Holzworth & Meng, GRL 2, p. 377, 1975.
;
; Originally written by J. Clemmons, June 1993.
; Corrected by J.Rauchleiba under the direction of Mike Temerin Apr 1997.
;-

function auroral_zone, mlt, q, $
	poleward=pole, $
	latitude=lat, $
	south=south, $
	degrees=deg

; Arrays are indexed (An, Q) where An is the nth
; best-fit constant and Q is the activity index
; n =	  0    1    2     3    4    5    6
apole=[[15.22,2.41,3.34,-0.85,1.01,0.32,0.90], $
       [15.85,2.70,3.32,-0.67,1.15,0.49,1.00], $
       [16.09,2.51,3.27,-0.56,1.30,0.42,0.94], $
       [16.16,1.92,3.14,-0.46,1.43,0.32,0.96], $
       [16.29,1.41,3.06,-0.09,1.35,0.40,1.03], $
       [16.44,0.81,2.99,0.14,1.25,0.48,1.05], $
       [16.71,0.37,2.90,0.63,1.59,0.60,1.00]]
aequa=[[17.36,3.03,3.46,0.42,2.11,-0.25,1.13], $
       [18.66,3.90,3.37,0.16,2.55,-0.13,0.96], $
       [19.73,4.69,3.34,-0.57,-1.41,-0.07,0.75], $
       [20.63,4.95,3.31,-0.66,-1.28,0.30,-0.58], $
       [21.56,4.93,3.31,-0.44,-0.81,-0.07,-0.75], $
       [22.32,4.96,3.29,-0.39,-0.72,-0.16,-0.52], $
       [23.18,4.85,3.34,-0.38,-0.62,-0.53,-0.16]]

; Columns 0, 1, 3, 5 in above arrays are in degrees.
; Columns 2, 4, 6 (args to cos function) are already in radians.

Andeg = [0,1,3,5]	; The An that are in degrees
apole(Andeg,*) = !dtor * apole(Andeg,*)
aequa(Andeg,*) = !dtor * aequa(Andeg,*)

; Phi relative to Sun pointer, in terms of MLT
; Since Holzworth mentions a different definition of phi(MLT)

p=(mlt-12.)/12.0*!pi

; Force i=q to within interval [0,6]

i=(fix(q)>0)<6

; Calculate corrected geomagnetic colattitude

if keyword_set(pole) then begin
  t=apole(0,i) + $
    apole(1,i)*cos(p+apole(2,i)) + $
    apole(3,i)*cos(2*(p+apole(4,i))) + $
    apole(5,i)*cos(3*(p+apole(6,i)))
endif else begin
  t=aequa(0,i) + $
    aequa(1,i)*cos(p+aequa(2,i)) + $
    aequa(3,i)*cos(2*(p+aequa(4,i))) + $
    aequa(5,i)*cos(3*(p+aequa(6,i)))
endelse

; Convert to magnetic lattitude

if keyword_set(lat) then t=!pi/2.0-t

; Southern auroral zone is merely a reflection

if keyword_set(south) then t=-t

; Conver to degrees if desired

if keyword_set(deg) then t=t*!radeg

return,t
end
