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

require "numru/ggraph"
include NumRu
include GGraph

velx_file = "../time_000000000-000172000/BS1998-all_VelX.nc"
velx_var ="VelX"
vely_file = "../time_000000000-000172000/BS1998-all_VelY.nc"
vely_var ="VelY"
velz_file = "../time_000000000-000172000/BS1998-all_VelZ.nc"
velz_var ="VelZ"
dens_file = "./restart_long-DensBZ.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

#dens = dens.mean('x').mean('y')

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


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

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