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

# 2014/09/14 Ishiwatari  modified

# 角運動量の式は
#   l = a u \cos \theta + \Omega a^2 \cos^2 \theta

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

wsn = 4

RPlanet = 6.371 * 10 ** 6  # [m]
Grav = 9.8   # [m/s2]

PI=3.141592

title = ""

# vapor の描画範囲
val_min = 0
val_max = 3e9
#val_interval = (val_max - val_min) / 30.0
val_interval = 1e8

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

LON_CNT=125  # 球のときに使用

time_start = 400
time_end   = 500

omega = 7.292106590880652e-5

dir = nil

# 簡易
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 '--range_time' then
    ARGV.shift
    tmp = ARGV[0].split(':')
    time_start = tmp[0].to_f
    time_end   = tmp[1].to_f
    ARGV.shift
  when '--range_val' then
    ARGV.shift
    tmp = ARGV[0].split(':')
    val_min = tmp[0].to_f
    val_max = tmp[1].to_f
    ARGV.shift
  when '--cint' then
    ARGV.shift
    val_interval = ARGV[0]
    ARGV.shift
  else
    raise "ERROR: invalid option: #{ARGV[0]}"
  end
end

dir='.'

# データ読み込み
path = File.join(dir, 'U.nc')

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

lat = GPhys::IO.open(path, 'lat')
nlat = lat.val.size

# 角運動量の計算
lat_rad = lat*PI/180.0


cos_lat_rad = lat_rad.cos
na_cos_lat_rad = cos_lat_rad.val

na_cos_lat_rad_4d = na_cos_lat_rad.reshape(1,nlat,1,1)

#p omega*RPlanet*RPlanet*na_cos_lat_rad_4d*na_cos_lat_rad_4d

gp_angularmom = RPlanet*gp_u*lat_rad.cos+omega*RPlanet*RPlanet*na_cos_lat_rad_4d*na_cos_lat_rad_4d

gp_angularmom_mean = gp_angularmom.mean('time').mean('lon')

# p gp_angularmom_mean

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)

# 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::uxmttl('T', '', 0.0)  # title (large character)

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

DCL.ugpset('INDEX', 9)

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

GGraph.tone_and_contour( gp_angularmom_mean, true, opt_contour)
GGraph.color_bar('landscape'=>true)

DCL.grcls