undef("demod")
function demod(y[*]:numeric, frq[1]:numeric)
;
; simple complex demodulation
;
local ny, pi, arg, amp, pha
begin

  ny = dimsizes(y)
  if (typeof(y).eq."double" .or. typeof(frq)) then
      pi = 4d0*atan(1d0)
  else
      pi = 4.0*atan(1.0)
  end if

  arg = ispan(0,ny-1,1)*frq*2*pi
  amp =  y*cos(arg)*2             ; real
  pha = -y*sin(arg)*2             ; imaginary (primary phase)

  amp@long_name = "amplitude"
  pha@long_name = "phase"         
  if (isatt(y,"units")) then
      amp@units = y@units
  end if

  copy_VarCoords(y, amp)
  copy_VarCoords(y, pha)

  return([/ amp, pha /])           ; return as a list
end