衛星の方位と仰角

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

観測点から衛星は何処に見えるか

 すでにケプラーの方程式を解いて衛星の軌道面内の座標を計算する方法。 また、観測点の位置を経緯度からG系地心直交座標に変換する方法。 さらには、衛星の軌道面の座標をG系地心直交座標に変換するために必要なグリニッジ平均恒星時を求める方法についても理解が進んだところです。 最終的に知りたい観測点からみた衛星の方位仰角は、これまでの知識と座標の回転から計算で求めることができます。 

 衛星の位置(空間に固定された座標)と観測点の位置(地球に固定された座標)を結びつけるには、お互いを共通の座標系に変換する必要があります。 ここでは共通の座標系としてG系地心直交座標にお互いの位置を変換します。 すでに観測点の位置をG系地心直交座標に変換する方法については具体的に示していますので、ここでは衛星軌道面のXY座標をG系地心直交座標に変換する手順について示します。

 

 行うことは座標の回転です。 はじめに衛星の軌道面内のXY座標 Xk、Yk を軌道要素から与えられる昇交点赤経 Omg近地点引数 w軌道傾斜角 Inc を使って座標を回転して、地心赤道直交座標 Xs、Ys、Zs に変換します。 これらの関係はつぎのようになります。

Xs = Xk * (cos Omg * cos W - sin Omg * cos Inc * sin W) - Yk * (cos Omg * sin W + sin Omg * cos Inc * cos W)
Ys = Xk * (sin Omg * cos W +cos Omg * cos Inc * sin W) - Yk * (sin Omg * sin W  - cos Omg * cos Inc  * cos W)
Zs = Xk * sin Inc * sin W + Yk * sin Inc *cos W

 

 さらに、グリニッジ平均恒星時 Tg を使って座標を回転して、G系地心直交座標 Xg、Yg、Zg に変換します。 これらの関係はつぎのようになります。

Xg = Xs * cos Tg + Ys * sin Tg
Yg = -Xs * sin Tg + Ys * cos Tg
Zg = Zs

 

 いろいろな座標系が出てきましたが、衛星と観測点の位置関係を知るためには、どれかの座標系でなくてはならいと言うことはないのですが、ここでは最終的に観測点から見た衛星の動きを求めることが目的なので、そのために便利なG系地心直交座標系で、相互の位置関係を計算することにしました。

 観測点を基準にした衛星の位置と距離は、G系地心直交座標に変換した衛星の位置 Xg、Yg、Zg と観測点の位置 Xmg、Ymg、Zmg から、観測点から見た衛星の位置 Xr、Yr、Zr と距離 Er とすれば、これらの関係はつぎのようになります。

Xr = Xg - Xmg
Yr = Yg - Ymg
Zr = Zg - Zmg
Er = (Xr ^ 2 + Yr ^ 2 + Zr ^ 2) ^ 0.5

 

 観測点からみた衛星の方向ベクトルと距離は簡単な計算で求めることができます。

 これまでの座標回転で距離 Er については具体的な値を求めることができましたが、方向についてはG系地心直交座標を観測点に移動したG系測心直交座標系(Xr、Yr、Zr)なので取り扱いに不便です。 一般的には衛星にアンテナを向ける場合には方位角 Eaz と仰角 Eel で表す地平座標系が便利です。 ここで、観測点の緯度 Ido と経度 Keido が与えられれば、これらの関係はつぎのようになります。

 

経度に対する座標回転(回転後の座標を X1、Y1、Z1 とする)

X1 = Xr * cos Keido + Yr * sin Keido
Y1 = Xr * -sin Keido + Yr * cos Keido
Z1 = Zr

 

緯度に対する座標回転(回転後の座標を X2、Y2、Z2 とする)

X2 = X1 * sin Ido - Z1 * cos Ido
Y2 = Y1
Z2 = X1 * cos Ido + Z1 * sin Ido

 

方位角(南が0で時計回りに2πまで)

Eaz = arctan (-Y2 / X2)

 ここで、X2 >= 0 であれば、第1象限、第4象限

      X2 < 0 であれば、第2象限、第3象限

 なお、北を方位0度とする場合には180度の補正が必要。


仰角(地平線が0で天頂のπ/2まで)

E1 = Z2 / Er (つぎの式がゴチャゴチャするので)
Eel = arctan (E1 / (-E1 * E1 + 1) ^ 0.5)

 

から、方位角仰角を求めることが できます。

 それでは、衛星の位置を設定して、観測点を移動しながら方位と仰角を計算して、フォームに表示するプログラムを使ってみましょう。

 

Calsat32AzEl.zipをダウンロードして、解凍してください。

 

 解凍したフォルダに含まれる"Calsat32AzEl.exe"を起動して、スライダーを左右に適当に動かすとつぎのような フォームが表示されます。

 

 

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


Private Sub mnuExit_Click()
    End
End Sub

Private Sub Slider1_Change()
    Call Disp
End Sub

Private Sub Slider2_Change()
    Call Disp
End Sub

 

 

Sub Disp()

 

    Dim dblEido As Double, dblEkeido As Double, dblEtakasa As Double
    Dim dblMido As Double, dblMkeido As Double, dblMtakasa As Double
    Dim dblEx As Double, dblEy As Double, dblEz As Double
    Dim dblMx As Double, dblMy As Double, dblMz As Double

    Dim dblX0 As Double, dblY0 As Double, dblZ0 As Double
    Dim dblX1 As Double, dblY1 As Double, dblZ1 As Double
    Dim dblX2 As Double, dblY2 As Double, dblZ2 As Double
    Dim dblEr As Double

    Dim dblTanEaz As Double, dblSinEel As Double
    Dim dblEaz As Double, dblEel As Double

    Dim dblPi As Double, dblRad As Double, dblDeg As Double
    Dim dblEa As Double


'エラー・トラップの設定
    On Error GoTo Atama_No_Maue

'πの値など
    dblPi = 4# * Atn(1#)
    dblRad = dblPi / 180#
    dblDeg = 180# / dblPi

'地球の半径
    dblEa = 6377#


'衛星直下点の緯度経度と高さを設定
    dblEido = Val(Form1.Text1.Text) * dblRad
    dblEkeido = Val(Form1.Text2.Text) * dblRad
    dblEtakasa = Val(Form1.Text3.Text) + dblEa

'衛星の位置を地心直交座標系に変換
    dblEx = dblEtakasa * Cos(dblEido) * Cos(dblEkeido)
    dblEy = dblEtakasa * Cos(dblEido) * Sin(dblEkeido)
    dblEz = dblEtakasa * Sin(dblEido)

'観測点の緯度経度を設定
    dblMido = -Val(Form1.Slider2.Value) * dblRad
    dblMkeido = Val(Form1.Slider1.Value) * dblRad
    dblMtakasa = dblEa
    Form1.Label11.Caption = Format(dblMido * dblDeg, "###0")
    Form1.Label12.Caption = Format(dblMkeido * dblDeg, "###0")

'観測点の位置を地球を完全な球として地心直交座標系に変換
    dblMx = dblMtakasa * Cos(dblMido) * Cos(dblMkeido)
    dblMy = dblMtakasa * Cos(dblMido) * Sin(dblMkeido)
    dblMz = dblMtakasa * Sin(dblMido)

'観測点からみた衛星の方向ベクトルと距離
    dblX0 = dblEx - dblMx
    dblY0 = dblEy - dblMy
    dblZ0 = dblEz - dblMz
    dblEr = (dblX0 ^ 2# + dblY0 ^ 2# + dblZ0 ^ 2#) ^ 0.5

'地平座標への変換
    dblX1 = dblX0 * Cos(dblMkeido) + dblY0 * Sin(dblMkeido)
    dblY1 = dblX0 * -Sin(dblMkeido) + dblY0 * Cos(dblMkeido)
    dblZ1 = dblZ0
    dblX2 = dblX1 * Sin(dblMido) - dblZ1 * Cos(dblMido)
    dblY2 = dblY1
    dblZ2 = dblX1 * Cos(dblMido) + dblZ1 * Sin(dblMido)

'方位角(南が0で時計回りに2πまで)
    dblTanEaz = -dblY2 / dblX2
    dblEaz = Atn(dblTanEaz)
    If dblX2 < 0# Then
        dblEaz = dblEaz + dblPi
    End If
    If dblEaz < 0# Then
        dblEaz = 2# * dblPi + dblEaz
    End If

'仰角
    dblSinEel = dblZ2 / dblEr
    dblEel = Atn(dblSinEel / (-dblSinEel * dblSinEel + 1#) ^ 0.5)

'計算結果を表示
    dblEaz = dblEaz + dblPi
    If dblEaz >= 2# * dblPi Then
        dblEaz = dblEaz - 2# * dblPi
    End If
    Form1.Text4.Text = Format(dblEaz * dblDeg, "###0.#")
    Form1.Text5.Text = Format(dblEel * dblDeg, "##0.#")

'正常終了
    Exit Sub

'エラー処理
    Atama_No_Maue:
    MsgBox ("観測点の真上に衛星が来たので、アークタンジェントの計算が出来ません。")

End Sub

 

 衛星の直下点の経緯度と地上からの高さは、テキストボックスに設定することができます。 スライドバーで観測点の位置を変更することができ、観測点から見た衛星の方位と仰角を表示します。 実際の運用では非常にまれなことですが、衛星が観測点の真上に位置した場合にはプログラム中のアークタンジェントの処理でエラーが生じるため、その旨の警告を表示します。 さまざまな条件を設定して衛星と観測点の位置から決まる方位と仰角の変化を確認してみましょう。

 

 

 観測点から見た衛星の方位と仰角を求める手順を知ることができました。 衛星の位置計算は特別に複雑な計算を必要とするものではありません。 手順を整理し、その意味を理解すれば、あなたにもできることです。 これまでに示した Visual BASIC のコードを組み合わせることで、軌道要素で与えられた元期の衛星位置を求めることが可能 なレベルに到達しています。