;+
;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
;    spurious 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 
; fail = if set to a named variable, returns 1 if an error occurs, 0 otherwise
;
;$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, fail=fail,_extra = _extra
;
  fail = 1

  tdps = -1 & fdps = -1 & dps = -1 ;init output variables
  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'
    return
  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)
  if(nspectra le 0) then begin
    message, /info, 'Not enough points for a calculation'
    return
  endif
  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
;
    if(nend le totalpoints-1) then begin
      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 = linfit(t, x, yfit = line) ;8-oct-2008, jmm
        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...', /info
        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[*]
    endif
;
  endfor
;
  fdps = freqcenter
  
  fail = 0
;
  return
end