【地球時・暦表時プログラム】

atime src download
このソースには

astro.h ------------ 天文定数、構造体
atime.c,atime.h ----- 計算本体
dt_tbl.c,dt_tbl.h ----- ΔTテーブル（内部定数用）
main.c ------------ main
dAT.tbl ----------- ΔAT(TAI - UTC)テーブル
Makefile ---------- gcc用MakeFile

以上が含まれています。

コンパイルは、
gccではmakeでatimeを作成。
VisualStudioは「win32コンソールプロジェクト」で*.h, *.cを読み込みコンパイル。『ｆopenやgmtimeは安全ではない』ワーニングが出ますが、無視します。


【使い方】
ΔAT(TAI-UTC)は外部ファイルです。atimeと同じフォルダーに入れて下さい。フォルダーや名前を変えたい場合は、main.cにある
#define        REAPSEC_FILE    "dAT.tbl"
を変えてください。

起動パラメータなしの場合、パソコンのシステム日時(UTC)で計算

> atime dTT(TT-UTC) 68.184000 TT 2016/10/09 03:31:08.1840 UTC 2016/10/09 03:30:00.0000

年月日時分秒を指定する場合、初期値を01/01 00:00:00にしているので年以降は省ける

> atime 2001 dTT(TT-UTC) 64.184000 TT 2001/01/01 00:01:04.1840 UTC 2001/01/01 00:00:00.0000 > atime 2001 12 dTT(TT-UTC) 64.184000 TT 2001/12/01 00:01:04.1840 UTC 2001/12/01 00:00:00.0000 > atime 2001 12 31 23 dTT(TT-UTC) 64.184000 TT 2001/12/31 23:01:04.1840 UTC 2001/12/31 23:00:00.0000

【ΔTの計算部分】
本プログラムでは、ΔTの計算を

1972年～
地球時の計算のΔATテーブル

1860年～1971年
TIMESCALESのデータテーブル・補間

～1859年
NASA POLYNOMIAL EXPRESSIONS FOR DELTA T
の近似式

以上の３種類になります。振り分けは
atime.cのget_dt()関数で行います。

double get_dt(double dJD) {     /* cal. parameter table */     if (dJD >= 2441317.5){         /* (1972～) */         return get_dtt_1972(dJD);     }     if (dJD >= 2400410.5){         /* 1860～1971 */         return get_dt_table(dJD);     }     /* other year */     return get_dt_approx(dJD); }

それぞれの処理を説明します。


・1972年～の処理
ここは60sの問題があるので、基本はTAIの計算プログラムと同じです。

get_dtt_1972() ----- UTCからΔTT
get_dtt_tt_1972() --- TTからΔTT

get_dtt_tt_1972()は、60sフラグ付きです。

外部ΔAT読み込みは、言語により読み方が違うためmain.cに置きました。
read_leepsec_table()内において

ptbl->dTime = (sec + DTAI_TT) / DAYSEC;

この行でΔAT（TAI-UTC）をΔTT（TAI-UTC+32.184s）にします。

この期間の表示はΔTと違うため

$ ./atime 1976 12 31 23 59 45 dTT(TT-UTC) 47.184000 TT 1977/01/01 00:00:32.1840 UTC 1976/12/31 23:59:45.0000

dTT(TT-UTC)という表記にしています。


・1860年～1971年
この期間の処理は、結局



この式にしました。ΔTはTIMESCALESのデータを使います。

ΔTデータはdt_tbl.cにある
REAPSEC tTimesDT[]
１年毎のΔTテーブルです。
最後の1972年
    {2441317.5, 42.184}
正しいΔTは 42.23ですが、1972年の42.184にあわせるためです。

このテーブルからΔTを計算するのは、atime.cにあるget_dt_table()関数で一次（直線）補間でΔTを求めています。

なお、関数の最初に行っている処理

/* about jd to year */     i = (int)((dJD - 2400410.5) / GREGOR_YEARDAY) - 1;     if (i < 0){         i = 0;     }

これは、該当テーブルが何番目になるか目途をつける処理です。テーブルの頭からなめるよりは早いでしょう。


・～1859年
ここは
NASA POLYNOMIAL EXPRESSIONS FOR DELTA T
の近似計算でΔTを計算します。

関数はatime.cのget_dt_approx()関数で、計算パラメータは、dt_tbl.cの
dDT_1941[]～dDT_m500[]です。
各式の係数数が違うので、
double *pDT[]のポインターのポインターでテーブルの先頭番地を集めています。

プログラムでは、NASAオリジナル（年・月）と違ってユリウス日（日時）で計算するようにしました。当然、ΔT値は変わります。

この図は、NASA計算式を年・月（NASA指定のまま）で計算した物とユリウス日（日時）で計算した物の差です。
-500年で差の符号が反転してますが、それでも±1.0s以内なので良しとしました。1860年以降大きく暴れているのは、NASA近似式以外の計算のためです。

で計算した期間は

$ ./atime 1958 dT(ET-UT) 32.180000 ET 1958/01/01 00:00:32.1800 UT 1958/01/01 00:00:00.0000

dT(ET-UT),ET,UT表記にしています。


【UTC⇔TTの計算部分】
全ての年で共通した関数にしています。

・UTC→TT変換
utc_tt()関数を使います。UTCは60s問題があるので、年月日時分秒の構造体渡しです。

if(dJD >= 2441317.5){         ddAT = get_dtt_1972(dJD);     }else{         ddAT = get_dt(dJD + dTimeDay);     }

この部分で1972年とそれ以外を振り分け、UTCが1972年なら（60s対策で）日付だけ渡し、1972年以前なら、日付＋時刻渡しで計算です。

・TT→UTC変換
tt_utc（）関数です。こちらも1972年で処理を分けます。

if (dJD_tt >= (2441317.5 + (10 + DTAI_TT) / DAYSEC)){         ddAT = get_dtt_tt_1972(dJD_tt, &n60sec);         if (n60sec != 0){             dJD_tt -= (1.0/DAYSEC);         }     }else{         /* ～1971 */         ddAT = 0;         do{             dJD = ddAT;             ddAT = get_dt(dJD_tt - dJD);         }while(fabs(dJD - ddAT) > (0.1 / DAYSEC));     }

1972年以後は一発で計算できますが、60s問題対策処理があります。

1972年以前は、日時でΔT値が変わるため、一定のΔT値になるまで繰返しが必要です。プログラムでは0.1s以下になるまでやってます。


【使用例】
60s確認のため、2015/6/30～7/1

$ ./atime 2015 6 30 23 59 59.987 dTT(TT-UTC) 67.184000 TT 2015/07/01 00:01:07.1710 UTC 2015/06/30 23:59:59.9870 $ ./atime 2015 6 30 23 59 60.987 dTT(TT-UTC) 67.184000 TT 2015/07/01 00:01:08.1710 UTC 2015/06/30 23:59:60.9870 $ ./atime 2015 7 1 0 0 0.987 dTT(TT-UTC) 68.184000 TT 2015/07/01 00:01:09.1710 UTC 2015/07/01 00:00:00.9870

TAIとTTの定義
1977/01/01.0_TAI = 1977/01/01 0:0:32.184_TT
1976/12/31_UTCのΔATは15sなので、1977/01/01.0_TAI = 1976/12/31 23:59:45_UTC

$ ./atime 1976 12 31 23 59 45 dTT(TT-UTC) 47.184000 TT 1977/01/01 00:00:32.1840 UTC 1976/12/31 23:59:45.0000

TAIとUT2の定義
1958/1/1.0_TAI = 1958/1/1.0_UT2
その時のΔTは32.18s

$ ./atime 1958 dT(ET-UT) 32.180000 ET 1958/01/01 00:00:32.1800 UT 1958/01/01 00:00:00.0000

暦表時の計算に戻る　　次は地心座標時の計算