#!/usr/bin/env ruby
# coding: utf-8
#
# 表題  ダイナモベンチマーク Case0 用解析スクリプト
#
# 履歴  2024/04/24  竹広 真一
#
require "numru/gphys"
include NumRu

Files = ARGV[0]
File0  = Files.sub(/-ip[0-9*]+/,"-ip0000")

gtURL=File0+'@ek'
ek = GPhys::IO.open_gturl(gtURL)

gtURL=Files+'@temp'
temp = GPhys::IO.open_gturl(gtURL)

gtURL=Files+'@vlon'
vlon = GPhys::IO.open_gturl(gtURL)

gtURL=Files+'@vrad'
vrad = GPhys::IO.open_gturl(gtURL)

lon  = temp.coord(0)
lat  = temp.coord(1)
rad  = temp.coord(2)
time = temp.coord(3)

Ro = rad.max.val ; Ri = rad.min.val
Tcut = time.val[-2..-1]

LatCut = 0.0 ; RadCut = (Ro+Ri)/2

ek_val=Array[]
vlon_val=Array[]
temp_val=Array[]
lon_val=Array[]

for it in 0..1 do
  ek_val.push ek.cut('t'=>Tcut[it]).val[0]
  temp_na = temp.cut('lat'=>LatCut,'rad'=>RadCut,'t'=>Tcut[it]).val
  vlon_na = vlon.cut('lat'=>LatCut,'rad'=>RadCut,'t'=>Tcut[it]).val
  vrad_na = vrad.cut('lat'=>LatCut,'rad'=>RadCut,'t'=>Tcut[it]).val
  lon_na  = lon.val ; nlon = lon_na.size 

  for i in 0..nlon-2 do
    if vrad_na[i]*vrad_na[i+1] <=0.0 && vlon_na[i+1]-vlon_na[i] > 0.0 
      vlon_val.push (vlon_na[i]*vrad_na[i+1]-vlon_na[i+1]*vrad_na[i]) \
                    /(vrad_na[i+1]-vrad_na[i])
      temp_val.push (temp_na[i]*vrad_na[i+1]-temp_na[i+1]*vrad_na[i]) \
                    /(vrad_na[i+1]-vrad_na[i])
      lon_val.push (lon_na[i]*vrad_na[i+1]-lon_na[i+1]*vrad_na[i]) \
                   /(vrad_na[i+1]-vrad_na[i])
      break
    end
   end

   printf(" === time = %f, index = %d, lon = %f, ek = %e, temp = %e, vlon = %e \n", Tcut[it], i, lon_val[it], ek_val[it], temp_val[it], vlon_val[it])
end  

ekp   = (ek_val[0]+ek_val[1])/2
vlonp = (vlon_val[0]+vlon_val[1])/2
tempp = (temp_val[0]+temp_val[1])/2      
omega = (lon_val[1]-lon_val[0])/(Tcut[1]-Tcut[0])/180.0*PI

printf(" ### Result : ek = %e, temp = %e, vlon = %e, omega = %e\n", ekp, tempp, vlonp, omega)