;+
; FUNCTION:   auroral_zone_ssc
;
; PURPOSE:    Given an array of MLT values, returns latitudes of the
;             auroral zone boundaries.  The definition of the 
;             auroral zone is that used by the Satellite Situation
;             Center.  It is basically the area between two
;             sinusoidally-perturbed concentric circles centered on
;             the magnetic pole. 
;             Output is equatorward latitude in Geomagnetic
;             Coordinates, and in radians unles DEGREES is set.
; ARGUMENTS:
;
;   MLTINPUT  The input array of MLT values.  Elements should look
;             like [0, .1, .2, ...,  23.9, 24.0]
;
; KEYWORDS:
;
;   POLEWARD  Set this to a named variable to receive the poleward
;             latitude values.
;   SOUTH     Set this keyword to one if a reflection through the
;             magnetic equator is desired.  This gives
;             S. Hem. coordinates.
;   DEGREES   Set this keyword to get the coordinates in degrees,
;             otherwise output is in radians.
;
; NOTES:      The definition of the auroral zones used here is that
;             accepted by the SSC.  See also: auroral_zone.pro.
;
; CREATED:    By Joseph Rauchleiba
;             98/1/6
;-

function auroral_zone_ssc, mltinput, $
                poleward=apole, $
                south=south, $
                degrees=deg

mlt = 2.*!pi * mltinput/24.
npts = dimen1(mlt)

C = fltarr(npts)
S = fltarr(npts)
S(*) = 1.0
;dayside = where(mlt GE 1.57080 AND mlt LE 4.71239)
;S(dayside) = -1.0 ; Only needed if using atan(y/x) instead of atan(y,x)

C(*) = 0.392699
switch = where(mlt GT 0.398517)
C(switch) = 5.89049

aequa = 90. - 25.*cos($
                               atan($
                                     (sin(mlt + C)/5.),  $
                                     (S*sqrt(1. - (sin(mlt + C)/5.)^2)) $
                                   )$
                             ) + 5.*cos(mlt + C)

apole = 90. - 15.*cos($
                               atan($
                                     (sin(mlt)/3.),  $
                                     (S*sqrt(1. - (sin(mlt)/3.)^2)) $
                                   )$
                             ) +5.*cos(mlt)

if keyword_set(south) then begin
    aequa = -aequa
    apole = -apole
endif
if NOT keyword_set(deg) then begin
    aequa = aequa * !dtor
    apole = apole * !dtor
endif

return, aequa
end




