C++ - ISS 位置・速度（TEME 座標）の算出！
Updated:
C++ で、 NASA 提供の最新の TLE（2行軌道要素形式）から任意の時刻（UT1; 世界時1）の ISS の位置・速度（TEME 座標）を、 SGP4 アルゴリズムを用いて計算してみました。
過去には Ruby, Python, Fortran で実装しています。（但し、 Ruby, Python 版はブログ記事にはしていない）
0. 前提条件
- Debian GNU/Linux 10.9 (64bit) での作業を想定。
- GCC 10.2.0 (G++ 10.2.0) (C++17) でのコンパイル＆ビルドを想定。
- ここでは、各種座標系、 SGP4 アルゴリズム（Simplified General Perturbations Satellite Orbit Model 4; NASA, NORAD が使用している、近地球域の衛星の軌道計算用で、周回周期225分未満の衛星に使用すべきアルゴリズム）等についての詳細は説明しない。（TEME 座標については、次項で概要のみ説明）
1. TEME 座標について
- TEME 座標とは「真赤道面平均春分点」を基準にした座標のことで、 “True Equator, Mean Equinox” の略。
- 今回算出する座標が TEME 座標なのは、参考にした SGP4 アルゴリズムのソースコードが TEME 座標を算出するようになっていて、それをほぼそのまま参考にしたため。
2. C++ ソースコードの作成
ここでは、実行部分のみ掲載。（全てのコードは GitHub リポジトリとして公開している）
File: iss_sgp4_teme.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
/***********************************************************
TLE から ISS の位置・速度を計算
: 但し、座標系は TEME(True Equator, Mean Equinox; 真赤道面平均春分点)
: 指定 JST における BLH 座標の計算は iss_sgp4_blh
DATE AUTHOR VERSION
2021.05.13 mk-mode.com 1.00 新規作成
Copyright(C) 2021 mk-mode.com All Rights Reserved.
引数 : UT1（世界時1）
書式：最大23桁の数字
（先頭から、西暦年(4), 月(2), 日(2), 時(2), 分(2), 秒(2),
1秒未満(9)（小数点以下9桁（ナノ秒）まで））
無指定なら現在(システム日時)と判断。
***********************************************************/
#include "sgp4.hpp"
#include "time.hpp"
#include "tle.hpp"
#include <cstdlib> // for EXIT_XXXX
#include <ctime>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
int main(int argc, char* argv[]) {
namespace ns = iss_sgp4_teme;
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 tm t = {}; // for work
struct timespec ut1; // UT1
std::vector<std::string> tle; // TLE
ns::Satellite sat; // 衛星情報
ns::PvTeme teme; // 位置・速度
try {
// 現在日時(UT1) 取得
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;
}
}
// TLE 読み込み, gravconst 取得
ns::Tle o_t(ut1);
tle = o_t.get_tle();
// ISS 初期位置・速度の取得
ns::Sgp4 o_s(ut1, tle);
sat = o_s.twoline2rv();
// 指定 UT1 の ISS 位置・速度の取得
teme = o_s.propagate(sat);
// Calculation & display
std::cout << "[" << ns::gen_time_str(ut1) << " UT1]" << std::endl;
std::cout << "TLE:" << tle[0] << std::endl;
std::cout << " " << tle[1] << std::endl;
std::cout << std::fixed << std::setprecision(8);
std::cout << "TEME POS:["
<< std::setw(16) << teme.r.x << ", "
<< std::setw(16) << teme.r.y << ", "
<< std::setw(16) << teme.r.z
<< "]" << std::endl;
std::cout << " VEL:["
<< std::setw(16) << teme.v.x << ", "
<< std::setw(16) << teme.v.y << ", "
<< std::setw(16) << teme.v.z
<< "]" << std::endl;
} catch (...) {
std::cerr << "EXCEPTION!" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
3. ソースコードのビルド（コンパイル＆リンク）
- リポジトリに
Makefileがあるので、それを使用してmakeするだけ。（リビルドする際はmake cleanをしてから）
$ make
4. 準備
- TLE（Two-line elements; 2行軌道要素形式）データを こちら からダウンロード。実行プログラムと同じディレクトリ内に
tle.txtとリネームして配置する。
5. 動作確認
コマンドライン引数には UT1（世界時1） を最大23桁（年(4)・月(2)・日(2)・時(2)・分(2)・秒(2)・ナノ秒(9)）で指定し、実行。（UT1 を指定しない場合は、システム時刻 UT1 とみなす）
$ ./iss_sgp4_teme 20210616120000
[2021-06-16 12:00:00.000 UT1]
TLE:1 25544C 98067A 21166.50000000 .00096922 00000-0 17562-2 0 173
2 25544 51.6425 351.0748 0003520 101.0606 51.5923 15.48957948 13
TEME POS:[ 5411.17576869, -3324.05178603, -2436.40925557]
VEL:[ 4.40352859, 3.26461608, 5.34294836]
以上。
Comments