# -*- 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 = 'Kinenergy'

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
comin = 1.1e-14
comax = 300

puts "fileopen"

velx0_file = "../#{filedir}/#{filehead}VelX_rank000000.nc"
velx1_file = "../#{filedir}/#{filehead}VelX_rank000001.nc"
velx2_file = "../#{filedir}/#{filehead}VelX_rank000002.nc"
velx3_file = "../#{filedir}/#{filehead}VelX_rank000003.nc"
velx4_file = "../#{filedir}/#{filehead}VelX_rank000004.nc"
velx5_file = "../#{filedir}/#{filehead}VelX_rank000005.nc"

velz0_file = "../#{filedir}/#{filehead}VelZ_rank000000.nc"
velz1_file = "../#{filedir}/#{filehead}VelZ_rank000001.nc"
velz2_file = "../#{filedir}/#{filehead}VelZ_rank000002.nc"
velz3_file = "../#{filedir}/#{filehead}VelZ_rank000003.nc"
velz4_file = "../#{filedir}/#{filehead}VelZ_rank000004.nc"
velz5_file = "../#{filedir}/#{filehead}VelZ_rank000005.nc"

time = GPhys::IO.open("../#{filedir}/#{filehead}VelZ_rank000000.nc", 't').val
densBZ = GPhys::IO.open("./restart_long-DensBZ.nc", 'DensBZ')
puts "file open"

velx   = GPhys::IO.open([velx0_file, \
                         velx1_file, \
                         velx2_file, \
                         velx3_file, \
                         velx4_file, \
                         velx5_file ], \
                        'VelX')
velz   = GPhys::IO.open([velz0_file, \
                         velz1_file, \
                         velz2_file, \
                         velz3_file, \
                         velz4_file, \
                         velz5_file ], \
                        'VelZ')

#zz = GPhys::IO.open(filename, 'z').val

kinener = densBZ * (velx**2.0 + velz**2.0) / 2.0

x1 = kinener.axis("x").pos * 1e-3 
x1.units = Units["km"]
z1 = kinener.axis("z").pos * 1e-3 
z1.units = Units["km"]
kinener.axis("x").set_pos(x1)
kinener.axis("z").set_pos(z1)


puts "detrend"
pre = kinener #.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")
#power = sp#power.mean('t')
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'=>'E')
#    GGraph.tone( power.cut('k'=>cutk,'t'=>cuttime),true )
    GGraph.line( power[1..1199, true].mean('t'), true, 'title'=>'Power spectrum', 'max'=>comax, 'min'=>comin )#,true, 'title'=>'Power spectrum' )
#    GGraph.tone( power[1..1199,true],true, 'title'=>'Power spectrum', 'max'=>comax, 'min'=>comin )
#    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