#!/usr/bin/env ruby
=begin
= mkfig_25time_series.rb -- make scatter plot for Hadley vs.Ferrel 

=end

require "getopts"        # for option_parse
require "numru/ggraph"
include NumRu
include Misc::EMath

def cal_correlation(x, y) # 相関係数を求める
  xa = x - x.mean(0)           # x の偏差
  ya = y - y.mean(0)           # y の偏差
  xstddev = x.stddev           # x の標準偏差(分散の平方根)
  ystddev = y.stddev           # y   ::
  
  crc = (( xa*ya ).mean(0))/(xstddev*ystddev) # 相関係数 
  return crc
end


################################################################
#                        make gphys 
################################################################

year = [1979, 2003]
month = ["01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12"]   # 初期値 




############      stream function        #########################

strm_nh_hadley_long_name = "hadley circulation intencity (NH)"
strm_sh_hadley_long_name = "hadley circulation intencity (SH)"

## make compisit-GPhys object 
# make array which set the gphys objects.

axmonth = []
i = 0

strm_name = ""
strm_unit = ""


strm_nh_hadley_long_name = "hadley circulation intencity (NH)"
strm_sh_hadley_long_name = "hadley circulation intencity (SH)"
strm_nh_ferrel_long_name = "ferrel circulation intencity (NH)"
strm_sh_ferrel_long_name = "ferrel circulation intencity (SH)"

## make compisit-GPhys object 
# make array which set the gphys objects.

axmonth = []
i = 0

strm_name = ""
strm_unit = ""

strm_hadley_north = []
strm_hadley_south = []
strm_ferrel_north = []
strm_ferrel_south = []

year[0].upto(year[1]) do |y|
  month.each do |m|
    path = "../../../STRM.NCEP/STRM.#{y}.NCEP/STRM_#{y}-#{m}_NCEP.nc" # stream function
    gp_climat = GPhys::NetCDF_IO.open("../../25MEANS.NCEP/25MEANS.1979-2003.NCEP/100hPa_PHYS_1979-2003_#{m}_NCEP.nc", 'strm')
    gp_strm =   GPhys::NetCDF_IO.open(path, 'strm').cut('level'=>100..1000) 

    gp_strm = GPhys::NetCDF_IO.open(path, 'strm')
#    if m =="02"
      strm_hadley_north <<  gp_strm.cut( 'lat'=>  0..30  ).data.val.max - gp_climat.cut( 'lat'=>  0..30  ).data.val.max
      strm_ferrel_north <<  gp_strm.cut( 'lat'=> 30..60  ).data.val.min - gp_climat.cut( 'lat'=> 30..60  ).data.val.min
#    elsif m == "08"
      strm_hadley_south <<  gp_strm.cut( 'lat'=>  0..-30 ).data.val.min - gp_climat.cut( 'lat'=>  0..-30 ).data.val.min
      strm_ferrel_south <<  gp_strm.cut( 'lat'=>-30..-60 ).data.val.max - gp_climat.cut( 'lat'=>-30..-60 ).data.val.max
#    end
# p gp_strm.cut( 'lat'=>0..30  ).data.val.max if m == "02"
# p gp_strm.cut( 'lat'=>0..-30 ).data.val.min if m == "08"

    strm_name =      gp_strm.data.name
    strm_unit =      gp_strm.data.get_att("units")

    axmonth << DCL::dateg3(year[0],1,1,y.to_i,m.to_i,1)
    i += 1
  end
end


### for intensity

  na_strm_hadley_north = NArray.to_na(strm_hadley_north)
  na_strm_hadley_south = NArray.to_na(strm_hadley_south)
  na_strm_ferrel_north = NArray.to_na(strm_ferrel_north)
  na_strm_ferrel_south = NArray.to_na(strm_ferrel_south)
                   
################################################################
#                        描画ルーチン
################################################################



##
# 事前準備

DCL.uscset('cyspos', 'B' )              # y 軸の単位の位置を下方へ 
rsizel2 = DCL.uzrget('rsizel2')         # 現在のラベルサイズを取得
DCL.uzrset('rsizel2', rsizel2*0.42 )    # ラベルサイズをデフォルトの 0.5 倍に

##
# お絵描きメイン

DCL.swpset('lsep',  true)    # ページ毎別々のファイルに落す
DCL.gropn(-2)
DCL.sldiv('T', 1, 2)
DCL.sgpset('lcntl', false)   # 
DCL.sgpset('lfull',true)     # フルスクリーン
DCL.uzfact(0.55)             # 
DCL.sgpset('lcorner',false)  # コーナーを取っ払う 
DCL.sgpset('lfprop',true)    # 
DCL.udpset('lmsg',false)     # 
DCL.uscset('cyspos', 'B' )   # y 軸の単位の位置を下方へ 

## y 軸範囲

axx = NArray[-5e10, 5e10] # 
axy = NArray[-5e10, 5e10] # 


### 1 枚目 : 北半球

DCL::grfrm                   # ページ確定 

x = na_strm_hadley_north
y = na_strm_ferrel_north

DCL::usspnt(axx, axy)
DCL::uspfit
DCL::grstrf

DCL::ussttl('Deviation of Hadley intencity (NH)', strm_unit, 'Deviation of Ferrel intencity ', strm_unit )
DCL::usdaxs

DCL::uusmkt(2)
DCL::uumrk(x, y)

print "北半球ハドレー循環強度とフェレル循環強度", cal_correlation(x, y), "\n"


### 2 枚目 : 南半球

DCL::grfrm                   # ページ確定 

x = na_strm_hadley_south
y = na_strm_ferrel_south

DCL::usspnt(axx, axy)
DCL::uspfit
DCL::grstrf

DCL::ussttl('Deviation of mass stream function (SH)', strm_unit, 'Deviation of Ferrel intencity ', strm_unit )
DCL::usdaxs

DCL::uusmkt(2)
DCL::uumrk(x, y)

print "北半球ハドレー循環強度とフェレル循環強度", cal_correlation(x, y), "\n"

DCL.grcls

##################################################
# 画像ファイル名を変更
File.rename('dcl_001.ps', 'FERREL_1973-2003_TIME-SERIES_NCEP_DEVIATION_sanpuzu.ps')