#!/usr/bin/ruby
require "numru/ggraph"
include NumRu

dir_ac   = "../prog02_calc_ac/out-dwn1e0"
fn_ac    = "LW_Case25_Earth_Tropics_CO2-300ppmv_ac.nc"
dir_sort = "../prog04_sort_ac_check/out-demo"
fn_sort  = "ac_reordered-04-08.nc"
fn_sort_noncor = "ac_reordered-04-11.nc"

press1 = 1e4
press2 = 1e5
molnum = 2
wns =  50e2; wne = 200e2
wns = 500e2; wne = 630e2
wns = 630; wne = 700


################################################

path = dir_ac + "/" + fn_ac
vname = 'AbsCoef'
gp = GPhys::NetCDF_IO.open( path, vname )

namelev = "WaveNum"
z = gp.axis(namelev).pos.convert_units( Units['cm-1'] )
z.long_name = 'wavenumber'
gp.axis(namelev).set_pos(z)

gp1 = gp.cut("Press"=>press1).cut("MolNum"=>molnum).cut("WaveNum"=>wns..wne)
gp2 = gp.cut("Press"=>press2).cut("MolNum"=>molnum).cut("WaveNum"=>wns..wne)

gp1 = gp1.copy
gp2 = gp2.copy

unit_org = gp1.units.to_s
gp1 = gp1.log10
gp2 = gp2.log10
gp1.units = "log("+unit_org+")"
gp2.units = "log("+unit_org+")"

#gp.long_name = "absorption coeffieicnt"
gp1.long_name = "log k"
gp2.long_name = "log k"


################################################

path = dir_sort + "/" + fn_sort
vname = 'AbsCoef'
gp_sort = GPhys::NetCDF_IO.open( path, vname )
#
path = dir_sort + "/" + fn_sort
vname = 'VMR'
gp_vmr = GPhys::NetCDF_IO.open( path, vname )


namelev = "WaveNum"
z = gp_sort.axis(namelev).pos.convert_units( Units['cm-1'] )
z.long_name = 'wavenumber'
gp_sort.axis(namelev).set_pos(z)

gp_sort1 = gp_sort.cut("Press"=>press1).cut("MolNum"=>molnum).cut("WaveNum"=>wns..wne)
gp_sort2 = gp_sort.cut("Press"=>press2).cut("MolNum"=>molnum).cut("WaveNum"=>wns..wne)

gp_sort1 = gp_sort1.copy
gp_sort2 = gp_sort2.copy

gp_sort1 = gp_sort1 / gp_vmr.cut("Press"=>press1).cut("MolNum"=>molnum)
gp_sort2 = gp_sort2 / gp_vmr.cut("Press"=>press2).cut("MolNum"=>molnum)

unit_org = gp_sort1.units.to_s
gp_sort1 = gp_sort1.log10
gp_sort2 = gp_sort2.log10
#gp_sort1 = (gp_sort1+1e-27).log10
#gp_sort2 = (gp_sort2+1e-27).log10
gp_sort1.units = "log("+unit_org+")"
gp_sort2.units = "log("+unit_org+")"

#gp.long_name = "absorption coeffieicnt"
gp_sort1.long_name = "log k"
gp_sort2.long_name = "log k"

wn = gp_sort1.coord('WaveNum').val
g = (wn-wns)/(wne-wns)
vag = VArray.new( g, {"long_name"=>"g", "units"=>"1"}, "WaveNum" )
gp_sort1.axis('WaveNum').set_pos(vag)
gp_sort2.axis('WaveNum').set_pos(vag)

################################################

path = dir_sort + "/" + fn_sort_noncor
vname = 'AbsCoef'
gp_sort_noncor = GPhys::NetCDF_IO.open( path, vname )
#
path = dir_sort + "/" + fn_sort_noncor
vname = 'VMR'
gp_vmr = GPhys::NetCDF_IO.open( path, vname )


namelev = "WaveNum"
z = gp_sort_noncor.axis(namelev).pos.convert_units( Units['cm-1'] )
z.long_name = 'wavenumber'
gp_sort_noncor.axis(namelev).set_pos(z)

gp_sort_noncor1 = gp_sort_noncor.cut("Press"=>press1).cut("MolNum"=>molnum).cut("WaveNum"=>wns..wne)
gp_sort_noncor2 = gp_sort_noncor.cut("Press"=>press2).cut("MolNum"=>molnum).cut("WaveNum"=>wns..wne)

gp_sort_noncor1 = gp_sort_noncor1.copy
gp_sort_noncor2 = gp_sort_noncor2.copy

gp_sort_noncor1 = gp_sort_noncor1 / gp_vmr.cut("Press"=>press1).cut("MolNum"=>molnum)
gp_sort_noncor2 = gp_sort_noncor2 / gp_vmr.cut("Press"=>press2).cut("MolNum"=>molnum)

unit_org = gp_sort_noncor1.units.to_s
gp_sort_noncor1 = gp_sort_noncor1.log10
gp_sort_noncor2 = gp_sort_noncor2.log10
#gp_sort1 = (gp_sort1+1e-27).log10
#gp_sort2 = (gp_sort2+1e-27).log10
gp_sort_noncor1.units = "log("+unit_org+")"
gp_sort_noncor2.units = "log("+unit_org+")"

#gp.long_name = "absorption coeffieicnt"
gp_sort_noncor1.long_name = "log k"
gp_sort_noncor2.long_name = "log k"

wn = gp_sort_noncor1.coord('WaveNum').val
g = (wn-wns)/(wne-wns)
vag = VArray.new( g, {"long_name"=>"g", "units"=>"1"}, "WaveNum" )
gp_sort_noncor1.axis('WaveNum').set_pos(vag)
gp_sort_noncor2.axis('WaveNum').set_pos(vag)
################################################


DCL.gropn(2)
DCL.sldiv('y',2,2)
DCL.sgpset('lcntl',false)
DCL.sgpset('lfull',true)
#DCL.uzfact(1.5)
DCL.sgpset('lfprop',true)

svx1 = 0.15; svx2 = 0.9; svy1 = 0.2; svy2 = 0.6
x1 = wns; x2 = wne; y1 = -1.0e-24; y2 = 8.01e-22
x1 = wns; x2 = wne; y1 = 0; y2 = 8e-22
x1 = wns; x2 = wne; y1 = 8e-25; y2 = 10e-22
x1 = wns; x2 = wne; y1 = 8e-25; y2 = 10e-22
x1 = wns; x2 = wne; y1 = -26; y2 = -20


itr = 1
#itr = 2

GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[x1,x2,y1,y2]
GGraph.line( gp1, true, "title"=>"", "annotate"=>false, "index"=>10 )

GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[0.0,1.0,y1,y2]
GGraph.line( gp_sort_noncor1, true, "title"=>"", "annotate"=>false, "index"=>310, "type"=>1 )
GGraph.line( gp_sort1, false, "title"=>"", "annotate"=>false, "index"=>10 )

GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[x1,x2,y1,y2]
GGraph.line( gp2, true, "title"=>"", "annotate"=>false, "index"=>10 )

GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[0.0,1.0,y1,y2]
GGraph.line( gp_sort2, true, "title"=>"", "annotate"=>false, "index"=>10 )

DCL.grcls