「中心差にみる重修大明暦と授時暦の類似性」のためのPCA処理スクリプト

anomaly_lines.rb

# -*- coding: utf-8 -*-
=begin
  Copyright (C) 2015-2019 Takashi SUGA

  You may use and/or modify this file according to the license
  described in the LICENSE.txt file included in https://github.com/suchowan/when_exe.
=end

# 
# 中心差にみる重修大明暦と授時暦の類似性
# The Similarity of the equation of center between the Revised Great Enlightenment Astronomical System and the Season-Granting System
#

require 'pp'
require 'matrix'
require 'when_exe'

if self.class.const_defined?(:Encoding)
  Encoding.default_external = 'UTF-8'
  Encoding.default_internal = 'UTF-8'
end

#
# 中心差を解析するクラス
#
class Anomaly

  include Math
  include When

  class << self
    #
    # 中心差の解析
    #
    def analyse
      calendars = [Rintoku, Almagest, Daien, Senmyo, Juji, Aryabhata, Surya, Kitai, Fuifui, Jokyo, Kepler, Daimei].map { |klass|
        calendar = klass.new
        pp [ calendar.class::Name,                                                  # 暦法または理論の名前
            [calendar.coef[0...5].map {|c| '%6.4f' % c}, '%5.2f' % calendar.error], # 正弦級数の係数/度と最大誤差/秒
             calendar.max_anomaly.map {|c| '%5.3f' % c},                            # 中心差の最大値/度, 真近点角/度, 平均近点角/度, (/日)
             calendar.interval.map    {|c| '%6.3f' % c},                            # 隣り合う節気と中気の間の平均近点角の差/度
           # calendar.interval_by_day.map   {|c| '%7.4f' % c},
           # calendar.max_anomaly_local.map {|c| '%5.3f' % c},
           ]
        calendar
      }
      pca(calendars)
      pp @@index
      pp coefficient(@@coef)
    # pp curve(Kepler, [0.01671])
    # pp curve(Almagest, [1/24.0])
      pp curve(Kepler,   (0..50).to_a.map {|i| 0.0005*i})
      pp curve(Almagest, (0..50).to_a.map {|i| 0.0010*i})
    end

    private

    #
    # 主成分分析
    #
    def pca(calendars)
      @@coef = Matrix[*calendars.map {|c| c.coef[0..10]}]
      mean = (0..10).to_a.map {|c| @@coef.column(c).inject(:+) / calendars.size}
      cov  = Matrix[*(0..10).to_a.map {|c1|
               (0..10).to_a.map {|c2|
                 (0...calendars.size).to_a.map {|i|
                   (@@coef[i,c1] - mean[c1]) * (@@coef[i,c2] - mean[c2])
                 }.inject(:+) / (calendars.size-1)
               }
             }]
      @@ev, @@ed = cov.eigensystem
      @@index    = (0..10).to_a.map {|i| [i,@@ed[i,i]]}.sort_by {|x| -x[1].abs}[0..2]
    end

    #
    # 理論曲線
    #
    def curve(klass, param)
      coef = Matrix[*param.map {|p| klass.new(p).coef[0..10]}]
      coefficient(coef)
    end

    #
    # 係数
    #
    def coefficient(coef)
      prod = coef * @@ev
      map  = (0...coef.row_size).to_a.map {|i| [prod[i, @@index[0][0]], prod[i, @@index[1][0]], prod[i, @@index[2][0]]]}
    end
  end

  #
  # 正弦級数の係数 / 度
  #
  # @return [Array<Numeric>]
  #
  attr_reader :coef

  #
  # 正弦級数の最大残差 / 秒
  #
  # @return [Numeric]
  #
  attr_reader :error

  #
  # 平均近点角　→　真近点角
  #
  # @param [Numeric] d 平均近点角 / 度
  #
  # @return [Numeric] 真近点角 / 度
  #
  def mean_to_true(d)
    d + anomaly(d)
  end

  #
  # 真近点角　→　平均近点角
  #
  # @param [Numeric] d 真近点角 / 度
  #
  # @return [Numeric] 平均近点角 / 度
  #
  def true_to_mean(d)
    Ephemeris.root(d, d) {|x| mean_to_true(x)}
  end

  #
  # 中心差の最大値 / 度
  #
  # @return [Array<Numeric>] 中心差の最大値/度, 真近点角/度, 平均近点角/度
  #
  def max_anomaly
    d = Ephemeris.root(90, nil) {|x| anomaly(x)}
    [anomaly(d), mean_to_true(d), d]
  end

  #
  # 隣り合う節気と中気の間の平均近点角の差 / 度
  #
  # @return [Array<Numeric>]
  #
  def interval
    (1..12).to_a.map {|n| true_to_mean(15*n) - true_to_mean(15*(n-1))}
  end

  private

  #
  # 正弦級数の係数と最大残差
  #
  # @param [Array<Numeric>] x 中心差の表 / 度
  # @param [Integer]        n 記憶する正弦級数の係数の数
  #
  # @return [Array<Array<Numeric>, Numeric>] 正弦級数の係数/度と最大残差/秒
  #
  def sins(x,n=(x.size==12) ? 11: 30)
    m = ft(x)
    c = (0...n).to_a.map {|k| m[k,0]}
    [c, (1...x.size).to_a.map {|d|
      (x[d] - (1..n).to_a.inject(0) {|s,k| s += c[k-1] * sin(k*d*PI/x.size)}).abs
    }.max*3600]
  end

  #
  # 高速でないフーリエ変換
  #
  # @param [Array<Numeric>] x 中心差の表 / 度
  #
  # @return [Matrix] 正弦級数の係数
  #
  def ft(x)
    Matrix.rows((1...x.size).to_a.map {|d|
      (1...x.size).to_a.map {|n| sin(n*d*PI/x.size)}
    }).inverse * Matrix.columns([x[1..-1]])
  end
end

#
# ケプラー
#
class Kepler < Anomaly

  Name = 'ケプラー'

  #
  # オブジェクトの生成
  #
  # @param [Numeric] eccentricity 離心率
  #
  def initialize(eccentricity=0.01671)
    @e = eccentricity
    @c = sqrt(1-@e*@e)
    @coef, @error = sins((0...180).to_a.map {|d| anomaly(d)})
  end

  #
  # 中心差
  #
  # @param [Numeric] 平均近点角 / 度
  #
  # @return [Numeric] 中心差 / 度
  #
  def anomaly(d)
    _M = d * PI / 180                                 # 平均近点角 / radian
    _E = Ephemeris.root(_M, _M) {|x| x - @e * sin(x)} # 離心近点角 / radian  E - e sin E = M
    atan2(@c*sin(_E), cos(_E) - @e) *                 # 真近点角 / radian    r cos υ = a(cos E - e)
    180 / PI - d                                      # 中心差 / 度          r sin υ = a √(1-e^2) sin E
  end
end

#
# アルマゲスト
#
class Almagest < Anomaly

  Name = 'アルマゲスト'

  #
  # オブジェクトの生成
  #
  # @param [Numeric] eccentricity 周転円の大きさ
  #
  def initialize(eccentricity=1/24.0)
    @e = eccentricity
    @coef, @error = sins((0...180).to_a.map {|d| anomaly(d)})
  end

  #
  # 中心差
  #
  # @param [Numeric] 平均近点角 / 度
  #
  # @return [Numeric] 中心差 / 度
  #
  def anomaly(d)
    theta = d * PI / 180                              # 平均近点角 / radian
    atan2(sin(theta), cos(theta) - @e) *              # 真近点角 / radian
    180 / PI - d                                      # 中心差 / 度
  end
end

#
# インド
#
class Hindu < Anomaly

  #
  # オブジェクトの生成
  #
  # @param [Numeric] circumm マンダ周転円の大きさ / 度
  #
  def initialize(circumm=self.class::Circumm)
    @circumm = circumm.map {|c| c / 360}
    @coef, @error = sins((0...180).to_a.map {|d| anomaly(d)})
  end

  #
  # 中心差
  #
  # @param [Numeric] 平均近点角 / 度
  #
  # @return [Numeric] 中心差 / 度
  #
  def anomaly(d)
    w       = (d - 90).abs / 90.0                     # 大きいマンダ周転円の重み
    circumm = w * @circumm[0] + (1-w) * @circumm[1]   # マンダ周転円の半径
    theta   = d * PI / 180                            # 平均近点角 / radian
    asin(circumm * sin(theta)) * 180 / PI             # マンダ補正 / 度
  end
end

#
# アールヤバティーヤ
#
class Aryabhata < Hindu
  Name    = 'アールヤバティーヤ'
  Circumm = [13+30/60.0] * 2
end

#
# スールヤ・シッダーンタ
#
class Surya < Hindu
  Name    = 'スールヤ・シッダーンタ'
  Circumm = [14+00/60.0, 13+40/60.0]
end

#
# 立成　―　表引き
#
class Table < Anomaly

  #
  # オブジェクトの生成
  #
  def initialize(table=self.class::Diff)
    @coef, @error = sins(table)
  end

  #
  # 中心差 ― フーリエ級数の計算
  #
  # @param [Numeric] 近点角 / 度
  #
  # @return [Numeric] 中心差 / 度
  #
  def anomaly(d)
    (1.. @coef.size).to_a.inject(0) {|s,k| s + @coef[k-1] * sin(k*d*PI/180)}
  end
end

#
# 立成　―　真近点角→中心差の表引き
#
class TrueAngleTable < Table

  #
  # オブジェクトの生成
  #
  def initialize
    @table = Table.new(self.class::Diff)
    @coef, @error = sins((0...180).to_a.map {|d| anomaly(d)})
  end

  #
  # 平均近点角　→　真近点角
  #
  # @param [Numeric] d 平均近点角 / 度
  #
  # @return [Numeric] 真近点角 / 度
  #
  def mean_to_true(d)
    Ephemeris.root(d, d) {|x| true_to_mean(x)}
  end

  #
  # 真近点角　→　平均近点角
  #
  # @param [Numeric] d 真近点角 / 度
  #
  # @return [Numeric] 平均近点角 / 度
  #
  def true_to_mean(d)
    d - @table.anomaly(d)
  end

  #
  # 中心差
  #
  # @param [Numeric] 平均近点角 / 度
  #
  # @return [Numeric] 中心差 / 度
  #
  def anomaly(d)
    @table.anomaly(mean_to_true(d))
  end

  #
  # 隣り合う節気と中気の間の平均近点角の差 / 日
  #
  # @return [Array<Numeric>]
  #
  def interval_by_day
    interval.map {|d| d * self.class::YEAR / 360}
  end
end

#
# 招差法　―　多項式近似
#
class Poly < Anomaly

  #
  # オブジェクトの生成
  #
  def initialize
    @diff = self.class::Diff
    @coef, @error = sins((0...180).to_a.map {|d| anomaly(d)})
  end

  #
  # 中心差
  #
  # @param [Numeric] 平均近点角 / 度
  #
  # @return [Numeric] 中心差 / 度
  #
  def anomaly(d)
    day = @diff[:year] / 360 * d                      # 平均近点角 / radian
    if day < @diff[:threshold]                        # 近日点側
      poly = @diff[:near]
    else                                              # 遠日点側
      poly = @diff[:far]
      day  = @diff[:year] / 2 - day
    end
    (day*(poly[0] + day*(poly[1] + day*poly[2])) /    # 中心差 / 一億分の一日
    100000000) /                                      # 中心差 / 日
    @diff[:year] * 360                                # 中心差 / 度
  end

  #
  # 中心差の最大値 / 度
  #
  # @return [Array<Numeric>] 中心差の最大値/度, 真近点角/度, 平均近点角/度, 平均近点角/日
  #
  def max_anomaly
    list = super
    list + [list.last * @diff[:year] / 360]
  end

  #
  # 多項式が極値となる位置
  #
  # @return [Numeric] 位置 / 日
  #
  def max_anomaly_local
    [:near, :far].map { |side|
      poly = @diff[side]
      root = sqrt(poly[1]*poly[1]-3*poly[2]*poly[0])
      (-poly[1]-root) /(3*poly[2])
    }
  end
end

#
# 麟徳暦(皇極暦,儀鳳暦)
#
class Rintoku < TrueAngleTable
  Name = '麟徳暦'
  YEAR = 489428.0 / 1340
  Diff = [   0,   722,  1340,  1854,  2368,  2986,  3708,  2986,  2368,  1854,  1340,   722].map {|d|
    d / 489428.0 * 360
  }
end

#
# 大衍暦
#
class Daien < TrueAngleTable
  Name = '大衍暦'
  YEAR = 1110343.0 / 3040
  Diff = [   0,  2353,  4198,  5588,  6564,  7152,  7366,  7152,  6564,  5588,  4198,  2353].map {|d|
    d / 1110343.0 * 360
  }
end

#
# 宣明暦
#
class Senmyo < TrueAngleTable
  Name = '宣明暦'
  YEAR = 3068055.0 / 8400
  Diff = [   0,  6000, 11000, 15000, 18000, 19800, 20400, 19800, 18000, 15000, 11000,  6000].map {|d|
    d / 3068055.0 * 360
  }
end

#
# 重修大明暦
#
class Daimei < TrueAngleTable
  Name = '重修大明暦'
  YEAR = 365 + 1274 / 5230.0
  Diff = [   0,  7059, 12979, 17697, 21150, 23276, 24015, 23276, 21150, 17697, 12979,  7059].map {|d|
    d / (YEAR * 10000) * 360
  }
end

#
# 授時暦
#
class Juji < Poly
  Name = '授時暦'
  Diff = {
    :year      => 365.2425,
    :threshold =>  88.9092250,
    :near      => [5133200, -24600, -31],
    :far       => [4870600, -22100, -27]
  }

  #
  # 中心差の最大値 / 度　―　近日点側と遠日点側の境界で最大
  #
  # @return [Array<Numeric>] 中心差の最大値/度, 真近点角/度, 平均近点角/度, 平均近点角/日
  #
  def max_anomaly
    d = @diff[:threshold] / @diff[:year] * 360
    [anomaly(d), mean_to_true(d), d, @diff[:threshold]]
  end
end

#
# 貞享暦
#
class Jokyo < Poly
  Name = '貞享暦'
  Diff = {
    :year      => 365.241696,
    :threshold =>  89.2539200,
    :near      => [4360000, -20000, -34],
    :far       => [4119800, -17640, -31]
  }
end

#
# Chinese-Uighur Calendar
#
class Kitai < Poly
  Name = 'キタイ暦'
  Amp  = 20000.0 / 9 * (365.2436 / 29.5306 + 1)
  Diff = {
    :year      => 364.0000000,
    :threshold =>  90.0000000,
    :near      => [182*Amp, -Amp, 0],
    :far       => [182*Amp, -Amp, 0]
  }
end

#
# 回回暦
#
class Fuifui < Table
  Name = '回回暦'
  Diff = [
    [0, 0, 0], [0, 2,11], [0, 4,22], [0, 6,33], [0, 8,44],  [0,10,55], [0,13, 5], [0,15,15], [0,17,25], [0,19,35],
    [0,21,44], [0,23,53], [0,26, 1], [0,28, 8], [0,30,15],  [0,32,21], [0,34,27], [0,36,32], [0,38,36], [0,40,40],
    [0,42,43], [0,44,45], [0,46,46], [0,48,46], [0,50,45],  [0,52,43], [0,54,40], [0,56,36], [0,58,31], [1, 0,24], # [1, 2,16]

    [1, 2,16], [1, 4, 7], [1, 5,57], [1, 7,46], [1, 9,33],  [1,11,18], [1,13, 2], [1,14,45], [1,16,27], [1,18, 7],
    [1,19,45], [1,21,22], [1,22,57], [1,24,31], [1,26, 3],  [1,27,33], [1,29, 1], [1,30,28], [1,31,53], [1,33,16],
    [1,34,38], [1,35,56], [1,37,13], [1,38,29], [1,39,43],  [1,40,55], [1,42, 5], [1,43,13], [1,44,19], [1,45,23], # [1,46,25]

    [1,46,25], [1,47,25], [1,48,23], [1,49,19], [1,50,13],  [1,51, 4], [1,51,53], [1,52,40], [1,53,25], [1,54, 8],
    [1,54,49], [1,55,28], [1,56, 5], [1,56,39], [1,57,11],  [1,57,41], [1,58, 8], [1,58,33], [1,58,56], [1,59,17],
    [1,59,36], [1,59,53], [2, 0, 8], [2, 0,20], [2, 0,30],  [2, 0,38], [2, 0,43], [2, 0,46], [2, 0,47], [2, 0,46], # [2, 0,43]

    [2, 0,43], [2, 0,38], [2, 0,30], [2, 0,20], [2, 0, 8],  [1,59,54], [1,59,37], [1,59,18], [1,58,57], [1,58,34],
    [1,58, 9], [1,57,42], [1,57,13], [1,56,42], [1,56, 9],  [1,55,34], [1,54,56], [1,54,16], [1,53,34], [1,52,50],
    [1,52, 4], [1,51,17], [1,50,28], [1,49,37], [1,48,44],  [1,47,49], [1,46,52], [1,45,53], [1,44,52], [1,43,49], # [1,42,45]

    [1,42,45], [1,41,39], [1,40,31], [1,39,21], [1,38,10],  [1,36,57], [1,35,42], [1,34,26], [1,33, 8], [1,31,48],
    [1,30,27], [1,29, 4], [1,27,40], [1,26,14], [1,24,47],  [1,23,19], [1,21,49], [1,20,17], [1,18,44], [1,17,10],
    [1,15,35], [1,13,58], [1,12,20], [1,10,41], [1, 9, 1],  [1, 7,20], [1, 5,37], [1, 3,53], [1, 2, 8], [1, 0,22], # [0,58,35]

    [0,58,35], [0,56,47], [0,54,58], [0,53, 8], [0,51,18],  [0,49,27], [0,47,35], [0,45,42], [0,43,48], [0,41,54],
    [0,39,59], [0,38, 3], [0,36, 7], [0,34,10], [0,32,13],  [0,30,15], [0,28,15], [0,26,16], [0,24,16], [0,22,17],
    [0,20,17], [0,18,16], [0,16,15], [0,14,13], [0,12,12],  [0,10,10], [0, 8, 8], [0, 6, 6], [0, 4, 4], [0, 2, 2], # [0, 0, 0]
  ].map {|d,m,s| d + m/60.0 + s/3600.0}
end

#
# メイン
#
Anomaly.analyse