module MkSunmoon::Compute

Public Instance Methods

compute_angle_ecliptic(jy, t, lambda, beta) click to toggle source
時刻(t)における黄経、黄緯(λ(jy),β(jy))の天体の方位角(ang)計算

@param:  lambda (天体の黄経(λ(T)(度)))
@param:  beta   (天体の黄緯(β(T)(度)))
@param:  jy     (ユリウス年)
@param:  t      (時刻(0.xxxx日))
@return: angle  (角度(xx.x度))
# File lib/mk_sunmoon/compute.rb, line 513
def compute_angle_ecliptic(jy, t, lambda, beta)
  res = eclip2equat(jy, lambda, beta)
  return compute_angle_equator(jy, t, *res)
end
compute_angle_equator(jy, t, alpha, delta) click to toggle source
時刻(t)における赤経、赤緯(α(jy),δ(jy))(度)の天体の方位角(ang)計算

@param:  jy    (ユリウス年)
@param:  t     (時刻 (0.xxxx日))
@param:  alpha (天体の赤経(α(jy)(度)))
@param:  delta (天体の赤緯(δ(jy)(度)))
@return: angle (角度(xx.x度))
# File lib/mk_sunmoon/compute.rb, line 527
def compute_angle_equator(jy, t, alpha, delta)
  time_sidereal = compute_sidereal_time(jy, t)
  hour_angle = time_sidereal - alpha
  a_0  = -1.0 * Math.cos(Const::PI_180 * delta) \
              * Math.sin(Const::PI_180 * hour_angle)
  a_1  = Math.sin(Const::PI_180 * delta) \
       * Math.cos(Const::PI_180 * @lat)  \
       - Math.cos(Const::PI_180 * delta) \
       * Math.sin(Const::PI_180 * @lat)  \
       * Math.cos(Const::PI_180 * hour_angle)
  angle  = Math.atan(a_0 / a_1) / Const::PI_180
  angle += 360.0 if a_1 > 0.0 && angle < 0.0
  angle += 180.0 if a_1 < 0.0
  return angle
end
compute_beta_moon(jy) click to toggle source
月視黄緯の計算

@param:  jy (ユリウス年(JST))
@return: lambda
# File lib/mk_sunmoon/compute.rb, line 313
def compute_beta_moon(jy)
  bm  =  0.0005 * Math.sin(Const::PI_180 * norm_angle(307.0   +    19.4    * jy))
  bm +=  0.0026 * Math.sin(Const::PI_180 * norm_angle( 55.0   +    19.34   * jy))
  bm +=  0.0040 * Math.sin(Const::PI_180 * norm_angle(119.5   +     1.33   * jy))
  bm +=  0.0043 * Math.sin(Const::PI_180 * norm_angle(322.1   +    19.36   * jy))
  bm +=  0.0267 * Math.sin(Const::PI_180 * norm_angle(234.95  +    19.341  * jy))
  bt  =  0.0003 * Math.sin(Const::PI_180 * norm_angle(234.0   + 19268.0    * jy))
  bt +=  0.0003 * Math.sin(Const::PI_180 * norm_angle(146.0   +  3353.3    * jy))
  bt +=  0.0003 * Math.sin(Const::PI_180 * norm_angle(107.0   + 18149.4    * jy))
  bt +=  0.0003 * Math.sin(Const::PI_180 * norm_angle(205.0   + 22642.7    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(147.0   + 14097.4    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle( 13.0   +  9325.4    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle( 81.0   + 10242.6    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(238.0   + 23281.3    * jy))
  bt +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(311.0   +  9483.9    * jy))
  bt +=  0.0005 * Math.sin(Const::PI_180 * norm_angle(239.0   +  4193.4    * jy))
  bt +=  0.0005 * Math.sin(Const::PI_180 * norm_angle(280.0   +  8485.3    * jy))
  bt +=  0.0006 * Math.sin(Const::PI_180 * norm_angle( 52.0   + 13617.3    * jy))
  bt +=  0.0006 * Math.sin(Const::PI_180 * norm_angle(224.0   +  5590.7    * jy))
  bt +=  0.0007 * Math.sin(Const::PI_180 * norm_angle(294.0   + 13098.7    * jy))
  bt +=  0.0008 * Math.sin(Const::PI_180 * norm_angle(326.0   +  9724.1    * jy))
  bt +=  0.0008 * Math.sin(Const::PI_180 * norm_angle( 70.0   + 17870.7    * jy))
  bt +=  0.0010 * Math.sin(Const::PI_180 * norm_angle( 18.0   + 12978.66   * jy))
  bt +=  0.0011 * Math.sin(Const::PI_180 * norm_angle(138.3   + 19147.99   * jy))
  bt +=  0.0012 * Math.sin(Const::PI_180 * norm_angle(148.2   +  4851.36   * jy))
  bt +=  0.0012 * Math.sin(Const::PI_180 * norm_angle( 38.4   +  4812.68   * jy))
  bt +=  0.0013 * Math.sin(Const::PI_180 * norm_angle(155.4   +   379.35   * jy))
  bt +=  0.0013 * Math.sin(Const::PI_180 * norm_angle( 95.8   +  4472.03   * jy))
  bt +=  0.0014 * Math.sin(Const::PI_180 * norm_angle(219.2   +   299.96   * jy))
  bt +=  0.0015 * Math.sin(Const::PI_180 * norm_angle( 45.8   +  9964.00   * jy))
  bt +=  0.0015 * Math.sin(Const::PI_180 * norm_angle(211.1   +  9284.69   * jy))
  bt +=  0.0016 * Math.sin(Const::PI_180 * norm_angle(135.7   +   420.02   * jy))
  bt +=  0.0017 * Math.sin(Const::PI_180 * norm_angle( 99.8   + 14496.06   * jy))
  bt +=  0.0018 * Math.sin(Const::PI_180 * norm_angle(270.8   +  5192.01   * jy))
  bt +=  0.0018 * Math.sin(Const::PI_180 * norm_angle(243.3   +  8206.68   * jy))
  bt +=  0.0019 * Math.sin(Const::PI_180 * norm_angle(230.7   +  9244.02   * jy))
  bt +=  0.0021 * Math.sin(Const::PI_180 * norm_angle(170.1   +  1058.66   * jy))
  bt +=  0.0022 * Math.sin(Const::PI_180 * norm_angle(331.4   + 13377.37   * jy))
  bt +=  0.0025 * Math.sin(Const::PI_180 * norm_angle(196.5   +  8605.38   * jy))
  bt +=  0.0034 * Math.sin(Const::PI_180 * norm_angle(319.9   +  4433.31   * jy))
  bt +=  0.0042 * Math.sin(Const::PI_180 * norm_angle(103.9   + 18509.35   * jy))
  bt +=  0.0043 * Math.sin(Const::PI_180 * norm_angle(307.6   +  5470.66   * jy))
  bt +=  0.0082 * Math.sin(Const::PI_180 * norm_angle(144.9   +  3713.33   * jy))
  bt +=  0.0088 * Math.sin(Const::PI_180 * norm_angle(176.7   +  4711.96   * jy))
  bt +=  0.0093 * Math.sin(Const::PI_180 * norm_angle(277.4   +  8845.31   * jy))
  bt +=  0.0172 * Math.sin(Const::PI_180 * norm_angle(  3.18  + 14375.997  * jy))
  bt +=  0.0326 * Math.sin(Const::PI_180 * norm_angle(328.96  + 13737.362  * jy))
  bt +=  0.0463 * Math.sin(Const::PI_180 * norm_angle(172.55  +   698.667  * jy))
  bt +=  0.0554 * Math.sin(Const::PI_180 * norm_angle(194.01  +  8965.374  * jy))
  bt +=  0.1732 * Math.sin(Const::PI_180 * norm_angle(142.427 +  4073.3220 * jy))
  bt +=  0.2777 * Math.sin(Const::PI_180 * norm_angle(138.311 +    60.0316 * jy))
  bt +=  0.2806 * Math.sin(Const::PI_180 * norm_angle(228.235 +  9604.0088 * jy))
  bt +=  5.1282 * Math.sin(Const::PI_180 * norm_angle( 93.273 +  4832.0202 * jy + bm))
  return bt
end
compute_diff_moon(jy) click to toggle source
月の視差の計算

@param:  jy (ユリウス年)
@return: t  (出入時刻(0.xxxx日))
# File lib/mk_sunmoon/compute.rb, line 375
def compute_diff_moon(jy)
  t  =  0.0003 * Math.sin(Const::PI_180 * norm_angle(227.0  +  4412.0   * jy))
  t +=  0.0004 * Math.sin(Const::PI_180 * norm_angle(194.0  +  3773.4   * jy))
  t +=  0.0005 * Math.sin(Const::PI_180 * norm_angle(329.0  +  8545.4   * jy))
  t +=  0.0009 * Math.sin(Const::PI_180 * norm_angle(100.0  + 13677.3   * jy))
  t +=  0.0028 * Math.sin(Const::PI_180 * norm_angle(  0.0  +  9543.98  * jy))
  t +=  0.0078 * Math.sin(Const::PI_180 * norm_angle(325.7  +  8905.34  * jy))
  t +=  0.0095 * Math.sin(Const::PI_180 * norm_angle(190.7  +  4133.35  * jy))
  t +=  0.0518 * Math.sin(Const::PI_180 * norm_angle(224.98 +  4771.989 * jy))
  t +=  0.9507 * Math.sin(Const::PI_180 * norm_angle(90.0))
  return t
end
compute_dist_sun(jy) click to toggle source
太陽の距離の計算

@param:  jy (ユリウス年(JST))
@return: distance
# File lib/mk_sunmoon/compute.rb, line 218
def compute_dist_sun(jy)
  r_sun  = 0.000007 * Math.sin(Const::PI_180 * norm_angle(156.0 +  329.6  * jy))
  r_sun += 0.000007 * Math.sin(Const::PI_180 * norm_angle(254.0 +  450.4  * jy))
  r_sun += 0.000013 * Math.sin(Const::PI_180 * norm_angle( 27.8 + 4452.67 * jy))
  r_sun += 0.000030 * Math.sin(Const::PI_180 * norm_angle( 90.0))
  r_sun += 0.000091 * Math.sin(Const::PI_180 * norm_angle(265.1 +  719.98 * jy))
  r_sun += (0.007256 - 0.0000002 * jy) * Math.sin(Const::PI_180 * norm_angle(267.54 + 359.991 * jy))
  r_sun  = 10.0 ** r_sun
  return r_sun
end
compute_dt(year, month, day) click to toggle source
ΔT の計算

* 1972-01-01 以降、うるう秒挿入済みの年+αまでは、以下で算出
    TT - UTC = ΔT + DUT1 = TAI + 32.184 - UTC = ΔAT + 32.184
  [うるう秒実施日一覧](http://jjy.nict.go.jp/QandA/data/leapsec.html)

@param:  year
@param:  month
@param:  day
@return: dt
# File lib/mk_sunmoon/compute.rb, line 61
def compute_dt(year, month, day)
  ymd = sprintf("%04d%02d%02d", year, month, day)
  tm = MkTime.new(ymd)
  return tm.dt
end
compute_height_ecliptic(jy, t, lambda, beta) click to toggle source
時刻(t)における黄経、黄緯(λ(jy),β(jy))の天体の高度(height)計算

@param:  jy     (ユリウス年)
@param:  t      (時刻(0.xxxx日))
@param:  lambda (天体の黄経(λ(T)(度)))
@param:  beta   (天体の黄緯(β(T)(度)))
@return: height (高度(xx.x度))
# File lib/mk_sunmoon/compute.rb, line 552
def compute_height_ecliptic(jy, t, lambda, beta)
  res = eclip2equat(jy, lambda, beta)
  return compute_height_equator(jy, t, *res)
end
compute_height_equator(jy, t, alpha, delta) click to toggle source
時刻(t)における赤経、赤緯(α(jy),δ(jy))(度)の天体の高度(height)計算

@param:  jy     (ユリウス年)
@param:  t      (時刻 (0.xxxx日))
@param:  alpha  (天体の赤経α(jy)(度))
@param:  delta  (天体の赤緯δ(jy)(度))
@return: height (高度(xx.x度))
# File lib/mk_sunmoon/compute.rb, line 566
def compute_height_equator(jy, t, alpha, delta)
  time_sidereal = compute_sidereal_time(jy, t)
  sidereal = time_sidereal - alpha
  height  = Math.sin(Const::PI_180 * delta) \
          * Math.sin(Const::PI_180 * @lat)  \
          + Math.cos(Const::PI_180 * delta) \
          * Math.cos(Const::PI_180 * @lat)  \
          * Math.cos(Const::PI_180 * sidereal)
  height  = Math.asin(height) / Const::PI_180

  # 大気差補正
  # [ 以下の内、3-2の計算式を採用 ]
  # # 1. 日月出没計算 by「菊池さん」による計算式
  # #   [ http://kikuchisan.net/ ]
  # h = 0.0167 / Math.tan( Const::PI_180 * ( height + 8.6 / ( height + 4.4 ) ) )

  # # 2. 中川用語集による計算式 ( 5度 - 85度用 )
  # #   [ http://www.es.ris.ac.jp/~nakagawa/term_collection/yogoshu/ll/ni.htm ]
  # h  = 58.1      / Math.tan( height )
  # h -=  0.07     / Math.tan( height ) ** 3
  # h +=  0.000086 / Math.tan( height ) ** 5
  # h *= 1 / 3600.0

  # # 3-1. フランスの天文学者ラドー(R.Radau)の平均大気差と1秒程度の差で大気差を求めることが可能
  # # ( 標準的大気(気温10゚C,気圧1013.25hPa)の場合 )
  # # ( 視高度30゚以上 )
  # h  = ( 58.294  / 3600.0 ) * Math.tan( Const::PI_180 * ( 90.0 - height ) )
  # h -= (  0.0668 / 3600.0 ) * Math.tan( Const::PI_180 * ( 90.0 - height ) ) ** 3

  # 3-2. フランスの天文学者ラドー(R.Radau)の平均大気差と1秒程度の差で大気差を求めることが可能
  # ( 標準的大気(気温10゚C,気圧1013.25hPa)の場合 )
  # ( 視高度 4゚以上 )
  a = Math.tan(Const::PI_180 * (90.0 - height))
  h  = (58.76   + \
       (-0.406  + \
       (- 0.0192) \
       * a) * a) * a * 1 / 3600.0

  # # 3-3. さらに、上記の大気差(3-1,3-2)を気温、気圧を考慮する
  # # ( しかし、気温・気圧を考慮してもさほど変わりはない )
  # pres = 1013.25 # <= 変更
  # temp = 30.0    # <= 変更
  # h *= pres / 1013.25
  # h *= 283.25 / ( 273.15 + temp )

  return height + h
end
compute_hour_angle_diff(alpha, delta, time_sidereal, height, div) click to toggle source
出入点(k)の時角(tk)と天体の時角(t)との差(dt=tk-t)を計算する

@param:  alpha(R.A.)  (赤経)
@param:  delta(Decl.) (赤緯)
@param:  time_sidereal(恒星時Θ(度))
@param:  height       (出没高度(度))
@param:  div          (0: 出, 1: 入, 2: 南中)
@return: hour angle difference (時角の差 dt)
# File lib/mk_sunmoon/compute.rb, line 458
def compute_hour_angle_diff(alpha, delta, time_sidereal, height, div)
  if div == 2
    tk = 0
  else
    tk  = Math.sin(Const::PI_180 * height)
    tk -= Math.sin(Const::PI_180 * delta) * Math.sin(Const::PI_180 * @lat)
    tk /= Math.cos(Const::PI_180 * delta) * Math.cos(Const::PI_180 * @lat)
    # 出没点の時角
    return 0.0 if tk.abs > 1
    tk = Math.acos(tk) / Const::PI_180
    # tkは出のときマイナス、入のときプラス
    tk = -tk if div == 0 && tk > 0
    tk = -tk if div == 1 && tk < 0
  end
  # 天体の時角
  t = time_sidereal - alpha
  dt = tk - t
  # dtの絶対値を180°以下に調整
  while dt >  180.0; dt -= 360.0; end
  while dt < -180.0; dt += 360.0; end
  return dt
end
compute_lambda_moon(jy) click to toggle source
月視黄経の計算

@param:  jy (ユリウス年(JST))
@return: lambda
# File lib/mk_sunmoon/compute.rb, line 235
def compute_lambda_moon(jy)
  am  = 0.0006 * Math.sin(Const::PI_180 * norm_angle( 54.0   +    19.3    * jy))
  am += 0.0006 * Math.sin(Const::PI_180 * norm_angle( 71.0   +     0.2    * jy))
  am += 0.0020 * Math.sin(Const::PI_180 * norm_angle( 55.0   +    19.34   * jy))
  am += 0.0040 * Math.sin(Const::PI_180 * norm_angle(119.5   +     1.33   * jy))
  rm  = 0.0003 * Math.sin(Const::PI_180 * norm_angle(280.0   + 23221.3    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle(161.0   +    40.7    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle(311.0   +  5492.0    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle(147.0   + 18089.3    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle( 66.0   +  3494.7    * jy))
  rm += 0.0003 * Math.sin(Const::PI_180 * norm_angle( 83.0   +  3814.0    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle( 20.0   +   720.0    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle( 71.0   +  9584.7    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle(278.0   +   120.1    * jy))
  rm += 0.0004 * Math.sin(Const::PI_180 * norm_angle(313.0   +   398.7    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(332.0   +  5091.3    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(114.0   + 17450.7    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(181.0   + 19088.0    * jy))
  rm += 0.0005 * Math.sin(Const::PI_180 * norm_angle(247.0   + 22582.7    * jy))
  rm += 0.0006 * Math.sin(Const::PI_180 * norm_angle(128.0   +  1118.7    * jy))
  rm += 0.0007 * Math.sin(Const::PI_180 * norm_angle(216.0   +   278.6    * jy))
  rm += 0.0007 * Math.sin(Const::PI_180 * norm_angle(275.0   +  4853.3    * jy))
  rm += 0.0007 * Math.sin(Const::PI_180 * norm_angle(140.0   +  4052.0    * jy))
  rm += 0.0008 * Math.sin(Const::PI_180 * norm_angle(204.0   +  7906.7    * jy))
  rm += 0.0008 * Math.sin(Const::PI_180 * norm_angle(188.0   + 14037.3    * jy))
  rm += 0.0009 * Math.sin(Const::PI_180 * norm_angle(218.0   +  8586.0    * jy))
  rm += 0.0011 * Math.sin(Const::PI_180 * norm_angle(276.5   + 19208.02   * jy))
  rm += 0.0012 * Math.sin(Const::PI_180 * norm_angle(339.0   + 12678.71   * jy))
  rm += 0.0016 * Math.sin(Const::PI_180 * norm_angle(242.2   + 18569.38   * jy))
  rm += 0.0018 * Math.sin(Const::PI_180 * norm_angle(  4.1   +  4013.29   * jy))
  rm += 0.0020 * Math.sin(Const::PI_180 * norm_angle( 55.0   +    19.34   * jy))
  rm += 0.0021 * Math.sin(Const::PI_180 * norm_angle(105.6   +  3413.37   * jy))
  rm += 0.0021 * Math.sin(Const::PI_180 * norm_angle(175.1   +   719.98   * jy))
  rm += 0.0021 * Math.sin(Const::PI_180 * norm_angle( 87.5   +  9903.97   * jy))
  rm += 0.0022 * Math.sin(Const::PI_180 * norm_angle(240.6   +  8185.36   * jy))
  rm += 0.0024 * Math.sin(Const::PI_180 * norm_angle(252.8   +  9224.66   * jy))
  rm += 0.0024 * Math.sin(Const::PI_180 * norm_angle(211.9   +   988.63   * jy))
  rm += 0.0026 * Math.sin(Const::PI_180 * norm_angle(107.2   + 13797.39   * jy))
  rm += 0.0027 * Math.sin(Const::PI_180 * norm_angle(272.5   +  9183.99   * jy))
  rm += 0.0037 * Math.sin(Const::PI_180 * norm_angle(349.1   +  5410.62   * jy))
  rm += 0.0039 * Math.sin(Const::PI_180 * norm_angle(111.3   + 17810.68   * jy))
  rm += 0.0040 * Math.sin(Const::PI_180 * norm_angle(119.5   +     1.33   * jy))
  rm += 0.0040 * Math.sin(Const::PI_180 * norm_angle(145.6   + 18449.32   * jy))
  rm += 0.0040 * Math.sin(Const::PI_180 * norm_angle( 13.2   + 13317.34   * jy))
  rm += 0.0048 * Math.sin(Const::PI_180 * norm_angle(235.0   +    19.34   * jy))
  rm += 0.0050 * Math.sin(Const::PI_180 * norm_angle(295.4   +  4812.66   * jy))
  rm += 0.0052 * Math.sin(Const::PI_180 * norm_angle(197.2   +   319.32   * jy))
  rm += 0.0068 * Math.sin(Const::PI_180 * norm_angle( 53.2   +  9265.33   * jy))
  rm += 0.0079 * Math.sin(Const::PI_180 * norm_angle(278.2   +  4493.34   * jy))
  rm += 0.0085 * Math.sin(Const::PI_180 * norm_angle(201.5   +  8266.71   * jy))
  rm += 0.0100 * Math.sin(Const::PI_180 * norm_angle( 44.89  + 14315.966  * jy))
  rm += 0.0107 * Math.sin(Const::PI_180 * norm_angle(336.44  + 13038.696  * jy))
  rm += 0.0110 * Math.sin(Const::PI_180 * norm_angle(231.59  +  4892.052  * jy))
  rm += 0.0125 * Math.sin(Const::PI_180 * norm_angle(141.51  + 14436.029  * jy))
  rm += 0.0153 * Math.sin(Const::PI_180 * norm_angle(130.84  +   758.698  * jy))
  rm += 0.0305 * Math.sin(Const::PI_180 * norm_angle(312.49  +  5131.979  * jy))
  rm += 0.0348 * Math.sin(Const::PI_180 * norm_angle(117.84  +  4452.671  * jy))
  rm += 0.0410 * Math.sin(Const::PI_180 * norm_angle(137.43  +  4411.998  * jy))
  rm += 0.0459 * Math.sin(Const::PI_180 * norm_angle(238.18  +  8545.352  * jy))
  rm += 0.0533 * Math.sin(Const::PI_180 * norm_angle( 10.66  + 13677.331  * jy))
  rm += 0.0572 * Math.sin(Const::PI_180 * norm_angle(103.21  +  3773.363  * jy))
  rm += 0.0588 * Math.sin(Const::PI_180 * norm_angle(214.22  +   638.635  * jy))
  rm += 0.1143 * Math.sin(Const::PI_180 * norm_angle(  6.546 +  9664.0404 * jy))
  rm += 0.1856 * Math.sin(Const::PI_180 * norm_angle(177.525 +   359.9905 * jy))
  rm += 0.2136 * Math.sin(Const::PI_180 * norm_angle(269.926 +  9543.9773 * jy))
  rm += 0.6583 * Math.sin(Const::PI_180 * norm_angle(235.700 +  8905.3422 * jy))
  rm += 1.2740 * Math.sin(Const::PI_180 * norm_angle(100.738 +  4133.3536 * jy))
  rm += 6.2887 * Math.sin(Const::PI_180 * norm_angle(134.961 +  4771.9886 * jy + am))
  rm += norm_angle(218.3161 + 4812.67881 * jy)
  return rm
end
compute_lambda_sun(jy) click to toggle source
太陽視黄経の計算

@param:  jy (ユリウス年(JST))
@return: lambda
# File lib/mk_sunmoon/compute.rb, line 188
def compute_lambda_sun(jy)
  rm  = 0.0003 * Math.sin(Const::K * norm_angle(329.7  +   44.43  * jy))
  rm += 0.0003 * Math.sin(Const::K * norm_angle(352.5  + 1079.97  * jy))
  rm += 0.0004 * Math.sin(Const::K * norm_angle( 21.1  +  720.02  * jy))
  rm += 0.0004 * Math.sin(Const::K * norm_angle(157.3  +  299.30  * jy))
  rm += 0.0004 * Math.sin(Const::K * norm_angle(234.9  +  315.56  * jy))
  rm += 0.0005 * Math.sin(Const::K * norm_angle(291.2  +   22.81  * jy))
  rm += 0.0005 * Math.sin(Const::K * norm_angle(207.4  +    1.50  * jy))
  rm += 0.0006 * Math.sin(Const::K * norm_angle( 29.8  +  337.18  * jy))
  rm += 0.0007 * Math.sin(Const::K * norm_angle(206.8  +   30.35  * jy))
  rm += 0.0007 * Math.sin(Const::K * norm_angle(153.3  +   90.38  * jy))
  rm += 0.0008 * Math.sin(Const::K * norm_angle(132.5  +  659.29  * jy))
  rm += 0.0013 * Math.sin(Const::K * norm_angle( 81.4  +  225.18  * jy))
  rm += 0.0015 * Math.sin(Const::K * norm_angle(343.2  +  450.37  * jy))
  rm += 0.0018 * Math.sin(Const::K * norm_angle(251.3  +    0.20  * jy))
  rm += 0.0018 * Math.sin(Const::K * norm_angle(297.8  + 4452.67  * jy))
  rm += 0.0020 * Math.sin(Const::K * norm_angle(247.1  +  329.64  * jy))
  rm += 0.0048 * Math.sin(Const::K * norm_angle(234.95 +   19.341 * jy))
  rm += 0.0200 * Math.sin(Const::K * norm_angle(355.05 +  719.981 * jy))
  rm += (1.9146 - 0.00005 * jy) * Math.sin(Const::K * norm_angle(357.538 + 359.991 * jy))
  rm += norm_angle(280.4603 + 360.00769 * jy)
  return norm_angle(rm)
end
compute_moon(div) click to toggle source
月の出・入・南中の計算

@param:  div(0: 出, 1: 入, 2: 南中)
@return: [time_str, hour_angle]
# File lib/mk_sunmoon/compute.rb, line 97
def compute_moon(div)
  time_val = compute_time_moon(div)
  if time_val == 0.0
    time_str = "--:--:--"
    angle    = "---"
  else
    time_str = val2hhmmss(time_val * 24.0)
    jy = (@jd + time_val + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_moon(jy)
    beta   = compute_beta_moon(jy)
    if div == 2
      angle = compute_height_ecliptic(jy, time_val, lambda, beta)
    else
      angle = compute_angle_ecliptic(jy, time_val, lambda, beta)
    end
  end
  return [time_str, angle]
end
compute_sidereal_time(jy, t) click to toggle source
恒星時Θ(度)の計算

@param:  jy (ユリウス年(JST))
@param:  t  (時刻(0.xxxx日))
@param:  longitude
@return: sidereal time
# File lib/mk_sunmoon/compute.rb, line 439
def compute_sidereal_time(jy, t)
  val  = 325.4606       + \
        (360.007700536  + \
        (  0.00000003879) \
        * jy) * jy        \
        + 360.0 * t + @lon
  return norm_angle(val)
end
compute_sun(div) click to toggle source
日の出・入・南中の計算

@param:  div(0: 出, 1: 入, 2: 南中)
@return: [time_str, hour_angle]
# File lib/mk_sunmoon/compute.rb, line 73
def compute_sun(div)
  time_val = compute_time_sun(div)
  if time_val == 0.0
    time_str = "--:--:--"
    angle    = "---"
  else
    time_str = val2hhmmss(time_val * 24.0)
    jy = (@jd + time_val + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_sun(jy)
    if div == 2
      angle = compute_height_ecliptic(jy, time_val, lambda, 0.0)
    else
      angle = compute_angle_ecliptic(jy, time_val, lambda, 0.0)
    end
  end
  return [time_str, angle]
end
compute_time_moon(div) click to toggle source
月の出/月の入/月の南中計算

@param:  div  (0: 出, 1: 入, 2: 南中)
@return: time (0.xxxx日)
# File lib/mk_sunmoon/compute.rb, line 156
def compute_time_moon(div)
  rev = 1    # 補正値初期値
  t   = 0.5  # 逐次計算時刻(日)初期設定

  while rev.abs > Const::CONVERGE
    jy = (@jd + t + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_moon(jy)
    beta   = compute_beta_moon(jy)
    alpha, delta = eclip2equat(jy, lambda, beta)
    unless div == 2  # 南中のときは計算しない
      diff = compute_diff_moon(jy)
      height = -1 * Const::REFRACTION - @incl + diff
    end
    time_sidereal = compute_sidereal_time(jy, t)
    hour_angle_diff = compute_hour_angle_diff(
      alpha, delta, time_sidereal, height, div
    )
    # 仮定時刻に対する補正値
    rev = hour_angle_diff / 347.8
    t += rev
  end
  # 月の出/月の入りがない場合は 0 とする
  t = 0.0 if t < 0.0 || t >= 1.0
  return t
end
compute_time_sun(div) click to toggle source
日の出/日の入/日の南中計算

@param:  div  (0: 出, 1: 入, 2: 南中)
@return: time (0.xxxx日)
# File lib/mk_sunmoon/compute.rb, line 122
def compute_time_sun(div)
  rev = 1    # 補正値初期値
  t   = 0.5  # 逐次計算時刻(日)初期設定

  while rev.abs > Const::CONVERGE
    jy = (@jd + t + @dt / 86400.0 - 2451545.0) / 365.25
    lambda = compute_lambda_sun(jy)
    dist   = compute_dist_sun(jy)
    alpha, delta = eclip2equat(jy, lambda, 0.0)
    unless div == 2  # 南中のときは計算しない
      r = 0.266994 / dist
      diff = 0.0024428 / dist
      height = -1 * r - Const::REFRACTION - @incl + diff
    end
    time_sidereal = compute_sidereal_time(jy, t)
    hour_angle_diff = compute_hour_angle_diff(
      alpha, delta, time_sidereal, height, div
    )
    return 0.0 if hour_angle_diff == 0.0
    # 仮定時刻に対する補正値
    rev = hour_angle_diff / 360.0
    t += rev
  end
  # 日の出・入がない場合は 0 とする
  t = 0.0 if t < 0.0 || t >= 1.0
  return t
end
eclip2equat(jy, lambda, beta) click to toggle source
黄道座標 -> 赤道座標変換

@param:  jy (ユリウス年(JST))
@param:  lambda (黄経)
@param:  beta   (黄緯)
@return: [alpha(R.A.), delta(Decl.)]
# File lib/mk_sunmoon/compute.rb, line 415
def eclip2equat(jy, lambda, beta)
  eps = (23.439291 - 0.000130042 * jy) * Const::PI_180
  lambda *= Const::PI_180
  beta   *= Const::PI_180
  a  =      Math.cos(beta) * Math.cos(lambda)
  b  = -1 * Math.sin(beta) * Math.sin(eps)
  b +=      Math.cos(beta) * Math.sin(lambda) * Math.cos(eps)
  c  =      Math.sin(beta) * Math.cos(eps)
  c +=      Math.cos(beta) * Math.sin(lambda) * Math.sin(eps)
  alpha  = b / a
  alpha  = Math.atan(alpha) / Const::PI_180
  alpha += 180 if a < 0
  delta  = Math.asin(c) / Const::PI_180
  return [alpha, delta]
end
gc2jd(year, month, day, hour = 0, min = 0, sec = 0) click to toggle source
Gregorian Calendar -> Julian Day

* フリーゲルの公式を使用する
  [ JD ] = int( 365.25 × year )
         + int( year / 400 )
         - int( year / 100 )
         + int( 30.59 ( month - 2 ) )
         + day
         + 1721088
  ※上記の int( x ) は厳密には、x を超えない最大の整数
    ( ちなみに、[ 準JD ]を求めるなら + 1721088.5 が - 678912 となる )

@param:  year
@param:  month
@param:  day
@param:  hour
@param:  minute
@param:  second
@return: jd ( ユリウス日 )
# File lib/mk_sunmoon/compute.rb, line 28
def gc2jd(year, month, day, hour = 0, min = 0, sec = 0)
  # 1月,2月は前年の13月,14月とする
  if month < 3
    year  -= 1
    month += 12
  end
  # 日付(整数)部分計算
  jd  = (365.25 * year).truncate
  jd += (year / 400.0).truncate
  jd -= (year / 100.0).truncate
  jd += (30.59 * (month - 2)).truncate
  jd += day
  jd += 1721088.125
  # 時間(小数)部分計算
  t  = sec / 3600.0
  t += min / 60.0
  t += hour
  t  = t / 24.0
  return jd + t
end
norm_angle(angle) click to toggle source
角度の正規化

@param:  angle
@return: angle
# File lib/mk_sunmoon/compute.rb, line 394
def norm_angle(angle)
  if angle < 0
    angle1  = angle * (-1)
    angle2  = (angle1 / 360.0).truncate
    angle1 -= 360 * angle2
    angle1  = 360 - angle1
  else
    angle1  = (angle / 360.0).truncate
    angle1  = angle - 360.0 * angle1
  end
  return angle1
end
val2hhmmss(val) click to toggle source
時刻:数値->時刻:時分変換 (xx.xxxx -> hh:mm:ss)

@param:  time value  (時刻, xx.xxxx日)
@return: time string (時刻, hh:mm:ss)
# File lib/mk_sunmoon/compute.rb, line 487
def val2hhmmss(val)
  val_h = val.truncate            # 整数部(時)
  val_2 = val - val_h             # 小数部
  val_m = (val_2 * 60).truncate   # (分)計算
  val_3 = val_2 - (val_m / 60.0)  # (秒)計算
  val_s = (val_3 * 60 * 60).round
  if val_s == 60
    val_s = 0
    val_m +=1
  end
  if val_m == 60
    val_m = 0
    val_h += 1
  end
  return sprintf("%02d:%02d:%02d", val_h, val_m, val_s)
end