#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; } |
/* ############################################# # # # 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 |
/* ############################################# # # # 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); } |
/* ############################################# # # # 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; } |
/* ============================================= = 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; } |