C++ で、海上保安庁・海洋情報部の「コンピュータによる天体の位置計算式」を利用して、太陽や月の視位置等を計算してみました。

過去には Ruby, Python, Fortran95 で行っています。

0. 前提条件

  • Debian GNU/Linux 10.8 (64bit) での作業を想定。
  • GCC 9.2.0 (G++ 9.2.0) (C++17) でのコンパイルを想定。

1. 海保略算式について

  • コンピュータによる天体の位置計算式」内の PDF ドキュメントや、(当記事冒頭に記載の)当ブログ過去記事を参照のこと。
  • 当然ながら、用意されている係数データファイルの年しか値を計算できない。

2. ソースコードの作成

ここでは、実行部分のみ掲載。(全てのコードは GitHub リポジトリとして公開している)

ephemeris_jcg.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
/***********************************************************
  海上保安庁の天測暦より太陽・月の視位置を計算
  (視黄経・視黄緯の計算を追加したもの)

    DATE          AUTHOR          VERSION
    2020.12.07    mk-mode.com     1.00 新規作成

  Copyright(C) 2020 mk-mode.com All Rights Reserved.

  引数 : UT1(世界時1)
           書式:最大23桁の数字
                 (先頭から、西暦年(4), 月(2), 日(2), 時(2), 分(2), 秒(2),
                             1秒未満(9)(小数点以下9桁(ナノ秒)まで))
                 無指定なら現在(システム日時)と判断。
***********************************************************/
#include "common.hpp"
#include "eph_jcg.hpp"

#include <cstdlib>   // for EXIT_XXXX
#include <ctime>
#include <iomanip>
#include <iostream>
#include <string>

int main(int argc, char* argv[]) {
  std::string tm_str;   // time string
  unsigned int s_tm;    // size of time string
  int s_nsec;           // size of nsec string
  int ret;              // return of functions
  struct timespec ut1;  // UTC
  struct tm t = {};     // for work
  namespace ns = ephemeris_jcg;

  try {
    // 日付取得
    if (argc > 1) {
      // コマンドライン引数より取得
      tm_str = argv[1];
      s_tm = tm_str.size();
      if (s_tm > 23) {
        std::cout << "[ERROR] Over 23-digits!" << std::endl;
        return EXIT_FAILURE;
      }
      s_nsec = s_tm - 14;
      std::istringstream is(tm_str);
      is >> std::get_time(&t, "%Y%m%d%H%M%S");
      ut1.tv_sec  = mktime(&t);
      ut1.tv_nsec = 0;
      if (s_tm > 14) {
        ut1.tv_nsec = std::stod(
            tm_str.substr(14, s_nsec) + std::string(9 - s_nsec, '0'));
      }
    } else {
      // 現在日時の取得
      ret = std::timespec_get(&ut1, TIME_UTC);
      if (ret != 1) {
        std::cout << "[ERROR] Could not get now time!" << std::endl;
        return EXIT_FAILURE;
      }
    }

    // Calculation & display
    ns::EphJcg o_e(ut1);
    std::cout << "[ UT1: " << ns::gen_time_str(ut1) << " ]" << std::endl;
    std::cout << std::fixed << std::setprecision(8);
    std::cout << "SUN     R.A. =  "
              << std::setw(12) << o_e.sun_ra   << " h"
              << " (= " << ns::hour2hms(o_e.sun_ra) << ")" << std::endl;
    std::cout << "        Dec. =  "
              << std::setw(12) << o_e.sun_dec  << " °"
              << " (= " << ns::deg2dms(o_e.sun_dec) << ")" << std::endl;
    std::cout << "       Dist. =  "
              << std::setw(12) << o_e.sun_dist << " AU" << std::endl;
    std::cout << "          hG =  "
              << std::setw(12) << o_e.sun_hg   << " h"
              << " (= " << ns::hour2hms(o_e.sun_hg) << ")" << std::endl;
    std::cout << "        S.D. =  "
              << std::setw(12) << o_e.sun_sd   << " ′"
              << " (= " << ns::deg2dms(o_e.sun_sd / 60.0) << ")" << std::endl;
    std::cout << "VENUS   R.A. =  "
              << std::setw(12) << o_e.vns_ra   << " h"
              << " (= " << ns::hour2hms(o_e.vns_ra) << ")" << std::endl;
    std::cout << "        Dec. =  "
              << std::setw(12) << o_e.vns_dec  << " °"
              << " (= " << ns::deg2dms(o_e.vns_dec) << ")" << std::endl;
    std::cout << "       Dist. =  "
              << std::setw(12) << o_e.vns_dist << " AU" << std::endl;
    std::cout << "          hG =  "
              << std::setw(12) << o_e.vns_hg   << " h"
              << " (= " << ns::hour2hms(o_e.vns_hg) << ")" << std::endl;
    std::cout << "        S.D. =  "
              << std::setw(12) << o_e.vns_sd   << " ″"
              << " (= " << ns::deg2dms(o_e.vns_sd / 3600.0) << ")" << std::endl;
    std::cout << "MARS    R.A. =  "
              << std::setw(12) << o_e.mrs_ra   << " h"
              << " (= " << ns::hour2hms(o_e.mrs_ra) << ")" << std::endl;
    std::cout << "        Dec. =  "
              << std::setw(12) << o_e.mrs_dec  << " °"
              << " (= " << ns::deg2dms(o_e.mrs_dec) << ")" << std::endl;
    std::cout << "       Dist. =  "
              << std::setw(12) << o_e.mrs_dist << " AU" << std::endl;
    std::cout << "          hG =  "
              << std::setw(12) << o_e.mrs_hg   << " h"
              << " (= " << ns::hour2hms(o_e.mrs_hg) << ")" << std::endl;
    std::cout << "        S.D. =  "
              << std::setw(12) << o_e.mrs_sd   << " ″"
              << " (= " << ns::deg2dms(o_e.mrs_sd / 3600.0) << ")" << std::endl;
    std::cout << "JUPITER R.A. =  "
              << std::setw(12) << o_e.jpt_ra   << " h"
              << " (= " << ns::hour2hms(o_e.jpt_ra) << ")" << std::endl;
    std::cout << "        Dec. =  "
              << std::setw(12) << o_e.jpt_dec  << " °"
              << " (= " << ns::deg2dms(o_e.jpt_dec) << ")" << std::endl;
    std::cout << "       Dist. =  "
              << std::setw(12) << o_e.jpt_dist << " AU" << std::endl;
    std::cout << "          hG =  "
              << std::setw(12) << o_e.jpt_hg   << " h"
              << " (= " << ns::hour2hms(o_e.jpt_hg) << ")" << std::endl;
    std::cout << "     S.D.(P) =  "
              << std::setw(12) << o_e.jpt_sd_p << " ″"
              << " (= " << ns::deg2dms(o_e.jpt_sd_p / 3600.0) << ")" << std::endl;
    std::cout << "     S.D.(E) =  "
              << std::setw(12) << o_e.jpt_sd_e << " ″"
              << " (= " << ns::deg2dms(o_e.jpt_sd_e / 3600.0) << ")" << std::endl;
    std::cout << "SATURN  R.A. =  "
              << std::setw(12) << o_e.sat_ra   << " h"
              << " (= " << ns::hour2hms(o_e.sat_ra) << ")" << std::endl;
    std::cout << "        Dec. =  "
              << std::setw(12) << o_e.sat_dec  << " °"
              << " (= " << ns::deg2dms(o_e.sat_dec) << ")" << std::endl;
    std::cout << "       Dist. =  "
              << std::setw(12) << o_e.sat_dist << " AU" << std::endl;
    std::cout << "          hG =  "
              << std::setw(12) << o_e.sat_hg   << " h"
              << " (= " << ns::hour2hms(o_e.sat_hg) << ")" << std::endl;
    std::cout << "     S.D.(P) =  "
              << std::setw(12) << o_e.sat_sd_p << " ″"
              << " (= " << ns::deg2dms(o_e.sat_sd_p / 3600.0) << ")" << std::endl;
    std::cout << "     S.D.(E) =  "
              << std::setw(12) << o_e.sat_sd_e << " ″"
              << " (= " << ns::deg2dms(o_e.sat_sd_e / 3600.0) << ")" << std::endl;
    std::cout << "R            =  "
              << std::setw(12) << o_e.r        << " h"
              << " (= " << ns::hour2hms(o_e.r) << ")" << std::endl;
    std::cout << "EPSILON      =  "
              << std::setw(12) << o_e.eps      << " °"
              << " (= " << ns::deg2dms(o_e.eps) << ")" << std::endl;
    std::cout << "MOON    R.A. =  "
              << std::setw(12) << o_e.mon_ra   << " h"
              << " (= " << ns::hour2hms(o_e.mon_ra) << ")" << std::endl;
    std::cout << "        Dec. =  "
              << std::setw(12) << o_e.mon_dec  << " °"
              << " (= " << ns::deg2dms(o_e.mon_dec) << ")" << std::endl;
    std::cout << "        H.P. =  "
              << std::setw(12) << o_e.mon_hp   << " °"
              << " (= " << ns::deg2dms(o_e.mon_hp) << ")" << std::endl;
    std::cout << "          hG =  "
              << std::setw(12) << o_e.mon_hg   << " h"
              << " (= " << ns::hour2hms(o_e.mon_hg) << ")" << std::endl;
    std::cout << "        S.D. =  "
              << std::setw(12) << o_e.mon_sd   << " ′"
              << " (= " << ns::deg2dms(o_e.mon_sd / 60.0) << ")" << std::endl;
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

3. ソースコードのビルド(コンパイル&リンク)

  • リポジトリMakefile があるので、それを使用して make するだけ。(リビルドする際は make clean をしてから)
  • 上記の Makefile 内では別途個別にインストールした c++92 コマンドを使用しているが、通常は c++ であるので注意。
1
$ make

4. 動作確認

  • 海保サイトからダウンロードした各年の係数ファイル(na99-data.txt という名称のファイル)を txt ディレクトリに配置する。(但し、環境によっては、文字コード(Character Set)を Shift-JIS から UTF-8 に変換する必要がある)
  • 説明書に記載の ΔT の値を一覧にしたファイル delta_t.txttxt ディレクトリに配置する。
  • そして、コマンドライン引数に UT1(世界時1)を「年・月・日・時・分・秒・ナノ秒」を最大23桁で指定して実行する。(UT1(世界時1)を先頭から部分的に指定した場合は、指定していない部分を 0 とみなす)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
$ ./ephemeris_jcg 20210225123456123456789
[ UT1: 2021-02-25 12:34:56.123 ]
SUN     R.A. =   22.58939081 h (=  22 h 35 m 21.807 s)
        Dec. =   -8.89332277 ° (= -  8 °53 ′35.962 ″)
       Dist. =    0.98992685 AU
          hG =    0.36567242 h (=   0 h 21 m 56.421 s)
        S.D. =   16.18301392 ′ (=    0 °16 ′10.981 ″)
VENUS   R.A. =   22.16846023 h (=  22 h 10 m  6.457 s)
        Dec. =  -12.69874078 ° (= - 12 °41 ′55.467 ″)
       Dist. =    1.69930817 AU
          hG =    0.78660301 h (=   0 h 47 m 11.771 s)
        S.D. =    4.88434065 ″ (=    0 ° 0 ′ 4.884 ″)
MARS    R.A. =    3.56720128 h (=   3 h 34 m  1.925 s)
        Dec. =   20.68476558 ° (=   20 °41 ′ 5.156 ″)
       Dist. =    1.43241090 AU
          hG =   19.38786195 h (=  19 h 23 m 16.303 s)
        S.D. =    3.28118141 ″ (=    0 ° 0 ′ 3.281 ″)
JUPITER R.A. =   21.23015770 h (=  21 h 13 m 48.568 s)
        Dec. =  -16.63843988 ° (= - 16 °38 ′18.384 ″)
       Dist. =    5.98670945 AU
          hG =    1.72490553 h (=   1 h 43 m 29.660 s)
     S.D.(P) =   15.38407715 ″ (=    0 ° 0 ′15.384 ″)
     S.D.(E) =   16.43640816 ″ (=    0 ° 0 ′16.436 ″)
SATURN  R.A. =   20.70890348 h (=  20 h 42 m 32.053 s)
        Dec. =  -18.68353162 ° (= - 18 °41 ′ 0.714 ″)
       Dist. =   10.83194503 AU
          hG =    2.24615975 h (=   2 h 14 m 46.175 s)
     S.D.(P) =    6.81318081 ″ (=    0 ° 0 ′ 6.813 ″)
     S.D.(E) =    7.63482456 ″ (=    0 ° 0 ′ 7.635 ″)
R            =   10.37280672 h (=  10 h 22 m 22.104 s)
EPSILON      =   23.43731626 ° (=   23 °26 ′14.339 ″)
MOON    R.A. =    9.15218237 h (=   9 h  9 m  7.857 s)
        Dec. =   20.79643783 ° (=   20 °47 ′47.176 ″)
        H.P. =    0.96526323 ° (=    0 °57 ′54.948 ″)
          hG =   13.80288087 h (=  13 h 48 m 10.371 s)
        S.D. =   15.78136265 ′ (=    0 °15 ′46.882 ″)

以上。