;+
;	Name:       THM_MAKE_AE
;
;	Purpose:    This routine calculates the "pseudo" AE, AL, and AU geomagnetic
;                   indices by using ground magnetometer data from THEMIS GBOs. The names
;                   of all stations used for calculation are printed on the screen.
;                   In future, it will be possible to include ground data from other
;                   magnetometer networks. Note that currently the calculation of
;                   the "pseudo" indices does not subtract quiet day variation but
;                   simply the median.
;
;       Syntax:     THM_MAKE_AE [, RES = float] [, SITES = string ] [, /NO_LOAD ]
;
;	Parameters: None.
;
;	Keywords:   res     = sampling interval (by default 60 sec)
;                   sites   = observatory name --> not implemented yet; the
;                             default is to load all available stations
;                   no_load = if set, use existing gmag (THEMIS) tplot variables which have
;                             already been loaded into the active TDAS environment
;                             if not set, load gmag data (either remotely or from computer disk)
;
;	Example:    see crib sheet "thm_crib_AE.pro"
;
;       Notes:      Written by Andreas Keiling, 15 May 2008
;
;       Modifications:
;                   Edited header, put subroutine before main body so that it would compile before
;                   being called by the main body, W.M.Feuerstein, 6/2/2008.
;
; $LastChangedBy: jimm $
; $LastChangedDate: 2009-02-06 12:58:33 -0800 (Fri, 06 Feb 2009) $
; $LastChangedRevision: 4908 $
; $URL $
;-


pro tdespike_AEALAU, lower, upper, quant=quant
; this routine removes artificial spikes
;
; variables:  lower = lower cutoff of spikes to be removed
;             upper = upper cutoff of spikes to be removed
;

get_data, quant, data=A_index
last=n_elements(A_index.y)-1
NaN=!values.f_nan

indices = where(A_index.y gt upper OR A_index.y lt lower , count)
if count ne 0 then begin
   for k=0,count-1 do begin
       i=indices[k]
       if (i eq 0 OR i eq 1 ) then begin
          A_index.y[0]=NaN
          A_index.y[1]=NaN
       endif else begin
          if (i eq last-1 OR i eq last) then begin
               A_index.y[last-1]=NaN
               A_index.y[last] = NaN
          endif else begin
               A_index.y[i-2]=NaN
               A_index.y[i-1]=NaN
               A_index.y[i] = NaN
               A_index.y[i+1]=NaN
               A_index.y[i+2]=NaN
          endelse
       endelse
   endfor
endif

store_data, quant+'_despike', data=A_index

end

;#######################################################

pro thm_make_AE, res=res, sites=sites, no_load=no_load


; set default time resolution
;----------------------------
if not keyword_set(res) then res = 60d   ; 60 sec resolution
res=double(res)


; load gmag data
;----------------
if not keyword_set(no_load) then thm_load_gmag,/subtract_median
split_vec,'thg_mag_????'
stations='thg_mag_????_x'


; calculate indices
;------------------
superpo_histo,stations,dif='thmAE',min='thmAL',max='thmAU',res=res


; thresholds can be changed if desired
;-------------------------------------
cut=1500
tdespike_AEALAU,-cut,cut, quant='thmAE'
tdespike_AEALAU,-cut,cut, quant='thmAL'
tdespike_AEALAU,-cut,cut, quant='thmAU'

clean_spikes, 'thmAE_despike', new_name = 'thmAE', thresh = 5
clean_spikes, 'thmAL_despike', new_name = 'thmAL', thresh = 5
clean_spikes, 'thmAU_despike', new_name = 'thmAU', thresh = 5

options,'thmAE',ytitle='THEMIS!CAE Index'
options,'thmAL',ytitle='THEMIS!CAL Index'
options,'thmAU',ytitle='THEMIS!CAU Index'


; delete obsolete data
;---------------------
del_data,'thmAE_despike'
del_data,'thmAU_despike'
del_data,'thmAL_despike'


; print station names used for calculation
;-----------------------------------------
print,'The following stations contributed to the index calculation:'
tplot_names, 'thg_mag_????_x'
print,'-----------------------------'

end