区間許容誤差　 ： 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


比較参照データ：情報処理学会研究報告：Vol.2018-HPC-164 No.1 ：平山　弘


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
!---------------------------------------------------------------------------