#
# draw the wavenumber-frequency diagram
#   like Fig.1 of Wheeler and Kiladis [1999] JAS
#   using GPhys::FFT
#
#-- create: 28 Feb. 2014 Eriko Nishimoto
#-- modified: 03 Jun. 2014 Eriko Nishimoto
#-- modified: 09 Mar. 2015 Shin-ya Murakami

require "numru/gphys"
include NumRu

year_len = 224.701 # days
nperiod = (year_len/4.0).floor.to_i
slide = (nperiod/3.0).floor.to_i

gp_s, gp_a = ["sym", "antisym"].map do |type|
#-----------------------------open data 
  dir="./"
  fname="daily_rmAC.nc"
  var="u_detrend"
  st, et = 70000, 80000
  
  slat, elat=-15, 15
  dum=GPhys::IO.open(dir+fname,var).cut(true,slat..elat,true)
  
  #-----------------------------make data
  #-- components of antisymmetric and symmetric about the equator
  nx,ny,nt=dum.shape
  mask=NArray.int(ny).indgen(-1,-1)+ny
  if type=="sym"
    #-----------------------------symmetric component
    #-- (a(x,y,t)+a(x,-y,t))/2
    tmp=(dum+dum[true,mask,true])/2.0
  elsif type=="antisym"
    #-----------------------------antisymmetric component
    #-- (a(x,y,t)-a(x,-y,t))/2
    tmp=(dum-dum[true,mask,true])/2.0
  end
  
  #-----------------------------FFT
  nx=tmp.shape[0]
  fc=[]
  power_sum = 0
  ntimes = ((et - st + 1 - slide * 2)/slide).floor.to_i
  (0..ntimes-1).each do |i|
    st0 = (st + i * slide).floor.to_i
    et0 = (st + i * slide + nperiod).floor.to_i
    #print "#{i}: st0..et0 = #{st0}..#{et0}\n"
    data=tmp.cut('time'=>st0..et0).detrend(2).cos_taper(2)
    nx,ny,nt=data.shape

    #----------------------foward FFT
    GPhys::fft_ignore_missing(true)
    fc=data.fft(false,0,2).mean(1)

    #----------------------change axes
    z0=fc.axis(0).pos
    fc.axis(0).set_pos(-z0)
    fc.coord(0).name = 'k'
    fc.coord(0).set_att('long_name', 'k-wavenumber')
    fc.coord(0).del_att('standard_name')
    fc.coord(0).del_att('axis')
    fc.coord(0).del_att('gt_calc_weight')

    z1=fc.axis(-1).pos.convert_units("day-1")/(2.0*Math::PI)
    fc.axis(-1).set_pos(z1)
    fc.coord(-1).name = 'freq'
    fc.coord(-1).set_att('long_name', 'frequency')
    fc.coord(-1).del_att('standard_name')
    fc.coord(-1).del_att('axis')
    fc.coord(-1).del_att('delta_t')
    
    #----------------------power spectral
    #-- |fc|^2
    fc.name='sd'
    fc.set_att('long_name', 'spectral density')

    power=(fc.abs**2).val * GPhys::COS_TAPER_SP_FACTOR
    power_sum += power/ntimes.to_f
  end
  fc.mean(1)
  fc.copy.replace_val(power_sum).rawspect2powerspect(0,-1)\
  .spect_zero_centering(0).spect_one_sided(-1)
end

type="sym"
outfile = NetCDF.create("sp-#{type}.nc")
GPhys::IO.write(outfile, gp_s)
outfile.close

type="antisym"
outfile = NetCDF.create("sp-#{type}.nc")
GPhys::IO.write(outfile, gp_a)
outfile.close