require "numru/ggraph"
include NumRu

wn1 =  400e2
wn2 = 2000e2

ncfn0 = "LBLRTM_Flux_Tendency_Case27_Spe.nc"
ncfn1 = "spe_run_ave/out/Case27_Spectral_Flux.nc"
ncfn0 = "LBLRTM_Flux_Tendency_Case20_Spe.nc"
ncfn1 = "spe_run_ave/out/Case20_Spectral_Flux.nc"

ncfn  = ncfn0
vname = "RadDwFlux"
gphys0 = GPhys::IO.open( ncfn, vname )
gphys0 = gphys0.cut("WaveNum"=>wn1..wn2).cut("Press"=>1e5)

ncfn  = ncfn1
vname = "RadDwFlux"
gphys1 = GPhys::IO.open( ncfn, vname )
gphys1 = gphys1.cut("Press"=>1e5).cut("WaveNum"=>wn1..wn2)
gphys1 = gphys1.copy



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.2,0.8,0.075,0.2], 'window'=>[wn1,wn2,0,5e-3]

gphys1.long_name = 'dn.flux'
GGraph.line gphys1, true, 'legend'=>'my calculation'
GGraph.line gphys0, false, 'index'=>20, 'legend'=>'LBLRTM', 'title'=>''
GGraph.line gphys1, false, 'legend'=>'my calculation'


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

gphyszero = gphys1.copy
gphyszero[true] = 0.0

gphyszero.long_name = 'diff.'

GGraph.set_fig 'itr'=> 1, 'viewport'=>[0.2,0.8,0.075,0.2], 'window'=>[wn1,wn2,-0.0001,0.0001]
GGraph.line gphyszero, true, 'type'=>4, 'title'=>''
GGraph.line gphys1-gphys0ip, false, 'legend'=>'difference'

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

GGraph.set_fig 'itr'=> 1, 'viewport'=>[0.2,0.8,0.075,0.2], 'window'=>[wn1,wn2,-0.1,0.1]
GGraph.line gphyszero, true, 'type'=>4, 'title'=>''
GGraph.line (gphys1-gphys0ip)/((gphys1+gphys0ip)/2), false, 'legend'=>'ratio'


DCL.grcls