N88-BASICで小惑星の軌道 (1回目 Apophis)

本日、2021.3.6(土)

小惑星アポフィスが地球に接近するそうです

と言っても月までの数十倍の距離だそうですが

詳しくはアポフィスで検索するといろいろ

出てきます

2068年に衝突する可能性が0ではないと

言われているそうです

1年前からの1年間の軌道を描画して見ました

青い軌跡が惑星、紫の軌跡が小惑星

水色の点が地球、赤色の点がアポフィスです

プログラムは

https://ameblo.jp/vlbasic/entry-12627629545.html

N88-BASICで惑星の軌道 (9回目)

を小惑星用に削って使用しています

式など詳しくは

https://ameblo.jp/vlbasic/entry-12623474486.html

N88-BASICで惑星の軌道 (1回目)

からの記事を見てください

アポフィスの軌道データは

JPL Small-Body Database Browser (nasa.gov)

を参照しました

データの扱いなどミスがあるかもしれませんので

自己責任で鑑賞して下さい

NL-BASICとnl30306.zip(ast001.bas)は

以下のリンクからダウンロードできます

NL-BASIC(N88-BASIC互換?)ホームページ

Readme.txtを読んで遊んで下さい

(動画もあります)

下記リストをマウスで選択しCtrl+cでコピーし、

NL-BASICの画面でAlt+v(Ctrl+vではないので注意)

でプログラムを読込めます。

ast001.bas

1000 '----------------------------------------------------------------------
1010 ' N88-BASICで小惑星の軌道 (1回目 Apophis) by ULproject 2021.3
1020 '
1030 ' 画面下が平均春分点の方向(春分の日の地球から太陽への方向)
1040 ' 画面上で地球が太陽の真下のとき秋分の日近くになります。
1050 ' 黄道(地球の軌道)面の真上から見た軌道を表示しています。
1060 ' 遠近法を使用していないので遠くも同じ大きさになります。
1070 '
1080 ' テキスト表示のあるプログラムはCtrl+sで一時停止できます。
1090 ' NL-BASICではAlt+6で多少速くできます。
1100 '----------------------------------------------------------------------
1110 DEFDBL A-H, J, L-S, U-Z: DEFINT I,K: DEFSNG T
1120 '--- 関数,定数定義
1130 DEF FNT(A) = A * PI / 180
1140 PI = ATN(1) * 4
1150 PI2 = PI * 2
1160 '--- 画面消去
1170 SCREEN 3: CONSOLE ,,, 1: COLOR 7,,,, 2: CLS 3
1180 '--- Select
1190 INPUT "見る角度(90など) "; X1
1200 X1 = FNT(X1)
1210 CLS
1220 '--- 基本データ読込み
1230 READ R, DN
1240 DIM JD0(DN-1)
1250 FOR I=0 TO DN-1
1260 READ YY, MM, DD: GOSUB *DATE2JD: JD0(I) = JD
1270 NEXT
1280 READ TL, TS
1290 TL = TL - 1
1300 READ KN, K3
1310 DIM C(KN), D(KN), A(KN), E(KN), I(KN), OO(KN), O(KN), M0(KN), P(KN)
1320 DIM X(KN), Y(KN), B(KN), T(15), M(15), MM(15, KN)
1330 '--- 軌道データ読込み
1340 GM = FNT(270) '--- 平均春分点の方向(画面下)
1350 FOR K=1 TO KN
1360 READ C(K), D(K), A(K), E(K), I(K), OO(K), O(K), M0(K), P(K)
1370 A(K) = A(K) * R
1380 I(K) = FNT(I(K) )
1390 OO(K) = FNT(OO(K)) + GM
1400 O(K) = FNT(O(K) ) + GM
1410 M0(K) = FNT(M0(K))
1420 X(K) = 640 '--- 画面外
1430 Y(K) = 400
1440 NEXT
1450 '---
1460 P = P(K3)^2 / A(K3)^3 '--- T^2/a^3 = const.
1470 FOR K=1 TO KN
1480 GOSUB *I2MATRIX '--- i
1490 B(K) = 1
1500 IF P(K) = 0 THEN P(K) = SQR(A(K)^3 * P): B(K) = 8+3
1510 P(K) = P(K) * 365.25 '--- ユリウス年を日に変換
1520 T = JD0(0) - JD0(D(K))
1530 N = PI2 / P(K) '--- n:平均運動(ラジアン/日)
1540 M0(K) = M0(K) + N * T '--- M:平均近点角
1550 NEXT
1560 '--- 画面表示
1570 X = 320
1580 Y = 200
1590 WINDOW(-X, -Y)-(X-1, Y-1)
1600 LINE(-X, 0)-(X-1, 0), 1
1610 LINE(0, -Y)-(0, Y-1), 1
1620 LOCATE 36, 0: PRINT "春分(平均)";
1630 LOCATE 0, 12: PRINT "夏至(平均)";
1640 LOCATE 36, 23: PRINT "秋分(平均)";
1650 LOCATE 70, 12: PRINT "冬至(平均)";
1660 '--- 方向表示
1670 COLOR 4: LOCATE 0, 0: PRINT "春分点方向(平均)"
1680 COLOR 2: LOCATE 0, 1: PRINT "近日点方向(赤:惑星)"
1690 COLOR 1: LOCATE 0, 2: PRINT "近日点方向(青:小惑星等)"
1700 COLOR 6: LOCATE 0, 3: PRINT "昇交点方向(各色)"
1710 COLOR 5: LOCATE 0, 4: PRINT "地球(水色)"
1720 COLOR 7
1730 R = 190
1740 LINE(0, 0)-(R * COS(GM ), -R * SIN(GM )), 4
1750 FOR K=KN TO 1 STEP -1
1760 B = 2
1770 IF B(K) > 7 THEN B = 1
1780 R = A(K)
1790 X = R * COS(O(K))
1800 Y = R * SIN(O(K))
1810 Z = 0
1820 GOSUB *PRODUCT
1830 LINE(0, 0)-(X1, -Y1), C(K)
1840 R = A(K) * (1 - E(K))
1850 X = R * COS(OO(K))
1860 Y = R * SIN(OO(K))
1870 Z = 0
1880 GOSUB *PRODUCT
1890 LINE(0, 0)-(X1, -Y1), B
1900 NEXT
1910 CIRCLE(0, 0), 3, 2,,,,F
1920 '--- Planets (T=0 → 元期0)
1930 FOR T = 0 TO TL STEP TS '--- ts日ずつ
1940 JD = JD0(0) + INT(T): GOSUB *JD2DATE
1950 WW = INT(JD + 0.5)
1960 WW = WW - INT(WW / 7) * 7
1970 LOCATE 64, 10
1980 PRINT USING "####/##/##("; YY, MM, DD;
1990 PRINT MID$("月火水木金土日", WW*2+1, 2); ")"
2000 FOR K=1 TO KN
2010 C = C(K) : A = A(K): E = E(K): OO = OO(K): O = O(K)
2020 M0 = M0(K): P = P(K): X = X(K): Y = Y(K) : B = B(K)
2030 GOSUB *ORBIT
2040 X(K) = X: Y(K) = Y
2050 NEXT
2060 NEXT
2070 LOCATE 0, 13
2080 END
2090 '--- Orbit
2100 *ORBIT
2110 N = PI2 / P '--- n:平均運動(ラジアン/日)
2120 M = M0 + N * T '--- M:平均近点角
2130 M = M - FIX(M / PI2) * PI2
2140 GOSUB *NEWTON '--- u:離心近点角
2150 F = ATN(SQR((1+E)/(1-E)) * TAN(U / 2)) * 2 '--- f:真近点角
2160 R = A * (1 - E*E) / (1 + E * COS(F)) '--- r:動径
2170 CIRCLE(X, -Y), 1, B,,,, F
2180 X = R * COS(F + OO)
2190 Y = R * SIN(F + OO)
2200 Z = 0
2210 GOSUB *PRODUCT
2220 X = X1: Y = Y1
2230 CIRCLE(X, -Y), 1, C,,,, F
2240 RETURN
2250 '--- Newton method newton(M, e) = u
2260 *NEWTON
2270 U = M + E * SIN(M)
2280 FOR I=20 TO 1 STEP -1 '--- 最大繰返し回数
2290 F = U - E * SIN(U) - M '--- f(u)
2300 FD = 1 - E * COS(U) '--- f'(u)
2310 DU = F / FD
2320 U = U - DU
2330 IF ABS(DU) <= 1E-5 THEN I = 1 '--- 終了条件(誤差)
2340 NEXT
2350 RETURN
2360 '--- Date(YY/MM/DD) → JD
2370 *DATE2JD
2380 IF MM <= 2 THEN JM = MM + 13: JY = YY - 1 ELSE JM = MM + 1: JY = YY
2390 JD = INT(JY * 365.25) + INT(JM * 30.601) + DD + 1721117.5 - 122 - 1
2400 IF JD >= 2299160.5 THEN JD = JD -INT(JY/100) +INT(JY/400) + 12 - 10
2410 RETURN
2420 '--- JD → Date(YY/MM/DD)
2430 *JD2DATE
2440 DD = JD + 0.5
2450 JY = INT(DD)
2460 DD = DD - JY
2470 IF JY < 2299161 THEN *JD2DATE1
2480 JY = JY - (12 - 10)
2490 YY = (JY - 1721118 + 1) / 365.2425
2500 JY = JY + INT(YY / 100) - INT(YY / 400)
2510 *JD2DATE1
2520 JY = JY - 31 - 28
2530 YY = INT((JY - 0.001) / 365.25)
2540 JY = JY - INT(YY * 365.25) + 122
2550 MM = INT(JY / 30.601)
2560 DD = DD + JY - INT(MM * 30.601)
2570 IF MM > 13 THEN MM = MM - 13: YY = YY + 1 ELSE MM = MM - 1
2580 YY = YY - 4712
2590 RETURN
2600 '--- 軌道傾斜角i(k),o(k) → MM(0-15, k)
2610 *I2MATRIX
2620 GOSUB *IDENTITY
2630 VIEWMODE = 0
2640 D = I(K ): O = O(K ): X = COS(O): Y = SIN(O): Z = 0
2650 GOSUB *ROTATE '--- iを傾ける
2660 VIEWMODE = 1
2670 D = -I(K3): O = O(K3): X = COS(O): Y = SIN(O): Z = 0
2680 GOSUB *ROTATE '--- 見る場所を黄道面真上にする
2690 D = -X1: X = 1: Y = 0: Z = 0
2700 GOSUB *ROTATE '--- x軸回転
2710 VIEWMODE = 0
2720 FOR I=0 TO 15
2730 MM(I, K) = M(I)
2740 NEXT
2750 RETURN
2760 '----------------------------------------------------------------------
2770 ' 行列とベクトルの積(p1 = Mp) ULproject
2780 '
2790 ' inp MM(0-15, k), x, y, z
2800 ' out x1, y1, z1
2810 '----------------------------------------------------------------------
2820 *PRODUCT
2830 X1 = MM(0, K) * X + MM(4, K) * Y + MM( 8, K) * Z + MM(12, K)
2840 Y1 = MM(1, K) * X + MM(5, K) * Y + MM( 9, K) * Z + MM(13, K)
2850 Z1 = MM(2, K) * X + MM(6, K) * Y + MM(10, K) * Z + MM(14, K)
2860 RETURN
2870 '----------------------------------------------------------------------
2880 ' 行列の積(M = MT model mode , M = TM view mode) ULproject
2890 '
2900 ' inp T,M
2910 ' out M
2920 ' tmp m1,m2,m3,m4
2930 '----------------------------------------------------------------------
2940 *MULT
2950 IF VIEWMODE THEN *MULT.TM
2960 FOR IY=0 TO 3
2970 M1 = M(IY): M2 = M(IY+4): M3 = M(IY+8): M4 = M(IY+12)
2980 FOR IX=0 TO 12 STEP 4
2990 M(IX+IY) = M1*T(IX) + M2*T(IX+1) + M3*T(IX+2) + M4*T(IX+3)
3000 NEXT
3010 NEXT
3020 RETURN
3030 *MULT.TM
3040 FOR IX=0 TO 12 STEP 4
3050 M1 = M(IX): M2 = M(IX+1): M3 = M(IX+2): M4 = M(IX+3)
3060 FOR IY=0 TO 3
3070 M(IX+IY) = T(IY)*M1 + T(IY+4)*M2 + T(IY+8)*M3 + T(IY+12)*M4
3080 NEXT
3090 NEXT
3100 RETURN
3110 '----------------------------------------------------------------------
3120 ' 単位行列 ULproject
3130 '
3140 ' out M
3150 '----------------------------------------------------------------------
3160 *IDENTITY
3170 M(0) = 1: M(4) = 0: M( 8) = 0: M(12) = 0
3180 M(1) = 0: M(5) = 1: M( 9) = 0: M(13) = 0
3190 M(2) = 0: M(6) = 0: M(10) = 1: M(14) = 0
3200 M(3) = 0: M(7) = 0: M(11) = 0: M(15) = 1
3210 RETURN
3220 '----------------------------------------------------------------------
3230 ' 並進行列 ULproject
3240 '
3250 ' inp x, y, z
3260 ' out T,M
3270 '----------------------------------------------------------------------
3280 *TRANSLATE
3290 T(0) = 1: T(4) = 0: T( 8) = 0: T(12) = X
3300 T(1) = 0: T(5) = 1: T( 9) = 0: T(13) = Y
3310 T(2) = 0: T(6) = 0: T(10) = 1: T(14) = Z
3320 T(3) = 0: T(7) = 0: T(11) = 0: T(15) = 1
3330 GOSUB *MULT
3340 RETURN
3350 '----------------------------------------------------------------------
3360 ' 回転行列 ULproject
3370 '
3380 ' inp d(rad), x, y, z
3390 ' out T,M
3400 ' tmp rr,tt,ss,nn,n1,n2,n3
3410 '----------------------------------------------------------------------
3420 *ROTATE
3430 RR = SQR(X*X + Y*Y + Z*Z)
3440 N1 = X / RR
3450 N2 = Y / RR
3460 N3 = Z / RR
3470 SS = COS(D)
3480 RR = 1.0 - SS
3490 T( 0) = N1 * N1 * RR + SS
3500 T( 5) = N2 * N2 * RR + SS
3510 T(10) = N3 * N3 * RR + SS
3520 SS = SIN(D)
3530 TT = N1 * N2 * RR: NN = N3 * SS: T(4) = TT - NN: T(1) = TT + NN
3540 TT = N1 * N3 * RR: NN = N2 * SS: T(2) = TT - NN: T(8) = TT + NN
3550 TT = N2 * N3 * RR: NN = N1 * SS: T(9) = TT - NN: T(6) = TT + NN
3560 T(12) = 0: T(13) = 0: T(14) = 0: T(15) = 1
3570 T( 3) = 0: T( 7) = 0: T(11) = 0
3580 GOSUB *MULT
3590 RETURN
3600 '----------------------------------------------------------------------
3610 '- 1AU = 1.495978707(億km)
3620 '- c:色(0黒,1青,2赤,3紫,4緑,5水,6黄,7白,+8+暗,8灰+暗)
3630 '- 軌道長半径a(AU), 離心率e
3640 '- 近日点黄経(Ω+ω近日点引数)(ﾟ), 昇交点黄経Ω(ﾟ)(平均春分点基準)
3650 '- 元期平均近点離角Mo(ﾟ)
3660 '- 周期P(ユリウス年)(=365.25日)
3670 '----------------------------------------------------------------------
3680 '- ～火星
3690 '----------------------------------------------------------------------
3700 DATA 100 '--- dot/AU(1AU = ?dot)
3710 DATA 4 '--- 元期(軌道要素を決めた日)数
3720 DATA 2020, 3, 6.0 '--- 元期d=0(表示開始日)
3730 DATA 2020,12,17.0 '--- 元期d=1(Mo惑星)
3740 DATA 2019, 4,27.0 '--- 元期d=2(Mo小惑星)
3750 DATA 2020,12,17.0 '--- 元期d=3(Mo小惑星)
3760 DATA 366, 1 '--- 表示期間(日),経過時間(日)
3770 DATA 7, 3 '--- 惑星数, 地球番号
3780 '--- c,d, a , e , i , Ω+ω , Ω, , Mo , P
3790 DATA 7,1, 0.3871, 0.2056, 7.004, 77.489, 48.305,183.661, 0.24085 '水
3800 DATA 6,1, 0.7233, 0.0068, 3.394,131.565, 76.622, 75.525, 0.61520 '金
3810 DATA 5,1, 1.0000, 0.0167, 0.003,103.006,174.823,342.791, 1.00002 '地
3820 DATA 3,1, 1.5237, 0.0934, 1.848,336.153, 49.496, 71.013, 1.88085 '火
3830 '--- c,d, a , e , i , Ω+ω , Ω , Mo , P
3840 DATA 14 ,2, 1.190, 0.190, 5.9,463.1 ,251.6,249.8 , 0 '竜宮
3850 DATA 12 ,2, 1.324, 0.280, 1.6,231.9 , 69.1,288.9 , 0 '糸川
3860 DATA 10 ,3, 0.923, 0.192, 3.3,330.6 ,204.0,110.6 , 0 'Apophis