#!/usr/bin/env ruby

#  * Drawing the value of integral for radiative heating. 

require "numru/ggraph"      # ggraph library
require "numru/dcl"      # dcl library
require "getoptlong" # A library for analysis of command line option
include NumRu

parser = GetoptLong.new
parser.set_options( ['--Time', GetoptLong::REQUIRED_ARGUMENT] )

begin
  parser.each_option do |name, arg|
    eval "$OPT_#{name.sub(/^--/, '').gsub(/-/, '_')} = '#{arg}'"
  end
    rescue
  exit(1)
end
# Refered to the web page 
# http://doc.loveruby.net/refm/api/view/library/getoptlong

T1 = $OPT_Time.to_f # time

# definition of constant value

Grav = 3.72
Psfc = 700.0 # surface pressure [Pa]
GasR = 188.913757257256
Cp   = 733.932446403617
Z1 = 0
Z2 = 1000
Z3 = 15000
Gint = 200
Xmin = 0
Xmax = 50000

#Min = 0.0
#Max = 1.0e5
Min = -2.0e-5
Max = 2.0e-5

GridNum = (Xmax - Xmin) / Gint

# file

file_1 = "MarsCond_BasicZ.nc"
file_2 = "MarsCond_PotTemp.nc"
file_3 = "MarsCond_Exner.nc"
file_4 = "MarsCond_Zprof.nc"

# GPhys object initialization

gp_thetaB = GPhys::IO.open(file_1, "PotTempBasicZ")
gp_thetaD = GPhys::IO.open(file_2, "PotTempDist")
gp_exnerB = GPhys::IO.open(file_1, "ExnerBasicZ")
gp_exnerD = GPhys::IO.open(file_3, "Exner")
gp_rad = GPhys::IO.open(file_4, "PotTempRad")

gp_rho = (Psfc / GasR) \
 * (gp_exnerB + gp_exnerD)**((Cp - GasR)/GasR) \
 / (gp_thetaB + gp_thetaD)
gp_rho = gp_rho.cut('x'=>0..50000)
gp_rho = gp_rho.mean('x') # zonal mean value of density
gp_rhoHeat = gp_rho.cut('z'=>Z1..Z2)
gp_rhoCool = gp_rho.cut('z'=>Z2..Z3)
gp_rhoHeatSum = gp_rhoHeat.sum('z') 
gp_rhoCoolSum = gp_rhoCool.sum('z')

#work1 = gp_rhoHeatSum.cut('t'=>300).val
#work2 = gp_rhoCoolSum.cut('t'=>300).val
#print(work1,"\n")
#print(work2,"\n")

gp_EB = gp_exnerB.mean('x')
EB1 = gp_EB.cut('z'=>100).val
EB2 = gp_EB.cut('z'=>1100).val
#EB3 = gp_EB.cut('z'=>10100).val
gp_radheat = gp_rad.cut('z'=>100) * EB1 / 2.0
gp_radcool = gp_rad.cut('z'=>1100) * EB2 / 2.0
#gp_radcool2 = gp_rad.cut('z'=>10100) * EB3 / 2.0

#print(EB1,"\n")
work = gp_radcool.cut('t'=>1.0).val
#work2 = gp_radcool2.cut('t'=>1.0).val
print(work,"\n")
#print(work2,"\n")

gp_Integral = \
   Cp * Gint * gp_radheat * gp_rhoHeatSum \
+  Cp * Gint * gp_radcool * gp_rhoCoolSum
#   Cp * GridNum * Gint * Gint * gp_radheat * gp_rhoHeatSum \
#+  Cp * GridNum * Gint * Gint * gp_radcool * gp_rhoCoolSum

# draw

DCL.gropn(4)
DCL.sgpset('lcntl', false)
DCL.uzfact(0.7)
DCL.sgpset('lfull', true)      # use full area in the window
DCL.sgpset('lfprop',true)      # use proportional font

GGraph::set_fig('view'=>[0.15,0.85,0.20,0.55])


GGraph.set_axes('xunits'=>'sec','yunits'=>'W/m^2','xtitle'=>'Time','ytitle'=>'Radiative Heating')

GGraph.line( gp_Integral, true, 'title'=>'Radiative Heating', 'max'=>Max, 'min'=>Min, 'exchange'=>false)
#GGraph.contour( gp_SqW, false, 'max'=>Max, 'min'=>Min, 'interval'=>Cint )
#GGraph.color_bar( 'vlength'=>0.5,'landscape'=>true, 'tickintv'=>1)

DCL.grcls