# -*- 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_000380000-000580000"
filehead = "BS1998_"
zcut = 10 # km で指定
cuttime = 380000..580000
#cutk = 2..3
ins = 2
itr = 1
puts "fileopen"
comax = 1.50e-2
comin = 0.0
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")
power = power#.cut('k'=>cutk)
#power = sp#.cut('k'=>cutk)

power_01 = power.cut('k'=>0.03490658..0.06981317)
#power_01 = power[1..2,0..-1]
power_01 = power_01.mean('k')


power_02 = power.cut('k'=>0.1047198..0.3141593)
#power_02 = power[3..9,0..-1]
power_02 = power_02.mean('k')

power_03 = power.cut('k'=>0.3490658..1.500983)
power_03 = power_03.mean('k')

power_04 = power.cut('k'=>1.53589..41.853)
power_04 = power_04.mean('k')

p power
DCL.gropn(ins)
DCL.uzfact(0.5)
DCL.sgpset('lfull',true)
#GGraph.set_fig( 'itr'=> 4)
DCL.swlset( 'ldump', true )
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('ytitle'=>'PTemp_PwSP')
#    GGraph.tone( power.cut('k'=>cutk,'t'=>cuttime),true )
GGraph.line( power_01, 
             true, 'exchange'=>false , 
             'index'=>21, 'type'=>1, #'label'=>'1',
             'max'=>comax, 'min'=>comin)   

GGraph.line( power_02, 
             false, 'exchange'=>false , 
             'index'=>31, 'type'=>1, #'label'=>'1',
             'max'=>comax, 'min'=>comin)   

GGraph.line( power_03, 
             false, 'exchange'=>false , 
             'index'=>41, 'type'=>1, #'label'=>'1',
             'max'=>comax, 'min'=>comin)   

GGraph.line( power_04, 
             false, 'exchange'=>false , 
             'index'=>61, 'type'=>1, #'label'=>'1',
             'max'=>comax, 'min'=>comin)   

#    GGraph.contour( power[1..1199,i], false )
#    GGraph.color_bar('vlength'=>0.5,"landscape"=>true,'tickintv'=>0)
=begin
for i in 2..3 do
#  if ( i.to_f % 20) == 0
#GGraph.set_fig( 'itr'=> itr, 'viewport'=>[0.25,0.7,0.15,0.6] )
  GGraph.line( power[i,1..500], 
               false, 'exchange'=>false , 
               'index'=>2, 'type'=>1, 'label'=>i.to_s,
               'max'=>comax, 'min'=>comin)   
end
=end

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

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