比較参照データ：情報処理学会研究報告：Vol.2018-HPC-164 No.1 ：平山 弘
区間許容誤差 ： Tol_Err = 10^(-19) で計算 小数点以下10桁 11桁目四捨五入 ◎赤字部は誤り
| step数 |
x1 |
y1 |
x2 |
y2 |
x3 |
y3 |
|
| 0 |
1.0000000000 | 3.0000000000 | -2.0000000000 | -1.0000000000 | 1.0000000000 | -1.0000000000 | |
| 10 | 390997 |
0.7784804109 |
0.1413922996 |
-2.0250924779 |
0.0972193841 |
1.1529857358 |
-0.1626108871 |
| 20 |
767693 |
3.0042926367 |
0.5119252350 |
-1.3886265371 |
-0.4704760499 |
-0.6916743523 |
0.0692256989 |
| 30 |
1185421 |
0.8563404988 |
2.2870936634 |
-0.8779838776 |
-0.8659638274 |
0.1885828028 |
-0.6794851361 |
| 40 |
1401411 |
-0.6220036904 |
1.8583181561 | 0.1735445568 |
-2.3684104427 |
0.2343665688 |
0.7797374606 |
| 50 |
1824156 |
-2.7014614093 |
-3.7972226828 |
1.5059381921 |
0.9608134251 |
0.4161262919 |
1.5096828696 |
| 60 |
2311751 |
0.7438074974 |
1.9399479437 |
0.2640103328 |
-0.7316243915 |
-0.6574927646 |
-0.5786692530 |
区間許容誤差 ： Tol_Err = 10^(-20) で計算 小数点以下10桁 11桁目四捨五入 ◎赤字部は誤り
| step数 |
x1 |
y1 |
x2 |
y2 |
x3 |
y3 |
|
| 0 |
1.0000000000 | 3.0000000000 | -2.0000000000 | -1.0000000000 | 1.0000000000 | -1.0000000000 | |
| 10 | 625754 |
0.7784804104 | 0.1413923001 | -2.0250924780 | 0.0972193841 |
1.1529857361 |
-0.1626108874 |
| 20 |
1220309 |
3.0042926367 |
0.5119252350 |
-1.3886265374 |
-0.4704760502 |
-0.6916743521 |
0.0692256991 |
| 30 |
1890422 |
0.8563404989 |
2.2870936636 |
-0.8779838787 |
-0.8659638277 |
0.1885828037 |
-0.6794851361 |
| 40 |
2234593 |
-0.6220036914 |
1.8583181573 |
0.1735445568 |
-2.3684104431 |
0.2343665694 |
0.7797374601 |
| 50 |
2912127 |
-2.7014614093 |
-3.7972226827 |
1.5059381922 |
0.9608134251 |
0.4161262918 |
1.5096828696 |
| 60 |
3682562 |
0.7438075001 |
1.9399479510 |
0.2640103344 |
-0.7316243948 | -0.6574927676 | -0.5786692548 |
区間許容誤差 ： Tol_Err = 10^(-21) で計算 小数点以下10桁 11桁目四捨五入 ◎赤字部は誤り (★工事中★)
| step数 |
x1 |
y1 |
x2 |
y2 |
x3 |
y3 |
|
| 0 |
1.0000000000 | 3.0000000000 | -2.0000000000 | -1.0000000000 | 1.0000000000 | -1.0000000000 | |
| 10 | 991999 |
0.7784804102 | 0.1413923003 | -2.0250924780 | 0.0972193841 |
1.1529857363 |
-0.1626108875 |
| 20 |
1930950 |
3.0042926367 |
0.5119252350 |
-1.3886265375 | -0.4704760502 |
-0.6916743520 |
0.0692256992 |
| 30 |
2996061 |
0.8563404989 |
2.2870936637 | -0.8779838789 | -0.8659638277 |
0.1885828038 |
-0.6794851361 |
| 40 |
3533400 |
-0.6220036917 |
1.8583181577 | 0.1735445568 |
-2.3684104432 | 0.2343665696 |
0.7797374599 |
| 50 |
4610183 |
-2.7014614093 |
-3.7972226827 |
1.5059381924 | 0.9608134250 | 0.4161262917 |
1.5096828697 |
| 60 |
5835471 |
0.7438075000 | 1.9399479508 | 0.2640103346 | -0.7316243947 | -0.6574927677 | -0.5786692547 |
区間許容誤差 ： Tol_Err = 10^(-21) で計算 小数点以下15桁 16桁目四捨五入 ◎60秒目のみ
| 60 |
x1 | y1 |
| 0.74380750002421 | 1.939947950844325 | |
| x2 | y2 | |
| 0.264010334550272 | -0.731624394743371 | |
| x3 | y3 | |
| -0.657492767654744 | ?0.578669254711899 |
OPTION BASE 1 ! 十進ＢＡＳＩＣ
PRINT "開始時刻： "; DATE$;" ";TIME$
CALL メイン_SUB
PRINT "終了時刻： "; DATE$;" ";TIME$
SUB メイン_SUB ! メインルーチン
!---------------------------------------------------------------------------
LET 終了時刻 = 75 !終了時刻 ０～７５秒
LET Tol_ERR = 10^(-12) !許容誤差 -19 / -20 / -21
LET Tol_ERR_KTA_cnt=Tol_ERR
FOR KKK=1 TO 50
LET Tol_ERR_KTA_cnt=Tol_ERR_KTA_cnt*10
LET Tol_ERR_桁 = KKK
IF Tol_ERR_KTA_cnt>=.999 THEN EXIT FOR
NEXT KKK
PRINT "Tol_ERR = -"; Tol_ERR_桁
LET dt = 0.0000001 !時間幅 dt 最適値=10^(-7)
PRINT "初期dt = "; dt
LET t10_60=1 !10,20,30,40,50,60,
!---------------------------------------------------------------------------
DECLARE EXTERNAL SUB テキスト印刷 ! 値引き渡し ByVal
DECLARE EXTERNAL SUB RK4 ! 値引き渡し ByVal
DECLARE EXTERNAL SUB RK4_A ! 値引き渡し ByVal
DECLARE EXTERNAL SUB RK4_B ! 値引き渡し ByVal
DECLARE EXTERNAL FUNCTION 適正dt ! 値引き渡し ByVal
DECLARE EXTERNAL SUB 図形出力 ! 値引き渡し ByVal
DECLARE EXTERNAL SUB 一次近似_出力
DIM rr(3, 3) ! rr（識別番号,軸）
DIM vv(3, 3) ! vv（識別番号,軸）
DIM m(3)
DIM pre_rr(3, 3) ! 一時保管 rr（識別番号,軸）
DIM RK4_ans(18)
DIM COM_ans(18)
!初期条件
LET m(1) = 3 ! 質量１
LET m(2) = 4 ! 質量２
LET m(3) = 5 ! 質量３
LET rr(1, 1) = 1 !（識別番号,軸）
LET rr(1, 2) = 3 !（識別番号,軸）
LET rr(1, 3) = 0 !（識別番号,軸）
LET rr(2, 1) = -2 !（識別番号,軸）
LET rr(2, 2) = -1 !（識別番号,軸）
LET rr(2, 3) = 0 !（識別番号,軸）
LET rr(3, 1) = 1 !（識別番号,軸）
LET rr(3, 2) = -1 !（識別番号,軸）
LET rr(3, 3) = 0 !（識別番号,軸）
LET vv(1, 1) = 0 !（識別番号,軸）
LET vv(1, 2) = 0 !（識別番号,軸）
LET vv(1, 3) = 0 !（識別番号,軸）
LET vv(2, 1) = 0 !（識別番号,軸）
LET vv(2, 2) = 0 !（識別番号,軸）
LET vv(2, 3) = 0 !（識別番号,軸）
LET vv(3, 1) = 0 !（識別番号,軸）
LET vv(3, 2) = 0 !（識別番号,軸）
LET vv(3, 3) = 0 !（識別番号,軸）
!---------------------------------------------------------------------------
LET cnt0 = 0 !初期化
LET cnt1 = 0 !初期化
! PRINT "dt" , "t" , "x1" , "y1" , "x2" , "y2" , "x3" , "y3" !先頭表記
!---------------------------------------------------------------------------
DO WHILE t <= 終了時刻
! 出力 dt t M1ｘ M1y M2ｘ M2y M3ｘ M3y
! PRINT dt,t,rr(1, 1),rr(1, 2),rr(2, 1),rr(2, 2),rr(3, 1),rr(3, 2) !（識別番号,軸）
! IF ROUND(t,2)=INT(t) OR ROUND(t,2)=CEIL(t) THEN ! 小数点下２桁が整数に近いとき 例) .99～1.01
IF t < t10_60 THEN ! 一時保管
LET pre_dt=dt
LET pre_t=t
MAT pre_rr=rr
END IF
IF t > t10_60 THEN
! CALL テキスト印刷((pre_dt),(pre_t),pre_rr) ! 値渡し () で括る
! CALL テキスト印刷((dt),(t),rr) ! 値渡し () で括る
CALL 一次近似_出力 ((cnt0),(pre_t),pre_rr,(t),rr)
LET t10_60=t10_60 + 1
END IF
! END IF
CALL 図形(rr)
LET cnt1 = cnt1 + 1
LET dt = dt * 2 !2倍にする
LET dt=適正dt(t, dt, m, rr, vv,RK4_ans,Tol_ERR)
! CALL RK4(t, dt, m, rr, vv,RK4_ans) ! ルンゲ＝クッタ RK4_ans()は行列
! CALL RK4_A(t, dt, m, rr, vv,RK4_ans) ! ルンゲ＝クッタ RK4_ans()は行列
CALL RK4_B(t, dt, m, rr, vv,RK4_ans) ! ルンゲ＝クッタ RK4_ans()は行列
LET rr(1, 1) = RK4_ans(1)
LET rr(1, 2) = RK4_ans(2)
LET rr(1, 3) = RK4_ans(3)
LET vv(1, 1) = RK4_ans(4)
LET vv(1, 2) = RK4_ans(5)
LET vv(1, 3) = RK4_ans(6)
LET rr(2, 1) = RK4_ans(7)
LET rr(2, 2) = RK4_ans(8)
LET rr(2, 3) = RK4_ans(9)
LET vv(2, 1) = RK4_ans(10)
LET vv(2, 2) = RK4_ans(11)
LET vv(2, 3) = RK4_ans(12)
LET rr(3, 1) = RK4_ans(13)
LET rr(3, 2) = RK4_ans(14)
LET rr(3, 3) = RK4_ans(15)
LET vv(3, 1) = RK4_ans(16)
LET vv(3, 2) = RK4_ans(17)
LET vv(3, 3) = RK4_ans(18)
LET cnt0 = cnt0 + 1
LET t = t + dt
LOOP
END SUB
END
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
EXTERNAL SUB 一次近似_出力(cnt0,pre_t,pre_rr(,),t,rr(,))
DIM 近似_rr(3,2)
LET NNN_t = INT(t) !整数部 1秒、2秒・・・、60秒
LET 差1 = NNN_t - pre_t
LET 差2 = t - NNN_t
LET 比1=差1/(差1+差2)
LET 比2=差2/(差1+差2)
LET 近似_rr(1,1) = 比1*rr(1,1) + 比2*pre_rr(1,1)
LET 近似_rr(1,2) = 比1*rr(1,2) + 比2*pre_rr(1,2)
LET 近似_rr(2,1) = 比1*rr(2,1) + 比2*pre_rr(2,1)
LET 近似_rr(2,2) = 比1*rr(2,2) + 比2*pre_rr(2,2)
LET 近似_rr(3,1) = 比1*rr(3,1) + 比2*pre_rr(3,1)
LET 近似_rr(3,2) = 比1*rr(3,2) + 比2*pre_rr(3,2)
LET 近似_rr(1,1) = ROUND(近似_rr(1,1),10) ! 10桁丸め
LET 近似_rr(1,2) = ROUND(近似_rr(1,2),10) ! 10桁丸め
LET 近似_rr(2,1) = ROUND(近似_rr(2,1),10) ! 10桁丸め
LET 近似_rr(2,2) = ROUND(近似_rr(2,2),10) ! 10桁丸め
LET 近似_rr(3,1) = ROUND(近似_rr(3,1),10) ! 10桁丸め
LET 近似_rr(3,2) = ROUND(近似_rr(3,2),10) ! 10桁丸め
PRINT cnt0;" ";NNN_t;" ";近似_rr(1, 1);" ";近似_rr(1, 2);" ";近似_rr(2, 1);" ";近似_rr(2, 2);" ";近似_rr(3, 1);" ";近似_rr(3, 2) !（識別番号,軸）
END SUB
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
EXTERNAL FUNCTION 適正dt(t,dt,m(),r0(,),v0(,),EX2_RK4_ans(),Tol_ERR)
DECLARE EXTERNAL SUB RK4_A ! 値引き渡し ByVal
DECLARE EXTERNAL SUB RK4_B ! 値引き渡し ByVal
DIM EX2A_RK4_ans(18)
DIM EX2B_RK4_ans(18)
LET dt00=dt
DO
MAT EX2A_RK4_ans = EX2_RK4_ans ! 行列の代入
MAT EX2B_RK4_ans = EX2_RK4_ans ! 行列の代入
CALL RK4_A(t,dt,m,r0,v0,EX2A_RK4_ans) ! ４次の計算
CALL RK4_B(t,dt,m,r0,v0,EX2B_RK4_ans) ! ５次の計算
IF 誤差_和(EX2A_RK4_ans,EX2B_RK4_ans) < Tol_ERR THEN
EXIT DO
ELSE
LET dt=dt/2
END IF
LOOP
LET 適正dt=dt
END FUNCTION
!---------------------------------------------------------------------------
EXTERNAL FUNCTION 誤差_和(ans_A(),ans_B())
FOR KK=1 TO 18
IF ans_A(KK)<>0 THEN ! ０の割り算回避
LET 差 = ans_A(KK) - ans_B(KK)
LET sum_aaa = sum_aaa + ABS(差/ans_A(KK))
END IF
NEXT KK
LET 誤差_和=sum_aaa
END FUNCTION
!---------------------------------------------------------------------------
EXTERNAL SUB RK4_B(t,dt,m(),r0(,),v0(,),EX_RK4_ans()) ! dt 1/2を２回
! OPTION ARITHMETIC DECIMAL_HIGH
DECLARE EXTERNAL SUB RK4 ! 値引き渡し ByVal
DIM EX_rr(3, 3) ! rr（識別番号,軸）
Dim EX_vv(3, 3) ! vv（識別番号,軸）
LET dt_2=dt/2
CALL RK4_A(t,dt_2,m,r0,v0,EX_RK4_ans) ! 1/2 1回目
LET EX_rr(1, 1) = EX_RK4_ans(1)
LET EX_rr(1, 2) = EX_RK4_ans(2)
LET EX_rr(1, 3) = EX_RK4_ans(3)
LET EX_vv(1, 1) = EX_RK4_ans(4)
LET EX_vv(1, 2) = EX_RK4_ans(5)
LET EX_vv(1, 3) = EX_RK4_ans(6)
LET EX_rr(2, 1) = EX_RK4_ans(7)
LET EX_rr(2, 2) = EX_RK4_ans(8)
LET EX_rr(2, 3) = EX_RK4_ans(9)
LET EX_vv(2, 1) = EX_RK4_ans(10)
LET EX_vv(2, 2) = EX_RK4_ans(11)
LET EX_vv(2, 3) = EX_RK4_ans(12)
LET EX_rr(3, 1) = EX_RK4_ans(13)
LET EX_rr(3, 2) = EX_RK4_ans(14)
LET EX_rr(3, 3) = EX_RK4_ans(15)
LET EX_vv(3, 1) = EX_RK4_ans(16)
LET EX_vv(3, 2) = EX_RK4_ans(17)
LET EX_vv(3, 3) = EX_RK4_ans(18)
CALL RK4_A(t,dt_2,m,EX_rr,EX_vv,EX_RK4_ans) ! 1/2 2回目
END SUB
!---------------------------------------------------------------------------
EXTERNAL SUB RK4_A(t,dt,m(),r0(,),v0(,),EX_RK4_ans())
DECLARE EXTERNAL SUB RK4 ! 値引き渡し ByVal
CALL RK4(t,dt,m,r0,v0,EX_RK4_ans)
END SUB
!---------------------------------------------------------------------------
EXTERNAL SUB RK4(t,dt,m(),r0(,),v0(,),RK4_ans()) ! ここでは、ｍが使われていない
DECLARE EXTERNAL FUNCTION 速度変位成分_func ! 値引き渡し ByVal
DECLARE EXTERNAL FUNCTION 位置変位成分_func ! 値引き渡し ByVal
!下準備-------
Dim ret(18)
Dim r1(3, 3), r2(3, 3), r3(3, 3), r4(3, 3) ! r(識別番号,軸）
Dim v1(3, 3), v2(3, 3), v3(3, 3), v4(3, 3) ! v(識別番号,軸）
Dim k1rr(3, 3), k2rr(3, 3), k3rr(3, 3), k4rr(3, 3) ! k1rr(識別番号,軸）
Dim k1vv(3, 3), k2vv(3, 3), k3vv(3, 3), k4vv(3, 3) ! k1vv(識別番号,軸）
!DIM 番号 As Integer, 軸 As Integer
!Dim cnt_00 As Integer
FOR 番号 = 1 TO 3
FOR 軸 = 1 TO 3
LET v1(番号, 軸) = v0(番号, 軸) !v(識別番号,軸) (3,3) 3体×3次元 前処理
LET r1(番号, 軸) = r0(番号, 軸) !r(識別番号,軸) (3,3) 3体×3次元 前処理
NEXT 軸
NEXT 番号
!下準備--------
!k1
For 番号 = 1 To 3 !識別番号
For 軸 = 1 To 3 !軸成分
LET k1vv(番号, 軸) = dt * 速度変位成分_func(番号, 軸, t, r1, v1, m) !速度変位_総和 m/s^2 ★dv/dt 3体×3次元＝9連立
LET k1rr(番号, 軸) = dt * 位置変位成分_func(番号, 軸, t, r1, v1, m) !位置変位成分 m/s ★dx/dt 3体×3次元＝9連立
Next 軸
Next 番号
!k2
For 番号 = 1 To 3
For 軸 = 1 To 3
LET v2(番号, 軸) = v1(番号, 軸) + k1vv(番号, 軸) / 2 !v(識別番号,軸) (3,3) 前処理
LET r2(番号, 軸) = r1(番号, 軸) + k1rr(番号, 軸) / 2 !r(識別番号,軸) (3,3) 前処理
NEXT 軸
NEXT 番号
For 番号 = 1 To 3 !識別番号
For 軸 = 1 To 3 !軸成分
LET k2vv(番号, 軸) = dt * 速度変位成分_func(番号, 軸, t + dt / 2, r2, v2, m) !速度変位_総和 m/s^2 ★dv/dt 3体×3次元＝9連立
LET k2rr(番号, 軸) = dt * 位置変位成分_func(番号, 軸, t + dt / 2, r2, v2, m) !位置変位成分 m/s ★dx/dt 3体×3次元＝9連立
Next 軸
Next 番号
!k3
For 番号 = 1 To 3
FOR 軸 = 1 TO 3
LET v3(番号, 軸) = v1(番号, 軸) + k2vv(番号, 軸) / 2 !v(識別番号,軸) (3,3) 前処理
LET r3(番号, 軸) = r1(番号, 軸) + k2rr(番号, 軸) / 2 !r(識別番号,軸) (3,3) 前処理
Next 軸
Next 番号
For 番号 = 1 To 3 !識別番号
For 軸 = 1 To 3 !軸成分
LET k3vv(番号, 軸) = dt * 速度変位成分_func(番号, 軸, t + dt / 2, r3, v3, m) !速度変位_総和 m/s^2 ★dv/dt 3体×3次元＝9連立
LET k3rr(番号, 軸) = dt * 位置変位成分_func(番号, 軸, t + dt / 2, r3, v3, m) !位置変位成分 m/s ★dx/dt 3体×3次元＝9連立
Next 軸
Next 番号
!k4
For 番号 = 1 To 3
For 軸 = 1 To 3
LET v4(番号, 軸) = v1(番号, 軸) + k3vv(番号, 軸) !v(識別番号,軸成分) (3,3) 前処理
LET r4(番号, 軸) = r1(番号, 軸) + k3rr(番号, 軸) !r(識別番号,軸成分) (3,3) 前処理
Next 軸
Next 番号
For 番号 = 1 To 3 !識別番号
For 軸 = 1 To 3 !軸成分
LET k4vv(番号, 軸) = dt * 速度変位成分_func(番号, 軸, t + dt, r4, v4, m) !速度変位_総和 m/s^2 ★dv/dt 3体×3次元＝9連立
LET k4rr(番号, 軸) = dt * 位置変位成分_func(番号, 軸, t + dt, r4, v4, m) !位置変位成分 m/s ★dx/dt 3体×3次元＝9連立
Next 軸
Next 番号
For 番号 = 1 To 3
For 軸 = 1 To 3
LET v1(番号, 軸) = v1(番号, 軸) + 1 / 6 * (k1vv(番号, 軸) + 2 * k2vv(番号, 軸) + 2 * k3vv(番号, 軸) + k4vv(番号, 軸)) ! m/s
LET r1(番号, 軸) = r1(番号, 軸) + 1 / 6 * (k1rr(番号, 軸) + 2 * k2rr(番号, 軸) + 2 * k3rr(番号, 軸) + k4rr(番号, 軸)) ! m
Next 軸
Next 番号
!便宜的に戻り値を行列 ret(18)でまとめ
LET RK4_ans(1) = r1(1, 1)
LET RK4_ans(2) = r1(1, 2)
LET RK4_ans(3) = r1(1, 3)
LET RK4_ans(4) = v1(1, 1)
LET RK4_ans(5) = v1(1, 2)
LET RK4_ans(6) = v1(1, 3)
LET RK4_ans(7) = r1(2, 1)
LET RK4_ans(8) = r1(2, 2)
LET RK4_ans(9) = r1(2, 3)
LET RK4_ans(10) = v1(2, 1)
LET RK4_ans(11) = v1(2, 2)
LET RK4_ans(12) = v1(2, 3)
LET RK4_ans(13) = r1(3, 1)
LET RK4_ans(14) = r1(3, 2)
LET RK4_ans(15) = r1(3, 3)
LET RK4_ans(16) = v1(3, 1)
LET RK4_ans(17) = v1(3, 2)
LET RK4_ans(18) = v1(3, 3)
END SUB
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
EXTERNAL SUB 図形(rr_ex(,))
! OPTION ARITHMETIC DECIMAL_HIGH
SET WINDOW -5 , 5 , -5 , 5
DRAW AXES(1,1)
SET POINT STYLE 1
! 1 ･ 2 + 3 ＊ 4 ○ 5 × 6 ■ 7 ●
! 0白， 1黒， 2青， 3緑， 4赤， 5水色， 6黄色， 7赤紫，
! 8 灰色，9 濃い青,10 濃い緑，11 青緑， 12 えび茶，13 オリーブ色，14 濃い紫，15 銀色，･･･
SET POINT COLOR 2 ! 濃い青 Ｍ３
PLOT POINTS: rr_ex(1, 1) , rr_ex(1, 2)
SET POINT COLOR 4 ! 赤 Ｍ４
PLOT POINTS: rr_ex(2, 1) , rr_ex(2, 2)
SET POINT COLOR 1 ! 黒 Ｍ５
PLOT POINTS: rr_ex(3, 1) , rr_ex(3, 2)
END SUB
!---------------------------------------------------------------------------
!---------------------------------------------------------------------------
EXTERNAL FUNCTION 速度変位成分_func(被番号, 軸, t, r(,), v(,), m()) ! （速度変位 m/s^2）を返す
! OPTION ARITHMETIC DECIMAL_HIGH
Dim 差異(3)
LET G = 1 ! 万有引力定数
LET rrr = 0 ! 初期化
LET sum_a = 0 ! 初期化
! dv/dt = -Ｇ・Ｍ・(→r)／r^3 ・・・ 万有引力の法則
For 基準天体 = 1 To 3 !番号
If 被番号 <> 基準天体 Then ! 相違する場合だけ計算する
LET 差異(1) = r(被番号, 1) - r(基準天体, 1) !r(識別番号,軸) x成分
LET 差異(2) = r(被番号, 2) - r(基準天体, 2) !r(識別番号,軸) y成分
LET 差異(3) = r(被番号, 3) - r(基準天体, 3) !r(識別番号,軸) z成分
LET rrr = SQR((差異(1) ^ 2 + 差異(2) ^ 2 + 差異(3) ^ 2) ) !r = √(x^2+y^2+z^2)
LET sum_a = -G * m(基準天体) * 差異(軸) / rrr / rrr / rrr + sum_a !相違し、且つ、計算する軸成分だけを積算する
END IF
Next 基準天体
LET 速度変位成分_func = sum_a !加速度 １成分のみ
End Function
!---------------------------------------------------------------------------
EXTERNAL FUNCTION 位置変位成分_func(被番号,軸, t, r(,), v(,), m()) ! （位置変位 m/s）を返す
! OPTION ARITHMETIC DECIMAL_HIGH
! dx/dt = v
LET 位置変位成分_func = v(被番号, 軸) !単に速度ベクトルを返す １成分のみ
END FUNCTION
!---------------------------------------------------------------------------