; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
;                         deconvolib
;
;     Library of routines for deconvolution of UBF waveforms
;
; Allows despinning and calibration in nT of waveforms to low
; frequency (typically  0.1-10 Hz)
; 
; Assumes existence of function which returns the 
; complex gain of the antennas in V/nT for a given frequency, 
; antenna,calibration file for a given satellite, and sample 
; frequency.  If given a negative frequency, it must return the
; complex conjugate of the gain. The form of the function must be 
; gainant(f, ix, isat, fe),
; 
; Based on deconvolib by P. Robert CNRS/CETP
;
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

pro thm_scm_apotrap, s, ns
; thm_scm_apotrap: apply a trapeziodal window to a waveform. waveform
;                  may be real or complex
; Based on P. Robert's apotrap in deconvolib.f90

na = ns/16

ap = dindgen(na)/double(na-1)

ap_ind = lindgen(na)
s[ap_ind] *= ap

ap_ind = ns -1 - ap_ind
s[ap_ind] *= ap

end

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

pro thm_scm_corgain, s, ns, df, fe, gainant, ix, pathfil,frq, init=init
; thm_scm_corgain: apply correction for antenna gain to complex 
;                  spectrum s
; Based on P. Robert's corgain in deconvolib.f90

  if keyword_set(init) then s = dcomplexarr(ns)+complex(1.0,0)
  frq = findgen(ns)*df
  frq[ns/2+1:*] -= float(ns)*df
  ;; IDL version of gainant handles negative frequencies.
  if keyword_set(gainant) then $
     c = call_function(gainant, frq, ix, pathfil, fe) $
  else c = complex(1.0, 0)
  smallc = where(abs(c) lt 1.e-6, nsmallc)
  if nsmallc gt 0 then c[smallc] = complex(1.e-6, 0.)
  s /= c
  ;; FFT of real function must be real at f=0
  s[0]=abs(s[0])*((real_part(s[0]) lt 0) ? -1 : 1) 
end

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

pro thm_scm_coretar, s, ns, fe, tau
; thm_scm_coretar: apply correction for sample delay of tau seconds 
;                  to complex spectrum s
; Based on P. Robert's coretar in deconvolib.f90

  t = double(ns)/fe
  ns2=ns/2
  
  n = lindgen(ns)
  n[ns2+1:*] -= ns
  teta = 2.0d*!dpi*n*double(tau)/t
  srot = exp(dcomplex(0, 1)*teta)
  s *= srot

end

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
     
pro thm_scm_filtspe, s, ns, df, f1, f2
; thm_scm_filtspe: filter a complex spectrum s below f1 and above f2
; Based on P. Robert's filtspe in deconvolib.f90

f = findgen(ns)*df
f[ns/2+1:*] -= float(ns)*df
f = abs(f)

; light smoothing below f1 [sic]

;arg = (f/(2*df))**2
;arg = arg < 37.
;af2 = 1.-exp(-arg)

; changed January 21, 2002: rectangular filter
af2 = 1.

filt = where(f lt f1 or f gt f2, nfilt)
if nfilt gt 0 then s[filt] = complex(0., 0.)

;filt = where(f ge f1 && f le f2, nfilt)
;if nfilt gt 0 then s[filt] *= af2

end

; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

pro thm_scm_deconvo_vec, xio, nbp, nk, f1, f2, fe, gainant, ix, pathfil, tau, $
                         edge_truncate=edget, edge_wrap=edgew, edge_zero=edgez, $
                         plotk=plotk

;+
;Name:
; thm_scm_deconvo_vec
;Purpose:
; continuous calibration accomplished by 
; convolveing a signal with a kernel.
; The kernel is derived by taking the fft of the 
; inverse of the frequency reponse.  The  kernel is also constructed
; to account for any sample delay.
;Inputs:
;    xio: input/output waveform
;    nbp: number of points in xio
;     nk: number of points in kernel: power of 2.
;     f1: low frequency
;     f2: high frequency
;     fe: sample frequency
;gainant: gain function, such as thm_scm_gainant
;     ix: antenna number (1,2,3)
;pathfil: calibration file name
;    tau: time correction
;-

  cx = ['','X','Y','Z']

  df = fe/float(nk)

  ;; pre-centering
  xavg = total(xio)/nbp
  xio -= xavg
  
  ;; generate kernel that compensates for gain and sample delay.
  ;; first get desired complex frequency reponse of filter

  s = dcomplexarr(nk)+dcomplex(1.0, 0)

  ;; corgain: inverse of antenna gain (s,nk,df,fe,gainant,ix,pathfil)
  thm_scm_corgain, s, nk, df, fe, gainant, ix, pathfil, frq

  ;; coretar: correct for sample delay (s,nk,fe,tau)
  thm_scm_coretar, s, nk, fe, tau
  
  ;; filtspe: bandpass filter between f1 and f2 (s,nbp,df,f1,f2)
  thm_scm_filtspe, s, nk, df, f1, f2

  ;; inverse fft to obtain kernel
  ks = fft(s, 1) 
  kernel = real_part(ks)

  ;; if we did everything right, the imaginary part truly is negligible
  pr = total(kernel*kernel)
  pi = total(imaginary(ks)*imaginary(ks))
  if pi/pr gt 1.e-5 then begin
     print, '*** thm_scm_deconvo: Imag/Real for impulse reponse for ', cx[ix], $
            ' =', pi/pr
  endif

  ;; zero phase of the kernel is at index 0. shift that to center, nk/2,
  ;; for use with IDL convol routine.
  kernel = shift(kernel, nk/2)
  
  ;; As this is a continuous calibration, the window must be applied to the 
  ;; kernel, rather than to the waveform.

  kernel_orig = kernel
  kernel *= hanning(nk, /double) 
  ;; note: application of window introduces a slight offset, which must be 
  ;; removed from the signal afterwards.  
  ;; Correcting for the offset in the kernel itself 
  ;; would nullify the benefit the of window.

  if keyword_set(plotk) then begin
     ;; visualize effects of window on frequency response:
     ;; (we expect to see an effect on the phase, when compared
     ;;  to ideal correction spectrum, due to centering of kernel)

     !p.multi=[0,0,4]
     ;; plot kernel - compare windowed to unwindowed
     plot, kernel_orig, /xstyle, color=1
     oplot, kernel, color=2
     re = max(kernel)
     oplot, [nk/2, nk/2], [-re, +re], color=3

     ;; plot amplitude response, compare shifted windowed and unwindowed to ideal.
     nk21 = nk/2+1
     plot, shift(frq,-nk21), shift(abs(s),-nk21), /ylog, $
           yrange=[0.01,max(abs(s))], /xstyle
     oplot, shift(frq,-nk21), shift(abs(fft(kernel_orig)),-nk21), color=1
     oplot, shift(frq,-nk21), shift(abs(fft(kernel)),-nk21), color=2

     ;; plot amplitude response, compare shifted windowed and unwindowed to ideal.
     plot, shift(frq,-nk21), shift(atan(s,/ph)*!radeg,-nk21), /xstyle
     oplot, shift(frq,-nk21),shift(atan(fft(kernel_orig),/ph)*!radeg,-nk21),color=1
     oplot, shift(frq,-nk21), shift(atan(fft(kernel),/ph)*!radeg,-nk21), color=2

     ;; plot group delay
     plot, shift(frq,-nk21), -(shift(atan(s, /ph),-nk21-1) - $
                               shift(atan(s,/ph),-nk21))/(df*2*!dpi), /xstyle
  endif
  ;; perform the convolution.
  ;;  notes on edge behavior:
  ;;  default: zero output when kernel overlaps edge
 ;;;  /edge_zero: usually good -- can emphasize low frequency trends, i.e.
 ;;;                             artifiacts of despin
 ;;;  /edge_wrap; similar to edge zero for cal signal.
 ;;;  /edge_truncate: usually bad

  kernel /= nk

  xio = convol(xio, kernel, /center, $
               edge_t=edget, edge_w=edgew, edge_z=edgez)

  ;; post-centering of waveform (necessary because of window)
  xavg = total(xio)/nbp
  xio -= xavg
  
end