#!/usr/bin/ruby

require "numru/ggraph"
include NumRu

# flag_show_all
#   true  : show all spectral range in 1 panel with logarithm x axis
#   false : show all spectral range in several panel with linear x axis
flag_show_all = true
#flag_show_all = false

# plev : pressure level
plev      = 1000e2
#plev      =  500e2
#plev      =  100e2
#plev      =   10e2

# path_ac   : path to netcdf file of absorption coefficient 
# path_ptcl : path to netcdf file of particle optical property
# path_pf   : (array of) path to netcdf file of spectrum of Planck function
#           : and incident steller/solar radiation
# label_pf  : (array of) variable name for spectrum specified in path_pf
path_ac   = "../../prog02.1_calc_ac/out/Earth_ICRCCM_LW_Case27_MLS_CO2-300ppmv_ac.nc"
path_ptcl = "../../prog02.2_calc_optprop_particle/out/Earth_particle_opt_prop.nc"
#path_pf   = []
#label_pf  = []
path_pf   = "../../prog02.3_calc_stellarspectrum/out/solar_flux_Gueymard2004_1366.1Wm-2.nc"
label_pf  = "InFlux"

# optimized 23 bands
band_bnds = [10.0, 100.0, 300.0, 400.0, 500.0, 600.0, 700.0, 1000.0, 1100.0, 1400.0, 2000.0, 2400.0, 3500.0, 4500.0, 5500.0, 6500.0, 12500.0, 14500.0, 18500.0, 25500.0, 30500.0, 32500.0, 34500.0, 50000.0]


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

# volume mixing ratio
gpvmr       = GPhys::NetCDF_IO.open( path_ac, 'VMR' )
# absorption coefficient
gpac        = GPhys::NetCDF_IO.open( path_ac, 'AbsCoef' )
chgwnunit( gpac )
# Rayleigh scattering coefficient
gprsc       = GPhys::NetCDF_IO.open( path_ptcl, 'RayScatCoef' )
chgwnunit( gprsc )
# Rayleigh scattering coefficient for non radiatively active specie
gprscnonact = GPhys::NetCDF_IO.open( path_ptcl, 'RayScatCoefNonRadAct' )
chgwnunit( gprscnonact )


#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-35
y2 = 1e-20

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

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


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

x1  = wn[0]
x2  = wn[-1]
if ( flag_show_all ) then
  wnrange = x2 - x1
  nloop   = 1
  itr     = 3
else
  wnrange = 10000.0
  nloop   = ( ( x2 - x1 ) / wnrange ).to_i + 1
  itr     = 1
end

for iloop in 0..(nloop-1)

  x1  = ( wn[0] + wnrange * iloop ) * 0.95
  x2  = ( x1 + wnrange            ) * 1.05

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

  flag_first = true

  # Absorption coefficient
  for i in 0..(molnum.size-1)
    gpout = gpac.cut('Press'=>plev).cut('MolNum'=>molnum[i])
    gpout = gpout * gpvmr.cut('Press'=>plev).cut('MolNum'=>molnum[i]).val[0]

    if i == 0 then
      gptot = gpout
    else
      gptot += gpout
    end

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

  # Total
  gpout = gptot
  gpout = gpout + y1*0.99
  gpout = gpout.log10
  gpout.long_name = 'log10(k*vmr)'
#  GGraph.line( gpout, flag_first, "exchange"=>false, "index"=>(0+1)*10, "title"=>"", 'legend'=>"Total", 'legend_vx'=>0.15, 'annotate'=>false )


  # Rayleigh scattering coefficient (total)
  vmrnonact = 1.0
  for i in 0..(molnum.size-1)
    vmrnonact -= gpvmr.cut('Press'=>plev).cut('MolNum'=>molnum[i]).val[0]
  end
  gptot = gprscnonact * vmrnonact
  for i in 0..(molnum.size-1)
    gpout = gprsc.cut('MolNum'=>molnum[i])
    gpout = gpout * gpvmr.cut('Press'=>plev).cut('MolNum'=>molnum[i]).val[0]
    gptot += gpout
  end
  gpout = gptot
  gpout = gpout + y1*0.99
  gpout = gpout.log10
  gpout.long_name = 'log10(k*vmr)'
  GGraph.line( gpout, flag_first, "exchange"=>false, "index"=>(0+1)*10, 'type'=>2, 
               "title"=>"",
               'legend'=>"Rayleigh scattering coefficient", 'legend_vx'=>0.15,
               'annotate'=>false )


  # Planck function
  if ( path_pf.instance_of?(Array) ) then
    a_path_pf = path_pf
  else
    a_path_pf = [path_pf]
  end
  if ( label_pf.instance_of?(Array) ) then
    a_label_pf = label_pf
  else
    a_label_pf = [label_pf]
  end
  for i in 0..(a_path_pf.size-1)
    path  = a_path_pf[i]
    label = a_label_pf[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'=>label, 'legend_vx'=>0.15, 'annotate'=>false )
  end

  band_bnds.each{ |wn|
    DCL.uulin([wn,wn],[log10y1,log10y2])
#    DCL.uulin([wn,wn],[log10y2-2,log10y2-1])
  }
  print band_bnds.size-1, " bands", "\n"
  print band_bnds, "\n"

end

DCL.grcls