衛星の地表高度の計算アルゴリズム



● (No.271) 衛星の地表高度の計算アルゴリズム (2001年 6月24日)
 -----------------------------------------------------------

ケプラー方程式ニュートン近似法

(1)  M = MA*2*3.14/256                      MAをラジアンに変換
(2)  to = M+e*sin(M)+0.5*e^2*sin(2M)        初期値
(3)  mo = to-e*sin(to)                      ケプラー方程式
(4)  dto = (M-mo)/(1-e*cos(to))             e=離心率
(5)  t1 = to+dto                            第一近似値
(6)  m1 = t1-e*sin(t1)                      小計算
(7)  dt1 = (M-m1)/(1-e*cos(t1))             小計算
(8)  t2 = t1+dt1                            第二近似値
(9)  m2 = t2-e*sin(t2)                      小計算
(10) dt2 = (M-m2)/(1-e*cos(t2))             小計算
(11) t3 = t2+dt2                            第三近似値
(12) t4 = tan(t3/2)                         小計算
(13) u = root((1+e)/(1-e))*t4               小計算
(14) s = 2*arctan(u)                        真近点離角
(15) p = 24*3600/n                          周期(単位秒)
(16) 4*3.14^2*(a^3/p^2) = G*Q               万有引力の法則
(17) a = (((G*Q*p^2)/(4*3.14^2))^(1/3))     軌道長半径
(18) r = a*(1-e^2)/(1+e*cos(s))             地球の中心と衛星との距離
(19) rh = r-637814000                       地表面と衛星との距離


この式(17)を補足します。

   a = (((G*Q*p^2)/(4*3.14^2))^(1/3))
     = (6.6732*10^(-8)*5.9732*10^27*(p*p)/(4*3.14*3.14))^(1/3)


ここで、項目の中の G は「万有引力定数」で、Q は「地球の質量」です。
cgs系で考えると、単位の [N](ニュートン) は、

   [N] = [(kg*m)/s^2] = 10^3 * 10^2 [(g*cm)/s^2]
                      = 10^5        [(g*cm)/s^2] なので、

   [m^2] = 10^4 [cm^2],  [kg^2] = 10^6 [g^2] にも注意して、

   G = 6.6732 * 10^(-11)                       [N*m^2/kg^2]
     = 6.6732 * 10^(-11) * 10^5 * 10^4 / 10^6  [cm^3/(g*(s^2))]
     = 6.6732 * 10^(-8)                        [cm^3/(g*(s^2))]

   Q = 5.9732 * 10^24  [kg]
     = 5.9732 * 10^27  [g]


衛星の質量を R [g],  軌道長半径を a [cm], 周期を p [s] として、
等速円運動と仮定して万有引力の法則とともに運動方程式をたてると、

   R * a * (2π/p)^2 = (G * Q * R) / a^2

   p^2 = (4π^2 / (G * Q)) * a^3

つまり、両辺を (1/3)乗 して変形し整理すると、

   a = ((G * Q) * p^2 / 4π^2 )^(1/3)   となります。

この式に cgs系の数値を代入したものが、上記の式(17)です。

a の単位は、 [(cm^3/(g*(s^2)) * g * s^2)^(1/3)] = [cm]  です。


--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
           平均近点角         離心率     初期値   ケプラー方程式    第一近似値                 第二近似値                 第三近似値                  真近点離角    平均運動    周期(総秒) (時)  (分)   軌道長半径    地心距離      地表距離
 MA        MA'     M          e          to       Mo       dto      t1       m1       dt1      t2       m2       dt2      t3          t4           u        s        n             p        p1    p2    a             r             rh
 
(sample1 : 16 Jan 1994, AO-13)
 22        30.9    0.53969    0.72099    1.139    0.484    0.079    1.218    0.541   -0.002    1.216    0.540    0.000    1.216       0.696        1.729    2.093    2.09721276    41198    11    27    2578169623    1933173674    12953.6
 
(sample2 : 28 May 2001, AO-40)
0.01        0.0    0.00025    0.81492    0.001    0.000    0.000    0.001    0.000    0.001    0.002    0.000    0.001    0.003       0.002        0.006    0.012    1.27026844    68017    18    54    3601408456     666581751      287.7
128       180.0    3.14000    0.81492    3.140    3.139    0.001    3.141    3.141   -0.001    3.140    3.139    0.001    3.141    3374.653    10567.542    3.141    1.27026844    68017    18    54    3601408456    6536251656    58984.4
255.99    360.0    6.27975    0.81492    6.275    6.282   -0.010    6.265    6.280   -0.001    6.264    6.280   -0.001    6.263      -0.010       -0.031   -0.062    1.27026844    68017    18    54    3601408456     667135754      293.2
 
(sample3 : 23 Jun 2001, AO-40 before arcjet)
0.01        0.0    0.00025    0.81508    0.001    0.000    0.000    0.001    0.000    0.001    0.002    0.000    0.001    0.003       0.002        0.006    0.012    1.2711484     67970    18    53    3599749208     665697946      278.8
 86       121.4    2.10969    0.81508    2.517    2.040    0.042    2.559    2.111   -0.001    2.558    2.109    0.000    2.558       3.329       10.430    2.950    1.2711484     67970    18    53    3599749208    6046194040    54083.8
128       180.0    3.14000    0.81508    3.140    3.139    0.001    3.141    3.141   -0.001    3.140    3.139    0.001    3.141    3374.653    10572.585    3.141    1.2711484     67970    18    53    3599749208    6533816936    58960.0
 
(sample4 : 23 Jun 2001, AO-40 after arcjet)
0.01        0.0    0.00025    0.81501    0.001    0.000    0.000    0.001    0.000    0.001    0.002    0.000    0.001    0.003       0.002        0.006    0.012    1.2710861     67973    18    53    3599855129     665944691      281.3
100       140.5    2.45313    0.81501    2.645    2.257    0.114    2.759    2.455   -0.001    2.758    2.453    0.000    2.758       5.150       16.132    3.018    1.2710861     67973    18    53    3599855129    6321346299    56835.3
128       180.0    3.14000    0.81501    3.140    3.139    0.001    3.141    3.141   -0.001    3.140    3.139    0.001    3.141    3374.653    10570.598    3.141    1.2710861     67973    18    53    3599855129    6533782042    58959.7
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------


The solution of this algorithm corrected by W3IWI and JE9PEL :

                     Before       After
  Perigee Height :   278.8        281.3  ....... +2.5 km
  Apogee Height  : 58960.0      58959.7  ....... -0.3 km


The report of Peter Guelzow, DB2OS :

                     Before       After
  Perigee Height :   279.754      282.2427  .... +2.5 km
  Apogee Height  : 58971.166    58971.024   .... -0.1 km


Excel sheet for Satellite Height
http://www.ne.jp/asahi/hamradio/je9pel/satrangj.zip



[補足]

2001年6月23日に実施された衛星AO-40 のアークジェットの噴射前と噴射後の
近地点, 遠地点(Perigee/Apogee) の「地表距離(hr)」は、上記 Excel sheet
内の sample3, sample4 の計算結果を見ると上のようになり、下記の AO-40 
コマンドチームの DB2OS, Peter Guelzow氏による amsat-bbへの広報の中の
数値と比較すると ほんの少し差違があるものの、その変移はほぼ同じ結果と
なりました。


> Subject: AO-40: more successful ARCJET operations
> From: Peter Guelzow
> To: amsat-bb
> Date: Sat, 23 Jun 2001 10:08:15 +0000
> 
> Dear All,
> 
> The arc-jet thruster has been invoked on Orbit 296, MA 118-135,
> which gave about 1 hour of thrusting gas only. 
> 
> The S2 TX is OFF from MA 100-180 to spare power for the ATOS. 
> The gas generator for the ammonia draws about 120 - 130 W of power
> when cycled on by the thermostat. 
> 
> The IHU-2 is running and logging telemetry into a circular buffer
> and hold about 2.5 days worth of data. Downloaded telemetry from
> this 1h burn indicated positive power budget and everything looks
> nominally. 
> 
> The thrust on orbit 296 started at MA 121.4 and lasted for 3618s.
> The acceleration is guestimated to be 54E-6 m/s^2 ,
> and the direction of acceleration is towards alon 274, alat -2 
> (the current attitude).
> 
> Give or take the unknowns, the expected outcome of this was:
> 
>                       Before             After
>   --------------------------------------------
>   Epoch year            2001              2001
>   Epoch time       173.12145         173.16312
>   Inclination         5.2833         5.2833592
>   R.A.A.N          180.71591         180.70361
>   Eccentricity      0.815077         0.8150139
>   Arg perigee      288.69088         288.71333
>   Mean Anomaly         121.4         140.46836
>   Mean motion MM   1.2711484         1.2710861
>   Revolution             296               296
>   SMA                36003.6         36004.773
>   Perigee Height     279.754          282.2427  .... +2.5 km
>   Apogee Height    58971.166         58971.024
>   --------------------------------------------
> 
> So the perigee was raised by this 1h burn by about +2.5 km, but only
> a thrust rate of 50% was used. 
> 
> Since everything went so well, the computer onboard AO-40 was 
> commanded to initiate 2h "burns" starting around apogee on orbit 297
> for the next three orbits. 
> 
> The first 2-hour burn stopped at Orbit 297 MA 142, which is 2001 
> Jun 23 0154 utc. Again, all telemetry is looking good and we may soon
> expand to 4h "burns" and possibly increase the thrust level as well.
> 
> While we haven't yet seen the latest NORAD data, some effects of the
> changing orbit should be soon noticeable. 
> 
> Indeed, we all are very happy with the successful results of the ATOS
> (Arcjet Thruster on OSCAR Satellite) so far.
> 
> Some more information about the ATOS system is available here:
> 
> http://www.irs.uni-stuttgart.de/RESEARCH/EL_PROP/PROJ/e_atos.html
> 
> best 73s
> Peter DB2OS for the AO-40 team


> Subject: Re: AO-40: more successful ARCJET operations
> From: Tom Clark
> To: Mineo Wakita, amsat-bb
> Date: Sun, 24 Jun 2001 18:08:22 -0400
> 
> Mineo -- Nice work. The perigee difference is probably due to your
> definition of Height, i.e.  Height from the surface of the earth.
> You said:
> 
> > (18) r = a(1-e^2)/(1+e*cos(s)) : Height from the center of the earth
> > (19) rh = r-6371               : Height from the face of the earth
> 
> The equatorial radius of the oblate earth is 6378.160 km, larger than
> the 6371 km you use. Since the inclination is only 5.3 deg, it is
> appropriate to use the equatorial value.
> 
> Therefore, using your formulation, I come up with:
> 
> > The solution of this algorithm [posted by JE9PEL]
> >                      Before       After
> >   Perigee Height :   286.0        288.4  ....... +2.4 km
> >
> > The report of Peter Guelzow, DB2OS
> >                      Before       After
> >   Perigee Height :   279.754      282.2427  .... +2.5 km
> 
>   The JE9PEL perigee corrected by W3IWI:
>                        Before       After
>     Perigee Height :   278.8        281.3  ....... +2.5 km
> 
> This agrees to ~1 km with the results that Peter reported. Since the
> object of the "burn" was to raise perigee (with very little change
> in apogee), and since a post-burn mean motion is not yet available,
> I haven't bothered to examine the apogee differences.
> 
> 73 de Tom, W3IWI


> Subject: Re: AO-40: more successful ARCJET operations
> From: Mineo Wakita
> To: jamsat-bb
> Date: Mon, 25 Jun 2001 22:26:26 +0900
> 
> 昨日「衛星高度の計算アルゴリズム」についての原稿の英語版を amsat-bb
> にも投稿したのですが、それに対してその道の専門家の一人の Tom Clark,
> W3IWI氏 から、直々に返信を受け取りました。
> 
> このアルゴリズムの最後の式の、「(19) r=r0-6371 」で使用している数値
> が間違っている、というものでした。  さっそく天文年鑑で確認すると、
> 確かに地球の赤道半径は、「6378.140 km」とありました。
> 
> 現在の衛星AO-40 の軌道傾斜角は わずか 5.3度なので、この数値で代用で
> きる、ということなので、この数値に差し替えて計算しなおしてみると、
> 近地点距離は、昨日の Peter Guelzow / DB2OS氏 の広報の Perigee と、
> 完璧に 1km 以内の誤差で一致しました。
> 
> 昨日の原稿の一部と 計算のための Excelシートを、下記のように差し替え
> ます。 (訂正済み)
> 
> 
> (18) r = a(1-e^2)/(1+e*cos(s)) : 地球の中心と衛星との距離
> (19) rh = r-637814000          : 地表面と衛星との距離
> 
> The solution of this algorithm corrected by W3IWI and JE9PEL :
> 
>                      Before       After
>   Perigee Height :   278.8        281.3  ....... +2.5 km
>   Apogee Height  : 58960.0      58959.7  ....... -0.3 km
> 
> The report of Peter Guelzow, DB2OS :
> 
>                      Before       After
>   Perigee Height :   279.754      282.2427  .... +2.5 km
>   Apogee Height  : 58971.166    58971.024   .... -0.1 km
> 
> Excel sheet for Satellite Height
> http://www.ne.jp/asahi/hamradio/je9pel/satrangj.zip



《参考》

なお、和製の衛星軌道計算ソフトとして有名な 『Calsat32』 の作者、JR1HUO
相田政則氏のホームページの中で、「CALSAT32のからくり」と題して、衛星の
位置計算に必要な事柄について、項目ごとに大変詳しく、なおかつ、わかりや
すく説明をしています。例題のプログラムも付いていて、必見の内容です。
http://jr1huo.my.coocan.jp/jr1huo_calsat32/Calsat32Karakuri.htm

暦歴 (ケプラー方程式)
http://www.sanmao.cocolog-nifty.com/reki/cat24015528/


 トップ へ戻る.
 前のページ へ戻る.
 次のページ へ移る.
 ホームページ(目次) へ戻る.