#!/usr/bin/ruby
# -*- coding: utf-8 -*-

#DATA_HOME = "/home/noda/data"
#DATA_HOME = "/home/noda/data/exps_101220_diff_T21L16_sr_Omega_sweep"

require "numru/ggraph"
include NumRu
require "pp"

wsn = 4

# ベクトルの描画間隔 (単位は格子点数)
xintv = 6
yintv = 2

#ITR=30   # 球
ITR=1

LON_CNT=125  # 球のときに使用

#xfact1 = 0.01
#yfact1 = 0.01
xfact1 = 1.0
yfact1 = 1.0

# 凡例の単位ベクトルの長さ
uxunit = 0.03
uyunit = 0.03

time_start = 1000
time_end   = 2000

#omegas = [0.0, 0.15, 0.5, 1.0]
#omegas = [1.0]
#omegas = [0.0]

omega = 1.0

type = 1  # 1: mean, 2: disturbance

cp = 1004.6

# 簡易
while ARGV.size > 0
  case ARGV[0]
  when '--omega' then
    ARGV.shift
    omega = ARGV[0]
    ARGV.shift
  when '--wsn' then
    ARGV.shift
    wsn = ARGV[0]
    ARGV.shift
  when '--xintv' then
    ARGV.shift
    xintv = ARGV[0]
    ARGV.shift
  when '--yintv' then
    ARGV.shift
    yintv = ARGV[0]
    ARGV.shift
  when '--xfact1' then
    ARGV.shift
    xfact1 = ARGV[0]
    ARGV.shift
  when '--yfact1' then
    ARGV.shift
    yfact1 = ARGV[0]
    ARGV.shift
  when '--uxunit' then
    ARGV.shift
    uxunit = ARGV[0]
    ARGV.shift
  when '--uyunit' then
    ARGV.shift
    uyunit = ARGV[0]
    ARGV.shift
  when '--type' then
    ARGV.shift
    type = ARGV[0].to_i
    ARGV.shift
  when '--range_time' then
    ARGV.shift
    tmp = ARGV[0].split(':')
    time_start = tmp[0].to_f
    time_end   = tmp[1].to_f
    ARGV.shift
  else
    raise "ERROR: invalid option: #{ARGV[0]}"
  end
end



puts omega

if type == 1 then
  puts "all"
elsif type == 2 then
  puts "disturbance"
else
  raise "Type #{type} is not supported."
end

#dir = File.join(DATA_HOME, 
#      "101220_diff-DryAir-Vapor_T21L16_sr_Omega" + omega.to_s + "_dt20m")

dir='.'


# geopotential 
temperature_file = dir + '/Temp.nc'
system("ruby /GFD_Dennou_Work3/momoko/SyncRotEarthRad/script/geopotential_time.rb --var Temp.nc@Temp --range_time #{time_start}:#{time_end}")


# input
path = File.join(dir, 'U.nc')
na_sig_weight = GPhys::IO.open(path, 'sig_weight').val
nsig = na_sig_weight.size

var = 'U'
path = File.join(dir, var+'.nc')
gp_u = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

var = 'V'
path = File.join(dir, var+'.nc')
gp_v = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)

var = 'Temp'
path = File.join(dir, var+'.nc')
gp_h = GPhys::IO.open(path, var).cut('time'=>time_start..time_end)*cp
gp_h = gp_h.convert_units('J kg-1')

# Phi
gp_phi = GPhys::IO.open('Phi.nc', 'Phi').cut('time'=>time_start..time_end)

# dry static energy
gp_q = gp_h + gp_phi


gp_u_q = gp_u * gp_q
gp_v_q = gp_v * gp_q

# 時間平均
gp_u_q = gp_u_q.mean('time')
gp_v_q = gp_v_q.mean('time')
if type == 2 then
  gp_u = gp_u.mean('time')
  gp_v = gp_v.mean('time')
  gp_q = gp_q.mean('time')
end

# 擾乱成分
if type == 2 then
  gp_ud_qd = gp_u_q - gp_u * gp_q
  gp_vd_qd = gp_v_q - gp_v * gp_q
end

# 鉛直積分
if type == 1 then
  gp_u_q = (gp_u_q * na_sig_weight.reshape(1,1,nsig)).sum('sig')
  gp_v_q = (gp_v_q * na_sig_weight.reshape(1,1,nsig)).sum('sig')
elsif type == 2 then
  gp_ud_qd = (gp_ud_qd * na_sig_weight.reshape(1,1,nsig)).sum('sig')
  gp_vd_qd = (gp_vd_qd * na_sig_weight.reshape(1,1,nsig)).sum('sig')
end

# for debug
#pp gp_u_q


# 後始末
#File.unlink('./Phi.nc')


DCL.gropn(wsn)

#DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7)
DCL.sgpset('lcntl', false) ; DCL.uzfact(1.0)

GGraph.set_axes('xtickint'=>10, 'xlabelint'=>30)
GGraph.set_axes('ytickint'=>10, 'ylabelint'=>30)

# 球の時
#GGraph.set_fig 'itr'=>ITR, 'viewport'=>[0.15,0.95,0.2,0.8], 'map_axis'=>[LON_CNT, 0, 0]
#GGraph.set_fig 'viewport'=>[0.15,0.95,0.2,0.8]

# same as gpview
GGraph.set_fig 'itr'=>ITR, 'viewport'=>[0.15,0.85,0.2,0.55], 'map_axis'=>[LON_CNT, 0, 0]


DCL.ugpset('LNRMAL', false)  # scale vector automatically

DCL.ugpset('XFACT1', xfact1)
DCL.ugpset('YFACT1', yfact1)

DCL.ugpset('UXUNIT', uxunit) # size of unit vector by dimentional value
DCL.ugpset('UYUNIT', uyunit)

DCL::ugrset('VXUNIT', 0.1)  # size of unit vector
DCL::ugrset('VYUNIT', 0.1)

DCL::uxmttl('T', '', 0.0)  # title (large character)

DCL.ugpset('LMSG', false)   # no message

DCL.ugpset('INDEX', 3)

opt_vector = {
  'xintv'=>xintv, 'yintv'=>yintv, 'title'=>'',
  'flow_vect'=>false, 'annotate'=>false, 'unit_vect'=>true
}

if type == 1 then
  # mean
  GGraph.vector( gp_u_q, gp_v_q, true, opt_vector)
elsif type == 2 then
  # disturbance
  GGraph.vector( gp_ud_qd, gp_vd_qd, true, opt_vector)
else
  raise "Type #{type} is not supported."
end

DCL.grcls