#!/usr/bin/ruby

require "numru/ggraph"
include NumRu

path = "out_hitemp_dwn1e1/Venus_VIRA+GasVMR+CldMMR+UVabsMMR_H2O_CO2_CO_SO2_HF_OCS_N2_ac.nc"
path_pfs = [ "../prog02.5_calc_stellarspectrum/out/solar_flux_Gueymard2004_1366.1Wm-2.nc", \
             "../prog02.5_calc_stellarspectrum/out/pf0150K.nc",
             "../prog02.5_calc_stellarspectrum/out/pf0750K.nc" ]
pf_labels = [ "solar flux", "150 K", "750 K" ]
path_ptcl   = "../prog02.1_calc_optprop_particle/out_dwn1e1/particle_opt_prop.nc"

band_bnds = [10.0, 350.0, 500.0, 630.0, 700.0, 820.0, 980.0, 1080.0, 1180.0, 1390.0, 1480.0, 1800.0, 2080.0, 2250.0, 2380.0, 2600.0, 3250.0, 4000.0, 4650.0, 5150.0, 6150.0, 7700.0, 8050.0, 12850.0, 16000.0, 22650.0, 29000.0, 38000.0, 50000.0]
#band_bnds = NArray.to_na(band_bnds)

def chgwnunit( gp )
  va_in  = gp.coord('WaveNum')
  wn     = va_in.val
  wn     = wn * 1e-2
  va_out = VArray.new( wn, { "long_name"=>va_in.long_name,
                             "units"=>"cm-1" }, "WaveNum" )
  gp.axis('WaveNum').set_pos(va_out)
end

gpac  = GPhys::NetCDF_IO.open( path, 'AbsCoef' )
chgwnunit( gpac )
gpvmr = GPhys::NetCDF_IO.open( path, 'VMR' )


itr = 3


#print "1: Display,  2: File\n"
#citr = gets
#citr = citr.chomp!
#DCL.gropn(citr.to_i)
##DCL.gropn(1)
iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL.gropn(iws)

#DCL.sldiv('y',1,1)
DCL.sgpset('lcntl',false)
DCL.sgpset('lfull',true)
DCL.uzfact(0.75)
DCL.sgpset('lfprop',true)
DCL.sglset('lclip',true)

y1 = 1e-40
y2 = 1e-20

log10y1 = log10(y1)
log10y2 = log10(y2)

wn = gpac.coord('WaveNum').val
molnum = gpac.coord('MolNum').val

x1 = wn[0]
x2 = wn[-1]

svx1 = 0.12; svx2 = 0.93; svy1 = 0.15; svy2 = 0.65

molname = ["H2O", "CO2", "CO", "SO2", "HF", "OCS", "N2"]

GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[x1,x2,log10y1,log10y2]

plev = 100e5

gptot = gpac.cut('Press'=>plev).cut('MolNum'=>molnum[0])
gptot = gptot * gpvmr.cut('Press'=>plev).cut('MolNum'=>molnum[0]).val[0]

for i in 0..(molnum.size-1)
  if i == 0 then
    flag_first = true
  else
    flag_first = false
  end
  gpout = gpac.cut('Press'=>plev).cut('MolNum'=>molnum[i])
#  gpout = gpout * gpvmr.cut('Press'=>plev).cut('MolNum'=>molnum[i])
  gpout = gpout * gpvmr.cut('Press'=>plev).cut('MolNum'=>molnum[i]).val[0]

  gptot = gptot + gpout

  gpout = gpout + y1*0.9
  gpout = gpout.log10
  gpout.long_name = 'log(k*vmr)'
  GGraph.line( gpout, flag_first, "exchange"=>false, "index"=>(i+1)*10, "title"=>"", 'legend'=>molname[i], 'legend_vx'=>0.2, 'annotate'=>false )
end

gpout = gptot
gpout = gpout + y1*0.9
gpout = gpout.log10
gpout.long_name = 'log(k*vmr)'
GGraph.line( gpout, flag_first, "exchange"=>false, "index"=>(i+1)*10, "title"=>"", 'legend'=>"Total", 'legend_vx'=>0.2, 'annotate'=>false )



for i in 0..(path_pfs.size-1)
  path = path_pfs[i]
  #
  gppf = GPhys::NetCDF_IO.open( path, 'InFlux' )
  chgwnunit( gppf )
  #
  gppf = gppf / gppf.max.val
  gppf = gppf * ( log10y2 - log10y1 ) + log10y1
  #
  GGraph.line( gppf, flag_first, "exchange"=>false, "index"=>10, "type"=>3, "title"=>"", 'legend'=>pf_labels[i], 'legend_vx'=>0.2, 'annotate'=>false )
end

band_bnds.each{ |wn|
  DCL.uulin([wn,wn],[log10y1,log10y2])
}
p band_bnds

DCL.grcls