mikeo_410

 人工衛星の軌道（アルマナック、2行軌道要素形式、エフェメリス）

  GPSの受信機が計測しているのは送信時間と受信時間の差で、光速から距離を計算するものだと知りました。時間はGPS衛星自体の時計によるメッセージへの打刻と、受信機の時計によって測られるようです。したがって、メッセージの内容自体を問題にするものではないようです。
  GPS受信機は計測された距離とGPS衛星の位置から、その球面の交点として、GPS受信機の位置を算出します。
  GPS衛星の位置は既定の軌道から計算されるもので、必ずしもGPS衛星自体から得る必要があるものではないようです。GPS衛星は、どの衛星も同じ内容の、全GPS衛星の軌道情報を放送するもののようです。
  GPS衛星がある時刻に、どこに在るか、あるいは在ったのかは、一義のもので、受信機の性能とは無関係な問題のようです。

  しかし、後述のように、アルマナックや2行軌道要素形式の情報から計算される衛星の位置は大まかなもののようです。パラメータは荒く、使用されている定数も3次元での位置を正確に表すには精度不足です。
  軌道情報は、GPS衛星に限ったものではなく、基本的には全人工衛星について提供されているもののようです。
  おそらく光学観測のために天球上の2次元的な位置を示すのには十分な精度なのだろうと思います。

  GPS衛星が個体の位置情報として放送しているのはエフェメリスと言う情報のようです。

  軌道情報にはドラグ（Drag）と呼ばれる項目があります。N1、大気抵抗、抗力、平均運動の変化率などと呼ばれるようです。
  近い軌道を回るほど、真空度が低く、影響を受けると言うことだと思います。
  半日周期で周回するGPS衛星は、遠い軌道を回る衛星です。

人工衛星の軌道

  人工衛星の軌道の概要は軌道半径で表されるようです。起動半径 r と角速度 ω は以下の関係にあるようです。
  r3 = g・R2 / ω2
  重力加速度 g は、g = 9.80665 m/s2
  地球の赤道半径 R は、R = 6378137 m

  静止衛星は地球の自転と同じ速さで回る必要があり、比較的遠い軌道を回ることになります。
  1日で一周するのは ω = 2・π / (24・60・60) と考えれば r = 42,253.112km と計算されます。
  地上から見て静止しているためには、赤道上空にあって、地球の自転と方向に進んでいなければなりません。また、地球の自転周期は、約23時間56分4.06秒のようです。3075m/s ぐらいの速さで進んでいることになります。

  GPS衛星は、当初24個の衛星で、地表の同じ点で常に4つのGPS衛星の電波が受信可能なように設計されたようです。単純には1つの衛星は6時間程度見えていることになります。したがって、GPS衛星も比較的遠いところを回っています。  GPS衛星の周期は11時間58分02秒のようです。r = 26,569.300km と計算されます。3875m/s ぐらいの速さです。

  これは概要で、実際の衛星は、地心を中心に真円を描いているのではないようです。

  1. 地心は準拠楕円体の自転軸と赤道面の交点ですが、質量の中心とは限らないものと思います。
  2. 地表での重力方向は一点を指さないものと思いますが、人工衛星についても同様だと思います。
  3. 人工衛星軌道も楕円で近似され、地球は焦点の1つと見るようです。

  大きく誇張して描くと衛星軌道は左図のようなイメージで捉えられているようです。
  地球は衛星軌道の焦点の一方の位置にあり、衛星の位置は近地点との角度で表されます。
  公転面は、地球の赤道面となす角と、交線の方向で表されます。

  軌道を楕円と考えるので、長半径 a と短半径 b があります。アルマナックでは、a と離心率 e で表されます。b2 = a2・（1-e2) の関係にあります。
  また、地球と軌道の中心との距離は a・e です。

平均近点角、真近点角、離心近点角

  アルマナックの軌道情報は、基準となる時刻の衛星の位置を平均近点角で与えています。平均近点角は図の γ です。∠NOU の 点U は、図では明瞭ではありません。
  衛星は点P にあり、α（∠NFP） は真近点角と呼ぶようです。
  また、楕円の媒介変数表記のβ（∠NOQ）は、離心近点角と呼ぶようです。

  点Uは、扇形NOUが扇形NFQと同じ面積になるように定めるもののようです。
  Fは楕円の焦点の1つだと考えられます。
  扇形NOUの面積は、円の面積の γ / 2π です。
  扇形NFQの面積は、扇形NOUの面積から、三角形QOFの面積を引いたものです。
  扇形NOUの面積は、円の面積の β / 2π です。三角形QOFの面積は、線分OFの長さがa・e、高さが a・sin β なので、 (a・e)(a・sin β) / 2 です。
  ここで、e は軌道となっている楕円の離心率で、アルマナックの軌道情報に含まれます。

  中心の重なった軌道の長半径 a を半径とする円の面積を S とすると、
  S(γ / 2π) = S(β / 2π) －(a・e)(a・sin β) / 2
  S = π・a2
  したがって、平均近点角 γ と離心近点角 βの間には γ = β - e・sin β の関係があります。
  γ に応じた β の値は、
  y = β
  y = γ + e・sin β
  の、交点となります。

  e・sin β は小さな値です。γ を初期値に漸化的に解くようです。

  1. > e <- 0.5140781403E-02  # 離心率
  2. > C <- 0.6236332231      # 平均近点角
  3. > d <- 1
  4. > B <- C
  5. > while(d>1e-16)
  6. + {
  7. +   d <- (C - B + e*sin(B)) / (1-e*cos(B))
  8. +   B <- B + d
  9. + }
  10. > B
  11. [1] 0.626647946455

  （a・cos β、b・sin β）は、楕円の媒介変数表示そのもので、点P（衛星の位置）を表します。
  b は、軌道の短半径で、b = √(a2・（1-e2) ) です。

アルマナック

  GPS衛星が放送したり、JAXAが提供する、衛星軌道を計算するための情報は、アルマナック（almanac）と呼ばれるようです。、

１つのGPS衛星のアルマナックの例
項目 値 備考
ID 02 GPS衛星の個体識別で、
PRN-02のように使われる
Health 015
Eccentricity 1.369953156E-02 離心率
Time of Applicability(s) 405504.0000 基準時刻（週初起点）
Orbital Inclination(rad) 0.9401468643 軌道傾斜角
Rate of Right Ascen(r/s) -8.1260527684E-09 昇交点赤経の時間変化率
SQRT(A)  (m 1/2) 5153.601074 軌道長半径の平方根
Right Ascen at Week(rad) -3.105959433E+00 週初昇交点経度
Argument of Perigee(rad) -2.397643430 近地点離角
Mean Anom(rad) 6.846340514E-01 平均近点角
Af0(s) 5.121231079E-04
Af1(s/s) 3.637978807E-12
week 781 1999年8月22日起点

  JAXAが提供しているアルマナックのファイルは、PRN-01 から PRN-32までのGPS衛星の情報が含まれています。
  どの衛星も同じ項目からなります。
  図は、1つのGPS衛星のアルマナックの例です。

  公転面上の軌道は、 長半径の平方根（SQRT(A)）と離心率（Eccentricity）で与えられます。

  衛星の公転面と地球の関係は、公転面と赤道面となす角「軌道傾斜角」、週始めにおける昇交点経度「Right Ascen at Week」、昇交点赤経の時間変化率（Rate of Right Ascen）で与えられます。

  軌道上の衛星の位置は、Time of Applicability の時刻の衛星の位置で与えられます。この時刻は、週初（日曜日の午前0時）を起点にした経過時間です。
  この時刻に衛星は軌道上の平均近点角（Mean Anom）の位置にありました。

軌道の計算

  軌道を計算すると言うことは、目的によって、いろいろなことを指すと思います。GPS受信機では、ある時点の衛星の座標をを特定することです。

計算に使用される定数
定数 値 備考
地球の自転速度 7.292115×10-5
7.2921151467E-05		 rad / s
(2π / 23時間56分4.098904秒)
地心重力定数
GM		 3.986005E+14
3.986004418E+14
3986004.36e8		 (m3/s2)、GRS80

  準拠楕円体上の経度、地理緯度を得ることがGPSの目的なので、地心を原点とし、自転軸をZ軸（北極方向）、直交する経度ゼロ方向をX軸、東経90°方向をY軸とした直交座標で考えます。

  地心は、GPS衛星の楕円軌道の焦点の一方だとして考えます。ただし、常時測位に使う衛星なので、その軌道は、ほぼ円だと思います。
  地心は重心ではありませんし、いろいろな場所での重力方向が1点に集まるものでもないと思いますが、地心が衛星軌道の焦点にあるとして考えます。

  まず、公転面上の軌道を考えます。
  軌道の長半径 a は、アルマナック SQRT(A) の2乗で、前述のPRN-02の例の場合は、a = (5153.601074)2 = 26559604.0299 (m) です。
  離心率 e は、e = 1.369953156E-02 なので、短半径 b は、b = 26557111.5974 (m) で、差は2.5km程度のようです。
  軌道である楕円の中心を原点、遠地点から近地点へ向かう直線をX軸、直交する軸をY軸と考えます。地心は（e・a、0）にあります。
  前述のように、平均近点角（Mean Anom）6.846340514E-01 から離心近点角 β を計算すると、0.693390081049 となります。
  基準時刻をTaとすると、Ta = 405504.0000 秒における衛星の位置は、（a・cos β、b・sin β）=（20426557.865、16973927.5775）となります。
  Ta と同じ表記で表した時刻 T の衛星の位置は、平均近点角に（T－Ta）・√(GM / a3)を加算して計算します。 

  衛星の公転面上の座標を、地心を中心としたXYZ直交座標で表すので、原点を軌道の中心から地心に相当する焦点に移動します。
  （20426557.865 - e・a、16973927.5775）=（20062703.7314 、16973927.5775）

  次に、衛星の公転面と、地球の赤道面の関係を考えます。
  公転面と赤道面の傾きは、軌道傾斜角（Orbital Inclination）i で、例では i = 0.9401468643 となっていて、約53.866° です。
  公転面のX軸と赤道面のなす角 ω は、近地点離角（Argument of Perigee）です。
  公転面と赤道面の交線の方向を示す Ω は、週初昇降点経度（Right Ascen at Week）のようです。

  i と ω は、公転面と赤道面の関係ですが、Ω は、XYZ直交座標のX軸からの角度で、地球の自転に従って変化します。また、長い周期の歳差運動があるものと思います。
  前者の変化率は「地球の自転速度」によって角速度で与えられます。これは、最後に調整します。
  後者は、昇交点赤経の時間変化率（Rate of Right Ascen）です。例では -8.1260527684E-09 と言う値で、その逆数が1周する時間になります。1424日以上の長い周期で自転方向とは逆に回っているようです。
  Ω  は、「週初」ではなく、他の値と同様、基準時刻（Time of Applicability）Ta における値として使われているようです。
  時刻 T の Ω は、週初昇降点経度（Right Ascen at Week）を R 、昇交点赤経の時間変化率（Rate of Right Ascen）を ΔR  とすると、
  Ω = R + ΔR(T - Ta)
  で計算されます。

  1. > x0 <- 20062703.7314 
  2. > y0 <- 16973927.5775
  3. > w <- -2.397643430 
  4. > i <- 0.9401468643
  5. > o <- -3.105959433
  6. > sn <- sin(w); cs <- cos(w)
  7. > mW <- matrix(c(cs,-sn,0,
  8. +                sn, cs,0,
  9. +                 0,  0,1),ncol=3,byrow=T)
  10. > sn <- sin(i); cs <- cos(i)
  11. > mI <- matrix(c(1, 0,  0,
  12. +                0,cs,-sn,
  13. +                0,sn, cs),ncol=3,byrow=T)
  14. > sn <- sin(o); cs <- cos(o)
  15. > mO <- matrix(c(cs,-sn,0,
  16. +                sn, cs,0,
  17. +                 0,  0,1),ncol=3,byrow=T)
  18. > x1 <- mO%*%(mI%*%(mW%*%matrix(c(x0,y0,0))))
  19. > x1
  20.                [,1]
  21. [1,]   2717548.6544
  22. [2,]  15482778.6040
  23. [3,] -21060028.3578
  24. > rE <- 7.2921151467E-05 # 地球自転角速度
  25. > rE <- 7.292115E-05
  26. > t <- rE * 405504       # 基準時刻
  27. > sn <- sin(-t); cs <- cos(-t)
  28. > mT <- matrix(c(cs,-sn,0,
  29. +                sn, cs,0,
  30. +                 0,  0,1),ncol=3,byrow=T)
  31. > mT %*% x1
  32.                 [,1]
  33. [1,] -15638461.46633
  34. [2,]  -1593746.20221
  35. [3,] -21060028.35775

  公転面での座標を、Z軸に ω 、X軸に i 、Z軸に Ω 回転します。
  例の基準時刻の値は、
  ω = -2.397643430 、
  i =  0.9401468643、
  Ω = -3.105959433 です。
  （20062703.7314 、16973927.5775）は、
  （2717548.654、15482778.604、-21060028.358）になります。

  最後に、地球の自転による影響を調整します。天球で考えると動いていない昇交点も、経度で表すと時点に従って回ります。
  その量は、（地球の自転の角速度 x 週初からの経過時間）で、自転と逆向きになります。 

  週初昇降点経度（Right Ascen at Week）の値も、この変換が行われます。週初昇降点経度の値は、地球の自転に対する補正が行われる前の、基準時刻の経度だと言うことになります。
  そのままで、基準時刻の経度ではないので、「週初からの」と言う名前になったようです。

  回転の結果は（-15638461.466、1593746.202、21060028.358）となりますが、これは地球の自転の角速度を 7.292115E-05 とした場合です。
  地球の自転周期は 23時間56分4.098 903 691秒 のようです。これで 2πを割った値を使うと、
  （-15638462.41445、-1593736.89885、-21060028.35775）
と、なります。
  大きな差がありますが、おそらく、前者が予定された値です。

2行軌道要素形式

  アルマナックはGPSを調べて知りましたが、衛星軌道情報は伝統的に「2行軌道要素形式（Two-Line Element Sets）」で提供されてきたようです。パンチカード2枚1組でデータを入力していたことから来ているようです。カードにカラム固定で数字を並べた形式です。
  現在は、24桁の衛星名の行を先行させた3行で1組のようです。

  1. GPS BIIR-13 (PRN 02)    
  2. 1 28474U 04045A   14227.20083466 -.00000016  00000-0  00000+0 0  9350
  3. 2 28474  53.8406 140.3425 0138399 223.2946  45.1386  2.00564675 71700

  軌道を表す要素は下表のようです。

2行軌道要素とアルマナック
項目 2行軌道要素 アルマナック
行 桁 長さ 内容
元期 1 19 14 2桁の年、12桁の日時。通日で、
時間も日が単位。分解能は0.864ms。		 week、Time of Applicability
軌道傾斜角 2 9 8 度。小数点以下4桁。小数点を含む。 Orbital Inclination
昇交点の赤経 18 8 度。小数点以下4桁。小数点を含む。
離心率 27 7 先頭が小数点以下1桁目の値 Eccentricity
近地点引数 35 8 度。小数点以下4桁。小数点を含む。
平均近点角 44 8 度。小数点以下4桁。小数点を含む。
平均運動 53 11 回 / 日。
通産周回数 54 5 昇交点通過回数

  元期は軌道要素がいつの時点のものかを示します。
  2行軌道要素の元期の例は、14228.82639742 のような14桁の文字列です。最初の2桁は 2014年を示します。続く12桁は年初からの日数です。例は228日と 0.82639742 日を示します。 228日目は8月16日です。0.82639742日は、一日の秒数を乗じて、
  0.82639742 ・.86400 = 71400.737088 （秒）
  で、19時間50分と解釈できます。
  アルマナックの場合は、1999年8月22日を起点とした週の数（week）と、その週初からの経過秒数（Time of Applicability）で元期が示されます。2行軌道要素の例と同じ元期は、 781 週と 589824秒と表されます。  781 週の週初は、2014年8月10日のようです。

  1. > T <- as.Date("1999-08-22")
  2. > T + (781*7)
  3. [1] "2014-08-10"

  経過秒数は、589824 = 6・86400 + 19*3600 + 50・60 + 24 で、週初から6日と19時間50分24秒であることが分かります。

 アルマナックのSQRT(A) のような軌道半径の直接的な表記がありません。これは平均運動から計算されるようです。
  平均運動 N は、衛星が日に何回昇交点を通過するかを示す値のようです。1日は24時間のようです。
  1秒間の角速度 n は、
  n = 2πN / 24時間 = (2π / 86400)・N
  長半径 a と n は、n = √（GM / a3）の関係にあります。地心重力定数 GM は、GM = 3.986005e+14 です。
  (2π / 86400)2 ・N2 = GM / a3 
  a3 =  (GM / (2π / 86400)2 )・(1 / N2 )= 7.5371227345e+22 / N2
  a = 42241097.73 / N2/3

  アルマナックでは、公転面（軌道面）と赤道面の交線の方向を、Right Ascen at Week で示します。2行軌道要素は昇交点の赤経で示します。前者は経度と説明され、後者は赤経です。
  衛星軌道を示す目的では赤経は妥当です。赤経のゼロ度は天球の春分点で、地球の自転によって移動するものではありません。
  経度を使って天球の決まった方向を表すと、時間の関数になります。アルマナックは、Right Ascen at Week （R）と 「地球の自転速度」（rE）の2つの値を使って時間の関数として表しているのだと思います。Ω(t) = R  + rE・t 
  1つの値で示す方が優れていますが、地表の場所と関連付ける必要があるなら大差はないことになります。
  2行軌道要素は赤経で与えられるので、経度との換算に関する情報は含まれません。この情報はグリニッジ平均恒星時を計算することで解決されるようです。グリニッジ恒星時

  グリニッジ恒星時がゼロと言うことは、赤経のゼロと、経度のゼロが一致していることを示し、天球と地球のXYZ直交座標が重なることです。恒星時は2つの座標系のZ軸の回りの回転量を表すことになります。

 エフェメリス

  エフェメリス（Ephemeris）は、GPS衛星が送信している、自衛星の、より正確な位置情報のようです。

エフェメリス・メッセージ
記号 名称 bits
IODE Ephemeris Issue of Data. Used to uniquely identify the current set of navigation data. 8
Crs Amplitude of the Sine Harmonic Correction Term to the Orbit Radius 16
Δn Mean motion difference from computed value 16
M0 Mean Anomaly at Reference Time 32
Cuc Amplitude of the Cosine Harmonic Correction Term to the Argument of Latitude 16
e 離心率 32
Cus Amplitude of the Sine Harmonic Correction Term to the Argument of Latitude 16
√A 長半径の平方根 32
Toe Reference time of Ephemeris 16
Fit interval flag 1
AODO 5
Cic Amplitude of the Cosine Harmonic Correction Term to the Angle of Inclination 16
Ω0 Longitude of Ascending Node of Orbit Plane at Weekly Epoch 32
i0 Inclination Angle at Reference Time 32
Crc Amplitude of the Cosine Harmonic Correction Term to the Orbit Radius 16
ω Argument of Perigee 32
Ωdot Rate of Right Ascension 24
IODE 8
IDOT Rate of inclination angle 14

  エフェメリスはGPS衛星のメッセージの2サブフレームの情報です。サブフレームは300ビットで、30ビット/ワードで扱われるようです。ワードの有効なビット数は24のようです。サブフレームの最初の2ワードはテリメトリ・ワード（TLM）、ハンドオーバ・ワード（HOW）に決められているので、エフェメリス情報は16ワード、384ビットの情報と言うことになります。含まれる数値情報は32ビットで格納されており、ワードを跨いで格納されているようです。
  左表の順に項目が割り当てられ、19項目、364ビットが使用されているようです。

  しかし、実際のエフェメリスの受信データを見る方法は分からないので、JAXAのサイトで配布されているファイルを見るよりありません。このファイルは、GPS衛星32個全ての分を格納しています。1つの衛星（PRN-22）の分は下図のようです。
  7個の有効桁の短い数値と、29個の有効桁の長い数値から成り立っています。

  最初の22はGPS衛星番号で PRN-22 のデータであることを表しています。
  続く6つのデータは日時を表しています。(20)14年8月17日23時59分44.0秒です。この エフェメリスの元期です。
  この情報は2つのサブフレームには含まれていない情報のようです。

  日時に続く数値は、順に p1 から番号を付けて指し示すことにします。
  p12 =  8.638400000000E+04 は、86384 で Toe だと思います。
  p11 = 5.153723068237E+03 は、√A で、軌道の長半径は 26560861.5 m です。
  p9 = 7.382828858681E-03 は、離心率 e で、短半径は 26560137.6 m と言うことになります。 

  1. 22 14  8 17 23 59 44.0 2.598208375275D-04 3.524291969370D-12 0.000000000000D+00
  2.     3.000000000000D+00 2.750000000000D+00 5.687022601604D-09 8.875469507845D-01
  3.     2.216547727585D-07 7.382828858681D-03 8.562579751015D-06 5.153723068237D+03
  4.     8.638400000000D+04-3.725290298462D-08-2.165875207385D+00 1.136213541031D-07
  5.     9.232751538875D-01 1.972500000000D+02-2.049857990253D+00-8.635716855127D-09
  6.    -2.996553389905D-10 1.000000000000D+00 1.806000000000D+03 0.000000000000D+00
  7.     2.000000000000D+00 0.000000000000D+00-1.722946763039D-08 3.000000000000D+00
  8.     7.923000000000D+04 4.000000000000D+00

  エフェメリスからの軌道計算の例は以下のR言語のスクリプトのように行われるようです。
  p5からp20までの15のパラメータが使用されています。

  1. >   TY <- 2014; TM <- 8; TD <- 17;  # ターゲット時刻の年、月、日
  2. >   Th <- 23; Tm <- 59; Ts <- 44;   # 時、分、秒
  3. >   # エフェメリス・ファイルのデータの例
  4. >   prn <- 22 # 衛星番号
  5. >   Y <- 14; M <- 8; D <- 17;  # 年の下2桁、月、日
  6. >   h <- 23; m <- 59; s <- 44; # 時、分、秒
  7. >   p5 <- 2.750000000000E+00   # Crs
  8. >   p6 <- 5.687022601604E-09   # Mean motion difference from computed value
  9. >   p7 <- 8.875469507845E-01   # M0. Mean Anomaly at Reference Time
  10. >   p8 <- 2.216547727585E-07   # Cuc
  11. >   p9 <- 7.382828858681E-03   # e. 離心率
  12. >   p10 <-  8.562579751015E-06 # Cus
  13. >   p11 <-  5.153723068237E+03 # Square Root of Semi-Major Axis
  14. >   p12 <-  8.638400000000E+04 # Toe. Reference time of Ephemeris
  15. >   p13 <- -3.725290298462E-08 # Cic
  16. >   p14 <- -2.165875207385E+00 # Ω0. Longitude of Ascending Node of 
  17. >                              #     Orbit Plane at Weekly Epoch
  18. >   p15 <-  1.136213541031E-07 # Cis
  19. >   p16 <-  9.232751538875E-01 # i0. Inclination Angle at Reference Time
  20. >   p17 <-  1.972500000000E+02 # Crc
  21. >   p18 <- -2.049857990253E+00 # ω. Argument of Perigee
  22. >   p19 <- -8.635716855127E-09 # Ωd. Rate of Right Ascension
  23. >   p20 <- -2.996553389905E-10 # Id. Rate of inclination angle
  24. >   # 定数
  25. >   jd0 <- julian(as.Date("1980-1-6"),-2440588) - 0.5 # GPS時間の起点
  26. >   GM <- 3986005e8  # 地心重力定数（m^3/s^2)
  27. >   W <- 7292115e-11 # 地球の自転の角速度（1/s）
  28. >   # ターゲット日時を時間を含んだユリウス通日にして、GPS時間起点からの通日を計算
  29. >   jd <- julian(as.Date(paste(TY,"-",TM,"-",TD,sep="")),-2440588)
  30. >   jd <- jd - 0.5 + (Th*3600 + Tm*60 + Ts) / 86400 - jd0
  31. >   weeks <- floor(jd[1]/7)     # GPS週
  32. >   seconds <- (jd[1]%%7)*86400 # 週初からの経過時間（秒）
  33. >   elapse <- seconds - p12 # データの時刻からターゲット時刻までの秒数
  34. >   a <- p11^2 # 軌道の長半径
  35. >   e <- p9    # 離心率
  36. >   # 平均近点角を初期値に離心近点角 E を計算 
  37. >   E <- p7 + elapse * (sqrt(GM/(a^3)) + p6) # P7:M0,P6:ΔN（平均運動の補正値）
  38. >   d <- 1
  39. >   while(d>1e-16)
  40. +   {
  41. +     d <- (p7 - E + e*sin(E)) / (1-e*cos(E))
  42. +     E <- E + d
  43. +   }
  44. >   # 真近点角 Vk を計算
  45. >   sVk <- sqrt(1 - (e^2)) * sin(E)
  46. >   cVk <- cos(E) -e
  47. >   Vk <- atan(sVk/cVk)
  48. >   # 地心-近地点と赤道面とのなす角ωを加算。Pkは赤道面から衛星への角度
  49. >   Pk <- Vk + p18 # p18:ω(近地点引数)
  50. >   # 緯度引数。 p10:Cus,p8:Cuc(緯度引数の補正量)
  51. >   Uk <- Pk + p10*sin(2*Pk)+p8*cos(2*Pk)
  52. >   # 地心-衛星間の距離。p5:Crs,p17:Crc(動径方向の補正量)
  53. >   Rk <- a*(1-e*cos(E)) + p5*sin(2*Pk)+p17*cos(2*Pk)
  54. >   # p16:i0(軌道傾斜角),p20:Id(変化率),p15:Cis,p13:Cic(軌道傾斜角の補正量)
  55. >   Ik <- p16 + p20*elapse + p15*sin(2*Pk)+p13*cos(2*Pk)
  56. >   # 昇交点経度。 p14:Ω0(昇交点赤経),p19:Ωd(時間変化率),p12:Toe
  57. >   Ok <- p14+(p19-W)*elapse-W*p12
  58. >   # 軌道面上での衛星の座標
  59. >   x0 <- Rk * cos(Uk)
  60. >   y0 <- Rk * sin(Uk)
  61. >   # X軸で回転行列
  62. >   sn <- sin(Ik); cs <- cos(Ik)
  63. >   mI <- matrix(c(1, 0,  0,
  64. +                  0,cs,-sn,
  65. +                  0,sn, cs),ncol=3,byrow=T)
  66. >   # Z熟で回転行列
  67. >   sn <- sin(Ok); cs <- cos(Ok)
  68. >   mO <- matrix(c(cs,-sn,0,
  69. +                  sn, cs,0,
  70. +                   0,  0,1),ncol=3,byrow=T)
  71. >   # 回転
  72. >   x1 <- mO%*%(mI%*%matrix(c(x0,y0,0)))
  73. >   x1
  74.               [,1]
  75. [1,] -18111544.953
  76. [2,]   -473944.479
  77. [3,] -19253699.856
mikeo_410@hotmail.com