地心座標時の計算

2016/10/19 改訂（公式・プログラム変更）

【固有時（Proper Time）】
地心座標時（TCG:Geocentric Coordinate Time）の計算をする前に、

相対論の考えでは速度が速い程、重力の影響が大きい程時間の進みが遅くなるらしい。一般に空間は（x, y, z）の三次元であるが、相対論では（t, x, y, z）の四次元で表す。これを時空（時空間）と呼ぶ。※私には理解できないので、ローレンツ変換等で調べてください。
その時空の時間を固有時と呼び、この固有時を用いれば、重力で歪んだ空間でも物理法則が成り立つ（らしい）。

地表に住む我々の固有時が地球時（TT）であり、（地表面は最大の地球重力下）地球中心座標の固有時が地心座標時（TCG）である。（地球中心では重力＝０となる）。
TCGはTTから地球重力分のみ引かれているためTCGの方がTTよりも速く進む。※自転の影響分は引かないのか？

【地心座標時（TCG : Geocentric Coordinate Time）】
地心座標時（TCG)は、地球重心を中心にする運動、人工衛星や歳差、章動、月に用いる。
具体的には、人工衛星の位置観測はTTで行い、その位置データから軌道計算はTCGに変換。その逆の人工衛星の位置計算はTCGで、観測用の予報位置はTTに変換。こんな感じになると思われます。

【TT⇔TCGの計算】
・TT⇒TCG

　・・・（式３－１－１）

出典Explanatory Supplement to the Astronomical Almanac 3rd Edition

この（式３－１－１）は、多くの文献で紹介されてますが近似式です。

　・・・（式３－１－２）

出典IERS Techical Note No.36

この（式３－１－２）が正式の様です。LGの定義は

１から比率を引いた物の様です。

この違いを計算してみました。

LG = 6.969290134000E-10 = 60.2146667578nsec
LG/(1-LG) = 6.969290138857E-10 = 60.2146667997nsec

（式３－１－１）と（式３－１－２）の差が秒に出てくるには億年先です。（式３－１－１）の近似式でも何ら問題はありません。
この式を見れば、TTとTCGは比例であり、比率がLGの値です。１日当り　0.000 060s = 60nsecTCGが速く進みます。ほんの僅かの差ですが５０年間で約１秒になります。

・TCG⇒TT
（式３－１－１）近似式で求めたTCGの場合

　・・・（式３－２－１）

（式３－１－２）で求めたTCGの場合

　・・・（式３－２－２）

となります。どちらの式をつかっても、その逆変換でLG/(1±LG)という計算を行う必要があります。

【TT⇔TCG変換プログラム】

program src download

プログラムは式の通りですが、doubleでは精度が足らない！LGは、小数以下１９桁になります。これにユリウス日の７桁で合計２６桁以上必要です。
２６桁以上の浮動小数点は、変数の有効桁で書いた16バイト以上のlong doubleになり、これが使えるのは64bit版Linuxのみ。前回は long double版を作りましたが、Windowsでも動かしたいので・・・
対策としてユリウス日はdouble変数２個使う方法にしました。

定義は astro.hの
typedef struct tag_ljd{
double dDate;
double dTime;
} LONGJD;

これを見たらわかる通り、日付用と時間用の２つの変数でユリウス日を表現しています。日付用は０時のユリウス日で常に xxx.5になる。時間用は日の単位で小数のみ。これでwindows上でもlong doubleの精度が保てます。ついでに、longjdの採用でユリウス日も±億年が使えるようにしました。（年をdouble扱いにした）

TCGの計算部分は日付と時間を分離して計算します。

double tt_tcg(LONGJD *pdJDtt, LONGJD *pdJDtcg)
{
    double dFrac;
    double dInt;
    double dTemp;

    /* cal. dTCG */
    dTemp = CONST_LGG * (pdJDtt->dDate - JD1977);
    dFrac = CONST_LGG * (pdJDtt->dTime - DTAI_TT / DAYSEC);

    dInt= (int)dTemp;
    dFrac += (dTemp - dInt);

    /* Calc. TCG */
    pdJDtcg->dDate = pdJDtt->dDate + dInt;
    pdJDtcg->dTime = pdJDtt->dTime + dFrac;
    ljd_normalize(pdJDtcg);

    return ((dFrac + dInt) * DAYSEC);
}

CONST_LGGが LG/(1-LG)です。

dInt= (int)dTemp;
ここは小数部も符号を残したいので、floor()ではなく切り捨ての(int)です。modf()関数でも構いません。

デバッグ込みでプログラムを動かしてみます。今回はWinXP上で動かしました。
パラメータなし＝PCのシステムタイム（UTC）

>atime
UTC 2016/10/27 17:33:00.000000
JD  2457689.231250

dTT(TT-UTC) 68.184000
TT  2016/10/27 17:34:08.184000
dTCG(TCG-TT) 0.875806
TCG 2016/10/27 17:34:09.059806

TT  2016/10/27 17:34:08.184000
UTC 2016/10/27 17:33:00.000000

dTCGが０になるのは1977/1/1.0 TAIですが、LGの単位から考えてTAIとUTCの差は無視。

>atime 1977
UTC 1977/01/01 00:00:00.000000
JD  2443144.500000

dTT(TT-UTC) 48.184000
TT  1977/01/01 00:00:48.184000
dTCG(TCG-TT) 0.000000
TCG 1977/01/01 00:00:48.184000

TT  1977/01/01 00:00:48.184000
UTC 1977/01/01 00:00:00.000000

dTCGにLGの値が出てくるのは、1977年から100億日後のJDで10002443144.5。
今回からJD入力可能にしたのでやってみます。

>atime jd 10002443144.5
UTC 27381047/01/27 00:00:00.000000
JD  10002443144.500000

dTT(TT-UTC) 69.184000
TT  27381047/01/27 00:01:09.184000
dTCG(TCG-TT) 602146.667997
TCG 27381047/02/02 23:16:55.851997

TT  27381047/01/27 00:01:09.184000
UTC 27381047/01/27 00:00:00.000000

dTCGの値が LG/(1-LG)の数字になりました。（式３－１－１）の場合は602146.667578です。

ｄTCGがマイナスのLG値になる日を試します。
1977年から100億日を引けば良い。2443144.5 - 10000000000 = -9997556855.5ですが、dTが74322.666517sと大きいので (74322.666517 - 32.184)/86400 引きます。
結果 -9997556856.359844。これでやってみます。

>atime jd -9997556856.359844
UT  -27376531/02/03 03:21:49.460449
JD  -9997556856.359844

dT(ET-UT) 74322.666517
ET  -27376531/02/04 00:00:32.126966
dTCG(TCG-TT) -602146.667997
TCG -27376531/01/28 00:44:45.458969

TT  -27376531/02/04 00:00:32.126966
UT  -27376531/02/03 03:21:49.460449

dTCGの値が -LG/(1-LG)の数字になりました。※このデバッグのためにJD入力を作りました。

TE,TTプログラムに戻る　　次は太陽系座標時