10 '************************************************************
15 '　　Runge-Kutta法のBASIC汎用サブルーチン(N88-BASIC)
20 ' 　モンキーハンティング（空気の抵抗がある場合と無い場合）の例
30 ' 　          　　　　　　　　１９８５年８月　加藤徳善
40 '************************************************************
60 '==========ルンゲ・クッタ法サブルーチン関係領域スタート===========
65 '　　　　　　　　　ユーザーは手を加える必要なし
70 DIM RKX(10,3),RKV(10,3),RKDX(10,10,3),RKR(10,10)
80 DIM RKIK(3,10,3),RKIM(3,10,3)
90 DIM RKAA(10,10),RKAL(10,10),RKAP(10,10),RKB(10),RKC(10),RKD(10,3)
100 DIM RKX0(10,3),RKMASS(10)
110 DIM RKX1(10,3),RKV1(10,3)
120 DIM RKV2(10),RKVN(10,3),RKCOEF(3),RKW(10,3),RKW2(10,3)
130 '
140 '................. R.K. COEFFICIENT .......
150 RKCOEF(0)=1 : RKCOEF(1)=2 : RKCOEF(2)=2 : RKCOEF(3)=1
160 '........................................
170 GOTO *RKMAIN
180 '--------------------------------------- 力の計算
190 *RKFORCE
200 IF RKSWA%=0 AND RKSWC%=0 GOTO *RKSKA
210 FOR RKI%=0 TO RKN%
220  IF RKSWA%=0 GOTO *RKSKA2
230  FOR RKJ%=0 TO RKI%
240   RKFW=0
250   FOR RKK%=0 TO RKDIM%
260    IF RKI%=RKJ% THEN *RKDIA
270 '                                       .......... difference
280    RKDX(RKI%,RKJ%,RKK%)=RKX1(RKI%,RKK%)-RKX1(RKJ%,RKK%)
290    RKDX(RKJ%,RKI%,RKK%)=RKX1(RKJ%,RKK%)-RKX1(RKI%,RKK%)
300    GOTO *RKNDIA
310   *RKDIA
320    RKDX(RKI%,RKJ%,RKK%)=RKX1(RKI%,RKK%)-RKX0(RKJ%,RKK%)
330   *RKNDIA
340    RKFW=RKFW+RKDX(RKI%,RKJ%,RKK%)*RKDX(RKI%,RKJ%,RKK%)
350   NEXT RKK%
360   RKR(RKI%,RKJ%)=SQR(RKFW)
370 '                                          ....... normalization
380   FOR RKK%=0 TO RKDIM%
390    IF RKR(RKI%,RKJ%)<>0 THEN RKNRM(RKI%,RKJ%,RKK%)=RKDX(RKI%,RKJ%,RKK%)/RKR(RKI%,RKJ%)
400    RKNRM(RKJ%,RKI%,RKK%)=-RKNRM(RKI%,RKJ%,RKK%)
410   NEXT RKK%
420   RKR(RKI%,RKJ%)=RKR(RKI%,RKJ%)-RKAL(RKI%,RKJ%)
430   RKR(RKJ%,RKI%)=RKR(RKI%,RKJ%)
440   IF RKAP(RKI%,RKJ%)=0 THEN RKPOW(RKI%,RKJ%)=1 : GOTO *RKPOW
450    IF RKAP(RKI%,RKJ%)=-2 THEN RKPOW(RKI%,RKJ%)=1/(RKR(RKI%,RKJ%)*RKR(RKI%,RKJ%)) : GOTO *RKPOW
460   IF RKAP(RKI%,RKJ%)=1 THEN RKPOW(RKI%,RKJ%)=RKR(RKI%,RKJ%):GOTO *RKPOW
470   RKPOW(RKI%,RKJ%)=EXP(RKAP(RKI%,RKJ%)*LOG(RKR(RKI%,RKJ%)))
480  *RKPOW
490   RKPOW(RKI%,RKJ%)=RKPOW(RKI%,RKJ%)*RKAA(RKI%,RKJ%)
500   RKPOW(RKJ%,RKI%)=RKPOW(RKI%,RKJ%)
510  NEXT RKJ%
520 *RKSKA2
530  RKFW1=0
540  FOR RKK%=0 TO RKDIM%
550   RKFW1=RKFW1+RKV1(RKI%,RKK%)*RKV1(RKI%,RKK%)
560  NEXT RKK%
570  RKVA(RKI%)=SQR(RKFW1)
580  FOR RKK%=0 TO RKDIM%
590   IF RKVA(RKI%)<>0 THEN RKVN(RKI%,RKK%)=RKV1(RKI%,RKK%)/RKVA(RKI%)
600  NEXT RKK%
610 NEXT RKI%
620 *RKSKA
630 '                                          ......... force accum.
640 FOR RKI%=0 TO RKN%
650  FOR RKK%=0 TO RKDIM%
660   RKF(RKI%,RKK%)=0
670  NEXT RKK%
680 NEXT RKI%
690 FOR RKI%=0 TO RKN%
700  FOR RKK%=0 TO RKDIM%
710   IF RKSWA%=0 GOTO *RKSWA1
720   FOR RKJ%=0 TO RKN%
730    RKF(RKI%,RKK%)=RKF(RKI%,RKK%)+RKPOW(RKI%,RKJ%)*RKNRM(RKI%,RKJ%,RKK%)
740   NEXT RKJ%
750  *RKSWA1
760   IF RKSWB%=1 THEN RKF(RKI%,RKK%)=RKF(RKI%,RKK%)+RKB(RKI%)*RKV1(RKI%,RKK%)
770   IF RKSWC%=1 THEN RKF(RKI%,RKK%)=RKF(RKI%,RKK%)+RKC(RKI%)*RKVA(RKI%)*RKVA(RKI%)*RKVN(RKI%,RKK%)
780   IF RKSWD%=1 THEN RKF(RKI%,RKK%)=RKF(RKI%,RKK%)+RKD(RKI%,RKK%)
790  RKF(RKI%,RKK%)=RKF(RKI%,RKK%)/RKMASS(RKI%)
800  NEXT RKK%
810 NEXT RKI%
820 RETURN
830 '----------------------------------------- RUNGE-KUTTA 法
840 *RK
850 FOR RKI%=0 TO RKN%
860  FOR RKK%=0 TO RKDIM%
870   RKX1(RKI%,RKK%)=RKX(RKI%,RKK%)
880   RKV1(RKI%,RKK%)=RKV(RKI%,RKK%)
890   RKW(RKI%,RKK%)=0
900   RKW1(RKI%,RKK%)=0
910  NEXT RKK%
920 NEXT RKI%
930 GOSUB *RKFORCE
940 RKW=0 : RKW1=0
950 FOR RKL%=0 TO 2
960  FOR RKI%=0 TO RKN%
970   FOR RKK%=0 TO RKDIM%
980    RKIK(RKL%,RKI%,RKK%)=RKDT*RKV1(RKI%,RKK%)
990    RKIM(RKL%,RKI%,RKK%)=RKDT*RKF(RKI%,RKK%)
1000    RKX1(RKI%,RKK%)=RKX(RKI%,RKK%)+RKIK(RKL%,RKI%,RKK%)/RKCOEF(RKL%+1)
1010    RKV1(RKI%,RKK%)=RKV(RKI%,RKK%)+RKIM(RKL%,RKI%,RKK%)/RKCOEF(RKL%+1)
1020    RKW(RKI%,RKK%)=RKW(RKI%,RKK%)+RKIK(RKL%,RKI%,RKK%)*RKCOEF(RKL%)
1030    RKW1(RKI%,RKK%)=RKW1(RKI%,RKK%)+RKIM(RKL%,RKI%,RKK%)*RKCOEF(RKL%)
1040   NEXT RKK%
1050  NEXT RKI%
1060  GOSUB *RKFORCE
1070 NEXT RKL%
1080 FOR RKI%=0 TO RKN%
1090  FOR RKK%=0 TO RKDIM%
1100   RKIK(RKL%,RKI%,RKK%)=RKDT*RKV1(RKI%,RKK%)
1110   RKIM(RKL%,RKI%,RKK%)=RKDT*RKF(RKI%,RKK%)
1120   RKX1(RKI%,RKK%)=RKX(RKI%,RKK%)+(RKW(RKI%,RKK%)+RKIK(3,RKI%,RKK%))/6
1130   RKV1(RKI%,RKK%)=RKV(RKI%,RKK%)+(RKW1(RKI%,RKK%)+RKIM(3,RKI%,RKK%))/6
1140  NEXT RKK%
1150 NEXT RKI%
1160 FOR RKI%=0 TO RKN%
1170  FOR RKK%=0 TO RKDIM%
1180   RKX(RKI%,RKK%)=RKX1(RKI%,RKK%)
1190   RKV(RKI%,RKK%)=RKV1(RKI%,RKK%)
1200  NEXT RKK%
1210 NEXT RKI%
1220 RETURN
1230 '---------------- END OF R.K. METHOD ------------------------
1240 '==================================== ユーザーズエリア　スタート
1250 '*  RKX(i,k) : 位置    ｉ：物体の番号
1260 '*  RKV(i,k) : 速度    ｋ：座標の番号（ｘ＝０，ｙ＝１，ｚ＝２）　
1270 '-------------------------
1280 *RKMAIN
1290 '-----------------------------以下のパラメーターを設定して下さい。
1300 '............... 初期設定 .......................
1310 '***** 次元 *****
1320         RKDIM%=2             : RKDIM%=RKDIM%-1
1330 '***** 時間間隔△ｔ *****
1340         RKDT=.1
1350 '***** 物体数 *****
1360         RKN%=3               : RKN%=RKN%-1
1370 '***** 位置 *****
1380 '         x              y            z
1390         RKX(0,0)=360 : RKX(0,1)=50  : RKX(0,2)=0
1400         RKX(1,0)=60  : RKX(1,1)=350 : RKX(1,2)=0
1410         RKX(2,0)=60  : RKX(2,1)=350 : RKX(2,2)=0
1420 '***** 固定点の位置（力を及ぼしあう点として各物体に１点定義できる。） *****
1430 '         x              y            z
1440         RKX0(0,0)=200 : RKX0(0,1)=0 : RKX0(0,2)=0
1450         RKX0(1,0)=300 : RKX0(1,1)=0 : RKX0(1,2)=0
1460         RKX0(2,0)=400 : RKX0(2,1)=0 : RKX0(2,2)=0
1470 '***** 速度 *****
1480 '         x              y            z
1490         RKV(0,0)=0 : RKV(0,1)=0 : RKV(0,2)=0
1500         RKV(1,0)=50: RKV(1,1)=-50 : RKV(1,2)=0
1510         RKV(2,0)=50: RKV(2,1)=-50 : RKV(2,2)=0
1520 '***** 質量 *****
1530         RKMASS(0)=3
1540         RKMASS(1)=1
1550         RKMASS(2)=1
1560 '............................................
1570 '
1580 '***** 力の設定 *****
1590 '
1600 ' 物体ｉと物体ｊの間に働く力 （重力・クーロン力・バネの力など）
1610     RKSWA%=1 :'..........スイッチ（０：切る　１：入る）
1620 '     強さ
1630         RKAA(0,0)=0
1640         RKAA(1,0)=100000! : RKAA(1,1)=0
1650         RKAA(2,0)=100000! : RKAA(2,1)=0       : RKAA(2,2)=0
1660 '     べき乗（重力，クーロン力は-2，ばねの弾性力は1，非整数にもできる）
1670         RKAP(0,0)=0
1680         RKAP(1,0)=-3      : RKAP(1,1)=0
1690         RKAP(2,0)=-3      : RKAP(2,1)=0       : RKAP(2,2)=0
1700 '     バネの長さ（通常は0）
1710         RKAL(0,0)=0
1720         RKAL(1,0)=0       : RKAL(1,1)=0
1730         RKAL(2,0)=0       : RKAL(2,1)=0       : RKAL(2,2)=0
1740 '      （物体ｉと物体ｉの間の設定は，物体ｉと固定点ｉとの間の力の設定を意味する。）
1750 '          ....................... 対称化
1760 '                                   （作用・反作用の法則）
1770             IF RKN=0 GOTO 1850
1780             FOR I%=1 TO RKN%
1790              FOR J%=0 TO I%-1
1800               RKAA(J%,I%)=RKAA(I%,J%)
1810               RKAP(J%,I%)=RKAP(I%,J%)
1820               RKAL(J%,I%)=RKAL(I%,J%)
1830              NEXT J%
1840             NEXT I%
1850 '          ......................................
1860 ' 速さに比例する力 （粘性抵抗：粘性が大きく速さの小さい場合）
1870     RKSWB%=0 :'..........スイッチ（０：切る　１：入る）
1880         RKB(0)=0
1890         RKB(1)=0
1900         RKB(2)=0
1910 ' 速さの２乗に比例する力  （慣性抵抗：粘性が小さく速さの大きい場合）
1920     RKSWC%=1 :'..........スイッチ（０：切る　１：入る）
1930         RKC(0)=0
1940         RKC(1)=0
1950         RKC(2)=-.001
1960 ' 一定の力　（地球上の重力など）
1970     RKSWD%=1 :'..........スイッチ（０：切る　１：入る）
1980 '         x                   y                   z
1990         RKD(0,0)=0  : RKD(0,1)=9.8*RKMASS(0) : RKD(0,2)=0
2000         RKD(1,0)=0  : RKD(1,1)=9.8*RKMASS(1) : RKD(1,2)=0
2010         RKD(2,0)=0  : RKD(2,1)=9.8*RKMASS(2) : RKD(2,2)=0
2020 '--------------------- パラメーター設定終了 -------------------
2030 '================ ユーザーズプログラム　スタート ==============
2040 CLS 3
2050 SCREEN 3
2060 CONSOLE ,,0
2070 VIEW (1,1)-(638,398),0,5
2080 TSTART=0
2090 TEND=9
2110 LINE (60,350)-(360,50),2
2120 FOR TT=TSTART TO TEND STEP RKDT  :'  loop start
2130   CIRCLE (RKX(0,0),RKX(0,1)),5,7
2140   CIRCLE (RKX(1,0),RKX(1,1)),2,6
2150   CIRCLE (RKX(2,0),RKX(2,1)),4,5
2190 '..................
2200   GOSUB*RK    :'   <---- ＊＊＊ Runge-Kutta サブルーチン呼び出し　＊＊＊
2205 '　（１回呼ぶごとに△ｔ後の新たな位置と速度がRKX,RKVに計算されて入ってくる。）
2210 '..................
2230 NEXT TT
2390 END