require "numru/ggraph"
include NumRu

cuts_hash = { 'lat'=>0, 'lon'=>-65..65 }
energy_analysis_fileNames=["energy_analysis.nc"]#, "energy_analysis2.nc"]

def extract_energyVars(energy_analysis_fileNames)
  totalEVariations = []
  kineticEVariations = []
  tpotentialEVariations = []

  for id in 0..energy_analysis_fileNames.size-1
    p "file=#{energy_analysis_fileNames[id]}"
    totalE = GPhys::NetCDF_IO.open(energy_analysis_fileNames[id], 'TotalE')
    kineticE = GPhys::NetCDF_IO.open(energy_analysis_fileNames[id], 'KineticE')
    internalE = GPhys::NetCDF_IO.open(energy_analysis_fileNames[id], 'InternalE')
    potentialE = GPhys::NetCDF_IO.open(energy_analysis_fileNames[id], 'PotentialE')

    tSize = potentialE.data.shape_current[0]
    tpotentialE = GPhys.new(potentialE.grid_copy, VArray.new(NArray.sfloat(tSize), {"units"=>"J/m^3"}, 'TotalPotentialE'))
    tpotentialE = internalE + potentialE


    totalE0 = totalE.cut('time'=>0)
    kineticE0 = kineticE.cut('time'=>0) 
    internalE0 = internalE.cut('time'=>0)
    potentialE0 = potentialE.cut('time'=>0)
    tpotentialE0 = tpotentialE.cut('time'=>0)

    totalEVariations[id] = totalE - totalE0
    kineticEVariations[id] = kineticE - kineticE0
    tpotentialEVariation = tpotentialE.copy
    tpotentialEVariation[true] = tpotentialE[true] - tpotentialE0.val
    tpotentialEVariations[id] = tpotentialEVariation
  end

  return totalEVariations, kineticEVariations, tpotentialEVariations
end 

def draw_graph(totalEVariation, kineticEVariation, tpotentialEVariation, overplotNum)
  DCL.gropn(1)
  DCL.sgpset('lcntl', false) ; DCL.uzfact(0.7)

  max = kineticEVariation[0].val.max
  min = tpotentialEVariation[0].val.min

  p "max=#{max}, min=#{min}"

  p kineticEVariation[0].val
  p tpotentialEVariation[0].val

  title="Kinetic, total and total potential energy"
  for id in 0..overplotNum-1
   lineIndex = id+1

   if id == 0
    GGraph.line( kineticEVariation[id], true, "min"=>min, "max"=>max, "title"=>title, "type"=>lineIndex )
   else
    GGraph.line( kineticEVariation[id], false, "min"=>min, "max"=>max, "type"=>lineIndex  )
   end
   
   GGraph.line( tpotentialEVariation[id], false,  "type"=>lineIndex )
   GGraph.line( totalEVariation[id], false,  "type"=>lineIndex  )
  end 

  DCL.grcls

end

#
#

teTVariations, keTVariations, tpeTVariations = extract_energyVars(energy_analysis_fileNames)

draw_graph(teTVariations, keTVariations, tpeTVariations, energy_analysis_fileNames.size)