Fortran - ISS 位置・速度(TEME 座標)の算出!

Updated:


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

0. 前提条件

  • LMDE 3 (Linuz Mint Debian Edition 3; 64bit) での作業を想定。
  • GCC 6.3.0 でのビルド(コンパイル&リンク)を想定。
  • ここでは、各種座標系、 SGP4 アルゴリズム(Simplified General Perturbations Satellite Orbit Model 4; NASA, NORAD が使用している、近地球域の衛星の軌道計算用で、周回周期225分未満の衛星に使用すべきアルゴリズム)等についての詳細は説明しない。

1. TEME 座標について

  • TEME 座標とは「真赤道面平均春分点」を基準にした座標のことで、 “True Equator, Mean Equinox” の略。
  • 今回算出する座標が TEME 座標なのは、参考にした SGP4 アルゴリズムのソースコードが TEME 座標を算出するようになっていて、それをほぼそのまま参考にしたため。

2. ソースコードの作成

以下は、実行部分。

File: tle2rv.f95

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
!*******************************************************************************
! TLE から ISS の位置・速度を計算
! : 但し、座標系は TEME(True Equator, Mean Equinox; 真赤道面平均春分点)
!
!   date          name            version
!   2018.11.21    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
! ---
! 引数 : UT1(世界時1)
!          YYYYMMDD[HHMMSS[MMM]]
!          無指定なら現在(システム日時)を UT1 とみなす。
!          (上記最後の MMM はミリ秒)
!*******************************************************************************
!
program get_iss_pv
  use const,       only : SP, DP, UID, F_TLE, FMT_DT_0, FMT_DT_1, &
                        & t_cst, getgravconst
  use ext,         only : t_time, add_day, date_fmt
  use model,       only : t_sat
  use propagation, only : propagate
  use io
  implicit none
  type(t_time)  :: ut1
  type(t_sat)   :: satrec
  type(t_cst)   :: gravconst
  character(69) :: tle(2)
  real(DP)      :: r(3), v(3)

  ! コマンドライン引数(現在日時(UT1))取得
  call get_arg(ut1)
  if (ut1%year == 0) stop

  ! TLE 読み込み
  call get_tle(ut1, tle)

  ! WGS84 用定数の取得
  gravconst = getgravconst("wgs84")

  ! ISS 初期位置・速度の取得
  call twoline2rv(tle, gravconst, satrec)

  ! 指定 UT1 の ISS 位置・速度の取得
  call propagate(gravconst, satrec, ut1, r, v)

  ! Error 時、終了
  !   1 - mean elements, ecc >= 1.0 or ecc < -0.001 or a < 0.95 er
  !   2 - mean motion less than 0.0
  !   3 - pert elements, ecc < 0.0  or  ecc > 1.0
  !   4 - semi-latus rectum < 0.0
  !   5 - epoch elements are sub-orbital
  !   6 - satellite has decayed
  if (satrec%error /= 0) then
    print '(A, I0, A)', "ERROR! [", satrec%error, "]"
    stop
  end if

  ! 結果出力
  print '("[", A, " UT1]")',  date_fmt(ut1)
  print '("TLE:", A)',               tle(1)
  print '("    ", A)',               tle(2)
  print '("TEME POS:[", 3F16.8"]")', r(1:3)
  print '("     VEL:[", 3F16.8"]")', v(1:3)

  stop
contains
  ! コマンドライン引数取得
  ! * YYYYMMDD[HHMMSS[MMM]] 形式
  ! * 17桁超入力された場合は、18桁目以降の部分は切り捨てる
  ! * コマンドライン引数がなければ、システム日時を UT1 とする
  ! * 日時の整合性チェックは行わない
  !
  ! :param(out) type(t_time) ut1
  subroutine get_arg(ut1)
    implicit none
    type(t_time), intent(out) :: ut1
    character(17) :: gc
    integer(SP)   :: dt(8)
    integer(SP)   :: len_gc

    if (iargc() == 0) then
      call date_and_time(values=dt)
      ut1 = t_time(dt(1), dt(2), dt(3), dt(5), dt(6), dt(7), dt(8))
    else
      call getarg(1, gc)
      len_gc = len(trim(gc))
      if (len_gc /= 8 .and. len_gc /= 14 .and. len_gc /= 17) then
        print *, "Format: YYYYMMDD[HHMMSS[MMM]"
        return
      end if
      read (gc, FMT_DT_0) ut1
    end if
  end subroutine get_arg

  ! TLE 取得
  !
  ! :param(in)  type(t_time)     ut1: UT1
  ! :param(out) character(69) tle(2): TLE(2行)
  subroutine get_tle(ut1, tle)
    implicit none
    type(t_time),  intent(in)  :: ut1
    character(69), intent(out) :: tle(2)
    character(69) :: buf, buf_p(2)
    integer(SP)   :: ios
    type(t_time)  :: utc  ! TLE 内の日時
    integer(SP)   :: i, l, y
    real(DP)      :: d
    character(17) :: s_ut1, s_utc

    write (s_ut1, FMT_DT_1) ut1

    ! ファイル OPEN
    open (unit   = UID,         &
        & iostat = ios,         &
        & file   = F_TLE,       &
        & action = "read",      &
        & form   = "formatted", &
        & status = "old")
    if (ios /= 0) then
      print '("[ERROR:", I0 ,"] Failed to open file: ", A)', ios, F_TLE
      stop
    end if

    buf_p = ""
    do
      read (UID, '(A)', iostat = ios) buf
      if (ios < 0) then
        exit
      else if (ios > 0) then
        print '("[ERROR:", I0 ,"] Failed to read file: ", A)', ios, F_TLE
      end if
      if (len(trim(buf)) == 0) cycle
      if (buf(1:1) /= "1" .and. buf(1:1) /= "2") cycle
      if (buf(1:1) == "1") then
        read (buf(19:20), '(I2)') y
        y = 2000 + y
        read (buf(21:32), '(F12.8)') d
        call add_day(t_time(y, 1, 1, 0, 0, 0, 0), d, utc)
        write (s_utc, FMT_DT_1) utc
        if (s_utc > s_ut1) then
          tle = buf_p
          exit
        end if
      end if
      read (buf(1:1), '(I1)') l
      buf_p(l) = buf
    end do

    ! 上記の処理で該当レコードが得られなかった場合は、最初の2行
    if (len(trim(tle(1))) == 0) then
      rewind(UID)
      do i = 1, 2
        read (UID, '(A)', iostat = ios) tle(i)
      end do
    end if

    ! ファイル CLOSE
    close(UID)
  end subroutine get_tle
end program get_iss_pv

モジュール等を含めた全てのファイルは GitHub リポジトリとして公開しているので、そちらを参照のこと。

3. 準備

4. ビルド

Makefile を使用するので、以下のようにすればよい。

$ make

5. 実行

以下は、世界時1(UT1) 2019年1月1日(0時0分0.000秒)の ISS の位置・速度(TEME 座標)を計算する例。

$ ./get_iss_pv 20190101
[2019-01-01 00:00:00.000 UT1]
TLE:1 25544U 98067A   18332.51648093  .00016717  00000-0  10270-3 0  9000
    2 25544  51.6394 279.8120 0005142  90.0655 270.1087 15.54020670 24063
TEME POS:[   3654.00994957   2141.32653133  -5300.47067891]
     VEL:[     -3.31607446      6.88762457      0.49238676]
  • UT1(世界時1)は「年・月・日・時・分・秒・ミリ秒」を最大17桁(時・分・秒(ミリ秒)は省略可)で指定する。
  • UT1(世界時1)を指定しない場合は、システム日時を UT1 とみなす。

以上。





 

Sponsored Link

 

Comments