# -*- coding: utf-8 -*-
# Danilov の不等式の計算結果を一覧で表示する
require "numru/ggraph"
include NumRu
include Math

# ----- コマンドラインオプション -----
iws = ( ARGV[0] || 1 ).to_i     # 装置番号 (1,2,or 4)

# ----- 使用データ -----
file_res1_ek = './case1/ek_analysis/ek-analysis.nc'
file_res1_ens = './case1/ens_analysis/ens-analysis.nc'
file_res2_ek = './case2/ek-analysis.nc'
file_res2_ens = './case2/ens-analysis.nc'
file_res3_ek = './case3/ek-analysis.nc'
file_res3_ens = './case3/ens-analysis.nc'
file_res4_ek = './case4/ek-analysis.nc'
file_res4_ens = './case4/ens-analysis.nc'

# t0 = 1800
# t1 = 2000
# t2 = 300
# t3 = 500

# case1 の flux
data_res1_ek1 =  (GPhys::IO.open file_res1_ek , 'tro_flux1').mean("t")
data_res1_ek2 =  (GPhys::IO.open file_res1_ek , 'tro_flux2').mean("t")
data_res1_ek3 =  (GPhys::IO.open file_res1_ek , 'cli_flux1').mean("t")
data_res1_ek4 =  (GPhys::IO.open file_res1_ek , 'cli_flux2').mean("t")
data_res1_ek5 =  (GPhys::IO.open file_res1_ek , 'cli_flux3').mean("t")
data_res1_ek = data_res1_ek1 + data_res1_ek2 + data_res1_ek3 + data_res1_ek4 + data_res1_ek5
data_res1_ens1 =  (GPhys::IO.open file_res1_ens , 'tro_flux1').mean("t")
data_res1_ens2 =  (GPhys::IO.open file_res1_ens , 'tro_flux2').mean("t")
data_res1_ens3 =  (GPhys::IO.open file_res1_ens , 'cli_flux1').mean("t")
data_res1_ens4 =  (GPhys::IO.open file_res1_ens , 'cli_flux2').mean("t")
data_res1_ens5 =  (GPhys::IO.open file_res1_ens , 'cli_flux3').mean("t")
data_res1_ens = data_res1_ens1 + data_res1_ens2 + data_res1_ens3 + data_res1_ens4 + data_res1_ens5
# case2 の flux
data_res2_ek1 =  (GPhys::IO.open file_res2_ek , 'tro_flux1').mean("t")
data_res2_ek2 =  (GPhys::IO.open file_res2_ek , 'tro_flux2').mean("t")
data_res2_ek3 =  (GPhys::IO.open file_res2_ek , 'cli_flux1').mean("t")
data_res2_ek4 =  (GPhys::IO.open file_res2_ek , 'cli_flux2').mean("t")
data_res2_ek5 =  (GPhys::IO.open file_res2_ek , 'cli_flux3').mean("t")
data_res2_ek = data_res2_ek1 + data_res2_ek2 + data_res2_ek3 + data_res2_ek4 + data_res2_ek5
data_res2_ens1 =  (GPhys::IO.open file_res2_ens , 'tro_flux1').mean("t")
data_res2_ens2 =  (GPhys::IO.open file_res2_ens , 'tro_flux2').mean("t")
data_res2_ens3 =  (GPhys::IO.open file_res2_ens , 'cli_flux1').mean("t")
data_res2_ens4 =  (GPhys::IO.open file_res2_ens , 'cli_flux2').mean("t")
data_res2_ens5 =  (GPhys::IO.open file_res2_ens , 'cli_flux3').mean("t")
data_res2_ens = data_res2_ens1 + data_res2_ens2 + data_res2_ens3 + data_res2_ens4 + data_res2_ens5
# case3 の flux
data_res3_ek1 =  (GPhys::IO.open file_res3_ek , 'tro_flux1').mean("t")
data_res3_ek2 =  (GPhys::IO.open file_res3_ek , 'tro_flux2').mean("t")
data_res3_ek3 =  (GPhys::IO.open file_res3_ek , 'cli_flux1').mean("t")
data_res3_ek4 =  (GPhys::IO.open file_res3_ek , 'cli_flux2').mean("t")
data_res3_ek5 =  (GPhys::IO.open file_res3_ek , 'cli_flux3').mean("t")
data_res3_ek = data_res3_ek1 + data_res3_ek2 + data_res3_ek3 + data_res3_ek4 + data_res3_ek5
data_res3_ens1 =  (GPhys::IO.open file_res3_ens , 'tro_flux1').mean("t")
data_res3_ens2 =  (GPhys::IO.open file_res3_ens , 'tro_flux2').mean("t")
data_res3_ens3 =  (GPhys::IO.open file_res3_ens , 'cli_flux1').mean("t")
data_res3_ens4 =  (GPhys::IO.open file_res3_ens , 'cli_flux2').mean("t")
data_res3_ens5 =  (GPhys::IO.open file_res3_ens , 'cli_flux3').mean("t")
data_res3_ens = data_res3_ens1 + data_res3_ens2 + data_res3_ens3 + data_res3_ens4 + data_res3_ens5
# case4 の flux
data_res4_ek1 =  (GPhys::IO.open file_res4_ek , 'tro_flux1').mean("t")
data_res4_ek2 =  (GPhys::IO.open file_res4_ek , 'tro_flux2').mean("t")
data_res4_ek3 =  (GPhys::IO.open file_res4_ek , 'cli_flux1').mean("t")
data_res4_ek4 =  (GPhys::IO.open file_res4_ek , 'cli_flux2').mean("t")
data_res4_ek5 =  (GPhys::IO.open file_res4_ek , 'cli_flux3').mean("t")
data_res4_ek = data_res4_ek1 + data_res4_ek2 + data_res4_ek3 + data_res4_ek4 + data_res4_ek5
data_res4_ens1 =  (GPhys::IO.open file_res4_ens , 'tro_flux1').mean("t")
data_res4_ens2 =  (GPhys::IO.open file_res4_ens , 'tro_flux2').mean("t")
data_res4_ens3 =  (GPhys::IO.open file_res4_ens , 'cli_flux1').mean("t")
data_res4_ens4 =  (GPhys::IO.open file_res4_ens , 'cli_flux2').mean("t")
data_res4_ens5 =  (GPhys::IO.open file_res4_ens , 'cli_flux3').mean("t")
data_res4_ens = data_res4_ens1 + data_res4_ens2 + data_res4_ens3 + data_res4_ens4 + data_res4_ens5

name1 = 'k|2"' + DCL.csgi(143) + '_E"'
name2 = DCL.csgi(143) + '_Q"'
name3 = 'k|2"' + DCL.csgi(143) + '_E"' + DCL.csgi(190) + DCL.csgi(143) + '_Q"'

# k^2\Pi_E のGphys オブジェクト作成
# case1
va_wn = data_res1_ek.coord("n")
axis_wn = Axis.new.set_pos(va_wn)
wn_ek_flux_tmp1 = NArray.float(341)
ek_flux  = data_res1_ek.val
for k in 0..339
  wn_ek_flux_tmp1[k] = (k+1)**2 * ek_flux[k]
end
wn_ek_flux_tmp2 = VArray.new(wn_ek_flux_tmp1, nil,name1 )
data_res1 = GPhys.new(Grid.new(axis_wn),wn_ek_flux_tmp2)
# case2
va_wn = data_res2_ek.coord("n")
axis_wn = Axis.new.set_pos(va_wn)
wn_ek_flux_tmp1 = NArray.float(170)
ek_flux  = data_res2_ek.val
for k in 0..168
  wn_ek_flux_tmp1[k] = (k+1)**2 * ek_flux[k]
end
wn_ek_flux_tmp2 = VArray.new(wn_ek_flux_tmp1, nil,name1 )
data_res2 = GPhys.new(Grid.new(axis_wn),wn_ek_flux_tmp2)
# case3
va_wn = data_res3_ek.coord("n")
axis_wn = Axis.new.set_pos(va_wn)
wn_ek_flux_tmp1 = NArray.float(170)
ek_flux  = data_res3_ek.val
for k in 0..168
  wn_ek_flux_tmp1[k] = (k+1)**2 * ek_flux[k]
end
wn_ek_flux_tmp2 = VArray.new(wn_ek_flux_tmp1, nil,name1 )
data_res3 = GPhys.new(Grid.new(axis_wn),wn_ek_flux_tmp2)
# case4
va_wn = data_res4_ek.coord("n")
axis_wn = Axis.new.set_pos(va_wn)
wn_ek_flux_tmp1 = NArray.float(170)
ek_flux  = data_res4_ek.val
for k in 0..168
  wn_ek_flux_tmp1[k] = (k+1)**2 * ek_flux[k]
end
wn_ek_flux_tmp2 = VArray.new(wn_ek_flux_tmp1, nil,name1 )
data_res4 = GPhys.new(Grid.new(axis_wn),wn_ek_flux_tmp2)

check1 = data_res1 - data_res1_ens
check2 = data_res2 - data_res2_ens
check3 = data_res3 - data_res3_ens
check4 = data_res4 - data_res4_ens

check1.name = name3
data_res1_ens.name = name2
data_res1.name = name1

x = NArray.sfloat(341)
y = NArray.sfloat(341)
y.fill!(0.0)
x[0] = 1.0
# y[0] = 5.0*10**(-3)
for k in 0..339
    x[k+1] = x[k] + 1.0
#   y[k+1] = y[k]*1.0*10**(0.205)
end

x1 = NArray.sfloat(170)
y1 = NArray.sfloat(170)
y1.fill!(0.0)
x1[0] = 1.0
# y[0] = 5.0*10**(-3)
for k in 0..168
    x1[k+1] = x[k] + 1.0
#   y[k+1] = y[k]*1.0*10**(0.205)
end

# ----- DCL設定 -----
def prep_dcl(iws=1) # iws : DCL出力デバイス．1,4:画面, 2:PS
  DCL.sgscmn(39)                # カラーマップ番号
  DCL.swpset("iwidth",1000)     # 画面の幅
  DCL.swpset("iheight",1000)    # 画面の高さ
  DCL.swpset("lalt",true)      # 裏で描画（パラパラアニメ用）   
   DCL.swpset("ifl",1)          # png で出力
#  DCL.swpset("lsysfnt",true)   # システムフォントを使用
#  DCL.swcset("fontname","Hiragino Kaku Gothic")
  DCL.sglset("lcntl",true)     # 上付き下付き文字を有効
  DCL.sgiset("ifont",2)
  DCL.gropn(iws)
#  DCL.sldiv("y",2,2)           # 画面分割 (描画順("y"oko/"t"ate),数:横,数:縦)
  DCL.sgpset('lfull',true)     # 全画面表示
#  DCL.sgpset('isub', 96)       # 下付き添字を表す制御文字を '_' から '`' に
  DCL.uzfact(0.5)             # 座標軸の文字列サイズを定数倍
  DCL.sgpset('lcorner',false)  # コーナーマークを書かない
  DCL.glpset('lmiss',true)     # DCLの欠損値処理を on に
end

max = 8.5*10**(0)
min = -8.5*10**(0)

case1 = DCL.csgi(131) + DCL.csgi(164) + "= 0"
case2 = DCL.csgi(131) + DCL.csgi(164) + "= -5.280" + DCL.csgi(194) + '10|-17"'
case3 = DCL.csgi(131) + DCL.csgi(164) + "= 5.280" + DCL.csgi(194) + '10|-17"'
case4 = DCL.csgi(131) + DCL.csgi(164) + "= 1.560" + DCL.csgi(194) + '10|-16"'
ytitle = name1 + ',' + name2 + ',' + name3

vpt = NArray[0.15, 0.5, 0.15, 0.5]             # ビューポートサイズ (1:1)
vpt00 = ( vpt + ([0.0]*2 + [0.4]*2) ).to_a    # x,y方向にずらしてArray化
vpt01 = ( vpt + ([0.374]*2 + [0.4]*2) ).to_a    # x,y方向にずらしてArray化
vpt10 = ( vpt + ([0.0]*2 + [-0.02]*2) ).to_a    # x,y方向にずらしてArray化
vpt11 = ( vpt + ([0.374]*2 + [-0.02]*2) ).to_a    # x,y方向にずらしてArray化

index1 = 853
index2 = 503
index3 = 283

prep_dcl(iws)
GGraph.set_fig("itr"=>3,'viewport'=>vpt00)
GGraph.set_axes('xunits'=>'','yunits'=>'','xtitle'=>'k','ytitle'=>ytitle,'xside'=>'b', 'yside'=>'l') 
DCL.uzpset('labelxb',true)
GGraph.line(check1,true,'annot'=>false, 'titl'=>' ', "max"=>max,"min"=>min,"legend"=>true,"type"=>1,"index"=>index1,"annotate"=>false,"legend_vx"=>0.3,"legend_vy"=>0.94)
GGraph.line(data_res1_ens,false,'annot'=>false,'titl'=>' ',"legend"=>true,"type"=>1,"index"=>index2,"annotate"=>false,"legend_vx"=>0.45,"legend_vy"=>0.94)
GGraph.line(data_res1,false,'annot'=>false,'titl'=>' ',"legend"=>true,"type"=>1,"index"=>index3,"annotate"=>false,"legend_vx"=>0.58,"legend_vy"=>0.94)
DCL.uzpset('pad1',-2.0) ; DCL.uxsttl('t',case1,-1) ; DCL.uzpset('pad1',0.7)
DCL.uulinz(x,y,3,3)
DCL.sgstxs(0.02)

GGraph.set_fig("itr"=>3,'viewport'=>vpt01,'new_frame'=>false)
DCL.uzpset('labelyl',false)
GGraph.line(check2,true,'annot'=>false, 'titl'=>' ', "max"=>max,"min"=>min,"legend"=>false,"type"=>1,"index"=>index1,"annotate"=>false)
GGraph.line(data_res2_ens,false,'annot'=>false,'titl'=>' ',"legend"=>false,"type"=>1,"index"=>index2,"annotate"=>false)
GGraph.line(data_res2,false,'annot'=>false,'titl'=>' ',"legend"=>false,"type"=>1,"index"=>index3,"annotate"=>false)
DCL.uzpset('pad1',-2.0) ; DCL.uxsttl('t',case2,-1) ; DCL.uzpset('pad1',0.7)
DCL.uulinz(x1,y1,3,3)
DCL.sgstxs(0.02)

GGraph.set_fig("itr"=>3,'viewport'=>vpt10,'new_frame'=>false)
DCL.uzpset('labelyl',true); DCL.uzpset('labelxb',true)
GGraph.line(check3,true,'annot'=>false, 'titl'=>' ', "max"=>max,"min"=>min,"legend"=>false,"type"=>1,"index"=>index1,"annotate"=>false)
GGraph.line(data_res3_ens,false,'annot'=>false,'titl'=>' ',"legend"=>false,"type"=>1,"index"=>index2,"annotate"=>false)
GGraph.line(data_res3,false,'annot'=>false,'titl'=>' ',"legend"=>false,"type"=>1,"index"=>index3,"annotate"=>false)
DCL.uzpset('pad1',-2.0) ; DCL.uxsttl('t',case3,-1) ; DCL.uzpset('pad1',0.7)
DCL.uulinz(x1,y1,3,3)
DCL.sgstxs(0.02)

GGraph.set_fig("itr"=>3,'viewport'=>vpt11,'new_frame'=>false)
DCL.uzpset('labelyl',false)
GGraph.line(check4,true,'annot'=>false, 'titl'=>' ', "max"=>max,"min"=>min,"legend"=>false,"type"=>1,"index"=>index1,"annotate"=>false)
GGraph.line(data_res4_ens,false,'annot'=>false,'titl'=>' ',"legend"=>false,"type"=>1,"index"=>index2,"annotate"=>false)
GGraph.line(data_res4,false,'annot'=>false,'titl'=>' ',"legend"=>false,"type"=>1,"index"=>index3,"annotate"=>false)
DCL.uzpset('pad1',-2.0) ; DCL.uxsttl('t',case4,-1) ; DCL.uzpset('pad1',0.7)
DCL.uulinz(x1,y1,3,3)
DCL.sgstxs(0.02)

DCL::sgtxzv(0.5,vpt00[3]+0.07,'check sign',
	    1.15*DCL.uzpget('rsizec2'),0,0,3)

DCL.grcls      # 描画終了処理