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

# =begin
# = Nakajima et al. (1992) の Tg vs OLR の(ような)図を描く
#   特定の実験ディレクトリ内で実行する.
# 
# == 概要
# 
# * 折れ線
#   * Nakajima et al. (1992) の Tg vs OLR の図念頭に, 
# 
# * 散布図
#   * 平均: 全球, 昼半球, 夜半球
#   * 各格子点 (少なくとも昼半球と夜半球の区別, 
#     恒星直下点と対蹠点がどこかは分かるようにする)
# 
# == 更新履歴
# 
# * 2014/10/04  石渡 正樹   単一ディレクトリ用に変更
# * 2011/04/25  納多 哲史   新規作成
# 
# =end

require "getoptlong"        # for option_parse
require "numru/ggraph"
require "pp"
include NumRu


## オプション解析
parser = GetoptLong.new
parser.set_options(
  ###    global option   ###
  ['--area',                      GetoptLong::REQUIRED_ARGUMENT],
  ['--time_start',                GetoptLong::REQUIRED_ARGUMENT],
  ['--time_end',                  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)

time_start = ($OPT_time_start||1).to_f
time_end = ($OPT_time_end||9999).to_f

area = ($OPT_area||'180').to_s

#RH_SCATTER_MAX = 1.0
#RH_SCATTER_MIN = 0.4
SCATTER_MAX = 3
SCATTER_MIN = 0

RPlanet = 6.371 * 10 ** 6  # [m]

TITLE = 'OLR'
YNAME = 'OLR'
#YUNIT = '10^15 W'

longname_rh='relative humidity'
name_rh= 'RH'
units_rh= '1'

P_0 = 1.4 * 1.0e11          # NHA92 で用いる基準気圧 [Pa]
Grav = 9.8                  # 重力加速度 [m s-2]
GasRUniv = 8.314            # 普遍気体定数 [J K-1 mol-1]
MolWtWet = 18.01528e-3      # 凝結成分の平均分子量 [kg mol-1]
MolWtDry = 28.964e-3        # 乾燥成分の平均分子量 [kg mol-1]
R_v = GasRUniv / MolWtWet   # 凝結成分の気体定数, 
EpsV = MolWtWet / MolWtDry  # 凝結成分と大気の分子量比
LatentHeat = 43655.0         # 凝結成分の潜熱 [J mol]

#XMIN = 170
#XMAX = 340
#YMIN = 90
#YMAX = 450

#XMIN = 250
#XMAX = 350
#YMIN = 200
#YMAX = 400

XMIN = 250
XMAX = 320
YMIN = 100
YMAX = 330

#MARK_SIZE = 0.01

vars = [
        {'name'=>'SurfTemp', 'min'=>240, 'max'=>320},
#        {'name'=>'Evap', 'min'=>0, 'max'=>600},
#        {'name'=>'Rain', 'min'=>0, 'max'=>500},
        {'name'=>'OLRA', 'min'=>270, 'max'=>400},
#        {'name'=>'Sens', 'min'=>-20, 'max'=>100},
#        {'name'=>'SLR', 'min'=>0, 'max'=>200}
       ]

# 2 次元以上のデータを 1 次元配列に落とす (GGraph.scatter 対策)
def flatten_gphys(gp)
  na_dat = gp.val.flatten
  n = na_dat.size
  axis = dummy_axis(n)
  data = VArray.new( na_dat,
                     {"long_name"=>gp.data.long_name, 
                       "units"=>gp.data.units.to_s},
                     gp.data.name )
  gp = GPhys.new( Grid.new(axis), data )
end

# ダミーの gphys を作成 (GGraph.scatter での図の範囲指定のため)
def dummy_gphys(na)
  n = na.size
  axis = dummy_axis(n)
  data = VArray.new( na,
                     {"long_name"=>"dummy", 
                       "units"=>"1"},
                     "dummy" )
  gp = GPhys.new( Grid.new(axis), data )
end

def dummy_axis(n)
  axis_a = VArray.new( NArray.sfloat(n).indgen,
                       {"long_name"=>"dummy","units"=>"1"},
                       "dummy" )
  axis = Axis.new.set_pos(axis_a)
end

# 領域平均. ただし lon のみ cut したものを使っても構わない
def ave_area(gp, na_lat_weight)
  gp = gp.mean('lon') * na_lat_weight
  gp = 0.5 * gp.sum('lat')
  gp.val = [gp.val]
end


# 描画準備
DCL.uzfact(0.6)

hash_gp_Tg = Hash.new
hash_gp_OLR = Hash.new

  ### 読み込み
  # OLR and basic value
  var = 'OLRA'
  file = var + '.nc'
  gp_OLR = GPhys::IO.open(file, var)

  sig = GPhys::IO.open(file, 'sig')
  sig_weight = GPhys::IO.open(file, 'sig_weight')
  na_lon = GPhys::IO.open(file, 'lon').val
  na_lat = GPhys::IO.open(file, 'lat').val
  na_lat_weight = GPhys::IO.open(file, 'lat_weight').val
  nlon = na_lon.size
  nlat = na_lat.size
  nsig = sig.val.size

# SurfTemp
  var = 'SurfTemp'
  file = var + '.nc'
  gp_Tg = GPhys::IO.open(file, var)

## Temp
#  var = 'Temp'
#  file = var + '.nc'
#  path = File.join(DATA_HOME, dir, file)
#  gp_t = GPhys::IO.open(path, var)
#
## QVap
#  var = 'QVap'
#  file = var + '.nc'
#  path = File.join(DATA_HOME, dir, file)
#  gp_q = GPhys::IO.open(path, var)
#
## Ps
#  var = 'Ps'
#  file = var + '.nc'
#  path = File.join(DATA_HOME, dir, file)
#  gp_ps = GPhys::IO.open(path, var)

if (area == '360') then
  hash_gp_OLR = gp_OLR.cut('time'=>time_start..time_end).mean('time')
  hash_gp_Tg  = gp_Tg.cut('time'=>time_start..time_end).mean('time')
elsif (area == '180') then # 直下点を中心とする 120度×120度領域
  hash_gp_OLR = gp_OLR.cut('time'=>time_start..time_end,'lon'=>0..180, 'lat'=>-90..90).mean('time')
  hash_gp_Tg  = gp_Tg.cut('time'=>time_start..time_end,'lon'=>0..180, 'lat'=>-90..90).mean('time')
elsif (area == '120') then # 直下点を中心とする 120度×120度領域
  hash_gp_OLR = gp_OLR.cut('time'=>time_start..time_end,'lon'=>30..150, 'lat'=>-60..60).mean('time')
  hash_gp_Tg  = gp_Tg.cut('time'=>time_start..time_end,'lon'=>30..150, 'lat'=>-60..60).mean('time')
elsif (area == '90') then # 直下点を中心とする 90度×90度領域
  hash_gp_OLR = gp_OLR.cut('time'=>time_start..time_end,'lon'=>45..135, 'lat'=>-45..45).mean('time')
  hash_gp_Tg  = gp_Tg.cut('time'=>time_start..time_end,'lon'=>45..135, 'lat'=>-45..45).mean('time')
else
  p 'area is not set properly: area=', area
  exit
end

# 範囲指定のためのダミーデータ
na_xrange = NArray.to_na([XMIN, XMAX])
na_yrange = NArray.to_na([YMIN, YMAX])
gp_xrange = dummy_gphys(na_xrange)
gp_yrange = dummy_gphys(na_yrange)
gp_xrange.data.long_name = gp_Tg.data.long_name
gp_yrange.data.long_name = gp_OLR.data.long_name
gp_xrange.data.units = gp_Tg.data.units
gp_yrange.data.units = gp_OLR.data.units


# 描画準備
DCL.gropn(wsn)

GGraph.set_fig('viewport'=>[0.3,0.6,0.3,0.5])
GGraph.set_axes('ytitle'=>YNAME)

DCL.sgpset('lfull', true)     # 全画面表示
DCL.sgpset('lcntl', false)

# 図の範囲を指定する方法が分からないのでダミーの点を打つ
GGraph.scatter( gp_xrange, gp_yrange, true, 'type'=>1)

opt = {
#   'max'=>VAL_MAX, 'min'=>VAL_MIN,  \
  'annotate'=>false,
  'title'=>TITLE,
  'index'=>1
}

opt['type'] = 1  # 4
opt['min'] = SCATTER_MIN
opt['max'] = SCATTER_MAX

gp_Tg = flatten_gphys( hash_gp_Tg.cut('lon'=>0..179.99) )
gp_OLR = flatten_gphys( hash_gp_OLR.cut('lon'=>0..179.99) )
gp_color = gp_Tg.copy
gp_color.val = gp_color.val.fill(0)
GGraph.color_scatter( gp_Tg, gp_OLR, gp_color, 
                      false, opt)

DCL.grcls

exit(0)