#!/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 = 1
time_end   = 9999

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

ofile = NetCDF.create(dir + 'AngMom.nc')

# 時間ループ
GPhys::NetCDF_IO.each_along_dims_write(gp_u, ofile, 'time') { 
  |gu|  

  gp_angularmom = gu.copy

#  p gp_angularmom

  # 角運動量の計算
  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.val = RPlanet*gu.val*lat_rad.cos+omega*RPlanet*RPlanet*na_cos_lat_rad_4d*na_cos_lat_rad_4d
  gp_angularmom.val = RPlanet*gu.val*na_cos_lat_rad_4d+omega*RPlanet*RPlanet*na_cos_lat_rad_4d*na_cos_lat_rad_4d

  gp_angularmom.units = 'm2 s-1'
  gp_angularmom.long_name = 'angular momentum'
  gp_angularmom.name = 'AngMom'

  [gp_angularmom]
}

ofile.close
exit


# 以下, 描画用. 
# このスクリプトはデータを作るだけにしたので
# 以下は使わない.


# 描画
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