;+
;NAME:
; dpwrspc
;PURPOSE:
;    Called with times time and data quantity, dpwrspc returns a dps
;    spectrum at frequencies fdps. A Hanning window is applied to
;    the input data, and its power is divided out of the returned
;    spectrum. A straight line is subtracted from the data to
;    reduce spurius power due to sawtooth behavior of a background.
;    UNITS ARE (UNITS)^2 WHERE UNITS ARE THE UNITS OF time. fdps is in Hz.
;    THUS THE OUTPUT REPRESENTS THE MEAN SQUARED AMPLITUDE OF THE SIGNAL
;       AT EACH SPECIFIC FREQUENCY. THE TOTAL (SUM) POWER UNDER THE CURVE IS
;       EQUAL TO THE MEAN (OVER TIME) POWER OF THE OSCILLATION IN TIME DOMAIN.
;    NOTE: IF KEYWORD notperhz IS SET, THEN POWER IS IN UNITS OF NT^2 ELSE
;                              IT IS IN UNITS OF NT^2/HZ
;CALLING SEQUENCE:
; dpwrspc, time, quantity, tdps, fdps, dps, nboxpoints = nboxpoints, $
;          nshiftpoints = nshiftpoints, bin = bin, tbegin = tbegin,$
;          tend = tend, noline = noline, nohanning = nohanning, $
;          notperhz = notperhz
;INPUT:
; time = the time array
; quantity = the function for which you want to obtain a power
;            spectrum
;OUTPUT:
; tps = the time array for the dynamic power spectrum
; fdps = the frequency array (units =1/time units)
; dps = the power spectrum, (units of quantity)^2/frequency_units
;KEYWORDS:
; nboxpoints = the number of points to use for the hanning window, the
;              default is 256
; nshiftpoints = the default is 128
; bin = a binsize for binning of the data, the default is 3
; tbegin = a start time, the default is time[0] 
; tend = an end time, the default is time[n_elements(time)-1]
; noline = if set, no straight line is subtracted
; nohanning = if set, then no hanning window is applied to the input
; notperhz = if set, the output units are simply the square of the
;            input units 
;
;$LastChangedBy$
;$LastChangedDate$
;$LastChangedRevision$
;$URL$
;
;-
pro dpwrspc, time, quantity, tdps, fdps, dps, nboxpoints = nboxpoints, $
             nshiftpoints = nshiftpoints, bin = bin, tbegin = tbegin, $
             tend = tend, noline = noline, nohanning = nohanning, $
             notperhz = notperhz, _extra = _extra
;
if keyword_set(tend) then begin
    t2 = tend
endif else begin
    t2 = time(n_elements(time)-1)
endelse
;
if keyword_set(tbegin) then begin
    t1 = tbegin
endif else begin
    t1 = time(0)
endelse
;
igood=where((time ge t1) and (time le t2),jgood)
if (jgood gt 0) then begin
  time2process=time(igood)
  quantity2process=quantity(igood)
endif else begin
  print,'tbegin or tend incompatible with time array'
endelse
;
if keyword_set(nboxpoints) then begin
    nboxpnts = nboxpoints
endif else begin
    nboxpnts = 256
endelse
;
if keyword_set(nshiftpoints) then begin
    nshiftpnts = nshiftpoints
endif else begin
    nshiftpnts = 128
endelse
;
if keyword_set(bin) then begin
    binsize = bin
endif else begin
    binsize = 3
endelse
;
if not(keyword_set(nohanning)) then  window = hanning(nboxpnts)
;
nboxpnts=long(nboxpnts)
nshiftpnts=long(nshiftpnts)
totalpoints=n_elements(time2process)
nspectra=long((totalpoints-nboxpnts/2l)/nshiftpnts)
tdps=dblarr(nspectra)
nfreqs=long((long(nboxpnts/2l))/binsize)
dps=fltarr(nspectra,nfreqs)
;
for nthspectrum=0,nspectra-1 do begin
;
nbegin=long(nthspectrum*nshiftpnts)
nend=nbegin+nboxpnts-1l
;ncenter=nbegin+(nboxpnts)/2l-1l
;
t = double(time2process(nbegin:nend))
t0=t(0)
t=t-t0
x = double(quantity2process(nbegin:nend))
;tdps(nthspectrum)=double(time2process(ncenter))
tdps(nthspectrum)=double(time2process(nbegin))
;
if not(keyword_set(noline)) then begin
    c=poly_fit(t,x,1,line,band,sigma)
    x=x-line
endif
;
if not(keyword_set(nohanning)) then  x = x * window
;
bign=nboxpnts
if (bign-(bign/2l)*2l ne 0) then begin
    message,'needs an even number of data points, dropping last point...',/continue
    t = t(0:bign-2)
    x = x(0:bign-2)
    bign = bign - 1
endif
;
; following Numerical recipes in Fortran, p. 421, sort of...
;
dbign = double(bign)
k = [0,dindgen(bign/2)+1]
tres=median(t(1:n_elements(t)-1)-t(0:n_elements(t)-2))
fk = k/(bign*tres)
;
xs2 = abs(fft(x,1))^2
;
pwr = dblarr(bign/2+1)
pwr(0) = xs2(0)/dbign^2
pwr(1:bign/2-1) = (xs2(1:bign/2-1) + xs2(bign - dindgen(bign/2)))/dbign^2
pwr(bign/2) = xs2(bign/2)/dbign^2
;
if not keyword_set(nohanning) then begin
    wss = dbign*total(window^2)
    pwr = dbign^2*pwr/wss
endif
;
; Now reduce variance by summing bin neighbors in frequency domain...
;
dfreq=binsize*(fk(1)-fk(0))
;
npwr = n_elements(pwr)-1
nfinal=long(npwr/binsize)
iarray=lindgen(nfinal)
power = dblarr(nfinal)
;
; Note: zeroth point includes zero freq. power.
freqcenter =(fk(iarray*binsize+1)+fk(iarray*binsize+binsize))/2.
;
for i=0l,binsize-1 do begin
    power = power+pwr(iarray*binsize+i+1)
endfor
;
if not(keyword_set(notperhz)) then begin
    power = power(*) / dfreq
    power0=pwr(0) / (fk(1)-fk(0)) ; must also include zero freq power
endif
;
dps(nthspectrum,*)=power(*)
;dps(nthspectrum,*)=power(*)/total([power0,power])
;
endfor
;
fdps=freqcenter
;
return
end