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

pfdir = "../prog02.5_calc_stellarspectrum/out-dwn1e2"

#sym = "LW_Case25_Earth_Tropics_CO2-300ppmv"
sym = "Test_Earth_Tropics_H2O_CO2_O3_N2O_CO_CH4_O2"
#sym = "Test_Earth_Tropics_O2"
spenc = "../prog03_calc_rte/out-dwn1e0-lw/"+sym+"_Spectral_Flux.nc"
atmnc = "../prog01.0_mkprofile/out/"+sym+".nc"
#spenc = "../prog03_calc_rte/out-dwn1e1-spe-lw23/Spectral_Flux.nc"
#atmnc = "../prog01.0_mkprofile/out/LW_Case25_Earth_Tropics_CO2-300ppmv.nc"
#
#spenc = "../prog03_calc_rte/out-dwn1e1-spe-lw33/Spectral_Flux.nc"
#atmnc = "../prog01.0_mkprofile/out/LW_Case33_Earth_SubArcWinter_CO2-300ppmv.nc"
#
#spenc = "../prog03_calc_rte/out-dwn1e1-lw/Test_Earth_Iso270K_Tropics_CO2-300ppmv_Spectral_Flux.nc"
#atmnc = "../prog01.0_mkprofile/out/Test_Earth_Iso270K_Tropics_CO2-300ppmv.nc"
#
#spenc = "../prog03_calc_rte/out-dwn1e1-lw/Test_Earth_LinInv_Tropics_CO2-300ppmv_Spectral_Flux.nc"
#atmnc = "../prog01.0_mkprofile-mod/out/Test_Earth_LinInv_Tropics_CO2-300ppmv.nc"
#
press = [1000e2, 300e2, 100e2, 30e2, 10e2, 3e2, 1e2, 0]


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

def prepare_vars( spenc, press, vname, flagpress = true )

  path = spenc
  gp = GPhys::NetCDF_IO.open( path, vname )

  if flagpress then
    gp0 = gp.cut("Press"=>press)
  else
    gp0 = gp
  end
  gp0 = gp0.copy

#  gp0 = gp0 * 1e2

#  namelev = "WaveNum"
#  wn = gp.coord(namelev)
#  gp0 = gp0 * wn

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

#  gp0.long_name = "flux"
  gp0.long_name = "optical depth"


  if flagpress then
    gpwl0 = gp.cut("Press"=>press)
  else
    gpwl0 = gp
  end
  gpwl0 = gpwl0.copy
  gpwl0 = mkwlobj( gpwl0 )

  return gp0
#  return gpwl0
end

def mkwlobj( gp )

  namelev = "WaveNum"
  wn = gp.coord(namelev)

  gp = gp.copy
  wl = 1.0/wn

#  gp = gp / wl**2
#  gp = gp * wl

#  wl = wl.convert_units( Units['micron'] )
  wl.long_name = "wavelength"
  vawl = wl.convert_units( Units['micron'] )
#  vawl = VArray.new( wl.val, { "long_name"=>'wavelength', "units"=>wl.units.to_s }, namelev )
  gp.axis(namelev).set_pos(vawl)

#  gp.long_name = "flux"
  gp.long_name = "optical depth"

  return gp
end
################################################


#DCL.gropn(-1)
DCL.gropn(-2)
#DCL.sldiv('y',1,4)
DCL.sgpset('lcntl',false)
DCL.sgpset('lfull',true)
DCL.uzfact(0.5)
DCL.sgpset('lfprop',true)
DCL.sglset('LCLIP', true)
DCL.sglset('LCNTL', true)

def draw_panel( spenc, press_target, svx1, svx2, svy1, svy2, pfdir, flagfirst, flaglast=false )

#  vname = "InFlux"
#  path = pfdir + "/" + "pf0300K.nc"
#  gppf300 = prepare_vars( path, 0.0, vname, false ) * 3.141592
#  path = pfdir + "/" + "pf0280K.nc"
#  gppf280 = prepare_vars( path, 0.0, vname, false ) * 3.141592
#  path = pfdir + "/" + "pf0260K.nc"
#  gppf260 = prepare_vars( path, 0.0, vname, false ) * 3.141592
#  path = pfdir + "/" + "pf0240K.nc"
#  gppf240 = prepare_vars( path, 0.0, vname, false ) * 3.141592
#  path = pfdir + "/" + "pf0220K.nc"
#  gppf220 = prepare_vars( path, 0.0, vname, false ) * 3.141592
#  path = pfdir + "/" + "pf0200K.nc"
#  gppf200 = prepare_vars( path, 0.0, vname, false ) * 3.141592

  wns = 10; wne = 50000
  wns = 400; wne = 1200
  wns = 400; wne = 1600

#  itr = 1
#itr = 2
#  itr = 3

  x1 = wns; x2 = wne
  # for log x axis
#  y1 = 0; y2 = 400
  # for linear x axis
  y1 = 0; y2 = 5e-3
  y1 = 0; y2 = 5e-1

  # for optical depth
  y1 = 1e-2; y2 = 1e2


#  x1 = 1.0/(x2*1e2)*1e6; x2 = 1.0/(x1*1e2)*1e6
                                          
#  GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[x1,x2,y1,y2]
#  GGraph.line( gp, true, "title"=>"", "annotate"=>false, "index"=>10 )

#  if flagfirst then
#    DCL.sgfrm
#  else
    DCL.grfig
#  end
  DCL.sgsvpt(svx1,svx2,svy1,svy2)
  DCL.sgswnd(x1,x2,y1,y2)
#  DCL.grstrn(1)
  DCL.grstrn(2)
  DCL.grstrf
  if !flagfirst then
    DCL.uzlset( 'labelxb', false )
  end
  DCL.uzlset( 'labelyl', false )
  DCL.uzlset( 'labelyr', true )
#  DCL.uxsttl( 'b', 'wavenumber (cm-1)', 0.0 )
  if flaglast then
#    DCL.ussttl( 'wavenumber', 'cm-1', 'flux', 'Wm-2' )
#    DCL.ussttl( 'wavenumber', 'cm-1', 'flux', 'Wm-2(cm-1)-1' )
    DCL.ussttl( 'wavenumber', 'cm|-1"', 'opt.depth', '1' )
#    DCL.uysttl( 'r', 'flux (Wm-2)', 0.0 )
  else
#    DCL.ussttl( 'wavenumber', 'cm-1', 'flux', '' )
    DCL.ussttl( 'wavenumber', 'cm|-1"', 'opt.depth', '' )
#    DCL.uysttl( 'r', 'flux (Wm-2)', 0.0 )
  end

#  DCL.uuslnt(3)
#  gppf = gppf300
##  DCL.usgrph(gppf.coord('WaveNum').val,gppf.val)
#  DCL.uulin(gppf.coord('WaveNum').val,gppf.val)
#  gppf = gppf280
##  DCL.usgrph(gppf.coord('WaveNum').val,gppf.val)
#  DCL.uulin(gppf.coord('WaveNum').val,gppf.val)
#  gppf = gppf260
##  DCL.usgrph(gppf.coord('WaveNum').val,gppf.val)
#  DCL.uulin(gppf.coord('WaveNum').val,gppf.val)
#  gppf = gppf240
##  DCL.usgrph(gppf.coord('WaveNum').val,gppf.val)
#  DCL.uulin(gppf.coord('WaveNum').val,gppf.val)
#  gppf = gppf220
##  DCL.usgrph(gppf.coord('WaveNum').val,gppf.val)
#  DCL.uulin(gppf.coord('WaveNum').val,gppf.val)
#  gppf = gppf200
##  DCL.usgrph(gppf.coord('WaveNum').val,gppf.val)
#  DCL.uulin(gppf.coord('WaveNum').val,gppf.val)
#
#  if flagfirst then
#    DCL.sgstxs( 0.01 )
#    DCL.sgtxu( 550, 440e-3, '300 K' )
#    DCL.sgtxu( 575, 360e-3, '280 K' )
#    DCL.sgtxu( 600, 280e-3, '260 K' )
#    DCL.sgtxu( 625, 200e-3, '240 K' )
#    DCL.sgtxu( 650, 120e-3, '220 K' )
#    DCL.sgtxu( 675,  40e-3, '200 K' )
#  end
#
#  DCL.uuslnt(1)
#  DCL.usgrph(gp.coord('WaveNum').val,gp.val)

  vname = 'RadUwFlux'
  vname = 'OptDep'
  path = spenc
  gp = GPhys::NetCDF_IO.open( path, vname )
  na_press = gp.coord("Press").val
  residual = 1.0e30
  for press1 in na_press
    if ( residual > ( press1-press_target ).abs )
      residual = ( press1-press_target ).abs
      press = press1
    end
  end
  gp = prepare_vars( spenc, press, vname )
  nrunmean = ( 1e0 / (gp.coord('WaveNum').val[1]-gp.coord('WaveNum').val[0]) ).to_i
#  gp = gp+1e-6
#  gp = (gp+1e-4).log
  gp = (-gp).exp
  gp = gp.running_mean('WaveNum',nrunmean)
#  gp = gp.exp
  gp = -gp.log+1e-6
  DCL.uulin(gp.coord('WaveNum').val,gp.val)

  DCL.usdaxs

  DCL.sglset('LCLIP', false)
  presslabel = sprintf("%6.2f", press*1e-2)
  presslabel = presslabel + ' hPa'
  DCL.sgstxc(-1)
#  DCL.sgtxu( 1250, 420e-3, presslabel )
#  DCL.sgtxu( x1+(x2-x1)*0.7, y2*(y2/y1)**0.1, presslabel )
  DCL.sgtxu( x1, y2*(y2/y1)**0.09, presslabel )
  DCL.sgstxc(0)
  DCL.sglset('LCLIP', true)

end


def draw_atmtemp( atmnc )

  path = atmnc
  vname = "Temp"
  gp = GPhys::NetCDF_IO.open( path, vname )

  namelev = "Press"
  z = gp.axis(namelev).pos.convert_units( Units['hPa'] )
  gp.axis(namelev).set_pos(z)

#  svx1 =  0.08; svx2 = 0.17
  svx1 =  0.19; svx2 = 0.30
  svy1 =  0.12; svy2 = 0.78

  x1 = 190; x2 = 310; y1 = 1020; y2 = 1

#  if flagfirst then
#    DCL.sgfrm
#  else
    DCL.grfig
#  end
  DCL.sgsvpt(svx1,svx2,svy1,svy2)
  DCL.sgswnd(x1,x2,y1,y2)
  DCL.grstrn(2)
  DCL.grstrf
  DCL.uzlset( 'labelyl', false )
  DCL.uspstx('dxl',50)
  DCL.ussttl( 'temperature', gp.units.to_s, 'pressure', gp.coord('Press').units.to_s )
#  DCL.ussttl( 'temperature', gp.units.to_s, '', '' )

  
  DCL.uuslnt(3)
  val = 300; DCL.uulin( [val,val], [y1,y2] )
  val = 280; DCL.uulin( [val,val], [y1,y2] )
  val = 260; DCL.uulin( [val,val], [y1,y2] )
  val = 240; DCL.uulin( [val,val], [y1,y2] )
  val = 220; DCL.uulin( [val,val], [y1,y2] )
  val = 200; DCL.uulin( [val,val], [y1,y2] )

  DCL.uuslnt(1)
  DCL.usgrph( gp.val, gp.coord('Press').val )

end
def draw_atmvmr( atmnc )

  path = atmnc
  vname = "VMR"
  gp = GPhys::NetCDF_IO.open( path, vname )

  namelev = "Press"
  z = gp.axis(namelev).pos.convert_units( Units['hPa'] )
  gp.axis(namelev).set_pos(z)

  svx1 =  0.07; svx2 = 0.17
  svx1 =  0.06; svx2 = 0.17
  svy1 =  0.12; svy2 = 0.78

  # x log axis
#  x1 = 1e-8; x2 = 1e-1; y1 = 1020; y2 = 1
  # x lin axis
  x1 = -8; x2 = 0; y1 = 1020; y2 = 1

#  if flagfirst then
#    DCL.sgfrm
#  else
    DCL.grfig
#  end
  DCL.sgsvpt(svx1,svx2,svy1,svy2)
  DCL.sgswnd(x1,x2,y1,y2)
  # x log axis
#  DCL.grstrn(4)
  # x lin axis
  DCL.grstrn(2)
  DCL.grstrf
  # x log axis
#  DCL.ussttl( 'mix. ratio', gp.units.to_s, 'pressure', gp.coord('Press').units.to_s )
  # x lin axis
  DCL.ussttl( 'log(mix.rat.)', gp.units.to_s, 'pressure', gp.coord('Press').units.to_s )


  DCL.sgstxs( 0.01 )
  DCL.sgstxc(-1)
  ylev = 8.0e-1 * 1.2
  # x log axis
#  x1 = 1e-8; x2 = 1e-3
  # x lin axis
  x1 = -8; x2 = -3
  #-----
#  molname = ["H2O", "CO2", "O3", "N2O", "CO", "CH4", "O2"]
  molname = ['H_2"O', 'CO_2"', 'O_3"', 'N_2"O', 'CO', 'CH_4"', 'O_2"']
  #-----
  (1..gp.coord('MolNum').val.size).reverse_each do |molnum|
    DCL.uuslnt(molnum)
    # x log axis
#    DCL.uulin( gp.cut('MolNum'=>molnum).val, gp.coord('Press').val )
    # x lin axis
    DCL.uulin( (gp+1e-30).log10.cut('MolNum'=>molnum).val, gp.coord('Press').val )
    DCL.sglset('LCLIP', false)
    ylev = ylev/1.2
    DCL.uulin( [x1,x2], [ylev,ylev] )
    # x log axis
#    DCL.sgtxu( x2*1.5, ylev, molname[molnum-1] )
    # x lin axis
    DCL.sgtxu( x2+1.5, ylev, molname[molnum-1] )
    DCL.sglset('LCLIP', true)
  end
  #-----
  DCL.sgstxc(0)
  DCL.uuslnt(1)

#  DCL.uxaxdv( 'b', 1.0, 2.0 )
#  DCL.uxaxdv( 't', 1.0, 2.0 )
#  DCL.ulylog( 'l', 3, 3 )
#  DCL.ulylog( 'r', 3, 3 )

  DCL.uspstx('dxl',2)
  DCL.usdaxs

end

draw_atmvmr( atmnc )
draw_atmtemp( atmnc )

#delsvy = 0.15
#panelheight = 0.12
#svy1 = -0.07

delsvy = 0.108
panelheight = 0.095
panelheight = 0.09
svy1 = -0.025

svx1 = 0.2; svx2 = 0.65
svx1 = 0.325; svx2 = 0.65


for i in 0..(press.size-1)
  flagfirst = false
  flaglast  = false
  if i == 0 then
    flagfirst = true
  elsif i == press.size-1 then
    flaglast = true
  end
  svy1 = svy1 + delsvy
  svy2 = svy1 + panelheight
  draw_panel( spenc, press[i], svx1, svx2, svy1, svy2, pfdir, flagfirst, flaglast )

end

DCL.grcls