require "qreki/version"
require "date"

# =========================================================================
# Based on 「旧暦計算サンプルプログラム」
#     Copyright (C) 1993,1994 by H.Takano
#     http://www.vector.co.jp/soft/dos/personal/se016093.html
# =========================================================================

module Qreki

  class Srk
    attr_accessor :year, :month, :day
  end

  class Qrk
    attr_accessor :year, :uruu, :month, :day, :rokuyou, :sekki
  end

  $k = Math::PI / 180.0  # Pi
  $tz = - (9.0 / 24.0)   # タイムゾーンオフセット
  $rm_sun0 = 0           # 太陽黄経

  #=========================================================================
  # 旧暦算出
  #
  # 引数　 .... 日付
  # 戻り値
  #    q.year  : 旧暦年
  #    q.uruu  : 平月／閏月 flag .... 平月:false 閏月:true
  #    q.month  : 旧暦月
  #    q.day : 旧暦日
  #    q.rokuyou : 六曜
  #    q.sekki : 二十四節気
  #=========================================================================
  def self.calc(year, month, day)
    calc_from_date(Date.new(year, month, day))
  end

  def self.calc_from_date(tm)
    array = calc_kyureki(tm.year, tm.month, tm.day)
    qreki = Qrk.new
    qreki.year    = array[0]
    qreki.uruu    = array[1]
    qreki.month   = array[2]
    qreki.day     = array[3]
    qreki.rokuyou = rokuyou(tm.year, tm.month, tm.day)
    qreki.sekki   = sekki(tm.year, tm.month, tm.day)
    qreki
  end

  # =========================================================================
  #  新暦に対応する、旧暦を求める。
  #
  #  呼び出し時にセットする変数：
  #    year : 計算する年（JSTユリウス日）
  #    mon  : 計算する月（JSTユリウス日）
  #    day  : 計算する日（JSTユリウス日）
  #  戻り値：
  #    this.year  : 旧暦年
  #    this.uruu  : 平月／閏月 flag .... 平月:false 閏月:true
  #    this.month : 旧暦月
  #    this.day  : 旧暦日
  # =========================================================================
  def self.calc_kyureki(year, mon, day)

    chu     = Array.new(5, [])
    m       = Array.new(5).map{Array.new(3,0)}
    saku    = Array.new
    kyureki = Array.new

    tm = ymdt2jd(year, mon, day, 0, 0, 0)

    # -----------------------------------------------------------------------
    #  計算対象の直前にあたる二分二至の時刻を求める
    # -----------------------------------------------------------------------
    chu[0] = calc_chu(tm, 90)

    # -----------------------------------------------------------------------
    #  上で求めた二分二至の時の太陽黄経をもとに朔日行列の先頭に月名をセット
    # -----------------------------------------------------------------------
    m[0][0] = ( $rm_sun0 / 30.0 ).to_i + 2

    # -----------------------------------------------------------------------
    #  中気の時刻を計算（３回計算する）
    #  chu[i]:中気の時刻
    # -----------------------------------------------------------------------
    (1..3).each do |i|
      chu[i] = calc_chu(chu[i - 1] + 32, 30)
    end
    # -----------------------------------------------------------------------
    #  計算対象の直前にあたる二分二至の直前の朔の時刻を求める
    # -----------------------------------------------------------------------
    saku[0] = calc_saku(chu[0])

    # -----------------------------------------------------------------------
    #  朔の時刻を求める
    # -----------------------------------------------------------------------
    (1..4).each do |i|
      saku[i] = calc_saku(saku[i - 1] + 30)
      # 前と同じ時刻を計算した場合（両者の差が26日以内）には、初期値を
      # +33日にして再実行させる。
      if (saku[i - 1].to_i.abs - saku[i].to_i) <= 26
        saku[i] = calc_saku(saku[i - 1] + 35)
      end
    end

    # -----------------------------------------------------------------------
    #  saku[1]が二分二至の時刻以前になってしまった場合には、朔をさかのぼり過ぎ
    #  たと考えて、朔の時刻を繰り下げて修正する。
    #  その際、計算もれ（saku[4]）になっている部分を補うため、朔の時刻を計算
    #  する。（近日点通過の近辺で朔があると起こる事があるようだ...？）
    # -----------------------------------------------------------------------
    if saku[1].to_i <= chu[0].to_i
      (0..4).each do |i|
        saku[i] = saku[i + 1]
      end
      saku[4] = calc_saku(saku[3] + 35)

    # -----------------------------------------------------------------------
    #  saku[0]が二分二至の時刻以後になってしまった場合には、朔をさかのぼり足
    #  りないと見て、朔の時刻を繰り上げて修正する。
    #  その際、計算もれ（saku[0]）になっている部分を補うため、朔の時刻を計算
    #  する。（春分点の近辺で朔があると起こる事があるようだ...？）
    # -----------------------------------------------------------------------
    elsif saku[0].to_i > chu[0].to_i
      [4,3,2,1].each do |i|
        saku[i] = saku[i - 1]
      end
      saku[0] = calc_saku(saku[0] - 27)
    end

    # -----------------------------------------------------------------------
    #  閏月検索Ｆｌａｇセット
    #  （節月で４ヶ月の間に朔が５回あると、閏月がある可能性がある。）
    #  lap=false:平月  lap=true:閏月
    # -----------------------------------------------------------------------
    lap = (saku[4].to_i <= chu[3].to_i)

    # -----------------------------------------------------------------------
    #  朔日行列の作成
    #  m[i][0] ... 月名（1:正月 2:２月 3:３月 ....）
    #  m[i][1] .... 閏フラグ（false:平月 true:閏月）
    #  m[i][2] ...... 朔日のjd
    # -----------------------------------------------------------------------
    m[0][1] = false
    m[0][2] = saku[0].to_i
    (1..4).each do |i|
      if lap && i > 1
        if (chu[i - 1] <= saku[i - 1].to_i) || (chu[i - 1] >= saku[i].to_i)
          m[i-1][0] = m[i-2][0]
          m[i-1][1] = true
          m[i-1][2] = saku[i - 1].to_i
          lap = false
        end
      end
      m[i][0] = m[i-1][0] + 1
      m[i][0] -= 12 if m[i][0] > 12
      m[i][2] = saku[i].to_i
      m[i][1] = false
    end

    # -----------------------------------------------------------------------
    #  朔日行列から旧暦を求める。
    # -----------------------------------------------------------------------
    state = 0
    i_tmp = 0
    (1..4).each do |i|
      i_tmp = i
      if tm.to_i < m[i][2].to_i
        state = 1
        break
      elsif tm.to_i == m[i][2].to_i
        state = 2
        break
      end
    end
    i_tmp -= 1 if state == 0 || state == 1
    i = i_tmp

    kyureki[1] = m[i][1]
    kyureki[2] = m[i][0]
    kyureki[3] = tm.to_i - m[i][2].to_i + 1

    # -----------------------------------------------------------------------
    #  旧暦年の計算
    #  （旧暦月が10以上でかつ新暦月より大きい場合には、
    #    まだ年を越していないはず...）
    # -----------------------------------------------------------------------
    a = jd2ymdt(tm)
    kyureki[0] = a[0]
    if kyureki[2] > 9 && kyureki[2] > a[1]
      kyureki[0] -= 1
    end

    [kyureki[0], kyureki[1], kyureki[2], kyureki[3]]
  end

  # =========================================================================
  #  直前の二分二至／中気の時刻を求める
  #
  #  パラメータ
  #    tm ............ 計算基準となる時刻（JSTユリウス日）
  #    longitude ..... 求める対象（90:二分二至,30:中気））
  #  戻り値
  #    求めた時刻（JSTユリウス日）を返す。
  #    グローバル変数 $rm_sun0 に、その時の太陽黄経をセットする。
  #
  #  ※ 引数、戻り値ともユリウス日で表し、時分秒は日の小数で表す。
  #     力学時とユリウス日との補正時刻=0.0secと仮定
  # =========================================================================
  def self.calc_chu(tm, longitude)
    # -----------------------------------------------------------------------
    #  時刻引数を小数部と整数部とに分解する（精度を上げるため）
    # -----------------------------------------------------------------------
    tm1 = tm.to_i
    tm2 = tm - tm1 + $tz  # JST -> UTC

    # -----------------------------------------------------------------------
    # 直前の二分二至の黄経 λsun0 を求める
    # -----------------------------------------------------------------------
    t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
    rm_sun  = longitude_sun(t)
    $rm_sun0 = longitude * (rm_sun / longitude.to_f).to_i

    # -----------------------------------------------------------------------
    #  繰り返し計算によって直前の二分二至の時刻を計算する
    #  （誤差が±1.0 sec以内になったら打ち切る。）
    # -----------------------------------------------------------------------
    delta_t1 ||= 0.0
    delta_t2 ||= 1.0

    while (delta_t1 + delta_t2).abs > (1.0 / 86400.0)
      # -------------------------------------------------------------------
      #  λsun(t) を計算
      #    t = (tm + .5 - 2451545) / 36525;
      # -------------------------------------------------------------------
      t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
      rm_sun = longitude_sun(t)

      # -------------------------------------------------------------------
      #  黄経差 Δλ＝λsun －λsun0
      # -------------------------------------------------------------------
      delta_rm = rm_sun - $rm_sun0

      # -------------------------------------------------------------------
      #  Δλの引き込み範囲（±180°）を逸脱した場合には、補正を行う
      # -------------------------------------------------------------------
      if delta_rm > 180.0
        delta_rm -= 360.0
      elsif delta_rm < -180.0
        delta_rm += 360.0
      end

      # -------------------------------------------------------------------
      #  時刻引数の補正値 Δt
      #  delta_t = delta_rm * 365.2 / 360;
      # -------------------------------------------------------------------
      delta_t1 = (delta_rm * 365.2 / 360.0).to_i
      delta_t2 = (delta_rm * 365.2 / 360.0) - delta_t1

      # -------------------------------------------------------------------
      #  時刻引数の補正
      #  tm -= delta_t;
      # -------------------------------------------------------------------
      tm1 = tm1 - delta_t1
      tm2 = tm2 - delta_t2
      if tm2 < 0
        tm1 -= 1.0
        tm2 += 1.0
      end
    end

    # -----------------------------------------------------------------------
    #  戻り値の作成
    #    時刻引数を合成し、戻り値（JSTユリウス日）とする
    # -----------------------------------------------------------------------
    tm2 + tm1 - $tz
  end

  # =========================================================================
  #  直前の朔の時刻を求める
  #
  #  呼び出し時にセットする変数
  #    tm ........ 計算基準の時刻（JSTユリウス日）
  #  戻り値
  #    朔の時刻
  #
  #  ※ 引数、戻り値ともJSTユリウス日で表し、時分秒は日の小数で表す。
  #     力学時とユリウス日との補正時刻=0.0secと仮定
  # =========================================================================
  def self.calc_saku(tm)
    # -----------------------------------------------------------------------
    #  ループカウンタのセット
    # -----------------------------------------------------------------------
    lc = 1

    # -----------------------------------------------------------------------
    #  時刻引数を小数部と整数部とに分解する（精度を上げるため）
    # -----------------------------------------------------------------------
    tm1 = tm.to_i
    tm2 = tm - tm1 + $tz  # JST -> UTC

    # -----------------------------------------------------------------------
    #  繰り返し計算によって朔の時刻を計算する
    #  （誤差が±1.0 sec以内になったら打ち切る。）
    # -----------------------------------------------------------------------
    delta_t1 ||= 0.0
    delta_t2 ||= 1.0

    while ( delta_t1 + delta_t2 ).abs > ( 1.0 / 86400.0 )
      lc += 1

      # -------------------------------------------------------------------
      #  太陽の黄経λsun(t) ,月の黄経λmoon(t) を計算
      #    t = (tm + .5 - 2451545) / 36525;
      # -------------------------------------------------------------------
      t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0
      rm_sun = longitude_sun(t)
      rm_moon = longitude_moon(t)

      # -------------------------------------------------------------------
      #  月と太陽の黄経差Δλ
      #  Δλ＝λmoon－λsun
      # -------------------------------------------------------------------
      delta_rm = rm_moon - rm_sun

      # -------------------------------------------------------------------
      #  ループの１回目（lc=1）で delta_rm < 0 の場合には引き込み範囲に
      #  入るように補正する
      # -------------------------------------------------------------------
      if lc==1 && delta_rm < 0
        delta_rm = normalization_angle(delta_rm)
      # -------------------------------------------------------------------
      #  春分の近くで朔がある場合（0 ≦λsun≦ 20）で、
      #  月の黄経λmoon≧300 の場合には、
      #  Δλ＝ 360 － Δλ と計算して補正する
      # -------------------------------------------------------------------
      elsif rm_sun >= 0 && rm_sun <= 20 && rm_moon >= 300
        delta_rm = normalization_angle(delta_rm)
        delta_rm = 360 - delta_rm
      # -------------------------------------------------------------------
      #  Δλの引き込み範囲（±40°）を逸脱した場合には、補正を行う
      # -------------------------------------------------------------------
      elsif delta_rm.abs > 40
        delta_rm = normalization_angle(delta_rm)
      end

      # -------------------------------------------------------------------
      #  時刻引数の補正値 Δt
      #  delta_t = delta_rm * 29.530589 / 360;
      # -------------------------------------------------------------------
      delta_t1 = ( delta_rm * 29.530589 / 360.0 ).to_i
      delta_t2 = ( delta_rm * 29.530589 / 360.0 ) - delta_t1

      # -------------------------------------------------------------------
      #  時刻引数の補正
      #  tm -= delta_t;
      # -------------------------------------------------------------------
      tm1 = tm1 - delta_t1
      tm2 = tm2 - delta_t2
      if( tm2 < 0 )
        tm1 -= 1.0
        tm2 += 1.0
      end

      # -------------------------------------------------------------------
      #  ループ回数が15回になったら、初期値 tm を tm-26 とする。
      # -------------------------------------------------------------------
      if( lc == 15 && (delta_t1 + delta_t2).abs > (1 / 86400.0) )
        tm1 = (tm - 26).to_i
        tm2 = 0
      # -------------------------------------------------------------------
      #  初期値を補正したにも関わらず、振動を続ける場合には初期値を答えとし
      #  て返して強制的にループを抜け出して異常終了させる。
      # -------------------------------------------------------------------
      elsif( lc > 30 && (delta_t1 + delta_t2).abs > (1 / 86400.0) )
        tm1 = tm
        tm2 = 0
        break
      end
    end

    # -----------------------------------------------------------------------
    #  戻り値の作成
    #    時刻引数を合成し、戻り値（ユリウス日）とする
    # -----------------------------------------------------------------------
    tm2 + tm1 - $tz
  end

  # =========================================================================
  #  角度の正規化を行う。すなわち引数の範囲を ０≦θ＜３６０ にする。
  # =========================================================================
  def self.normalization_angle(angle)
    if angle >= 0.0
      return         angle - 360.0 * ( angle / 360.0 ).to_i
    else
      return 360.0 + angle - 360.0 * ( angle / 360.0 ).to_i
    end
  end

  # =========================================================================
  #  太陽の黄経 λsun(t) を計算する（t は力学時）
  # =========================================================================
  def self.longitude_sun(t)
    # -----------------------------------------------------------------------
    #  摂動項の計算
    # -----------------------------------------------------------------------
    th  = 0.0004 * Math.cos( $k * normalization_angle(  31557.0  * t + 161.0  ) )
    th += 0.0004 * Math.cos( $k * normalization_angle(  29930.0  * t +  48.0  ) )
    th += 0.0005 * Math.cos( $k * normalization_angle(   2281.0  * t + 221.0  ) )
    th += 0.0005 * Math.cos( $k * normalization_angle(    155.0  * t + 118.0  ) )
    th += 0.0006 * Math.cos( $k * normalization_angle(  33718.0  * t + 316.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(   9038.0  * t +  64.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(   3035.0  * t + 110.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(  65929.0  * t +  45.0  ) )
    th += 0.0013 * Math.cos( $k * normalization_angle(  22519.0  * t + 352.0  ) )
    th += 0.0015 * Math.cos( $k * normalization_angle(  45038.0  * t + 254.0  ) )
    th += 0.0018 * Math.cos( $k * normalization_angle( 445267.0  * t + 208.0  ) )
    th += 0.0018 * Math.cos( $k * normalization_angle(     19.0  * t + 159.0  ) )
    th += 0.0020 * Math.cos( $k * normalization_angle(  32964.0  * t + 158.0  ) )
    th += 0.0200 * Math.cos( $k * normalization_angle(  71998.1  * t + 265.1  ) )
    th -= 0.0048 * Math.cos( $k * normalization_angle(  35999.05 * t + 267.52 ) ) * t
    th += 1.9147 * Math.cos( $k * normalization_angle(  35999.05 * t + 267.52 ) )

    # -----------------------------------------------------------------------
    #  比例項の計算
    # -----------------------------------------------------------------------
    ang = normalization_angle( 36000.7695 * t )
    ang = normalization_angle( ang + 280.4659 )
    th  = normalization_angle( th + ang )

    th
  end

  # =========================================================================
  #  月の黄経 λmoon(t) を計算する（t は力学時）
  # =========================================================================
  def self.longitude_moon(t)
    # -----------------------------------------------------------------------
    #  摂動項の計算
    # -----------------------------------------------------------------------
    th =  0.0003 * Math.cos( $k * normalization_angle( 2322131.0  * t + 191.0  ) )
    th += 0.0003 * Math.cos( $k * normalization_angle(    4067.0  * t +  70.0  ) )
    th += 0.0003 * Math.cos( $k * normalization_angle(  549197.0  * t + 220.0  ) )
    th += 0.0003 * Math.cos( $k * normalization_angle( 1808933.0  * t +  58.0  ) )
    th += 0.0003 * Math.cos( $k * normalization_angle(  349472.0  * t + 337.0  ) )
    th += 0.0003 * Math.cos( $k * normalization_angle(  381404.0  * t + 354.0  ) )
    th += 0.0003 * Math.cos( $k * normalization_angle(  958465.0  * t + 340.0  ) )
    th += 0.0004 * Math.cos( $k * normalization_angle(   12006.0  * t + 187.0  ) )
    th += 0.0004 * Math.cos( $k * normalization_angle(   39871.0  * t + 223.0  ) )
    th += 0.0005 * Math.cos( $k * normalization_angle(  509131.0  * t + 242.0  ) )
    th += 0.0005 * Math.cos( $k * normalization_angle( 1745069.0  * t +  24.0  ) )
    th += 0.0005 * Math.cos( $k * normalization_angle( 1908795.0  * t +  90.0  ) )
    th += 0.0006 * Math.cos( $k * normalization_angle( 2258267.0  * t + 156.0  ) )
    th += 0.0006 * Math.cos( $k * normalization_angle(  111869.0  * t +  38.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(   27864.0  * t + 127.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(  485333.0  * t + 186.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(  405201.0  * t +  50.0  ) )
    th += 0.0007 * Math.cos( $k * normalization_angle(  790672.0  * t + 114.0  ) )
    th += 0.0008 * Math.cos( $k * normalization_angle( 1403732.0  * t +  98.0  ) )
    th += 0.0009 * Math.cos( $k * normalization_angle(  858602.0  * t + 129.0  ) )
    th += 0.0011 * Math.cos( $k * normalization_angle( 1920802.0  * t + 186.0  ) )
    th += 0.0012 * Math.cos( $k * normalization_angle( 1267871.0  * t + 249.0  ) )
    th += 0.0016 * Math.cos( $k * normalization_angle( 1856938.0  * t + 152.0  ) )
    th += 0.0018 * Math.cos( $k * normalization_angle(  401329.0  * t + 274.0  ) )
    th += 0.0021 * Math.cos( $k * normalization_angle(  341337.0  * t +  16.0  ) )
    th += 0.0021 * Math.cos( $k * normalization_angle(   71998.0  * t +  85.0  ) )
    th += 0.0021 * Math.cos( $k * normalization_angle(  990397.0  * t + 357.0  ) )
    th += 0.0022 * Math.cos( $k * normalization_angle(  818536.0  * t + 151.0  ) )
    th += 0.0023 * Math.cos( $k * normalization_angle(  922466.0  * t + 163.0  ) )
    th += 0.0024 * Math.cos( $k * normalization_angle(   99863.0  * t + 122.0  ) )
    th += 0.0026 * Math.cos( $k * normalization_angle( 1379739.0  * t +  17.0  ) )
    th += 0.0027 * Math.cos( $k * normalization_angle(  918399.0  * t + 182.0  ) )
    th += 0.0028 * Math.cos( $k * normalization_angle(    1934.0  * t + 145.0  ) )
    th += 0.0037 * Math.cos( $k * normalization_angle(  541062.0  * t + 259.0  ) )
    th += 0.0038 * Math.cos( $k * normalization_angle( 1781068.0  * t +  21.0  ) )
    th += 0.0040 * Math.cos( $k * normalization_angle(     133.0  * t +  29.0  ) )
    th += 0.0040 * Math.cos( $k * normalization_angle( 1844932.0  * t +  56.0  ) )
    th += 0.0040 * Math.cos( $k * normalization_angle( 1331734.0  * t + 283.0  ) )
    th += 0.0050 * Math.cos( $k * normalization_angle(  481266.0  * t + 205.0  ) )
    th += 0.0052 * Math.cos( $k * normalization_angle(   31932.0  * t + 107.0  ) )
    th += 0.0068 * Math.cos( $k * normalization_angle(  926533.0  * t + 323.0  ) )
    th += 0.0079 * Math.cos( $k * normalization_angle(  449334.0  * t + 188.0  ) )
    th += 0.0085 * Math.cos( $k * normalization_angle(  826671.0  * t + 111.0  ) )
    th += 0.0100 * Math.cos( $k * normalization_angle( 1431597.0  * t + 315.0  ) )
    th += 0.0107 * Math.cos( $k * normalization_angle( 1303870.0  * t + 246.0  ) )
    th += 0.0110 * Math.cos( $k * normalization_angle(  489205.0  * t + 142.0  ) )
    th += 0.0125 * Math.cos( $k * normalization_angle( 1443603.0  * t +  52.0  ) )
    th += 0.0154 * Math.cos( $k * normalization_angle(   75870.0  * t +  41.0  ) )
    th += 0.0304 * Math.cos( $k * normalization_angle(  513197.9  * t + 222.5  ) )
    th += 0.0347 * Math.cos( $k * normalization_angle(  445267.1  * t +  27.9  ) )
    th += 0.0409 * Math.cos( $k * normalization_angle(  441199.8  * t +  47.4  ) )
    th += 0.0458 * Math.cos( $k * normalization_angle(  854535.2  * t + 148.2  ) )
    th += 0.0533 * Math.cos( $k * normalization_angle( 1367733.1  * t + 280.7  ) )
    th += 0.0571 * Math.cos( $k * normalization_angle(  377336.3  * t +  13.2  ) )
    th += 0.0588 * Math.cos( $k * normalization_angle(   63863.5  * t + 124.2  ) )
    th += 0.1144 * Math.cos( $k * normalization_angle(  966404.0  * t + 276.5  ) )
    th += 0.1851 * Math.cos( $k * normalization_angle(   35999.05 * t +  87.53 ) )
    th += 0.2136 * Math.cos( $k * normalization_angle(  954397.74 * t + 179.93 ) )
    th += 0.6583 * Math.cos( $k * normalization_angle(  890534.22 * t + 145.7  ) )
    th += 1.2740 * Math.cos( $k * normalization_angle(  413335.35 * t +  10.74 ) )
    th += 6.2888 * Math.cos( $k * normalization_angle( 477198.868 * t +  44.963) )

    #-----------------------------------------------------------------------
    # 比例項の計算
    #-----------------------------------------------------------------------
    ang = normalization_angle( 481267.8809 * t )
    ang = normalization_angle( ang + 218.3162 )
    th  = normalization_angle( th + ang )

    th
  end

  #=========================================================================
  # 年月日、時分秒（世界時）からユリウス日（JD）を計算する
  #
  # ※ この関数では、グレゴリオ暦法による年月日から求めるものである。
  #    （ユリウス暦法による年月日から求める場合には使用できない。）
  #=========================================================================
  def self.ymdt2jd(year, mon, day, hour, min, sec)
    if mon < 3
      year -= 1
      mon  += 12
    end

    jd  = ( 365.25 * year ).to_i
    jd += ( year / 400.0 ).to_i
    jd -= ( year / 100.0 ).to_i
    jd += ( 30.59 * ( mon - 2 ) ).to_i
    jd += 1721088
    jd += day

    t  = sec / 3600.0
    t += min / 60.0
    t += hour
    t  = t / 24.0

    jd += t

    jd
  end

  #=========================================================================
  # ユリウス日（JD）から年月日、時分秒（世界時）を計算する
  #
  # 戻り値の配列TIME[]の内訳
  # TIME[0] ... 年  TIME[1] ... 月  TIME[2] ... 日
  # TIME[3] ... 時  TIME[4] ... 分  TIME[5] ... 秒
  #
  # ※ この関数で求めた年月日は、グレゴリオ暦法によって表されている。
  #
  #=========================================================================
  def self.jd2ymdt(jd)
    time = []

    x0 = ( jd + 68570.0  ).to_i
    x1 = ( x0 / 36524.25 ).to_i
    x2 = x0 - ( 36524.25 * x1 + 0.75 ).to_i
    x3 = ( ( x2+1 ) / 365.2425 ).to_i
    x4 = x2 - ( 365.25 * x3 ).to_i + 31.0
    x5 = ( x4.to_i / 30.59 ).to_i
    x6 = ( x5.to_i / 11.0 ).to_i

    time[2] = x4 - ( 30.59 * x5 ).to_i
    time[1] = x5 - 12 * x6 + 2
    time[0] = 100 * ( x1 - 49 ) + x3 + x6

    if time[1] == 2 && time[2] > 28
      if time[0] % 100 == 0 && time[0] % 400 == 0
        time[2] = 29
      elsif time[0] % 4 ==0
        time[2]=29
      else
        time[2]=28
      end
    end

    tm = 86400.0 * ( jd -  jd.to_i )
    time[3] = (  tm / 3600.0 ).to_i
    time[4] = ( (tm - 3600.0 * time[3] ) / 60.0 ).to_i
    time[5] = (  tm - 3600.0 * time[3] - 60 * time[4] ).to_i

    time
  end

  #=========================================================================
  # 六曜算出関数
  #
  # 引数　 .... 計算対象となる年月日　$year $mon $day
  # 戻り値 .... 六曜  (大安 赤口 先勝 友引 先負 仏滅)
  #=========================================================================
  def self.rokuyou(year, mon, day)
    rokuyou = %w(大安 赤口 先勝 友引 先負 仏滅)

    q_yaer, uruu, q_mon, q_day = calc_kyureki(year, mon, day)

    rokuyou[ (q_mon + q_day) % 6 ]
  end

  #=========================================================================
  # 今日が２４節気かどうか調べる
  #
  # 引数　 .... 計算対象となる年月日　$year $mon $day
  # 戻り値 .... ２４節気の名称
  #=========================================================================
  def self.sekki(year, mon, day)
    #-----------------------------------------------------------------------
    # ２４節気の定義
    #-----------------------------------------------------------------------
    sekki24 = %w(春分 清明 穀雨 立夏 小満 芒種 夏至 小暑 大暑 立秋 処暑 白露
                 秋分 寒露 霜降 立冬 小雪 大雪 冬至 小寒 大寒 立春 雨水 啓蟄)

    tm = ymdt2jd(year, mon, day, 0, 0, 0)

    #-----------------------------------------------------------------------
    # 時刻引数を分解する
    #-----------------------------------------------------------------------
    tm1  = tm.to_i
    tm2  = tm - tm1
    tm2 -= 9.0 / 24.0
    t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0

    #今日の太陽の黄経
    rm_sun_today = longitude_sun(t)

    tm += 1
    tm1  = tm.to_i
    tm2  = tm - tm1
    tm2 -= 9.0 / 24.0
    t = (tm2 + 0.5) / 36525.0 + (tm1 - 2451545.0) / 36525.0

    # 明日の太陽の黄経
    rm_sun_tommorow = longitude_sun(t)

    #
    rm_sun_today0    = 15.0 * (rm_sun_today    / 15.0).to_i
    rm_sun_tommorow0 = 15.0 * (rm_sun_tommorow / 15.0).to_i

    if rm_sun_today0 != rm_sun_tommorow0
      return sekki24[rm_sun_tommorow0 / 15]
    else
      return ''
    end
  end

  def self.shunbun(year)
    day = (20.8431 + 0.242194 * ( year - 1980 ) - ( year - 1980 ) / 4).to_i

    sreki = Srk.new
    sreki.year  = year
    sreki.month = 3
    sreki.day   = day
    sreki
  end

  def self.shuubun(year)
    day = (23.2488 + 0.242194 * ( year - 1980 ) - ( year - 1980 ) / 4).to_i

    sreki = Srk.new
    sreki.year  = year
    sreki.month = 9
    sreki.day   = day
    sreki
  end
end
