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

# = anim_qvap-sigdot-uv.rb
# == 説明
#  QVap のトーン, SigDot のコンター, (u,v) のベクトルを重ねがきする.

#DATA_HOME = "/home/noda/data"
DATA_HOME = "/GFD_Dennou_Work3/momoko/SyncRotEarthRad"

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

RPlanet = 6.371 * 10 ** 6  # [m]
#Ps = 1e5  # [Pa]
Grav = 9.8   # [m/s2]

title = ""

#YUNIT_vapor = 'kg m-2'
#epsfn = "qvap_area_column_mass.eps"

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

# 球面プロットと平面プロットの切り替えは
# この部分と下の設定部分を変更する.

#ITR=30   # 球
ITR=1   # 平面プロット

#LON_CNT=125  # 球のときに使用 (野田くん)
#LON_CNT=90  # 球のときに使用

lon_cont = 125

dir = nil


## オプション解析 new version ここから
parser = GetoptLong.new
parser.set_options(
   ###    global option   ###
   ['--time_start',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--time_end',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--xfact1',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--yfact1',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--uxunit',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--uyunit',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--val_min',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--val_max',                      GetoptLong::REQUIRED_ARGUMENT],
   ['--wsn',                      GetoptLong::REQUIRED_ARGUMENT]
                    )
parser.each_option do |name, arg|
   eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'"  # strage option value to $OPT_val
end

wsn = ($OPT_wsn||4)

xfact1 = ($OPT_xfact1||1.0).to_f
yfact1 = ($OPT_yfact1||1.0).to_f

# 凡例の単位ベクトルの長さ (デフォルト値)
uxunit = ($OPT_uxunit||0.03).to_f
uyunit = ($OPT_uyunit||0.03).to_f

# 時間範囲
time_start = ($OPT_time_start||1).to_f
time_end = ($OPT_time_end||20).to_f

# SigDot の描画範囲
#val_min = -1e-7
#val_max = 1e-7
val_min = ($OPT_val_min||-1e-7).to_f
val_max = ($OPT_val_max||1e-7).to_f

val_interval = (val_max - val_min) / 30.0


# 簡易
while ARGV.size > 0
  case ARGV[0]
  when '--dir' then
    ARGV.shift
    dir = 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 '--loncnt' then
    ARGV.shift
    lon_cnt = ARGV[0]
    ARGV.shift
  else
    raise "ERROR: invalid option: #{ARGV[0]}"
  end
end


if dir == nil then
  dir = File.join(DATA_HOME,"ExpOcean")
end


# input
var = 'U'
path = var+'.nc'
gp_u = GPhys::IO.open(path, var).cut('time'=>time_start..time_end, 'sig'=>0.1)

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

var = 'QVap'
path = var+'.nc'
gp_qvap = GPhys::IO.open(path, var).cut('time'=>time_start..time_end, 'sig'=>0.1)

var = 'SigDot'
path = var+'.nc'
gp_sigdot = GPhys::IO.open(path, var).cut('time'=>time_start..time_end, 'sigm'=>0.1)


# 時間平均しない
#gp_u = gp_u.mean('time')
#gp_v = gp_v.mean('time')
#gp_sigdot = gp_sigdot.mean('time')

DCL.gropn(wsn)

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

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

# 球の時 (平面プロットの時は 下の 2 行をコメントアウトするだけ)
if ITR == '30' then
  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]
end

# same as gpview
vx0,vx1,vy0,vy1 = 0.15,0.85,0.2,0.55
GGraph.set_fig 'itr'=>ITR, 'viewport'=>[vx0,vx1,vy0,vy1], 'map_axis'=>[lon_cnt, 0, 0]

DCL.uscset('cyspos', 'B' )              # move unit y axis

DCL.udlset('LMSG', false)

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('LNRMAL', false)  # scale vector automatically
#DCL.ugpset('LMSG', false)   # no message
DCL.ugpset('LUMSG', false)   # no message

# 矢印の設定
DCL.ugpset('INDEX', 5)
#DCL.ugpset('INDEX', 959)

# オプション設定
opt_vector = {
  'xintv'=>xintv, 'yintv'=>yintv,
  'flow_vect'=>false, 'annotate'=>true, 'unit_vect'=>true
#  'flow_vect'=>true, 'annotate'=>false, 'unit_vect'=>true
}

opt_contour = {
  'min'=>val_min, 'max'=>val_max, 'interval'=>val_interval,
  'title'=>title, 'annotate'=>false
}

# カラーマップの設定
#DCL.sgscmn(9)   # colormap: white-blue-black
#DCL.sgscmn(8)   # colormap: white-yellow-red-black
#DCL.sgscmn(2)   # colormap: black-red-yellow-white
DCL.sgscmn(7)   # colormap: 

# 描画
for d in gp_u.coord('time').val
  GGraph.tone( gp_qvap.cut('time'=>d), true)
  GGraph.contour( gp_sigdot.cut('time'=>d), false)
  GGraph.vector( gp_u.cut('time'=>d), gp_v.cut('time'=>d), false, opt_vector)
  GGraph.color_bar('landscape'=>true, 'voff'=>0.1)
end


DCL.grcls

exit

# label of vector
# Ref.:
# * definition of unit_vect in dclext.rb
# * http://ruby.gfd-dennou.org/products/ruby-dcl/ruby-dcl-doc/grph1/node62.html

space = 0.0
vxuloc = vx1 + space
vyuloc = vy0 + space
rsizet = 0.017
index = 1

p xunit = (uxunit.to_f * 1e7).round / 1e7
p yunit = (uyunit.to_f * 1e7).round / 1e7

# mean
sfxunit = [xunit.to_s, "m s-1"].join(' ')
sfyunit = [yunit.to_s, "m s-1"].join(' ')

DCL.sgtxzv(vxuloc, vyuloc-1.2*rsizet,
           sfxunit, rsizet, 0, -1, index)

DCL.sgtxzv(vxuloc+2.5*rsizet, vyuloc+0.5*rsizet,
           sfyunit, rsizet, 90, -1, index)

DCL.grcls
