require "numru/ggraph"
include NumRu

wn1 = 852e2
wn2 = 853e2
#wn1 =   10e2
#wn2 = 2000e2
#wn1 =  476e2
#wn2 =  477e2
#wn1 =  800e2
#wn2 = 1200e2

ncfn  = "data1.nc"
#ncfn  = "data1-modified.nc"
ncfn  = "data1_delod1_852.nc"
vname = "value"
gphys0 = GPhys::IO.open( ncfn, vname )
#gphys0 = gphys0.cut("wn"=>wn1..wn2)
colden0 = 2.3166947E+24 * 1.0e4
colden0 = 2.3174016E+24 * 1.0e4
colden0 = 3.8138349E+22 * 1.0e4 # H2O Only
gphys0 = gphys0 / colden0

p1 = 101300
p2 =  90200
p  = 95750
delp = p1 - p2
grav = 9.80665
mmr = 2.8758911851455217E-002
avogn = 6.022140857e23
vmr1 = 0.0187514704586114
vmr2 = 0.0137988787636849
ncfn  = "out-dwn1e0-0010-2000/LW_Case27_Earth_MidLatSummer_CO2-300ppmv_Spectral_Flux.nc"
ncfn  = "out-dwn1e0-0010-2000/LW_Case20_Earth_MidLatSummer_H2OOnly_Spectral_Flux.nc"
vname = "DelOptDep"
gphys1 = GPhys::IO.open( ncfn, vname )
gphys1 = gphys1.cut("PressLC"=>p).cut("WaveNum"=>wn1..wn2)
colden = (delp)/grav / (mmr/avogn) * (vmr1+vmr2)/2
gphys1 = gphys1 / colden

colden1 = colden
p colden0
p colden1
p 1-colden0/colden1

ncfn  = "../prog02.0_calc_ac-selected-lines-852/out-dwn1e0/LW_Case20_Earth_MidLatSummer_H2OOnly_ac.nc"
vname = "AbsCoef"
gphys1 = GPhys::IO.open( ncfn, vname )
gphys1 = gphys1.cut('MolNum'=>1)
gphys1 = gphys1.cut("WaveNum"=>wn1..wn2)
gphys1 = ( gphys1.cut("Press"=>p1) + gphys1.cut("Press"=>p2) ) / 2

#gphys0 = gphys0 * 4e27
#gphys1 = gphys1 * 4e27

#ncfn  = "out-dwn1e0-lw-lor/LW_Case20_Earth_MidLatSummer_H2OOnly_Spectral_Flux.nc"
#ncfn  = "out-dwn1e0-lw-fs85/LW_Case20_Earth_MidLatSummer_H2OOnly_Spectral_Flux.nc"
#ncfn  = "out-dwn1e0-lw-width/LW_Case20_Earth_MidLatSummer_H2OOnly_Spectral_Flux.nc"
#ncfn  = "out-dwn1e0-lw/LW_Case20_Earth_MidLatSummer_H2OOnly_Spectral_Flux.nc"
#vname = "DelOptDep"
#gphys2 = GPhys::IO.open( ncfn, vname )
#gphys2 = gphys2.cut("PressLC"=>1e5).cut("WaveNum"=>wn1..wn2)


iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL.gropn(iws)

DCL.sldiv('y',1,3)
DCL.sgpset('lclip', true)
DCL.sgpset('isub', 96)   # control character of subscription: '_' --> '`'
DCL.sgpset('lfull',true)
DCL.uzfact(0.6)
GGraph.set_fig 'itr'=> 1, 'viewport'=>[0.25,0.7,0.1,0.2], 'window'=>[wn1,wn2,0,2e-26]

gphys1.long_name = 'abs.coef.'

GGraph.line gphys1, true , 'legend'=>'my calculation'
GGraph.line gphys0, false, 'index'=>20, 'legend'=>'LBLRTM'
GGraph.line gphys1, false, 'legend'=>'my calculation'
#GGraph.line gphys1*factor, false, 'index'=>30, 'legend'=>'(my calculation)*'+factor.to_s
#GGraph.line gphys2, false, 'index'=>40, 'legend'=>'test'

wn = gphys1.coord('WaveNum').val
wn = VArray.new(wn,{"units"=>"m-1"}, "wn")
gphys0ip = gphys0.interpolate(wn)
p gphys1-gphys0ip

gphyszero = gphys1.copy
gphyszero[true] = 0.0

GGraph.set_fig 'itr'=> 1, 'viewport'=>[0.25,0.7,0.1,0.2], 'window'=>[wn1,wn2,-1e-28,1e-28]

gphyszero.long_name = 'diff.'

GGraph.line gphyszero, true, 'type'=>4
GGraph.line gphys1-gphys0ip, false, 'legend'=>'difference'
#GGraph.line gphys1*factor-gphys0ip, false, 'index'=>30, 'legend'=>'(my calculation)*'+factor.to_s
#GGraph.line gphys2-gphys0ip, false, 'index'=>40, 'legend'=>'test'

GGraph.set_fig 'itr'=> 1, 'viewport'=>[0.25,0.7,0.1,0.2], 'window'=>[wn1,wn2,-0.05,0.05]

gphyszero.long_name = 'rel.err.'
gphyszero.units = '1'

GGraph.line gphyszero, true, 'type'=>4
GGraph.line (gphys1-gphys0ip)/((gphys1+gphys0ip)/2), false, 'legend'=>'ratio'
#GGraph.line (gphys1*factor-gphys0ip)/gphys0ip, false, 'index'=>30, 'legend'=>'(my calculation)*'+factor.to_s
#GGraph.line (gphys2-gphys0ip)/gphys0ip, false, 'index'=>40, 'legend'=>'test'

DCL.grcls