Module: When::Ephemeris

Includes:
Math
Included in:
Coordinates::Spatial, CelestialObject, Coords, Formula, Hindu::Graha, Moon, Sun
Defined in:
lib/when_exe/ephemeris.rb,
lib/when_exe/ephemeris/sun.rb,
lib/when_exe/region/indian.rb,
lib/when_exe/ephemeris/moon.rb,
lib/when_exe/region/chinese.rb,
lib/when_exe/ephemeris/planets.rb,
lib/when_exe/ephemeris/term_table.rb,
lib/when_exe/ephemeris/phase_table.rb

Overview

Copyright (C) 2016 Takashi SUGA

You may use and/or modify this file according to the license described in the LICENSE.txt file included in this archive.

Defined Under Namespace

Modules: V50 Classes: CelestialObject, ChineseTrueLunation, Coords, Datum, Earth, Formula, FormulaWithTable, Hindu, JGD2000, LunarFormulaWithTable, MeanLunation, Moon, Shadow, SolarFormulaWithTable, Star, Sun

Constant Summary collapse

DAY =

日 / 秒

86400.0
BCENT =

ベッセル世紀

36524.2194
JCENT =

ユリウス世紀

36525.0
JYEAR =

ユリウス年

JCENT/100
EPOCH1800 =

1800 I 0.5

2378496.0
EPOCH1900 =

1900 I 1.5

2415021.0
EPOCH1975 =

1975 I 0.0

2442412.5
EPOCH2000 =

2000 I 1.5

2451545.0
DEG =

度 / radian

PI / 180
CIRCLE =

全円周 / radian

2 * PI
AU =

天文単位距離 / km

1.49597870E8
C0 =

光速度 / (km/ユリウス世紀)

299792.458*86400*JCENT
PSEC =

パーセク / AU

360.0*60.0*60.0/(2*PI)
FARAWAY =

遠距離物体のための仮定義

10000.0 * PSEC
LIN =

演算の番号

-1
SIN = 
0
COS = 
1
SINT = 
2
COST = 
3
SINL = 
4
COSL = 
5
SINLT = 
6
COSLT = 
7
AcS = 
8
Mercury =

水星

Datum::Far.new(   2440.0,  0.00347  ,  1.16, {:phi=>P1L, :theta=>P1B,  :radius=>P1R})
Venus =

金星

Datum::Near.new(  5988.0,  0.00484  , -4.00, {:phi=>P2L, :dl=>P2dL, :theta=>P2B, :radius=>P2Q})
Mars =

火星

Datum::Near.new(  3397.0,  0.00700  , -1.30, {:phi=>P4L, :dl=>P4dL, :theta=>P4B, :radius=>P4Q})
Jupiter =

木星

Datum::Big.new(  71398.0,  0.01298  , -8.93, {:phi=>P5L, :nn=>P5dL, :theta=>P5B, :radius=>P5Q,
:jsn=>P5n, :jsl=>P5l, :jst=>P5t, :jsr=>P5r})
Saturn =

土星

Datum::Big.new(  60330.0,  0.01756  , -8.68, {:phi=>P6L, :nn=>P6dL, :theta=>P6B, :radius=>P6Q,
:jsn=>P6n, :jsl=>P6l, :jst=>P6t, :jsr=>P6r})
Uranus =

天王星

Datum::Far.new(  25400.0,  0.02490  , -6.85, 2433283, 2473460, {:phi=>P7L, :theta=>P7B, :radius=>P7R})
Neptune =

海王星

Datum::Far.new(  24300.0,  0.03121  , -7.05, 2433283, 2473460, {:phi=>P8L, :theta=>P8B, :radius=>P8R})
Pluto =

冥王星

Datum::Far.new(   1180.0,  0.03461  , -1.00, 2433283, 2473460, {:phi=>P9L, :theta=>P9B, :radius=>P9R})

Class Method Summary collapse

Class Method Details

._adjust(y1, y0, delta) ⇒ Object

y が周期量である場合の周期の補正



350
351
352
 
# File 'lib/when_exe/ephemeris.rb', line 350

def _adjust(y1, y0, delta)
  (-2..+2).to_a.map {|i| y1 + i * delta}.sort_by {|y| (y-y0).abs}[0]
end

._rot(x, y, t) ⇒ Array<Numeric>

回転(2次元)

Parameters:

Returns:



220
221
222
223
 
# File 'lib/when_exe/ephemeris.rb', line 220

def _rot(x, y, t)
  c, s = cosc(t), sinc(t)
  return [x*c - y*s, x*s + y*c]
end

._to_p2(x, y) ⇒ Array<Numeric>

直交座標→極座標(2次元)

Parameters:

Returns:

  • (Array<Numeric>) —

    ( phi, radius )

    phi - 経度 / CIRCLE
    radius - 距離 


204
205
206
207
 
# File 'lib/when_exe/ephemeris.rb', line 204

def _to_p2(x, y)
  return [0.0, 0.0] if x==0 && y==0
  return [atan2(y,x)/CIRCLE, sqrt(x*x+y*y)]
end

._to_p3(x, y, z) ⇒ Array<Numeric>

直交座標→極座標(3次元)

Parameters:

Returns:

  • (Array<Numeric>) —

    ( phi, theta, radius )

    phi - 経度 / CIRCLE
    theta - 緯度 / CIRCLE
    radius - 距離 


188
189
190
191
192
 
# File 'lib/when_exe/ephemeris.rb', line 188

def _to_p3(x, y, z)
  phi, radius = _to_p2(x,  y)
  theta, radius = _to_p2(radius, z)
  return [phi, theta, radius]
end

._to_r3(phi, theta, radius) ⇒ Array<Numeric>

極座標→直交座標(3次元)

Parameters:

Returns:

  • (Array<Numeric>) —

    ( x, y, z )

    x - x 座標
    y - y 座標
    z - z 座標 


171
172
173
174
 
# File 'lib/when_exe/ephemeris.rb', line 171

def _to_r3(phi, theta, radius)
  c, s = cosc(theta), sinc(theta)
  return [radius*c*cosc(phi), radius*c*sinc(phi), radius*s]
end

.acos(x) ⇒ Numeric

arc cos / radian

Parameters:

Returns:



113
114
115
 
# File 'lib/when_exe/ephemeris.rb', line 113

def acos(x)
  atan2(sqrt(1-x*x), x)
end

.asin(x) ⇒ Numeric

arc sin / radian

Parameters:

Returns:



106
107
108
 
# File 'lib/when_exe/ephemeris.rb', line 106

def asin(x)
  atan2(x, sqrt(1-x*x))
end

.cosc(x) ⇒ Numeric

円周単位のcos

Parameters:

Returns:



141
142
143
 
# File 'lib/when_exe/ephemeris.rb', line 141

def cosc(x)
  cos(x * CIRCLE)
end

.cosd(x) ⇒ Numeric

度のcos

Parameters:

Returns:



120
121
122
 
# File 'lib/when_exe/ephemeris.rb', line 120

def cosd(x)
  cos(x * DEG)
end

.delta_e(c) ⇒ Numeric Also known as: deltaE

Δε

Parameters:

  • c (Numeric) —

    2000年からの経過世紀

Returns:

  • (Numeric) —

    緯度の章動 / CIRCLE



251
252
253
254
255
 
# File 'lib/when_exe/ephemeris.rb', line 251

def delta_e(c)
  trigonometric(c, 
    [[COS   ,  125.04      ,   -1934.136    , +0.00256  ],
     [COS   ,  200.93      ,  +72001.539    , +0.00016  ]]) / 360
end

.delta_p(c) ⇒ Numeric Also known as: deltaP

Δφ

Parameters:

  • c (Numeric) —

    2000年からの経過世紀

Returns:

  • (Numeric) —

    経度の章動 / CIRCLE



264
265
266
267
268
 
# File 'lib/when_exe/ephemeris.rb', line 264

def delta_p(c)
  trigonometric(c, 
    [[SIN   ,  125.04      ,   -1934.136    , -0.00478  ],
     [SIN   ,  200.93      ,  +72001.539    , +0.00037  ]]) / 360
end

.julian_century_from_2000(tt) ⇒ Numeric

時間の単位の換算 - 2000年元期の dynamical_time / ユリウス世紀

Parameters:

  • tt (Numeric) —

    ユリウス日(Terrestrial Time)

Returns:



241
242
243
 
# File 'lib/when_exe/ephemeris.rb', line 241

def julian_century_from_2000(tt)
  return (tt - EPOCH2000) / JCENT
end

.julian_year_from_1975(tt) ⇒ Numeric

時間の単位の換算 - 1975年元期の dynamical_time / ユリウス年

Parameters:

  • tt (Numeric) —

    ユリウス日(Terrestrial Time)

Returns:



231
232
233
 
# File 'lib/when_exe/ephemeris.rb', line 231

def julian_year_from_1975(tt)
  return (tt - EPOCH1975) / JYEAR
end

.obl(c) ⇒ Numeric

黄道傾角

Parameters:

  • c (Numeric) —

    2000年からの経過世紀

Returns:

  • (Numeric) —

    黄道傾角 / CIRCLE



277
278
279
 
# File 'lib/when_exe/ephemeris.rb', line 277

def obl(c)
  return (23.43929 + -0.013004*c + 1.0*delta_e(c)) / 360
end

.polynomial(t, equ) ⇒ Numeric

多項式

Parameters:

Returns:



72
73
74
 
# File 'lib/when_exe/ephemeris.rb', line 72

def polynomial(t, equ)
  equ.reverse.inject(0) {|sum, v| sum * t + v}
end

.root(t0, y0 = nil, delta = 0, count = 10, error = 1E-6, &func) ⇒ Numeric

func の逆変換

Parameters:

  • t0 (Numeric) —

    独立変数の初期近似値

  • y0 (Numeric) (defaults to: nil) —

    逆変換される関数値(nil なら極値を求める)

  • delta (Numeric) (defaults to: 0) —

    戻り値が周期量である場合の周期

  • count (Numeric) (defaults to: 10) —

    収束までの最大繰り返し回数

  • error (Float) (defaults to: 1E-6) —

    収束と判断する誤差

  • func (Block) —

    逆変換される関数

Returns:

Raises:

  • (RangeError) 


292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
 
# File 'lib/when_exe/ephemeris.rb', line 292

def root(t0, y0=nil, delta=0, count=10, error=1E-6, &func)

  # 近似値0,1
  diff  = [0.01, error*1000].min
  t     = [t0-diff,         t0+diff        ]
  y     = [func.call(t[0]), func.call(t[1])]
  y.map! {|y1| _adjust(y1, y0, delta)} unless delta==0

  # 近似値2(1次関数による近似)
  t << (y0 ? (t[1]-t[0])/(y[1]-y[0])*(y0-y[0])+t[0] : t0)

  # 繰り返し
  i = count
  while ((t[2]-t[1]).abs > error) && (i > 0)
    # 予備計算
    y << func.call(t[2])
    y[2] = _adjust(y[2], y0, delta) unless delta==0
    break if y0 && (y[2]-y0).abs <= error / 10

    # printf("t=%20.7f,L=%20.7f\n",t[2],y[2])
    t01     = t[0]-t[1]
    t02,y02 = t[0]-t[2], y[0]-y[2]
    t12,y12 = t[1]-t[2], y[1]-y[2]

    # 2次関数の係数
    a = ( y02     / t02 - y12     / t12) / t01
    b = (-y02*t12 / t02 + y12*t02 / t12) / t01
    c = y[2]

    if y0
      # 判別式
      if (discriminant = b*b-4*a*(c-y0)) < 0
        # 近似値(1次関数による近似)
        if y12 == 0
          i = -1
          break
        end
        t << t12/y12*(y0-y[1]) + t[1]
      else
        # 近似値(2次関数による近似)
        sqrtd = Math.sqrt(discriminant)
        sqrtd = -sqrtd if b < 0
        t << (t[2] + 2*(y0-c)/(b+sqrtd)) # <-桁落ち回避
      end
    else
      t << (t[2] - b / (2*a))
    end

    t.shift
    y.shift
    i -= 1
  end
  return t[2] if i > 0
  return t[1] + (t[0]-t[1]) / (y[0]-y[1]) * (y0-y[1]) if y0 && count < 6
  raise RangeError, "The result does not converge - #{t}->#{y}."
end

.sinc(x) ⇒ Numeric

円周単位のsin

Parameters:

Returns:



148
149
150
 
# File 'lib/when_exe/ephemeris.rb', line 148

def sinc(x)
  sin(x * CIRCLE)
end

.sind(x) ⇒ Numeric

度のsin

Parameters:

Returns:



127
128
129
 
# File 'lib/when_exe/ephemeris.rb', line 127

def sind(x)
  sin(x * DEG)
end

.tanc(x) ⇒ Numeric

円周単位のtan

Parameters:

Returns:



155
156
157
 
# File 'lib/when_exe/ephemeris.rb', line 155

def tanc(x)
  tan(x * CIRCLE)
end

.tand(x) ⇒ Numeric

度のtan

Parameters:

Returns:



134
135
136
 
# File 'lib/when_exe/ephemeris.rb', line 134

def tand(x)
  tan(x * DEG)
end

.trigonometric(t, equ, dl = 0.0, count = 0) ⇒ Numeric

三角関数の和

Parameters:

  • t (Numeric) —

    独立変数

  • equ (Array<Numeric>) —

    係数の Array

  • dl (Numeric) (defaults to: 0.0) —

    位相補正値(デフォルト 補正なし)

  • count (Integer) (defaults to: 0) —

    打ち切り項数(デフォルト 打ち切りなし)

Returns:



86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
 
# File 'lib/when_exe/ephemeris.rb', line 86

def trigonometric(t, equ, dl=0.0, count=0)
  t2 = t * t
  equ[0..(count-1)].inject(0) do |sum, v|
    op, epoch, freq, amp, sqr = v
    ds = epoch + t * freq
    if (op < 0) # 直線
      ds += t2 * amp
    else        # 三角関数
      ds += t2 * sqr if sqr
      ds += dl  if (op[2]==1) # delta L is exist
      ds  = amp * ((op[0]==1) ? cosd(ds) : sind(ds))
      ds *= t   if (op[1]==1) # time proportional
    end
    sum += ds
  end
end