#!/usr/bin/ruby
# coding: utf-8

require "narray"
require './solar_system_param'
require './planetorbit'
require './insolation_fig'


NOrbitDiv = 181
JMax      =  64

if ARGV.size == 0 then

  iws           = 1
  FigFN         = "dcl"

  print "Select mode:\n"
  print "  1 : Seaonal variation of insolation\n"
  print "  2 : Annual mean\n"
  print "> "
  sMode = gets.to_i

  iPlanet = 3-1
  sSMAxis       = A_SMAxis[iPlanet]
  sEccen        = A_Eccen[iPlanet]
  sPerLon       = A_PerLon[iPlanet]
  #sSecInAYear = tryear[iPlanet] * 365.0 * ( 60.0 * 60.0 * 24.0 )
  sSecInAYear   = A_TrYear[iPlanet] * 365.242191 * ( 60.0 * 60.0 * 24.0 )
  sEpsOrb       = A_Obliq[iPlanet]
  sPerLs        = A_PerLs[iPlanet]
  sMLonAtEpoch  = A_MLonEp[iPlanet]

else

  iws          = 2
  FigFN        = ARGV[0]
  sMode        = ARGV[1].to_i
  sSMAxis      = ARGV[2].to_f
  sEccen       = ARGV[3].to_f
  sPerLon      = ARGV[4].to_f
  sSecInAYear  = ARGV[5].to_f
  sEpsOrb      = ARGV[6].to_f
  sPerLs       = ARGV[7].to_f
  sMLonAtEpoch = ARGV[8].to_f

end

y_Index = NArray.float(JMax     ).indgen!(1.0,1.0)
y_Lat   = -90.0 + 180.0/JMax/2 + 180.0/JMax*(y_Index-1)

a_Time, ay_Ins = ty_insolation( NOrbitDiv, y_Lat, sSMAxis, sEccen, sPerLon, sSecInAYear, sEpsOrb, sPerLs, sMLonAtEpoch )

draw_preparation( iws, FigFN )

case sMode
when 1 then
  ay_InsWithNeg = NArray.float(NOrbitDiv,JMax)
  ay_InsWithNeg[true,true] = ay_Ins[true,true]
  ay_InsWithNeg[ ay_Ins.le 0.0 ] = -1.0e-10

  draw_insolation( a_Time, y_Lat, ay_InsWithNeg )
  #draw_insolation( a_Ls, y_Lat, ay_Ins )
when 2 then
  y_CosLat = NMath.cos( y_Lat * PI/180.0 )
  y_Ins = ay_Ins.mean(0)
  y_Ins /= ( y_Ins * y_CosLat * PI / JMax ).sum
  draw_insolation_1D( y_Lat, y_Ins )
end

draw_finish