C++ - JPL 天文暦データから惑星間の距離を計算!

Updated:


NASA の機関である JPL(Jet Propulsion Laboratory) が惑星探査用に編纂・発行している太陽・月・惑星の暦の最新版 DE430 には太陽・月・惑星の位置(ICRS座標系)の情報が格納されています。

それらの値を使用して、太陽・月・その他惑星の任意の2天体間の距離を C++ で計算してみました。

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

0. 前提条件

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

1. 天文暦バイナリデータについて

当ブログ過去記事を参照のこと。

また、天文暦データには各種バージョンが存在するが、日本の国立天文台が現在使用している DE430 を当方も採用する。

2. ソースコードの作成

ここでは、実行部分のみ掲載。(全てのコードは 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. ソースコードのビルド(コンパイル&リンク)

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

4. 準備

5. 動作確認

第1引数に対象天体番号(113)、第2引数に基準天体番号(113)、第3引数にユリウス日を指定して実行。(第1、第2引数は必須。天体番号一覧は dist_jpl.cpp 内のコメントを参照)

$ ./dist_jpl 3 11 2459283.625
DISTANCE [ Earth <=> Sun ] (JD: 2459283.62500000)
= 0.99311846 AU

以上、





 

Sponsored Link

 

Comments