;  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 wavelet as defined by Sato.

function p_jk_xL, H, d, xh, j, k, xL

  psihatL=p_jk_wL(H, d, xh, j, k, wL)
  n = n_elements(wL)
  dn = d * n
  xL=[findgen(n/2)-n/2,findgen(n/2)] * d + (xh-H)
  psiL=shift(float(fft( shift(psihatL,n/2)/dn, /inverse)),n/2)
  
  return, psiL

end

pro p_jk_xL_test

  H = 32.
  d = 1.
  j=3
  k=1
  xh=0.
  psirefrL = refr_jk_xL(H, d, xh, 0, 0, xL)
  psirefrjkL = refr_jk_xL(H, d, xh, j, k, xL)
  psipL = phase_jk_xL(H, d, xh, 0, 0, xL)
  psipjkL = phase_jk_xL(H, d, xh, j, k, xL)
  window, 0, xs=1000,ys=500
  !p.multi=[0,2,0]
  plot, xL/H, psirefrjkL, xr=xh + [-5,5], /xs,/ys
  oplot, xL/H, psirefrL, li=2

  plot, xL/H, psipjkL, xr=xh + [-5,5], /xs,/ys
  oplot, xL/H, psipL, li=2

end

pro p_jk_xL_test2

  H = 32.
  d = 1.
  xh=32000.
  psirefrL = refr_jk_xL(H, d, xh, 0, 0, xL)
  c = 0.5 / max(psirefrL)
  refrL = exp(-(xL-xH)/H) * (1 + c * psirefrL)
  
  psipL = p_jk_xL(H, d, xh, 0, 0, xL)
  pL = exp(-(xL-xH)/H) * (1 + c * psipL*0.5)
  n=n_elements(xL)
  pL2 = fltarr(n)
  for i=0, n-1 do $
    pL2[i] = atmint(reverse(xL)*1e5, reverse(refrL), xL[i]*1e5, xL[n-1]*1e5, 1., EXPON=expon)

  pL2 = pL2 / (H*1e5)
  
  window, 0, xs=1000,ys=500
  !p.multi=[0,2,0]
  plot, (xL-xH)/H, refrL, xr=[-15,5], /xs,/ylog

  plot, (xL-xH)/H, pL, xr=[-15,5], /xs,/ylo
  oplot, (xL-xH)/H, pL2, /li

end

pro p_jk_xL_test3

  H = 32.
  d = 1.
  xh=32000.
  psirefrL = refr_jk_xL(H, d, xh, 0, 0, xL)
  c = 0.5 / max(psirefrL)
  refrL = exp(-(xL-xH)/H) * (1 + c * psirefrL)
  
  psipL = p_jk_xL(H, d, xh, 0, 0, xL)
  pL = exp(-(xL-xH)/H) * (1 + c * psipL)
  n=n_elements(xL)
  pL2 = fltarr(n)
  for i=0, n-1 do $
    pL2[i] = atmint(reverse(xL)*1e5, reverse(refrL), xL[i]*1e5, xL[n-1]*1e5, 1., EXPON=expon)

  pL2 = pL2 / (H*1e5)
  
  window, 0, xs=1000,ys=500
  !p.multi=[0,2,0]
  plot, (xL-xH)/H, refrL/exp(-(xL-xH)/H), xr=[-19,5], /xs

  plot, (xL-xH)/H, pL/exp(-(xL-xH)/H), xr=[-19,5], /xs
  oplot, (xL-xH)/H, pL2/exp(-(xL-xH)/H), /li, ps=-4

end