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

Updated:


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

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

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

2. ソースコードの作成

以下は、実行部分。

File: iss_sgp4.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
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
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
!*******************************************************************************
! Getting of ISS position.
!
!   date          name            version
!   2018.11.27    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
! ---
! 引数 : JST(日本標準時)
!          YYYYMMDD[HHMMSS[MMM]]
!          無指定なら現在(システム日時)を JST とみなす。
!          (上記最後の MMM はミリ秒)
! ---
! MEMO:
!   TEME: True Equator, Mean Equinox; 真赤道面平均春分点
!    PEF: Pseudo Earth Fixed; 擬地球固定座標系
!   ECEF: Earth Centered, Earth Fixed; 地球中心・地球固定直交座標系
!*******************************************************************************
!
program iss_sgp4
  use const,       only : SP, DP, JST_UTC, UID, F_TLE, SEC_D, TT_TAI, &
                        & DAY_JC, FMT_DT_0, FMT_DT_1, &
                        & getgravconst
  use propagation, only : propagate
  use io
  use model
  use ext
  use eph
  implicit none
  type(t_time)  :: jst, utc, ut1, tai, tt
  type(t_eop)   :: eop
  type(t_sat)   :: satrec
  type(t_cst)   :: gravconst
  character(69) :: tle(2)
  integer(SP)   :: dat
  real(DP)      :: jd_ut1, jd_tt, t_tt, gmst, om, gmst_g, r(3), v(3)
  real(DP)      :: mtx_z(3, 3), mtx_pm(3, 3)
  real(DP)      :: r_pef(3), r_ecef(3), om_e(3), v_pef(3), v_ecef(3)
  real(DP)      :: lat, lon, ht, vel_ecef

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

  ! 初期処理
  call init
  print '(A)', "[初期データ]"
  print '(A, " JST")',            date_fmt(jst)
  print '(A, " UTC")',            date_fmt(utc)
  print '(A, " UT1")',            date_fmt(ut1)
  print '("TAI - UTC = ", I0)',             dat
  print '(A, " TAI")',            date_fmt(tai)
  print '(A, "  TT")',            date_fmt(tt )
  print '("JD(UT1) = ", F11.3, " day")', jd_ut1
  print '("JD(TT ) = ", F16.8, " day")',  jd_tt
  print '(" T(TT ) = ", F10.8)',           t_tt
  print '("TLE = ", A)',                 tle(1)
  print '("      ", A)',                 tle(2)

  ! ISS 位置・速度取得
  call get_iss
  print '(/A)', "[途中経過]"
  print '("TEME POS = ", 3F16.8)', r
  print '("     VEL = ", 3F16.8)', v

  ! TEME -> BLH 変換
  call teme2blh
  print '("  GMST = ", F10.8, " rad.")',                 gmst
  print '("    om = ", F10.8, " rad.")',                   om
  print '("GMST_G = ", F10.8, " rad.")',               gmst_g
  print '("ROTATE MATRIX(for GMST) =", 3(/3(E20.8)))',  mtx_z
  print '("ROTATE MATRIX(for PM) =",   3(/3(E20.8)))', mtx_pm
  print '("POSITION(PEF) = ",             /3(F16.8))',  r_pef
  print '("POSITION(ECEF) = ",            /3(F16.8))', r_ecef
  print '("om_earth =",                   /3(E20.8))',   om_e
  print '("VELOCITY(PEF) = ",             /3(F16.8))',  v_pef
  print '("VELOCITY(ECEF) = ",            /3(F16.8))', v_ecef

  ! 最終結果出力
  print '(/A)', "[計算結果]"
  print '(A)',  "WGS84(BLH):"
  print '("  POSITION  LAT = ", F9.4, " °"  )',      lat
  print '("            LON = ", F9.4, " °"  )',      lon
  print '("         HEIGHT = ", F9.4, " km"  )',       ht
  print '("  VELOCITY      = ", F9.4, " km/s")', vel_ecef

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

    if (iargc() == 0) then
      call date_and_time(values=dt)
      jst = 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) jst
    end if
  end subroutine get_arg

  ! 初期処理
  subroutine init
    implicit none

    ! 各種時刻換算
    call add_day(jst, -JST_UTC / 24.0_DP, utc)
    call get_eop(utc, eop)
    call add_day(utc, eop%dut1/ SEC_D, ut1)
    call get_dat(utc, dat)
    call add_day(utc, dat / SEC_D, tai)
    call add_day(tai, TT_TAI / SEC_D, tt)
    call gc2jd(ut1, jd_ut1)
    call gc2jd(tt, jd_tt)
    t_tt = (jd_tt - 2451545.0_DP) / DAY_JC

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

  ! ISS 位置・速度取得
  subroutine get_iss
    implicit none

    ! 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
  end subroutine get_iss

  ! TEME -> BLH 変換
  subroutine teme2blh
    implicit none

    ! GMST(グリニッジ平均恒星時)計算
    gmst = calc_gmst(jd_ut1)
    ! Ω(月の平均昇交点黄経)計算(IAU1980章動理論)
    om = calc_om(t_tt)
    ! GMST に運動項を適用(1997年より新しい場合)
    gmst_g = apply_kinematic(gmst, jd_ut1, om)
    ! GMST 回転行列(z軸を中心とした回転)
    mtx_z = r_z(gmst_g)
    ! 極運動(Polar Motion)回転行列
    mtx_pm = r_pm(eop%pm_x, eop%pm_y, t_tt)
    ! PEF 座標の計算(GMST 回転行列の適用)
    r_pef = matmul(r, mtx_z)
    ! ECEF 座標(位置)の計算(極運動(Polar Motion)の適用)
    r_ecef = matmul(r_pef, mtx_pm)
    ! Ω_earth値の計算
    om_e = calc_om_e(eop%lod)
    ! PEF 座標(速度)の計算(GMST 回転行列の適用)
    v_pef = matmul(v, mtx_z) - v_cross(om_e, r_pef)
    ! ECEF 座標(速度)の計算(極運動(Polar Motion)の適用)
    v_ecef = matmul(v_pef, mtx_pm)
    ! ECEF 座標 => BLH(Beta, Lambda, Height) 変換
    call ecef2blh(r_ecef, lat, lon, ht)
    vel_ecef = sqrt(sum(v_ecef * v_ecef))
  end subroutine teme2blh

  ! 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

    tle   = ""
    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 iss_sgp4

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

3. 準備

4. ビルド

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

$ make

5. 実行

以下は、日本標準時(JST) 2019年1月1日(0時0分0.000秒)の ISS の位置・速度(WGS84(BLH) 座標)を計算する例。
計算の流れが分かるよう、途中経過も出力するようにしている。

$ ./iss_sgp4 20190101
[初期データ]
2019-01-01 00:00:00.000 JST
2018-12-31 15:00:00.000 UTC
2018-12-31 14:59:59.963 UT1
TAI - UTC = 37
2018-12-31 15:00:37.000 TAI
2018-12-31 15:01:09.184  TT
JD(UT1) = 2458484.125 day
JD(TT ) = 2458484.12580074 day
 T(TT ) = 0.18998291
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 =    -923.64052601   6327.37451246  -2273.32176709
     VEL =      -5.27233259      1.18198845      5.43472644
  GMST = 5.67215877 rad.
    om = 2.05236340 rad.
GMST_G = 5.67215878 rad.
ROTATE MATRIX(for GMST) =
      0.81905952E+00     -0.57370855E+00      0.00000000E+00
      0.57370855E+00      0.81905952E+00      0.00000000E+00
      0.00000000E+00      0.00000000E+00      0.10000000E+01
ROTATE MATRIX(for PM) =
      0.10000000E+01     -0.43289968E-10      0.48424645E-09
      0.43289968E-10      0.10000000E+01     -0.13400880E-08
     -0.48424645E-09      0.13400880E-08      0.10000000E+01
POSITION(PEF) =
  -4386.58540889   4652.59588753  -2273.32176709
POSITION(ECEF) =
  -4386.58541019   4652.59589039  -2273.32175873
om_earth =
      0.00000000E+00      0.00000000E+00      0.72921151E-04
VELOCITY(PEF) =
     -4.65719844     -1.73678852      5.43472644
VELOCITY(ECEF) =
     -4.65719844     -1.73678852      5.43472644

[計算結果]
WGS84(BLH):
  POSITION  LAT =  -19.6855 °
            LON =  133.3144 °
         HEIGHT =  410.7864 km
  VELOCITY      =    7.3649 km/s
  • JST(日本標準時)は「年・月・日・時・分・秒・ミリ秒」を最大17桁(時・分・秒(ミリ秒)は省略可)で指定する。
  • JST(日本標準時)を指定しない場合は、システム日時を JST とみなす。

当方、このプログラムを応用して 10秒間隔2日分の JSON データを作成するプログラムも作成しています。

さらに、 Python や GMT と連携し、 ISS が特定の地点に近付いたら自動でツイートするようにもしております。

以上。





 

Sponsored Link

 

Comments