# -*- coding: utf-8 -*-
# 運動エネルギー密度を計算し, 絵を描く

require "numru/ggraph"
include NumRu
include GGraph

velx_file = "./BS1998-all_VelX.nc"
velx_var ="VelX"
vely_file = "./BS1998-all_VelY.nc"
vely_var ="VelY"
velz_file = "./BS1998-all_VelZ.nc"
velz_var ="VelZ"
dens_file = "./BS1998_restart_rank000000.nc"
dens_var = "DensBZ"



DCL.gropn(2)

velx = GPhys::IO.open("#{velx_file}", "#{velx_var}")
vely = GPhys::IO.open("#{vely_file}", "#{vely_var}")
velz = GPhys::IO.open("#{velz_file}", "#{velz_var}")
dens = GPhys::IO.open("#{dens_file}", "#{dens_var}")
time = GPhys::IO.open("#{velx_file}", 't').val
x_dco = GPhys::IO.open("#{dens_file}", 'x').val
x_co = GPhys::IO.open("#{velx_file}", 'x').val
z_co = GPhys::IO.open("#{velz_file}", 'z').val

p dens

denscut = dens[5..104, 0..-1, 5..104]
densallay = NArray.sfloat(400,1,100).fill(dens.cut)
denscut = denscut.mean('x')

puts denscut

kener = 0.50 * ( velx**2.0 + velz**2.0)
kener = kener.mean('x') * denscut
kener = kener.mean('z').mean('y')

p kener
puts "denscut"
p denscut
kenerdens = kener 
#kenerdens = kener / denscut


#=begin
GGraph.set_axes('xtitle'=>'time','ytitle'=>'K')
GGraph.set_fig('viewport'=>[0.1,0.8,0.4,0.6])
GGraph.line( kenerdens, true, 'title'=>'Kinetic energy density', 'max'=>30.0)

#GGraph.line( kenerdens.cut( 'y'=>100 ), true)
#=end
DCL.grcls