module EphBpn::Compute
Public Instance Methods
comp_gamma_bp()
click to toggle source
comp_gamma_p()
click to toggle source
comp_phi_bp()
click to toggle source
comp_phi_p()
click to toggle source
comp_psi_bp()
click to toggle source
comp_psi_p()
click to toggle source
comp_r_bias()
click to toggle source
¶ ↑
Bias 変換行列(一般的な理論) * 赤道座標(J2000.0)の極は ICRS の極に対して12時(x軸のマイナス側)の方 向へ 17.3±0.2 mas、18時(y軸のマイナス側)の方向へ 5.1±0.2 mas ズレ ているので、変換する。 さらに、平均分点への変換はICRSでの赤経を78±10 mas、天の極を中心に回 転させる。 18時の方向の変換はx軸を-5.1mas回転 + 1 0 0 + R1(θ ) = | 0 cosθ sinθ | + 0 -sinθ cosθ + 12時の方向の変換はy軸を-17.3mas回転 + cosθ 0 -sinθ + R2(θ ) = | 0 1 0 | + sinθ 0 cosθ + 天の極を中心に78.0mas回転 + cosθ sinθ 0 + R3(θ ) = | -sinθ cosθ 0 | + 0 0 1 + @param: <none> @return: r (回転行列)
¶ ↑
# File lib/eph_bpn/compute.rb, line 110 def comp_r_bias r = r_x( -5.1 * Const::MAS2R) r = r_y(-17.3 * Const::MAS2R, r) r = r_z( 78.0 * Const::MAS2R, r) return r rescue => e raise end
comp_r_bias_prec()
click to toggle source
¶ ↑
Bias + Precession 変換行列 * IAU 2006 (Fukushima-Williams 4-angle formulation) 理論 @param: <none> @return: r (変換行列)
¶ ↑
# File lib/eph_bpn/compute.rb, line 127 def comp_r_bias_prec gamma = comp_gamma_bp phi = comp_phi_bp psi = comp_psi_bp r = r_z(gamma) r = r_x(phi, r) r = r_z(-psi, r) r = r_x(-@eps, r) return r rescue => e raise end
comp_r_bias_prec_nut()
click to toggle source
¶ ↑
Bias + Precession + Nutation 変換行列 * IAU 2006 (Fukushima-Williams 4-angle formulation) 理論 @param: <none> @return: r (変換行列)
¶ ↑
# File lib/eph_bpn/compute.rb, line 148 def comp_r_bias_prec_nut gamma = comp_gamma_bp phi = comp_phi_bp psi = comp_psi_bp fj2 = -2.7774e-6 * @jc dpsi_ls, deps_ls = compute_lunisolar dpsi_pl, deps_pl = compute_planetary dpsi, deps = dpsi_ls + dpsi_pl, deps_ls + deps_pl dpsi += dpsi * (0.4697e-6 + fj2) deps += deps * fj2 r = r_z(gamma) r = r_x(phi, r) r = r_z(-psi - dpsi, r) r = r_x(-@eps - deps, r) return r rescue => e raise end
comp_r_nut()
click to toggle source
¶ ↑
nutation(章動)変換行列 * IAU 2000A nutation with adjustments to match the IAU 2006 precession. @param: <none> @return: r (変換行列)
¶ ↑
# File lib/eph_bpn/compute.rb, line 240 def comp_r_nut fj2 = -2.7774e-6 * @jc dpsi_ls, deps_ls = compute_lunisolar dpsi_pl, deps_pl = compute_planetary dpsi, deps = dpsi_ls + dpsi_pl, deps_ls + deps_pl dpsi += dpsi * (0.4697e-6 + fj2) deps += deps * fj2 r = r_x(@eps) r = r_z(-dpsi, r) r = r_x(-@eps-deps, r) return r rescue => e raise end
comp_r_prec()
click to toggle source
¶ ↑
precession(歳差)変換行列(J2000.0 用) * 歳差の変換行列 P(ε , ψ , φ , γ ) = R1(-ε ) * R3(-ψ ) * R1(φ ) * R3(γ ) 但し、R1, R2, R3 は x, y, z 軸の回転。 + 1 0 0 + + cosθ sinθ 0 + R1(θ ) = | 0 cosθ sinθ | , R3(θ ) = | -sinθ cosθ 0 | + 0 -sinθ cosθ + + 0 0 1 + + P_11 P_12 P_13 + P(ε, ψ, φ, γ) = | P_21 P_22 P_23 | とすると、 + P_31 P_32 P_33 + P_11 = cosψ cosγ + sinψ cosφ sinγ P_12 = cosψ sinγ - sinψ cosφ ̄cosγ P_13 = -sinψ sinφ P_21 = cosε sinψ cosγ - (cosε cosψ cosφ + sinε sinφ )sinγ P_22 = cosε sinψ cosγ + (cosε cosψ cosφ + sinε sinφ )cosγ P_23 = cosε cosψ sinφ - sinε cosφ P_31 = sinε sinψ cosγ - (sinε cosψ cosφ - cosε sinφ)sinγ P_32 = sinε sinψ cosγ + (sinε cosψ cosφ - cosε sinφ)cosγ P_33 = sinε cosψ sinφ + cosε cosφ @param: <none> @return: r (変換行列)
¶ ↑
# File lib/eph_bpn/compute.rb, line 192 def comp_r_prec gamma = comp_gamma_p phi = comp_phi_p psi = comp_psi_p r = r_z(gamma) r = r_x(phi, r) r = r_z(-psi, r) r = r_x(-@eps, r) return r rescue => e raise end
comp_r_prec_nut()
click to toggle source
¶ ↑
Precession + Nutation 変換行列 * IAU 2000A nutation with adjustments to match the IAU 2006 precession. @param: <none> @return: r (変換行列)
¶ ↑
# File lib/eph_bpn/compute.rb, line 213 def comp_r_prec_nut gamma = comp_gamma_p phi = comp_phi_p psi = comp_psi_p fj2 = -2.7774e-6 * @jc dpsi_ls, deps_ls = compute_lunisolar dpsi_pl, deps_pl = compute_planetary dpsi, deps = dpsi_ls + dpsi_pl, deps_ls + deps_pl dpsi += dpsi * (0.4697e-6 + fj2) deps += deps * fj2 r = r_z(gamma) r = r_x(phi, r) r = r_z(-psi - dpsi, r) r = r_x(-@eps - deps, r) return r rescue => e raise end
compute_d_mhb2000(jc)
click to toggle source
¶ ↑
Mean elongation of the Moon from the Sun (MHB2000) @param: jc (ユリウス世紀) @return: longitude (rad)
¶ ↑
# File lib/eph_bpn/compute.rb, line 488 def compute_d_mhb2000(jc) return (( 1072260.70369 + \ (1602961601.2090 + \ ( -6.3706 + \ ( 0.006593 + \ ( -0.00003169) \ * jc) * jc) * jc) * jc) % Const::TURNAS) * Const::AS2R rescue => e raise end
compute_d_mhb2000_2(jc)
click to toggle source
compute_f_iers2003(jc)
click to toggle source
¶ ↑
Mean longitude of the Moon minus that of the ascending node (IERS 2003) @param: jc (ユリウス世紀) @return: longitude (rad)
¶ ↑
# File lib/eph_bpn/compute.rb, line 471 def compute_f_iers2003(jc) return (( 335779.526232 + \ (1739527262.8478 + \ ( -12.7512 + \ ( -0.001037 + \ ( 0.00000417) \ * jc) * jc) * jc) * jc) % Const::TURNAS) * Const::AS2R rescue => e raise end
compute_f_mhb2000(jc)
click to toggle source
compute_l_iers2003(jc)
click to toggle source
¶ ↑
Mean anomaly of the Moon (IERS 2003) @param: jc (ユリウス世紀) @return: longitude (rad)
¶ ↑
# File lib/eph_bpn/compute.rb, line 437 def compute_l_iers2003(jc) return (( 485868.249036 + \ (1717915923.2178 + \ ( 31.8792 + \ ( 0.051635 + \ ( -0.00024470) \ * jc) * jc) * jc) * jc) % Const::TURNAS) * Const::AS2R rescue => e raise end
compute_l_mhb2000(jc)
click to toggle source
compute_lea_iers2003(jc)
click to toggle source
compute_lju_iers2003(jc)
click to toggle source
compute_lma_iers2003(jc)
click to toggle source
compute_lme_iers2003(jc)
click to toggle source
compute_lne_mhb2000(jc)
click to toggle source
compute_lp_mhb2000(jc)
click to toggle source
¶ ↑
Mean anomaly of the Sun (MHB2000) @param: jc (ユリウス世紀) @return: longitude (rad)
¶ ↑
# File lib/eph_bpn/compute.rb, line 454 def compute_lp_mhb2000(jc) return (( 1287104.79305 + \ (129596581.0481 + \ ( -0.5532 + \ ( 0.000136 + \ ( -0.00001149) \ * jc) * jc) * jc) * jc) % Const::TURNAS) * Const::AS2R rescue => e raise end
compute_lsa_iers2003(jc)
click to toggle source
compute_lunisolar()
click to toggle source
¶ ↑
日月章動(luni-solar nutation)の計算 @param: <none> @return: [dpsi_ls, deps_ls]
¶ ↑
# File lib/eph_bpn/compute.rb, line 367 def compute_lunisolar dp, de = 0.0, 0.0 begin l = compute_l_iers2003(@jc) lp = compute_lp_mhb2000(@jc) f = compute_f_iers2003(@jc) d = compute_d_mhb2000(@jc) om = compute_om_iers2003(@jc) Const::NUT_LS.map do |x| x[0, 5] + x[5..-1].map { |y| y.to_s.sub(/\./, "").to_i } end.reverse.each do |x| arg = (x[0] * l + x[1] * lp + x[2] * f + x[3] * d + x[4] * om) % Const::PI2 sarg, carg = Math.sin(arg), Math.cos(arg) dp += (x[5] + x[6] * @jc) * sarg + x[ 7] * carg de += (x[8] + x[9] * @jc) * carg + x[10] * sarg end return [dp * Const::U2R, de * Const::U2R] rescue => e raise end end
compute_lur_iers2003(jc)
click to toggle source
compute_lve_iers2003(jc)
click to toggle source
compute_obliquity(jc)
click to toggle source
¶ ↑
黄道傾斜角計算 * 黄道傾斜角 ε (単位: rad)を計算する。 以下の計算式により求める。 ε = 84381.406 - 46.836769 * T - 0.0001831 T^2 + 0.00200340 T^3 - 5.76 * 10^(-7) * T^4 - 4.34 * 10^(-8) * T^5 ここで、 T は J2000.0 からの経過時間を 36525 日単位で表したユリウス 世紀数で、 T = (JD - 2451545) / 36525 である。 @param: jc (ユリウス世紀) @return: eps (平均黄道傾斜角)
¶ ↑
# File lib/eph_bpn/compute.rb, line 74 def compute_obliquity(jc) return (84381.406 + \ ( -46.836769 + \ ( -0.0001831 + \ ( 0.00200340 + \ ( -5.76 * 10e-7 + \ ( -4.34 * 10e-8 ) \ * jc) * jc) * jc) * jc) * jc) * Const::AS2R rescue => e raise end
compute_om_iers2003(jc)
click to toggle source
¶ ↑
Mean longitude of the ascending node of the Moon (IERS 2003) @param: jc (ユリウス世紀) @return: longitude (rad)
¶ ↑
# File lib/eph_bpn/compute.rb, line 505 def compute_om_iers2003(jc) return (( 450160.398036 + \ (-6962890.5431 + \ ( 7.4722 + \ ( 0.007702 + \ ( -0.00005939) \ * jc) * jc) * jc) * jc) % Const::TURNAS) * Const::AS2R rescue => e raise end
compute_om_mhb2000(jc)
click to toggle source
compute_pa_iers2003(jc)
click to toggle source
compute_planetary()
click to toggle source
¶ ↑
惑星章動(planetary nutation) @param: <none> @return: [dpsi_pl, deps_pl]
¶ ↑
# File lib/eph_bpn/compute.rb, line 397 def compute_planetary dp, de = 0.0, 0.0 begin l = compute_l_mhb2000(@jc) f = compute_f_mhb2000(@jc) d = compute_d_mhb2000_2(@jc) om = compute_om_mhb2000(@jc) pa = compute_pa_iers2003(@jc) lme = compute_lme_iers2003(@jc) lve = compute_lve_iers2003(@jc) lea = compute_lea_iers2003(@jc) lma = compute_lma_iers2003(@jc) lju = compute_lju_iers2003(@jc) lsa = compute_lsa_iers2003(@jc) lur = compute_lur_iers2003(@jc) lne = compute_lne_mhb2000(@jc) Const::NUT_PL.map do |x| x[0, 14] + x[14..-1].map { |y| y.to_s.sub(/\./, "").to_i } end.reverse.each do |x| arg = (x[ 0] * l + x[ 2] * f + x[ 3] * d + x[ 4] * om + x[ 5] * lme + x[ 6] * lve + x[ 7] * lea + x[ 8] * lma + x[ 9] * lju + x[10] * lsa + x[11] * lur + x[12] * lne + x[13] * pa) % Const::PI2 sarg, carg = Math.sin(arg), Math.cos(arg) dp += x[14] * sarg + x[15] * carg de += x[16] * sarg + x[17] * carg end return [dp * Const::U2R, de * Const::U2R] rescue => e raise end end
gc2jd(t)
click to toggle source
¶ ↑
年月日(グレゴリオ暦)からユリウス日(JD)を計算する * フリーゲルの公式を使用する JD = int(365.25 * year) + int(year / 400) - int(year / 100) + int(30.59 (month - 2)) + day + 1721088 * 上記の int(x) は厳密には、 x を超えない最大の整数 * 「ユリウス日」でなく「準ユリウス日」を求めるなら、 `+ 1721088` を `- 678912` とする。 @param: t (Time Object) @return: jd (ユリウス日)
¶ ↑
# File lib/eph_bpn/compute.rb, line 22 def gc2jd(t) year, month, day = t.year, t.month, t.day hour, min, sec = t.hour, t.min, t.sec begin # 1月,2月は前年の13月,14月とする if month < 3 year -= 1 month += 12 end # 日付(整数)部分計算 jd = (365.25 * year).floor \ + (year / 400.0).floor \ - (year / 100.0).floor \ + (30.59 * (month - 2)).floor \ + day \ + 1721088.5 # 時間(小数)部分計算 t = (sec / 3600.0 + min / 60.0 + hour) / 24.0 return jd + t rescue => e raise end end
jd2jc(jd)
click to toggle source
r_x(phi, r = Const::R_UNIT)
click to toggle source
¶ ↑
回転行列(x 軸中心) @param: phi (回転量(rad)) @param: r (回転行列) @return: r_a (回転後)
¶ ↑
# File lib/eph_bpn/compute.rb, line 679 def r_x(phi, r = Const::R_UNIT) r_a = Array.new begin s, c = Math.sin(phi), Math.cos(phi) a10 = c * r[1][0] + s * r[2][0] a11 = c * r[1][1] + s * r[2][1] a12 = c * r[1][2] + s * r[2][2] a20 = - s * r[1][0] + c * r[2][0] a21 = - s * r[1][1] + c * r[2][1] a22 = - s * r[1][2] + c * r[2][2] r_a << r[0] r_a << [a10, a11, a12] r_a << [a20, a21, a22] return r_a rescue => e raise end end
r_y(theta, r = Const::R_UNIT)
click to toggle source
¶ ↑
回転行列(y 軸中心) @param: theta (回転量(rad)) @param: r (回転行列) @return: r_a (回転後)
¶ ↑
# File lib/eph_bpn/compute.rb, line 706 def r_y(theta, r = Const::R_UNIT) r_a = Array.new begin s, c = Math.sin(theta), Math.cos(theta) a00 = c * r[0][0] - s * r[2][0] a01 = c * r[0][1] - s * r[2][1] a02 = c * r[0][2] - s * r[2][2] a20 = s * r[0][0] + c * r[2][0] a21 = s * r[0][1] + c * r[2][1] a22 = s * r[0][2] + c * r[2][2] r_a << [a00, a01, a02] r_a << r[1] r_a << [a20, a21, a22] return r_a rescue => e raise end end
r_z(psi, r = Const::R_UNIT)
click to toggle source
¶ ↑
回転行列(z 軸中心) @param: psi (回転量(rad)) @param: r (回転行列) @return: r_a (回転後)
¶ ↑
# File lib/eph_bpn/compute.rb, line 733 def r_z(psi, r = Const::R_UNIT) r_a = Array.new begin s, c = Math.sin(psi), Math.cos(psi) a00 = c * r[0][0] + s * r[1][0]; a01 = c * r[0][1] + s * r[1][1]; a02 = c * r[0][2] + s * r[1][2]; a10 = - s * r[0][0] + c * r[1][0]; a11 = - s * r[0][1] + c * r[1][1]; a12 = - s * r[0][2] + c * r[1][2]; r_a << [a00, a01, a02] r_a << [a10, a11, a12] r_a << r[2] return r_a rescue => e raise end end