#!usr/bin/env ruby
#
# 軌道要素計算
# Berger (1978) の方法と係数を用いる
# 係数は外部ファイルから読み込む
# 
# 2013/09/19  Satoshi NODA   作成
# 

require 'narray'
include Math
include NMath

DATAFILE = './CoefsForOrbElem_Earth_B1978.dat'

# 計算する年の始点, 終点, 間隔. 西暦で指定.
t0 =  -1000000.0
t1 =    100000.0
interval = 100.0

ny = (1.0 + (t1-t0)/interval).to_i

if t0 > t1 then
  raise "ERROR: t0 > t1 is invalid."
end

if not File.exist?(DATAFILE) then
  raise "ERROR: file #{DATAFILE} does not exist."
end

### データ読み込み
file = open(DATAFILE)

# 項の数の設定
tmp = file.gets.split(' ').map{|s| s.to_i}
N_Obl  = tmp[0]
N_EPi  = tmp[1]
N_GPre = tmp[2]

obl_ampl  = NArray.float(N_Obl)  # Amplitude
obl_mrate = NArray.float(N_Obl)  # Mean rate
obl_phase = NArray.float(N_Obl)  # Phase
epi_ampl  = NArray.float(N_EPi)  # Amplitude
epi_mrate = NArray.float(N_EPi)  # Mean rate
epi_phase = NArray.float(N_EPi)  # Phase
gpre_ampl  = NArray.float(N_GPre)  # Amplitude
gpre_mrate = NArray.float(N_GPre)  # Mean rate
gpre_phase = NArray.float(N_GPre)  # Phase
obl0        = NArray.float(1)
gpre_mrate0 = NArray.float(1)
gpre0       = NArray.float(1)

[*0..N_Obl-1].each{|i|
  tmp = file.gets.split(' ').map{|s| s.to_f}
  obl_ampl[i] = tmp[0]
  obl_mrate[i] = tmp[1]
  obl_phase[i] = tmp[2]
}
[*0..N_EPi-1].each{|i|
  tmp = file.gets.split(' ').map{|s| s.to_f}
  epi_ampl[i] = tmp[0]
  epi_mrate[i] = tmp[1]
  epi_phase[i] = tmp[2]
}
[*0..N_GPre-1].each{|i|
  tmp = file.gets.split(' ').map{|s| s.to_f}
  gpre_ampl[i] = tmp[0]
  gpre_mrate[i] = tmp[1]
  gpre_phase[i] = tmp[2]
}
# 定数
tmp = file.gets.split(' ').map{|s| s.to_f}
obl0[0]        = tmp[0]
gpre_mrate0[0] = tmp[1]
gpre0[0]       = tmp[2]

file.close

# 角度の単位を radian に直す
obl_ampl   *= PI / (180.0 * 3600.0)
gpre_ampl  *= PI / (180.0 * 3600.0)

obl_mrate  *= PI / (180.0 * 3600.0)
epi_mrate  *= PI / (180.0 * 3600.0)
gpre_mrate *= PI / (180.0 * 3600.0)

obl_phase  *= PI / 180.0
epi_phase  *= PI / 180.0
gpre_phase *= PI / 180.0

obl0        *= PI / 180.0
gpre_mrate0 *= PI / (180.0 * 3600.0)
gpre0       *= PI / 180.0


### 計算
yearAD   = NArray.float(ny).indgen(t0,interval)
year   = yearAD - 1950.0  # 計算に使用する年. 西暦 1950 年起点.
obl    = NArray.float(ny).fill(0.0)
episin = NArray.float(ny).fill(0.0)
epicos = NArray.float(ny).fill(0.0)
gpre   = NArray.float(ny).fill(0.0)

# 桁落ちを防ぐため, 振幅の小さい項から足す
[*0..N_Obl-1].reverse.each{|i|
  obl += obl_ampl[i] * cos(obl_mrate[i]*year + obl_phase[i])
}
obl[true] = obl0 + obl[true]

[*0..N_EPi-1].reverse.each{|i|
  episin += epi_ampl[i] * sin(epi_mrate[i]*year + epi_phase[i])
  epicos += epi_ampl[i] * cos(epi_mrate[i]*year + epi_phase[i])
}

[*0..N_GPre-1].reverse.each{|i|
  gpre += gpre_ampl[i] * sin(gpre_mrate[i]*year + gpre_phase[i])
}
gpre += gpre_mrate0 * year
gpre[true] = gpre0 + gpre[true]

# Eccentricity
ecc = sqrt(episin**2 + epicos**2)

# Longitude of the Perihelion
fixedEquinox = acos(epicos / ecc)
mask = asin(episin) .lt 0.0
fixedEquinox[mask] = - fixedEquinox[mask]

lonOfPeri = (fixedEquinox + gpre) % (2.0*PI)
mask = lonOfPeri .lt 0.0
lonOfPeri[mask] = 2.0*PI + lonOfPeri[mask]

### output
# radian は degree に変換
w = 180.0 / PI
[*0..ny-1].each{|i|
  puts [yearAD[i], w*obl[i], ecc[i], w*lonOfPeri[i]].map{|f| f.to_s}.join(' ')
}

__END__