C++ - ISS 位置・速度(BLH(WGS84)座標)の算出!

Updated:


前回、 C++ で NASA 提供の最新の TLE(2行軌道要素形式)から任意の時刻(UT1; 世界時1)の ISS の位置・速度(TEME 座標)を、 SGP4 アルゴリズムを用いて計算しました。

今回は、これの応用で、取得した TEME 座標を WGS84 座標(いわゆる、緯度・経度・高度(BLH)という座標)に変換します。

過去には Ruby, Python, Fortran で実装しています。(但し、 Ruby 版はブログ記事にはしていない)

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分未満の衛星に使用すべきアルゴリズム)等についての詳細は説明しない。(各種座標系については次項で概要のみ説明)

1. 各種座標系について(概要)

  • TEME 座標とは「真赤道面平均春分点」のことで、 “True Equator, Mean Equinox” の略。
  • PEF 座標とは「擬地球固定座標系」のことで、 “Pseudo Earth Fixed” の略。
  • ECEF 座標とは「地球中心・地球固定直交座標系」のことで、 “Earth Centered, Earth Fixed” の略。
  • WGS84 座標とは「世界測地系1984」のことで、 “World Geodetic System 1984” の略。 GPS で使用する測地系。

2. C++ ソースコードの作成

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

File: iss_sgp4_blh.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
/***********************************************************
  TLE から ISS の位置・速度を計算
  : 但し、座標系は BLH(Beta(Latitude), Lambda(Longitude), Height)
  : 指定 UT1 からの TEME 座標計算は iss_sgp4_teme

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

  Copyright(C) 2021 mk-mode.com All Rights Reserved.
  ---
  引数 : JST(日本標準時)
           書式:最大23桁の数字
                 (先頭から、西暦年(4), 月(2), 日(2), 時(2), 分(2), 秒(2),
                             1秒未満(9)(小数点以下9桁(ナノ秒)まで))
                 無指定なら現在(システム日時)と判断。
  ---
  MEMO:
    TEME: True Equator, Mean Equinox; 真赤道面平均春分点
     PEF: Pseudo Earth Fixed; 擬地球固定座標系
    ECEF: Earth Centered, Earth Fixed; 地球中心・地球固定直交座標系
***********************************************************/
#include "blh.hpp"
#include "eop.hpp"
#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_blh;
  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 jst;           // JST
  struct timespec utc;           // UTC
  struct timespec ut1;           // UT1
  struct timespec tai;           // TAI
  double pm_x;                   // 極運動(x)
  double pm_y;                   // 極運動(y)
  std::string lod_str;           // LOD 文字列
  double lod;                    // LOD
  std::vector<std::string> tle;  // TLE
  ns::Satellite sat;             // 衛星情報
  ns::PvTeme    teme;            // 位置・速度(TEME)
  ns::PvBlh     blh;             // 位置・速度(BLH)

  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");
      jst.tv_sec  = mktime(&t);
      jst.tv_nsec = 0;
      if (s_tm > 14) {
        jst.tv_nsec = std::stod(
            tm_str.substr(14, s_nsec) + std::string(9 - s_nsec, '0'));
      }
    } else {
      // 現在日時の取得
      ret = std::timespec_get(&jst, TIME_UTC);
      if (ret != 1) {
        std::cout << "[ERROR] Could not get now time!" << std::endl;
        return EXIT_FAILURE;
      }
    }

    // EOP データ取得
    utc = ns::jst2utc(jst);
    ns::Eop o_e(utc);
    pm_x = o_e.pm_x;
    pm_y = o_e.pm_y;
    lod  = o_e.lod;
    ut1  = ns::utc2ut1(utc);
    tai  = ns::utc2tai(utc);
    std::cout << ns::gen_time_str(jst) << " JST" << std::endl;
    std::cout << ns::gen_time_str(utc) << " UTC" << std::endl;
    std::cout << ns::gen_time_str(ut1) << " UT1" << std::endl;
    std::cout << ns::gen_time_str(tai) << " TAI" << std::endl;

    // 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);

    // TEME -> BLH 変換
    ns::Blh o_b(ut1, tai, pm_x, pm_y, lod);
    blh = o_b.teme2blh(teme);

    // 結果出力
    std::cout << "---" << std::endl;
    std::cout << std::fixed << std::setprecision(8);
    std::cout << "EOP: " << o_e.eop << std::endl;
    std::cout << "---" << std::endl;
    std::cout << "TLE: " << tle[0] << std::endl;
    std::cout << "     " << tle[1] << std::endl;
    std::cout << "---" << 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;
    std::cout << "---" << std::endl;
    std::cout << std::fixed << std::setprecision(4);
    std::cout << "WGS84(BLH):" << std::endl
              << "     BETA(Latitude) = " << std::setw(9) << blh.r.b
              << " °" << std::endl
              << "  LAMBDA(Longitude) = " << std::setw(9) << blh.r.l
              << " °" << std::endl
              << "             HEIGHT = " << std::setw(9) << blh.r.h / 1000.0
              << " km" << std::endl
              << "           VELOCITY = " << std::setw(9) << blh.v
              << " km/s" << std::endl;
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

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

  • リポジトリMakefile があるので、それを使用して make するだけ。(リビルドする際は make clean をしてから)
$ make

4. 準備

  • うるう秒ファイルを こちら からダウンロードし、実行プログラムと同じディレクトリ内に Leap_Second.dat として配置する。
  • 地球姿勢(回転)パラメータ用のデータファイルを、当ブログ過去記事「C++ - EOP(地球姿勢パラメータ)データファイル 生成!」を参考に生成し、実行プログラムと同じディレクトリ内に eop.txt として配置する。(計算対象日の EOP データが存在しないと、計算できないので注意)
  • TLE(Two-line elements; 2行軌道要素形式)データを こちら からダウンロードし、実行プログラムと同じディレクトリ内に tle.txt として配置する。

5. 動作確認

コマンドライン引数には JST(日本標準時) を最大23桁(年(4)・月(2)・日(2)・時(2)・分(2)・秒(2)・ナノ秒(9))で指定し、実行。(JST を指定しない場合は、システム時刻を JST とみなす)

$ ./iss_sgp4_blh 20210629120000
2021-06-29 12:00:00.000 JST
2021-06-29 03:00:00.000 UTC
2021-06-29 02:59:59.855 UT1
2021-06-29 03:00:37.000 TAI
2021-06-29 03:01:09.184 TT
JD(UT1): 2459394.62499832 day
JD(TT) : 2459394.62580074 day
JCN(TT): 0.21491104
---
EOP: 2021-06-29 59394.00 P  0.201594 0.007940  0.416775 0.010428 P -0.1449344 0.0085261
---
TLE: 1 25544C 98067A   21162.50000000  .00044577  00000-0  80837-3 0    17
     2 25544  51.6415  10.8602 0003442  86.5921  66.0554 15.48948920    12
---
TEME: POS = [  -4126.99595304,    -889.83934472,   -5329.21689869]
      VEL = [      1.72112275,      -7.45393134,      -0.09112474]
---
WGS84(BLH):
     BETA(Latitude) =  -51.7888 °
  LAMBDA(Longitude) = -130.2547 °
             HEIGHT =  433.8826 km
           VELOCITY =    7.6506 km/s

以上。





 

Sponsored Link

 

Comments