衛星の位置計算では必修科目、久しぶりの サイン と コサイン

極座標系と直交座標系

　ベクトルとは方向と大きさをもった量で、衛星の位置も地球の重心からの方向と距離、すなわちベクトルで表すことができます。　地図などで北の方角を基準にすれば、「東回りに１２０度の方向に３００ｍ進んだ場所」のように、角度と距離で位置を表現することができ、この方法は極座標系と呼ばれています。　一方、北の方向と、これから９０度ずれた東の方向に２つの直交する軸をとり、「北へ−２００ｍ、東へ１００ｍ進んだ場所」のように、２つの軸の大きさで方向と距離を表現する方法もあり、直交座標系と呼ばれています。

衛星軌道面内での位置の表し方

　衛星の位置計算でも何らかの方法で位置を示す必要があります。　衛星の軌道面内でケプラーの方程式を解くと、衛星が地球の重心からみて近地点方向から衛星が衛星の進む方向にどのくらいの角度にいるのかを示す真近点角と地球の重心から衛星までの距離が得られます。　この表現方法は極座標系によるものです。　衛星の位置計算では、最終的には観測点からみた方位、仰角、距離、距離変化などが必要になります。　この変換には３次元の座標変換が必要で、座標変換を行うためには極座標系ではなく、直交座標系で取り扱うことが便利です。　具体的には衛星軌道面での極座標系の表現を、地心から近地点方向をＸ軸、近地点から反時計回りに９０度方向をＹ軸として軌道面内の直交座標系に変換します。　ちなみに、軌道面内でみると衛星は平面上の楕円の上を運動するために、Ｚ軸の値は０になります。

座標変換

　軌道面内の衛星の位置は、軌道要素の近地点引数、軌道傾斜角、昇交点赤経を使って地心直交赤道座標に変換されます。　さらに、グリニッジ平均恒星時を使って地球に固定した座標（Ｇ系地心直交座標）に変換します。　以上の座標変換を行うと、経緯度・高さと直交座標で示した観測点位置のＧ系地心直交座標と共通の座標になります。

方向余弦

　直交座標系 Ｘ，Ｙ，Ｚ で原点からある方向に向いた直線が、

　　Ｘ軸の正の向きと作る角を ａ

　　Ｙ軸の正の向きと作る角を ｂ

　　Ｚ軸の正の向きと作る角を ｃ

とするとき、

　　Ｌ＝ｃｏｓ ａ

　　Ｍ＝ｃｏｓ ｂ

　　Ｎ＝ｃｏｓ c

で、表される３つの値を方向余弦（Ｌ，Ｍ，Ｎ）と呼んでいます。

　また、Ｌ，Ｍ，Ｎ の間には、

　　Ｌ＾２ ＋ Ｍ＾２ ＋ Ｎ＾２ ＝ １

の関係があります。

Ｘ軸を回転軸とした場合

　Ｘ軸の方向を固定して、Ｘ軸の正の方向からみて反時計回りに角度 θ で回転したときの新しい直交座標系を Ｘｘ，Ｙｘ，Ｚｘ とするとき、直線の向きの方向余弦 Ｌｘ，Ｍｘ，Ｎｘ は、回転前の方向余弦を Ｌ，Ｍ，Ｎ とすれば、

　　Ｌｘ ＝ Ｌ

　　Ｍｘ ＝ Ｍ ｃｏｓθ ＋ Ｎ ｓｉｎθ

　　Ｎｘ＝ −Ｍ ｓｉｎθ ＋ Ｎ ｃｏｓθ

から計算することができます。

Ｙ軸を回転軸とした場合

　Ｙ軸の方向を固定して、Ｙ軸の正の方向からみて反時計回りに角度 θ で回転したときの新しい直交座標系を Ｘｙ，Ｙｙ，Ｚｙ とするとき、直線の向きの方向余弦 Ｌｙ，Ｍｙ，Ｎｙ は、回転前の方向余弦を Ｌ，Ｍ，Ｎ とすれば、

　　Ｌｙ ＝ Ｌ ｃｏｓθ − Ｎ ｓｉｎθ

　　Ｍｙ ＝ Ｍ

　　Ｎｙ＝ Ｌ ｓｉｎθ ＋ Ｎ ｃｏｓθ

から計算することができます。

Ｚ軸を回転軸とした場合

　Ｚ軸の方向を固定して、Ｚ軸の正の方向からみて反時計回りに角度 θ で回転したときの新しい直交座標系を Ｘｚ，Ｙｚ，Ｚｚ とするとき、直線の向きの方向余弦 Ｌｚ，Ｍｚ，Ｎｚ は、回転前の方向余弦を Ｌ，Ｍ，Ｎ とすれば、

　　Ｌｚ ＝ Ｌ ｃｏｓθ ＋ Ｍ ｓｉｎθ

　　Ｍｚ ＝ −Ｌ ｓｉｎθ ＋ Ｍ ｃｏｓθ

　　Ｎｚ＝ Ｎ

から計算することができます。

　それでは、地球の重心から見た「ここは何処」を計算して、フォームに表示するプログラムを使ってみましょう。

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

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

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

Private Sub Form_Load()

'起動時の計算
    Call Conv

End Sub
　

Private Sub mnuExit_Click()

'プログラムの終了
    End

End Sub

　

Private Sub Slider1_Change()

'Ｘ軸の角度が変わったら再計算
    Call Conv

End Sub

Private Sub Slider2_Change()

'Ｙ軸の角度が変わったら再計算
    Call Conv

End Sub

Private Sub Slider3_Change()

'Ｚ軸の角度が変わったら再計算
    Call Conv

End Sub

　

Sub Conv()

Dim dblL As Double, dblM As Double, dblN As Double
Dim dblLx As Double, dblMx As Double, dblNx As Double
Dim dblLy As Double, dblMy As Double, dblNy As Double
Dim dblLz As Double, dblMz As Double, dblNz As Double
Dim dblXa As Double, dblYa As Double, dblZa As Double
Dim dblPi As Double

dblPi = 4# * Atn(1#)

'方向余弦の初期値
    dblL = 1#
    dblM = 0#
    dblN = 0#

'Ｘ軸を回転
    dblXa = Val(Form1.Slider1.Value) * dblPi / 180#
    dblLx = dblL
    dblMx = dblM * Cos(dblXa) + dblN * Sin(dblXa)
    dblNx = dblM * -Sin(dblXa) + dblN * Cos(dblXa)

'Ｙ軸を回転
    dblYa = Val(Form1.Slider2.Value) * dblPi / 180#
    dblLy = dblLx * Cos(dblYa) + dblNx * -Sin(dblYa)
    dblMy = dblMx
    dblNy = dblLx * Sin(dblYa) + dblNx * Cos(dblYa)

'Ｚ軸を回転
    dblZa = Val(Form1.Slider3.Value) * dblPi / 180#
    dblLz = dblLy * Cos(dblZa) + dblMy * Sin(dblZa)
    dblMz = dblLy * -Sin(dblZa) + dblMy * Cos(dblZa)
    dblNz = dblNy

'計算結果の表示
    Form1.Label11.Caption = Format(dblLz, "#0.########")
    Form1.Label13.Caption = Format(dblMz, "#0.########")
    Form1.Label15.Caption = Format(dblNz, "#0.########")

End Sub
　

　方向余弦の初期値としては、Ｌ＝１、Ｍ＝０、Ｎ＝０ が与えられています。　それぞれの軸を回転させるためのスライドバーを動かして、変換後の数値を確認してみましょう。　このプログラムでは、Ｘ、Ｙ、Ｚ軸の３軸についてスライドバーで与えられた角度の座標回転を行っています。　スライドバーはマウスで直接動かすか、マウスでクリックしてフォーカスされているスライドバーに対して左右の矢印ボタンで動かすことができます。

　衛星の位置計算で頻繁に出てくるのは座標の回転です。　３次元での座標回転なので一度に３軸を考えるとやっかいなことになります。　もちろん、一度に３軸の回転角を与えて座標回転する式を作ることはできますが、プログラミングでは１軸ごとに計算すると間違いは少なくなるかと思います。　また、座標回転では回転行列と呼ばれるものを使って行列式で書かれることが多いのですが、Ｖｉｓｕａｌ　ＢＡＳＩＣ では行列演算を直接サポートしていないので、ここで示した行列式を展開した形にするか、自分でサブルーチンを作成することになります。　衛星の位置計算は、座標の回転と、回転量を求める作業であると言っても言い過ぎではないと思います。　私は フレミングの右手、左手の法則のように、親指、人差し指、中指を直角に開いて手首を回転させ、理解を深めました。　端から見るとかなり変ですけれど。