C++ - 太陽・月の視位置計算（海保略算式版）！
Updated:
C++ で、海上保安庁・海洋情報部の「コンピュータによる天体の位置計算式」を利用して、太陽や月の視位置等を計算してみました。
過去には Ruby, Python, Fortran95 で行っています。
- Ruby - 太陽・月の視赤経・視赤緯等の計算（海保略算式版）！
- Ruby - 太陽・月の視黄経・視黄緯の計算（海保略算式版）！
- Python - 太陽・月の視赤経・視赤緯等の計算（海保略算式版）！
- Python - 太陽・月の視黄経・視黄緯等の計算（海保略算式版）！
- Fortran - 太陽・月の視位置計算（海保略算式版）！
0. 前提条件
- Debian GNU/Linux 10.8 (64bit) での作業を想定。
- GCC 9.2.0 (G++ 9.2.0) (C++17) でのコンパイルを想定。
1. 海保略算式について
- 「コンピュータによる天体の位置計算式」内の PDF ドキュメントや、（当記事冒頭に記載の）当ブログ過去記事を参照のこと。
- 当然ながら、用意されている係数データファイルの年しか値を計算できない。
2. ソースコードの作成
ここでは、実行部分のみ掲載。（全てのコードは GitHub リポジトリとして公開している）
File: 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++であるので注意。
$ make
4. 動作確認
- 海保サイトからダウンロードした各年の係数ファイル(
na99-data.txtという名称のファイル)をtxtディレクトリに配置する。（但し、環境によっては、文字コード(Character Set)を Shift-JIS から UTF-8 に変換する必要がある） - 説明書に記載の ΔT の値を一覧にしたファイル
delta_t.txtをtxtディレクトリに配置する。 - そして、コマンドライン引数に UT1（世界時1）を「年・月・日・時・分・秒・ナノ秒」を最大23桁で指定して実行する。（UT1（世界時1）を先頭から部分的に指定した場合は、指定していない部分を
0とみなす）
$ ./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 ″)
以上。
Comments