XHFによる太陽系座標時の計算

【XHFによる太陽系座標時の計算】
太陽系座標時(TCB)の計算で一番簡単なのがXHFを使う方法です。

・XHFのダンロード
IERSのサイトからダウンロードします。
ftp://tai.bipm.org/iers/conv2010/chapter10/software

このサイトには以下のファイルがあり、全部ダウンロードします。

・fb, HF, XHFの違い
前に書いたTCBの計算式をTCB-TCGの形で表します。

IERS Technical Note No. 36によれば

fb は、上式の

の部分と思われ、
HFは、

の部分、
XHFは、HFを改良した物

と考えられますが・・・英語がよく解らないので違うかもしれません・・・でも、XHFを使うのが簡単そうです。XHFの精度は、TE405と比較して0.453ns(rms)、最大2.248nsと書かれています。

ところで、上式の残り

すなわち、地球上の観測者の項ですが、今はまだ地球の座標や速度等を調べてないので無視とします。

・xhf20020.fをCに移植
プログラムの参考例は xhf2002.f のフォートランで書かれており、Linuxならば gfortran でコンパイル・実行できる。xhf2002.fの実行結果は xhf2002.out になる。
まずは、xhf2002.fをCに移植してTESTします。

xhf2002.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// HF Data file
#define    HFDATA_PATH    "HF2002.DAT"
#define    HFDAT_NE        463
#define    HFDAT_NX        36
#define    HFDAT_TOTAL    (HFDAT_NE + HFDAT_NX + 1)

// Const
#define    HF_START_JDT    2305450.5
#define    HF_N1           218997.0

// HF data
double    dHf_a1[HFDAT_TOTAL];
double    dHf_a2[HFDAT_TOTAL];
double    dHf_f[HFDAT_TOTAL];
// HF const
double    dPI2;
double    dNd;

int read_hfdata(void){
    FILE *fp;
    int i;

    fp = fopen(HFDATA_PATH, "r");
    if (fp == NULL){
        return -1;
    }

    for(i=0; i<HFDAT_TOTAL; i++){
        fscanf(fp, "%lf %lf %lf", &dHf_a1[i], &dHf_a2[i], &dHf_f[i]);
    }
    fclose(fp);

    // makeconst
    dPI2 = 8.0 * atan(1.0);
    dNd = 4.0 / HF_N1;

    return 0;
}

double tcg_tcb(double jdt)
{
    double t;
    double dPI2t;
    double dXi;
    double s0, s1;
    double arg;
    int i, j;

    t = (jdt - HF_START_JDT) - 0.5 * HF_N1;
    dPI2t = dPI2 * t;
    dXi = dNd * t;

    s0 = dHf_a1[0] + (dHf_a2[0] + dHf_f[0] * dXi) * dXi;

    j = HFDAT_NE + HFDAT_NX;
    s1 = 0.0;
    for(i=0; i<HFDAT_NX; i++){
        arg = dPI2t * dHf_f[j];
        s1 += (dXi * (dHf_a1[j] * sin(arg) + dHf_a2[j] * cos(arg)));
        j--;
    }

    for (i=0; i<HFDAT_NE; i++){
        arg = dPI2t * dHf_f[j];
        s1 += (dHf_a1[j] * sin(arg) + dHf_a2[j] * cos(arg));
        j--;
    }

    return (s1 + s0);
}

int main(void) {
    double jdt;
    double dtcb;
    int i;

    if (read_hfdata() != 0){
        puts("Can not open HFDATA file");
        return -1;
    }

    for(i=1600; i<=2200; i+=100){
        jdt = 2451545.0 + (i - 2000) * 365.25;
        dtcb = tcg_tcb(jdt);
        printf("JDT = %.1f TCB-TCG = %.13f\n", jdt, dtcb);
    }

    return 0;
}

これをコンパイル・実行すると

$ ./xhf2002

JDT = 2305445.0 TCB-TCG = -176.1773514612386
JDT = 2341970.0 TCB-TCG = -129.4460734593428
JDT = 2378495.0 TCB-TCG = -82.7147423663766
JDT = 2415020.0 TCB-TCG = -35.9834211705210
JDT = 2451545.0 TCB-TCG = 10.7478544289339
JDT = 2488070.0 TCB-TCG = 57.4792196952081
JDT = 2524595.0 TCB-TCG = 104.2104780798795

xhf2002.out と同じ結果になり、プログラムが正しい事がわかる。

・xhfを用いたTCBプログラム
XHFは、TTからTCB-TCGを計算します。よって

TCB = TCG + XHF(TT)

という形になり、TCGとTTの２つの引数が必要になります。
他のTCB計算もあるので、今回のファイルは xhf_tcb.cとして独立させておきます。

xhf_tcb.h

/*
#############################################
#                                            #
#        TCB functions for XHF2002        #
#                                            #
#############################################
*/
#ifndef        _XHF_TCB_H
#define        _XHF_TCB_H

#include "longjd.h"

/*
---------------------------------------------
            proto type
---------------------------------------------
*/
int read_hfdata(void);
double calc_tcb_tcg(LONGJD *pLjd);

#endif

ヘッダーファイルはプロトタイプ宣言だけです。

xhf_tcb.c

/*
#############################################
#                                           #
#       TCB functions for XHF2002          #
#                                           #
#############################################
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "astro.h"
#include "xfh_tcb.h"

// HF Data file
#define   HFDATA_PATH   "HF2002.DAT"
#define   HFDAT_NE       463
#define   HFDAT_NX       36
#define   HFDAT_TOTAL   (HFDAT_NE + HFDAT_NX + 1)

// Const
#define   HF_START_JDT   2305450.5
#define   HF_N1           218997.0

// HF data
double   dHf_a1[HFDAT_TOTAL];
double   dHf_a2[HFDAT_TOTAL];
double   dHf_f[HFDAT_TOTAL];
// HF const
double   dPI2;
double    dNd;

/*
=============================================
=       Read XHF2002.DAT                   =
=exit:   success retun 0                   =
=============================================
*/
int read_hfdata(void)
{
   FILE *fp;
   int i;

   fp = fopen(HFDATA_PATH, "r");
   if (fp == NULL){
       return -1;
   }

   for(i=0; i<HFDAT_TOTAL; i++){
       fscanf(fp, "%lf %lf %lf", &dHf_a1[i], &dHf_a2[i], &dHf_f[i]);
   }
   fclose(fp);

   // makeconst
   dPI2 = 8.0 * atan(1.0);
   dNd = 4.0 / HF_N1;

   return 0;
}

/*
=============================================
=       calc. TCG-TCB                       =
=ent:   LONGJD *pLjd = TT Julianday       =
=exit:   retun TCG-TCB(sec)               =
=============================================
*/
double calc_tcb_tcg(LONGJD *pLjd)
{
   LONGJD lgTmp;
   double t;
   double dPI2t;
   double dXi;
   double s0, s1;
   double arg;
   int i, j;

   lgTmp.dDate = (pLjd->dDate - HF_START_JDT) - 0.5 * HF_N1;
   lgTmp.dTime = pLjd->dTime;

   t = ljdtojd(&lgTmp);
   dPI2t = dPI2 * t;
   dXi = dNd * t;

   s0 = dHf_a1[0] + (dHf_a2[0] + dHf_f[0] * dXi) * dXi;

   j = HFDAT_NE + HFDAT_NX;
   s1 = 0.0;
   for(i=0; i<HFDAT_NX; i++){
       arg = dPI2t * dHf_f[j];
       s1 += (dXi * (dHf_a1[j] * sin(arg) + dHf_a2[j] * cos(arg)));
       j--;
   }

   for (i=0; i<HFDAT_NE; i++){
       arg = dPI2t * dHf_f[j];
       s1 += (dHf_a1[j] * sin(arg) + dHf_a2[j] * cos(arg));
       j--;
   }

   return (s1 + s0);
}

xhf2002.cとほとんど同じ物で、HF2002.DATの読み込みと、TCB-TCG計算部分の２つです。
TCBの計算は、atime.c側で行います。

atime.c TCB計算部分

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

    // Get TCB - TCG
    dTCB = calc_tcb_tcg(pdJDtt);

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

    return dTCB;
}

引数はTTとTCGの２つが必要で、TTで計算したTCB-TCGをTCGに加えてTCBを求めます。

・TCBからTCG or TTの計算
次は、TCBからTCGを求めるプログラムを考えます。

XHFでは、TCBからTCGを一発で計算する方法はありません。
それで最初はTCBをTTとしてTCGを求め、そのTCGからTTを求め、そのTTから再びTCGを求める、逐次近似法を使うことにしました。

一回目 TCG = TCB - XHF(TCB)　→ 求めたTCGからTTを計算
二回目 TCG = TCB - XHF(TT) → 求めたTCGからTTを計算
・・・以下繰り返し
実際にやってみた所、２〜３回で収束しました。

atime.c TCB_TT関数

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

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

    return dTCB;
}

do_while間に、強制的に抜ける処理が入ってないのでバグになります。必要ならループ回数で抜けるように作ってください。

次は、太陽系力学時の計算です。

前に戻る次に進む