ＳＧＰ計算モデル

ホーム 上へ 衛星の軌道要素 経緯度・高さと直交座標 天空の時計と原子時計 ちょこっと（逆三角関数） 衛星軌道面での動き ベクトルと座標回転 衛星の方位と仰角 ちょこっと（ＮＴＰ） 可視範囲と等距離円 衛星の速度とＳＡ 地球の大きさと形 ちょこっと（方位磁石） ＳＧＰ計算モデル アンテナの自動制御 リグのコントロール ちょこっと（ドップラ補正） DDEインターフェース 位置計算の数値例 マクロを使った自動運転 ちょこっと（軌道コレクション）

純正の計算モデルを使ってみる

　教科書的な衛星の位置計算について、これまで述べてきました。　座標の回転などは教科書どおりの処理で問題ありませんが、より正確な衛星の位置を推定するには摂動について実際的な補正をする必要があります。

　NASA２ラインフォーマットには衛星毎の摂動補正項が含まれていて、これを使うことで衛星位置の推定精度が高まることが期待されます。　SGPと呼ばれる衛星の位置計算モデルは、これらの摂動項を利用するものです。　ＣＡＬＳＡＴ３２には計算モデルとして、通常(永年項のみ補正),SGP,SGP4/SDP4,SGP8/SDP8 が組み込まれていて、メニューバーから選択することができます。

　SGPでは、NASA２ラインフォーマットに含まれる軌道要素と元期からの経過時間を与えると、その時点における衛星の位置と速度ベクトルを求めることができます。　座標は地心赤道座標系なので、計算時点のグリニッジ恒星時を求めて座標回転を行えば地球に固定した地心直交座標系に変換することができます。

　それでは、パソコンの内部時計を使って１秒ごとに衛星の位置と速度のベクトルを求め、表示するプログラムを使ってみましょう。

Calsat32SGP.zipをダウンロードして、解凍して下さい。

　解凍したフォルダに含まれる"Calsat32SGP.exe"を起動すると、つぎのようなフォームが表示されます。





　

　また、"Calsat32SGP.vbp"をクリックすると Visual BASIC が起動して、つぎのようなコードを見ることができます。


　

Private Sub cmdExit_Click()


    End
　

End Sub
　


Private Sub cmdStart_Click()
　

    Timer1.Enabled = True
　

End Sub
　


Private Sub Timer1_Timer()

    Static vntDateBefor As Variant, dblSmjd As Double
    Dim dblDateTime As Double, vntDateNow As Variant, vntDateTime As Variant
    Dim dblCmjd As Double, dblCelem(12) As Double, dblTsince As Double
    Dim dblX As Double, dblY As Double, dblZ As Double
    Dim dblXdot As Double, dblYdot As Double, dblZdot As Double


'現在時刻を取得し、時刻が変化していなければ、つぎのタイマー割り込みまで待つ
    vntDateNow = Now
    If vntDateBefor = vntDateNow Then Exit Sub
    vntDateBefor = vntDateNow


'シリアルデイトをＭＪＤに変換
    dblDateTime = vntDateNow
    dblCmjd = dblDateTime + 15018#
    If dblSmjd = 0 Then dblSmjd = dblCmjd


    'dblCelem(0)　元期のＭＪＤ
    'dblCelem(1)　西暦
    'dblCelem(2)　通算日
    'dblCelem(3)　軌道傾斜角(rad)
    'dblCelem(4)　昇交点赤経
    'dblCelem(5)　離心率
    'dblCelem(6)　近地点引数
    'dblCelem(7)　平均近点角
    'dblCelem(8)　平均運動
    'dblCelem(9)　平均運動の１次項
    'dblCelem(10) 平均運動の２次項（乗数を含む）
    'dblCelem(11) Ｂ* drag term（乗数を含む）

'テスト用軌道要素（元期は計算開始時点のパソコンの時刻）
    dblCelem(0) = dblSmjd
    dblCelem(3) = 72.8435
    dblCelem(4) = 115.9689
    dblCelem(5) = 0.0086731
    dblCelem(6) = 52.6988
    dblCelem(7) = 110.5714
    dblCelem(8) = 16.05824518
    dblCelem(9) = 0.00073094
    dblCelem(10) = 0.00013844
    dblCelem(11) = 0.000066816


'dblTsinceは元期からの経過時間（分）
'dblCmjdは計算したい時点のMJD
'    dblCelem(0)は軌道要素の元期のMJD
    dblTsince = (dblCmjd - dblCelem(0)) * 1440#


'ＳＧＰ計算モデルの呼び出し（赤道座標系）
'戻り値 X,Y,Zは位置ベクトル(km)、Xdot,Ydot,Zdotは速度ベクトル(km/s)
    Call SGP(dblTsince, dblCelem(), dblX, dblY, dblZ, dblXdot, dblYdot, dblZdot)


'計算結果の表示
    Cls
    Print
    Print " 元期からの経過秒数："; Round((dblCmjd - dblSmjd) * 86400#)
    Print
    Print " 数値は赤道座標系です。"
    Print
    Print " 位置ベクトル(km)"
    Print " X:"; dblX
    Print " Y:"; dblY
    Print " Z:"; dblZ
    Print
    Print " 速度ベクトル(km/s)"
    Print " XDOT:"; dblXdot
    Print " YDOT:"; dblYdot
    Print " ZDOT:"; dblZdot
    Print
    Print " ｸﾞﾘﾆｯｼﾞ恒星時を求めて地球に固定した座標に変換します。"


End Sub
　

　

Function ACTAN(SINX As Double, COSX As Double) As Double

    ACTAN = 0#
    If COSX = 0# Then
    If SINX = 0# Then Exit Function
        If SINX > 0# Then
            ACTAN = 1.57079633
        Else
            ACTAN = 4.71238898
        End If
        Exit Function
    End If

    If COSX > 0# Then
        If SINX = 0# Then Exit Function
        If SINX < 0# Then
            ACTAN = 6.2831853
        End If
    Else
        ACTAN = 3.14159265
    End If

    ACTAN = ACTAN + Atn(SINX / COSX)

End Function
　


Function FMOD2P(X As Double) As Double

    FMOD2P = X
    FMOD2P = FMOD2P - Int(FMOD2P / 6.2831853) * 6.2831853
    I f FMOD2P < 0# Then
        FMOD2P = FMOD2P + 6.2831853
    End If

End Function
　


Sub SGP(TSINCE As Double, ELEM() As Double, X As Double, Y As Double, Z As Double, XDOT As Double, YDOT As Double, ZDOT As Double)


'ＳＧＰで使う変数の宣言
    Dim XMO As Double, XNODEO As Double, OMEGAO As Double, EO As Double, XINCL As Double
    Dim XNO As Double, XNDT2O As Double, XNDD6O As Double, BSTAR As Double
    Dim CK2 As Double, CK4 As Double, E6A As Double, QOMS2T As Double, S As Double, TOTHRD As Double
    Dim XJ3 As Double, XKE As Double, XKMPER As Double, XMNPDA As Double, AE As Double
    Dim DE2RA As Double, PI As Double, PIO2 As Double, TWOPI As Double, X3PIO2 As Double
    Dim QO As Double, SO As Double, XJ2 As Double, XJ4 As Double
    Dim A1 As Double, DEL1 As Double, AO As Double, DELO As Double, XNODP As Double
    Dim C1 As Double, C2 As Double, C3 As Double, C4 As Double, C5 As Double, C6 As Double
    Dim COSIO As Double, SINIO As Double
    Dim XLO As Double, PO2NO As Double, OMGDT As Double, XNODOT As Double
    Dim D1 As Double, PO As Double
    Dim D1O As Double, D2O As Double, D3O As Double, D4O As Double
    Dim A As Double, E As Double, P As Double
    Dim XNODES As Double, OMGAS As Double, XLS As Double
    Dim AXNSL As Double, AYNSL As Double, XL As Double
    Dim ITEM3 As Integer
    Dim U As Double, EO1 As Double
    Dim TEM2 As Double, TEM5 As Double, SINEO1 As Double, COSEO1 As Double
    Dim ECOSE As Double, ESINE As Double, EL2 As Double, PL As Double, PL2 As Double
    Dim R As Double, RDOT As Double, RVDOT As Double
    Dim SINU As Double, COSU As Double, SU As Double
    Dim SINUK As Double, COSUK As Double, SINNOK As Double, COSNOK As Double
    Dim SINIK As Double, COSIK As Double, XMX As Double, XMY As Double
    Dim UX As Double, UY As Double, UZ As Double, VX As Double, VY As Double, VZ As Double
    Dim SIN2U As Double, COS2U As Double, RK As Double, UK As Double
    Dim XNODEK As Double, XINCK As Double
    Dim TEMP As Double


'定数
        DE2RA = 0.0174532925
        E6A = 0.000001
        PI = 3.14159265
        PIO2 = 1.57079633
        TOTHRD = 0.66666667
        TWOPI = 6.2831853
        X3PIO2 = 4.71238898
        XJ3 = -0.00000253881
        XKE = 0.0743669161
        XKMPER = 6378.135
        XMNPDA = 1440#
        AE = 1#
        QO = 120#
        SO = 78#
        XJ2 = 0.001082616
        XJ4 = -0.00000165597
        CK2 = 0.5 * XJ2 * AE ^ 2#
        CK4 = -0.375 * XJ4 * AE ^ 4#
        QOMS2T = ((QO - SO) * AE / XKMPER) ^ 4#
        S = AE * (1# + SO / XKMPER)


'軌道要素を配列から変数にコピーしてラジアンに変換
    'ELEM(0)　元期のＭＪＤ
    'ELEM(1)　西暦
    'ELEM(2)　通算日
    'ELEM(3)　軌道傾斜角(rad)
        XINCL = ELEM(3)
    'ELEM(4)　昇交点赤経
        XNODEO = ELEM(4)
    'ELEM(5)　離心率
        EO = ELEM(5)
    'ELEM(6)　近地点引数
        OMEGAO = ELEM(6)
    'ELEM(7)　平均近点角
        XMO = ELEM(7)
    'ELEM(8)　平均運動
        XNO = ELEM(8)
    'ELEM(9)　平均運動の１次項
        XNDT2O = ELEM(9)
    'ELEM(10) 平均運動の２次項（乗数を含む）
        XNDD6O = ELEM(10)
    'ELEM(11) Ｂ* drag term（乗数を含む）
        BSTAR = ELEM(11)

        XNODEO = XNODEO * DE2RA
        OMEGAO = OMEGAO * DE2RA
        XMO = XMO * DE2RA
        XINCL = XINCL * DE2RA
        TEMP = TWOPI / XMNPDA / XMNPDA
        XNO = XNO * TEMP * XMNPDA
        XNDT2O = XNDT2O * TEMP
        XNDD6O = XNDD6O * TEMP / XMNPDA


'下ごしらえ
        A1 = (XKE / XNO) ^ TOTHRD
        TEMP = 1.5 * CK2 * (3# * Cos(XINCL) ^ 2# - 1#) / (1# - EO * EO) ^ 1.5
        DEL1 = TEMP / (A1 * A1)
        AO = A1 * (1# - DEL1 * (0.5 * TOTHRD + DEL1 * (1# + 134# / 81# * DEL1)))
        DELO = TEMP / (AO * AO)
        XNODP = XNO / (1# + DELO)
        BSTAR = BSTAR / AE
        C1 = CK2 * 1.5
        C2 = CK2 / 4#
        C3 = CK2 / 2#
        C4 = XJ3 * AE ^ 3# / (4# * CK2)
        COSIO = Cos(XINCL)
        SINIO = Sin(XINCL)
        A1 = (XKE / XNO) ^ TOTHRD
        D1 = C1 / A1 / A1 * (3# * COSIO * COSIO - 1#) / (1# - EO * EO) ^ 1.5
        AO = A1 * (1# - 1# / 3# * D1 - D1 * D1 - 134# / 81# * D1 * D1 * D1)
        PO = AO * (1# - EO * EO)
        QO = AO * (1# - EO)
        XLO = XMO + OMEGAO + XNODEO
        D1O = C3 * SINIO * SINIO
        D2O = C2 * (7# * COSIO * COSIO - 1#)
        D3O = C1 * COSIO
        D4O = D3O * SINIO
        PO2NO = XNO / (PO * PO)
        OMGDT = C1 * PO2NO * (5# * COSIO * COSIO - 1#)
        XNODOT = -2# * D3O * PO2NO
        C5 = 0.5 * C4 * SINIO * (3# + 0.5 * COSIO) / (1# + COSIO)
        C6 = C4 * SINIO


'永年項の補正と大気の抗力
        A = AO * (XNO / (XNO + (2# * XNDT2O + 3# * XNDD6O * TSINCE) * TSINCE)) ^ TOTHRD
        E = E6A
        If A > QO Then E = 1# - QO / A
        P = A * (1# - E * E)
        XNODES = XNODEO + XNODOT * TSINCE
        OMGAS = OMEGAO + OMGDT * TSINCE
        XLS = FMOD2P(XLO + (XNO + OMGDT + XNODOT + (XNDT2O + XNDD6O * TSINCE) * TSINCE) * TSINCE)


'長期項
        AXNSL = E * Cos(OMGAS)
        AYNSL = E * Sin(OMGAS) - C6 / P
        XL = FMOD2P(XLS - C5 / P * AXNSL)


'ケプラーの方程式
        U = FMOD2P(XL - XNODES)
        ITEM3 = 0
        EO1 = U
        TEM5 = 1#
KepLoop:
        SINEO1 = Sin(EO1)
        COSEO1 = Cos(EO1)
        If Abs(TEM5) < E6A Then GoTo KepExit
        If (ITEM3 >= 10) Then GoTo KepExit
        ITEM3 = ITEM3 + 1
        TEM5 = 1# - COSEO1 * AXNSL - SINEO1 * AYNSL
        TEM5 = (U - AYNSL * COSEO1 + AXNSL * SINEO1 - EO1) / TEM5
        TEM2 = Abs(TEM5)
        If TEM2 > 1# Then TEM5 = TEM2 / TEM5
        EO1 = EO1 + TEM5
        GoTo KepLoop
KepExit:


'短期項の大きさ
        ECOSE = AXNSL * COSEO1 + AYNSL * SINEO1
        ESINE = AXNSL * SINEO1 - AYNSL * COSEO1
        EL2 = AXNSL * AXNSL + AYNSL * AYNSL
        PL = A * (1# - EL2)
        PL2 = PL * PL
        R = A * (1# - ECOSE)
        RDOT = XKE * A ^ 0.5 / R * ESINE
        RVDOT = XKE * PL ^ 0.5 / R
        TEMP = ESINE / (1# + (1# - EL2) ^ 0.5)
        SINU = A / R * (SINEO1 - AYNSL - AXNSL * TEMP)
        COSU = A / R * (COSEO1 - AXNSL + AYNSL * TEMP)
        SU = ACTAN(SINU, COSU)


'短期項補正
        SIN2U = (COSU + COSU) * SINU
        COS2U = 1# - 2# * SINU * SINU
        RK = R + D1O / PL * COS2U
        UK = SU - D2O / PL2 * SIN2U
        XNODEK = XNODES + D3O * SIN2U / PL2
        XINCK = XINCL + D4O / PL2 * COS2U


'ベクトルの座標回転
        SINUK = Sin(UK)
        COSUK = Cos(UK)
        SINNOK = Sin(XNODEK)
        COSNOK = Cos(XNODEK)
        SINIK = Sin(XINCK)
        COSIK = Cos(XINCK)
        XMX = -SINNOK * COSIK
        XMY = COSNOK * COSIK
        UX = XMX * SINUK + COSNOK * COSUK
        UY = XMY * SINUK + SINNOK * COSUK
        UZ = SINIK * SINUK
        VX = XMX * COSUK - COSNOK * SINUK
        VY = XMY * COSUK - SINNOK * SINUK
        VZ = SINIK * COSUK


'位置ベクトル（km）
        X = (RK * UX * XKMPER / AE)
        Y = (RK * UY * XKMPER / AE)
        Z = (RK * UZ * XKMPER / AE)


'速度ベクトル（km/s）
        XDOT = (RVDOT * VX + RDOT * UX) * XKMPER / AE * XMNPDA / 86400#
        YDOT = (RVDOT * VY + RDOT * UY) * XKMPER / AE * XMNPDA / 86400#
        ZDOT = (RVDOT * VZ + RDOT * UZ) * XKMPER / AE * XMNPDA / 86400#


End Sub


　上のコードには、SGP計算モデル(SGP)、アークタンジェント求める関数(ACTAN)、角度ラジアンの２π以下の部分を求める関数(FMOD2P)と１秒ごとに衛星の位置を計算して数値を表示するメインフォームが含まれています。　経過時間は計算開始ボタンを押した時点からの経過秒数で、計算に使っている軌道要素は計算モデルの評価用のものです。　

　


　すでに、グリニッジ恒星時の求め方と座標の回転、観測点の緯度経度から地心直交座標への変換、観測点から見た方位、仰角、距離など基本的なもの、軌道要素のファイルからの読み込み、衛星の可視範囲の表示など関連的なものについて具体的なプログラム例を示していますので、上のSGPを使ったプログラム例と組み合わせると時々刻々変化する衛星の位置を求めることができます。

　これまでの説明で、アマチュア的な実用精度での位置計算について必要なことはすべて網羅されています。　すでにプログラムを作り始めている方や完成までにいま一歩の方も「ＣＡＬＳＡＴ３２のからくり」をご覧いただいたことをきっかけに、ＣＡＬＳＡＴ３２を越えるプログラムを開発して頂くことを期待します。