【INPOPによるTCBの計算(Linux only)】
INPOPという惑星座標のファイルがあり、この座標ファイルの中にTT-TDBかTCG-TCBの値が収められています。このINPOPを使って、TCB(太陽系座標時）を求めてみます。
今回は、INPOPから提供されているライブラリーを使う関係でLinux版のみとなります。

・INPOPの場所
INPOPは、以下のサイトからダウンロードする事ができる。

https://www.imcce.fr/fr/presentation/equipes/ASD/inpop/

2018年1月には、INPOP17aのバージョンになっていました。
INPOP17aの本体は、「DOWNLOAD 〜」 を押せば飛んでいきます。

https://www.imcce.fr/recherche/equipes/asd/inpop/download17a


・INPOP本体
INPOPファイルには、

  Binary files（TDB or TCB）
  Text files(TDB or TCB)
  SPICE data files(TDB only)
  JPL DExx compatible files(No Time transformation)

の４種類ありますが、サイズを考えればBinary Fileが良いです。
Binary filesには、TDB base とTCB baseのに種類あり、XHFとの比較も行いたいので TCB base で行います。
さらに 期間がJ2000±100年とJ2000±1000年の二種類あり、これはどちらでも構いません。

・ライブラリー（CALEPH library）
INPOPファイルを読むためには、ファイル構造を調べてプログラムを作る必要がありますが、今回は提供されているライブラリーを使うことにします。
CALEPHのダウンロードは、DOWNLOADページのUsageにある「CALEPH Library」のLINKで飛んで行きます。

https://www.imcce.fr/recherche/equipes/asd/calceph/

ここからライブラリーをDOWNLOADします。2018年1月の時点で、バージョンは 2.3.2になっています（calceph-2.3.2.tar.gz）。
同じページに、ドキュメント（Documentation）もありますのでDOWNLOADしておきます。

・ライブラリーのインストール
ライブラリーのインストールは、ドキュメント内にある「Quick instruction〜Linux」で行いました。
インストールにはgfortranが必要になりますので、gfortranをインストールしておきます（もしかしたら不要かもしれませんがTESTしていません）。インストールするフォルダーは何処でも構いません。

Windows-SDKでもインストール可能と書かれていますが、fortran77/90/95が必要になるため試していません。そもそもWindows-SDK(無料で配布されている)も現在Windows10用になっており、未だにWinXPを使っている私はインストール・実行できませんので・・・今回はLinuxのみとなりました。

インストールが成功すると、指定したフォルダーに bin include lib share の４つのフォルダーが作成されています。Cで使うのはこの内
 include/calceph.h
 lib/libcalceph.a
２つのファイルのみです。

・INPOP用プログラム
calceph-2.3.2.tar.gzを展開したフォルダー内にサンプル集（examples）があります。これを参考にプログラムを作ります。
calcephは、一つのINPOP-Fileを操作するシングル用と、複数のINPOP-Fileを操作できるマルチ用の２種類の関数があります。時間だけを求めるのでシングルを使います。

INPOP用のプログラムは簡単で

// File open
calceph_sopen("INPOP-File Name");

// Get TT-TDB(16) or TCG-TCB(17)
calceph_scompute(jdtt-integer, jdtt-fraction, 16 or 17, 0, Arrey);

// Close Inpop file
calceph_scompute();

この３つの関数のみでOKです。ただし、
TDB baseのファイルではTCG-TCBをリクエストするとエラーになり、
TCB baseのファイルではTT-TDBをリクエストするとエラーになります。
前もって、TDB baseかTCB baseか知っていた方が良いです。
ドキュメントによれば、Time baseのフラグが定数取得でわかります。

// check time base
calceph_sgetconstant("TIMESC", &dFlag)

dFlag = 0ならばTDB base
dFlag = 1ならばTCB base

になります。

さらに、ファイルの期間外を指定するとエラーになります。

CALCEPH ERROR : time 2.15936e+06 isn't in the interval [2414105 , 2488985]

期間を知るにはINPOP binary fileのヘッダーを読むしか方法は無い様なので、一度エラーさせて範囲を知った方が早いです。

INPOP.h
#ifndef		_INPOP_H
#define		_INPOP_H

#include "longjd.h"

// TDB or TCB function number
#define		FNUM_TT_TDB		16
#define		FNUM_TCG_TCB		17

#define		INPOP_NO_ERR		0
#define		INPOP_READ_ERR		-1
#define		INPOP_FILE_ERR		-2

#define		INPOP_TDB_BASE	FNUM_TT_TDB
#define		INPOP_TCB_BASE	FNUM_TCG_TCB

/*
 ------------------------------------------------
 	 	 	 proto type
 ------------------------------------------------
 */
int init_inpop(void);
void close_inpop(void);

#endif


INPOP.C
#include 
#include 
#include 
#include "astro.h"
#include "inpop.h"
#include "calceph.h"

#define	INPOP_FILE	"/home/kk/astro/inpop17a_TCB_m100_p100_tcg.dat"
#define	INPOP_SJD	2414105.5
#define	INPOP_LJD	2488985.5

// function number
int nFNum;

/*
=================================================
=		INPOP file open			=
=================================================
*/
int init_inpop(void)
{
	double dTmp;
	int nRes = INPOP_FILE_ERR;

	if (calceph_sopen(INPOP_FILE)){
		nRes = INPOP_NO_ERR;
		if (calceph_sgetconstant("TIMESC", &dTmp)){
			if ((int)dTmp == 0){
				nFNum = FNUM_TT_TDB;
			}else{
				nFNum = FNUM_TCG_TCB;
			}
		}
	}

	return nRes;
}

/*
=================================================
=		close INPOP file		=
=================================================
*/
void close_inpop(void)
{
	calceph_sclose();
}

/*
=================================================
=		Get INPOP time data		=
=================================================
*/
int inpop(LONGJD *pJD_tt, double *pdT)
{
	double dPV[6];
	int nRes = INPOP_READ_ERR;

	// heck interval
	if ((pJD_tt->dDate < INPOP_SJD) || (pJD_tt->dDate > INPOP_LJD)){
		return nRes;
	}

	// Read data
	if (calceph_scompute(pJD_tt->dDate, pJD_tt->dTime, nFNum, 0, dPV)){
		nRes = INPOP_NO_ERR;
		*pdT = dPV[0];
	}

	return nRes;
}

時間の取得は、inpop()関数です。今までと違いreturn値はフラグにしました。
期間外だった場合、低精度の計算に移行させるためです。

・計算部分
TCG-TCBの場合    TCB = TCG - (TCG - TCB)
TT-TDBの場合    TDB = TT - (TT - TDB)
と引き算になります。

逆のTCBやTDBからTT、TCGを求めるのは、XHFと同じく逐次近似法を使うしか方法はありません。

/*
=============================================
=		TT,TCG to TCB			=
=						=
=ent:	LONGJD *pdJDtt = TT JulianDay		=
=		LONGJD *pdJDtcg = TCG JulianDay	=
=exit:	return double TCB-TCG			=
=		LONGJD *pdJDtcb = TCB JulianDay	=
=============================================
*/
double inpop_tcg_tcb(LONGJD *pdJDtt, LONGJD *pdJDtcg, LONGJD *pdJDtcb)
{
	double	dTCB;

	// Get TCB - TCG
	inpop(pdJDtt, &dTCB);

	// Calc. TCB
	pdJDtcb->dDate = pdJDtcg->dDate;
	pdJDtcb->dTime = pdJDtcg->dTime - (dTCB / DAYSEC);
	ljd_normalize(pdJDtcb);

	return dTCB;
}

/*
=============================================
=		TCB to TT,TCG			=
=						=
=ent:	LONGJD *pdJDtcb = TCB JulianDay		=
=exit:	return double TCB-TCG			=
=		LONGJD *pdJDtt = TT JulianDay	=
=		LONGJD *pdJDtcg = TCG JulianDay	=
=============================================
*/
double inpop_tcb_tt(LONGJD *pdJDtcb, LONGJD *pdJDtt, LONGJD *pdJDtcg)
{
	double	dTCB;
	double nOldDay;

	*pdJDtt = *pdJDtcb;
	do{
		nOldDay = pdJDtt->dTime;
		inpop(pdJDtt, &dTCB);
		pdJDtcg->dDate = pdJDtcb->dDate;
		pdJDtcg->dTime = pdJDtcb->dTime + (dTCB / DAYSEC);
		tcg_tt(pdJDtcg, pdJDtt);
	}while(fabs(nOldDay - pdJDtt->dTime) > 1e-15);

	return dTCB;
}

まだ inpop()関数の戻り値での判断は行っていません。
低精度の計算式が出来てから飛ばします。

ソースDOWNLOAD(Linux only)
makeコマンドでコンパイルできます。
calcephのライブラリーも含んでいます。


さて、前回使ったXHFと値を比較してみます。
XCBがXHFで、TCBがINPOPです。
$ ./atime
UTC 2018/01/10 05:42:29.000000 JD = 2458128.737836
TT  2018/01/10 05:43:38.184000
TCG 2018/01/10 05:43:39.086271
XCB 2018/01/10 05:43:58.257891
TCB 2018/01/10 05:43:58.257891
TDB 2018/01/10 05:43:38.184204

TCB 2018/01/10 05:43:58.257891
TCG 2018/01/10 05:43:39.086271
TT  2018/01/10 05:43:38.184000
UTC 2018/01/10 05:42:29.000000

$ ./atime 2100 12 31 23 59
UTC 2100/12/31 23:59:00.000000 JD = 2488434.499306
TT  2101/01/01 00:00:09.184000
TCG 2101/01/01 00:00:11.911122
XCB 2101/01/01 00:01:09.856674
TCB 2101/01/01 00:01:09.856674
TDB 2101/01/01 00:00:09.183901

$ ./atime 1950
UT  1950/01/01 00:00:00.000000 JD = 2433282.500000
ET  1950/01/01 00:00:29.150000
TCG 1950/01/01 00:00:28.556163
XCB 1950/01/01 00:00:15.938376
TCB 1950/01/01 00:00:15.938376
TDB 1950/01/01 00:00:29.149930
μsec単位では、XHFとINPOPの差はみえません。

次回は、精度を落とした式です。


前に戻る        次に進む

inserted by FC2 system