O - C
昨日、
T(N) = a*N^3 + b*N^2 + c*N + d (1)
の係数を、データからRのnls関数で直接求めようとしたら、エラーで出来なかったと書きました。
そこで、P0を元期での周期、 E0を元期のHJDとして、この式を少し変形します。
T(N) = a*N^3 + b*N^2 + (c' + P0) * N + (d'+ E0)
= a*N^3 + b*N^2 + c' *N + P0 * N + d'+ E0
それから、P0*N + E0を左辺に移行します。
T(N) - (P0 * N + E0) = a*N^3 + b*N^2 + c' * N + d' (2)
すると左辺は有名なO - C になります。
これの右辺の係数はnls関数で出せました。具体的にはどの係数も大変に小さな値です。
a = 1.875027e-12
b = 1.940218e-08
c'= -8.900277e-06
d'= -7.395936e-04
最初の式(1)の係数は、この様にまずO - C のグラフを(2)式の3次曲線で近似し、その係数から求めました。
これらの係数と、P0 = 2.40599, E0 = 2456788.050、上記の(1)の式によってT(N)が計算できまるようになりました。(煩雑なので、得られた係数値を入れた式はいま書きません。)
T(N) = a*N^3 + b*N^2 + (c' + 2.40599) * N + (d'+ 2456788.050)
昨日書いた様にP(N)はこの式を微分して
P(N) = 3*a*N^2 + 2*b*N + (c' + P0)
となります。
N=0(ケプラーの観測の真ん中あたり)では
P(0) = c' + P0 = 2.40599 - 0.0000089 日となります。
周期の変化率PdotはP(N)をさらに微分して
Pdot = 6*a*N +2*b
となります。これらをグラフにしたものは先日載せました。
ところで
T(N) = N * P0 + E0
で描いた直線と観測で得られた曲線との差がO-Cになるわけで、最初のホフマイスターのデータ付近だとこの様に見えます。