C++ - JPL 天文暦データから惑星間の距離を計算!
Updated:
NASA の機関である JPL(Jet Propulsion Laboratory) が惑星探査用に編纂・発行している太陽・月・惑星の暦の最新版 DE430 には太陽・月・惑星の位置(ICRS座標系)の情報が格納されています。
それらの値を使用して、太陽・月・その他惑星の任意の2天体間の距離を C++ で計算してみました。
過去には Ruby で行っています。
0. 前提条件Permalink
- Debian GNU/Linux 10.8 (64bit) での作業を想定。
- GCC 10.2.0 (G++ 10.2.0) (C++17) でのコンパイルを想定。
1. 天文暦バイナリデータについてPermalink
当ブログ過去記事を参照のこと。
また、天文暦データには各種バージョンが存在するが、日本の国立天文台が現在使用している DE430 を当方も採用する。
2. ソースコードの作成Permalink
ここでは、実行部分のみ掲載。(全てのコードは GitHub リポジトリとして公開している)
File: dist_jpl.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
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
/***********************************************************
太陽/月/地球/その他惑星の任意の2星間の距離を計算
DATE AUTHOR VERSION
2021.04.02 mk-mode.com 1.00 新規作成
Copyright(C) 2021 mk-mode.com All Rights Reserved.
----------------------------------------------------------
引数 : [第1] 対象天体番号(必須。 1 - 13)
[第2] 基準天体番号(必須。 1 - 13)
1: 水星 (Mercury)
2: 金星 (Venus)
3: 地球 (Earth)
4: 火星 (Mars)
5: 木星 (Jupiter)
6: 土星 (Saturn)
7: 天王星 (Uranus)
8: 海王星 (Neptune)
9: 冥王星 (Pluto)
10: 月 (Moon)
11: 太陽 (Sun)
12: 太陽系重心 (Solar system Barycenter)
13: 地球 - 月の重心 (Earth-Moon Barycenter)
[第3] ユリウス日(省略可。省略時は現在日時のユリウス日)
***********************************************************/
#include "jpl.hpp"
#include <cstdlib>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
namespace dist_jpl {
// 定数
static constexpr char kAstrs[13][24] = {
"Mercury", "Venus", "Earth", "Mars", "Jupiter",
"Saturn", "Uranus", "Neptune", "Pluto", "Moon", "Sun",
"Solar system Barycenter", "Earth-Moon barycenter"
}; // 天体名称
static constexpr unsigned int kJstOffset = 32400;
static constexpr bool kFlgKm = false; // 単位フラグ
// true: km
// false: AU
static constexpr bool kFlgBary = true; // 基準フラグ
// true: 太陽系重心が基準
// false: 太陽が基準
static constexpr char kKm[] = "km"; // 単位文字列
static constexpr char kAu[] = "AU"; // 〃
/*
* @brief 変換: JST -> UTC
*
* @param[in] JST (timespec)
* @return UTC (timespec)
*/
struct timespec jst2utc(struct timespec ts_jst) {
struct timespec ts;
try {
ts.tv_sec = ts_jst.tv_sec - kJstOffset;
ts.tv_nsec = ts_jst.tv_nsec;
} catch (...) {
throw;
}
return ts;
}
/*
* @brief GC (グレゴリオ暦) -> JD (ユリウス日)
*
* @param[in] GC (timespec)
* @return JD (double)
*/
double gc2jd(struct timespec ts) {
struct tm t;
unsigned int year;
unsigned int month;
unsigned int day;
unsigned int hour;
unsigned int min;
unsigned int sec;
double jd;
try {
localtime_r(&ts.tv_sec, &t);
year = t.tm_year + 1900;
month = t.tm_mon + 1;
day = t.tm_mday;
hour = t.tm_hour;
min = t.tm_min;
sec = t.tm_sec;
// 1月,2月は前年の13月,14月とする
if (month < 3) {
--year;
month += 12;
}
// 日付(整数)部分
jd = static_cast<int>(365.25 * year)
+ static_cast<int>(year / 400.0)
- static_cast<int>(year / 100.0)
+ static_cast<int>(30.59 * (month - 2))
+ day
+ 1721088.5;
// 時間(小数)部分
jd += (sec / 3600.0 + min / 60.0 + hour) / 24.0;
// 時間(ナノ秒)部分
jd += ts.tv_nsec / 1000000000.0 / 3600.0 / 24.0;
} catch (...) {
throw;
}
return jd;
}
/*
* @brief 日時文字列生成
*
* @param[in] 日時 (timespec)
* @return 日時文字列 (string)
*/
std::string gen_time_str(struct timespec ts) {
struct tm t;
std::stringstream ss;
std::string str_tm;
try {
localtime_r(&ts.tv_sec, &t);
ss << std::setfill('0')
<< std::setw(4) << t.tm_year + 1900 << "-"
<< std::setw(2) << t.tm_mon + 1 << "-"
<< std::setw(2) << t.tm_mday << " "
<< std::setw(2) << t.tm_hour << ":"
<< std::setw(2) << t.tm_min << ":"
<< std::setw(2) << t.tm_sec << "."
<< std::setw(3) << ts.tv_nsec / 1000000;
return ss.str();
} catch (...) {
throw;
}
}
/*
* @brief 距離計算
*
* @param 基準天体から見た対象天体の座標
* @return 距離 (double)
*/
double calc_dist(double ps[3]) {
double d = 0.0;
unsigned int i;
try {
for (i = 0; i < 3; ++i) { d += ps[i] * ps[i]; }
d = std::sqrt(d);
} catch (...) {
throw;
}
return d;
}
} // namespace dist_jpl
int main(int argc, char* argv[]) {
namespace ns = dist_jpl;
unsigned int astr_t; // 天体番号(対象)
unsigned int astr_c; // 天体番号(基準)
double jd; // ユリウス日
int ret; // 関数実行結果
struct timespec jst; // 時刻オブジェクト
std::string unit; // 単位(位置/角位置)
double d; // 距離
try {
if (argc < 3) {
std::cout << "[USAGE] ./jpl_calc_430 TARGET_NO CENTER_NO [JULIAN_DAY]"
<< std::endl;
return EXIT_FAILURE;
}
// 天体番号(対象)取得
astr_t = atoi(argv[1]);
if (astr_t < 1 || astr_t > 15) {
std::cout << "[ERROR} !!! TARGET_NO must be between 1 and 15 !!!"
<< std::endl;
return EXIT_FAILURE;
}
// 天体番号(基準)取得
astr_c = atoi(argv[2]);
if (astr_c < 0 || astr_c > 13) {
std::cout << "[ERROR} !!! CENTER_NO must be between 0 and 13 !!!"
<< std::endl;
return EXIT_FAILURE;
}
// ユリウス日取得
if (argc > 3) {
jd = atof(argv[3]);
} else {
// 現在日時の取得
ret = timespec_get(&jst, TIME_UTC);
if (ret != 1) {
std::cout << "[ERROR] Could not get now time!" << std::endl;
return EXIT_FAILURE;
}
std::cout << jst.tv_sec << "." << jst.tv_nsec << std::endl;
jd = ns::gc2jd(ns::jst2utc(jst));
}
// バイナリファイル読み込み
ns::Jpl o_jpl(jd, ns::kFlgKm, ns::kFlgBary);
o_jpl.read_bin();
// 位置・速度(Positions(Radian), Velocities(Radian/Day)) 計算
o_jpl.calc_pv(astr_t, astr_c);
// 単位文字列
unit = ns::kAu;
if (ns::kFlgKm) { unit = ns::kKm; }
// 距離計算
d = ns::calc_dist(o_jpl.pos);
// 結果出力
std::cout << "DISTANCE [ "
<< ns::kAstrs[astr_t - 1] << " <=> "
<< ns::kAstrs[astr_c - 1] << " ] (JD: "
<< std::fixed << std::setprecision(8)
<< jd << ")" << std::endl;
std::cout << "= "
<< d << " " << unit << std::endl;
} catch (...) {
std::cerr << "EXCEPTION!" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
3. ソースコードのビルド(コンパイル&リンク)Permalink
- リポジトリに
Makefile
があるので、それを使用してmake
するだけ。(リビルドする際はmake clean
をしてから) - 上記の
Makefile
内では別途個別にインストールしたc++102
コマンドを使用しているが、通常はc++
であるので注意。
$ make
4. 準備Permalink
- JPL 天文暦バイナリデータ
JPLEPH
を実行ファイルと同じディレクトリ内に配置。
(参照「JPL 天文暦データのバイナリ化!」)
5. 動作確認Permalink
第1引数に対象天体番号(1
〜13
)、第2引数に基準天体番号(1
〜13
)、第3引数にユリウス日を指定して実行。(第1、第2引数は必須。天体番号一覧は dist_jpl.cpp
内のコメントを参照)
$ ./dist_jpl 3 11 2459283.625
DISTANCE [ Earth <=> Sun ] (JD: 2459283.62500000)
= 0.99311846 AU
以上、
Comments