;  Given a number of points and a spacing, return an
;  array of the wavelet fourier series and the
;  array of frequencies (2 pi/wavelength) for the
;  j, k=0 wavelet as defined by Sato.

function sato_wL, n, d, diln, posn, wL

  dn = d * n
  wL=[findgen(n/2)-n/2,findgen(n/2)] * 2. * !pi/dn
  wabsL = diln * abs(wL)
  wtimesL = wabsL * 3./!pi
  
  ; initialize the return array with zeros.
  psihatL = complexarr(n)
  
  ; fill in case # 1
  gd = where (wtimesL gt 2. and wtimesL lt 4., ngd)
  if ngd gt 0 then begin
     psihatL[gd] = 1.-sato_g(wabsL[gd])
  end

  ; fill in case # 2
  gd = where (wtimesL eq 4., ngd)
  if ngd gt 0 then begin
     psihatL[gd] = 1.
  end

  ; fill in case # 3
  gd = where (wtimesL gt 4. and wtimesL lt 8., ngd)
  if ngd gt 0 then begin
     psihatL[gd] = sato_g(wabsL[gd]/2.)
  end
  
  ; multiply by exp(-iw/2) * exp(-i posn w)
  gd = where (wtimesL gt 2. and wtimesL lt 8., ngd)
  ishift=(0.5 + posn) * complex(0,1)
  if ngd gt 0 then begin
     psihatL[gd] = sqrt(diln*psihatL[gd]) * exp(-ishift*wL[gd])
  end
  
  return, psihatL

end

pro sato_wL_test

  n = 256
  d = 8./n
  psihatL = sato_wL(n, d, 1, 0, wL)
  window, 0, xs=1000,ys=500
  !p.multi=[0,2,0]
  plot, wL, abs(psihatL), xr=[-10,10], yr=[0,1.2], /ys, ps=-4

end