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

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

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

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

• Fortran - ISS 位置・速度（WGS84(BLH) 座標）の算出！
• Python - ISS 位置／速度計算（SGP4 アルゴリズム）！

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

/***********************************************************
  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;
}

• GitHub - C++ source code to calculate the BLH position of the ISS.

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

---

以上。