# -*- coding: utf-8 -*-
require "numru/gphys"
require "numru/ggraph"
include NumRu

def usage 
  "\n USAGE: $ ruby #{$0} filename variablename \n"
end

#filename = ARGV[0] || raise(usage)
#varname = ARGV[1] || raise(usage)

varname = ARGV[0] || raise(usage)

puts "ignore"
GPhys::fft_ignore_missing(true, 0.0)

filedir = "time_000180000-000280000"
filehead = "BS1998_"
zcut = 10 # km で指定
cuttime = 180000..280000
cutk = 0.001..1
ins =2
itr = 3
puts "fileopen"

var0_file = "../#{filedir}/#{filehead}#{varname}_rank000000.nc"
var1_file = "../#{filedir}/#{filehead}#{varname}_rank000001.nc"
var2_file = "../#{filedir}/#{filehead}#{varname}_rank000002.nc"
var3_file = "../#{filedir}/#{filehead}#{varname}_rank000003.nc"
var4_file = "../#{filedir}/#{filehead}#{varname}_rank000004.nc"
var5_file = "../#{filedir}/#{filehead}#{varname}_rank000005.nc"
time = GPhys::IO.open("../#{filedir}/#{filehead}#{varname}_rank000000.nc", 't').val
puts "file open"

gphys = GPhys::IO.open([var0_file, \
                        var1_file, \
                        var2_file, \
                        var3_file, \
                        var4_file, \
                        var5_file], \
                        varname)


#zz = GPhys::IO.open(filename, 'z').val
x1 = gphys.axis("x").pos * 1e-3 
x1.units = Units["km"]
z1 = gphys.axis("z").pos * 1e-3 
z1.units = Units["km"]
gphys.axis("x").set_pos(x1)
gphys.axis("z").set_pos(z1)


puts "detrend"
pre = gphys #.detrend(3).cos_taper(3)
pre = pre.cut('z'=>zcut).mean('y')

puts "fft"
sp = pre.fft(false, 0).abs ** 2.0

puts "power"
power = sp.rawspect2powerspect(0)#.spect_one_sided(3).spect_zero_centering(0)#.cut("t"=>9000).mean("y")

DCL.gropn(ins)
DCL.uzfact(0.5)
DCL.sgpset('lfull',true)
#GGraph.set_fig( 'itr'=> 4)
GGraph.set_fig( 'itr'=> itr, 'viewport'=>[0.25,0.7,0.15,0.6] )
#GGraph.line( power[1..1199] )#, true, 'exchange'=>true )

#for i in 0...time.length do
#  if (i.to_i % 5) == 0 
    GGraph.set_axes('xtitle'=>'k', 'ytitle'=>'t')
#    GGraph.tone( power.cut('k'=>cutk,'t'=>cuttime),true )
    GGraph.tone( power[1..1199,true],true, 'title'=>'Power spectrum' )
#    GGraph.contour( power[1..1199,i], false )
    GGraph.color_bar('vlength'=>0.5,"landscape"=>true,'tickintv'=>0)

#  end
#end



outfile = NetCDF.create("line_PwSpect_#{varname}.nc")

#GPhys::IO.write(outfile, sp)
GPhys::IO.write(outfile, power)
outfile.close
DCL.grcls