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)式の３次曲線で近似し、その係数から求めました。

これらの係数と、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になるわけで、最初のホフマイスターのデータ付近だとこの様に見えます。