●円柱周りの流れ
        円柱周りの流れ考察
    ・円柱表面の流れの速度と圧力分布
        圧力分布について考察
        レイノルズ数←新出
●高度による気圧の変化
●陰線処理を施した投影図
●ライフゲーム(シミュレーション)
●LCR直列電気回路過渡応答
● スプライン補間曲線
● ビリヤードシミュレーション
● 平面パズル 3×3
　

¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨
（※このカッコ内は、ActiveBASIC  2.5x についての話です。現在のバージョンとは互換性がありません。少し悲しい…。
◆リストとコンパイラについて

このサブページのソースリストは、BASICコンパイラ用です。
コンパイラは、フリーソフト
Active  BASIC です。
ベクター　Vector から入手できます。

　掲載ソースはN88日本語BASIC(DOS版)インタープリター用記述からの転用です。
ActiveBASICは命令語と関数名については小文字と大文字の区別は無用です。
変数名と行ラベル名については区別されます。）



――――――――――――

これは C言語欄に載せている「 ★流線形物体表面圧力分布」の基礎部分の一部に当てはまります。

■複素速度ポテンシャル
複素数 z は実部を x ,虚部を y として、
z=x + i y
と表す事ができる。ここで i は虚数 √(-1)。
流体が非圧縮(密度がどこでも一定)で、渦が無い(流れが物体から剥がれない)時、流れが２次元(奥行き方向には変化がない)なら、流れの様子は複素速度ポテンシャルで表すことが出来る。
複素速度ポテンシャル f(z) は実部 Re(f)をφ、虚部 Im(f)をψ、虚数を i とすると、
f(z)=φ(x,y) + i ψ(x,y)
表現を簡略化して
f=φ+ψi
ψの値の等しい位置を結ぶとそれが流線になる。

[flow03.bas]
100 '一様流中に置かれた静止円柱周りの流れ
110 '流れの関数Ψを使って流線を描く
120 'Ψ＝constで x,c:given  の時、ニュートン法を使ってy を得る
130 'ｆ＝Ｕ(ｚ＋ｒ^2／ｚ)　（ｚ＝ｘ＋iｙ　ｆ＝Φ＋iΨ）
140 CLEAR:SCREEN 3:CLS 3
150 SX=32:SY=20:R=5:U=10
160 N=50:DX=2*SX/N
180 CIRCLE (10*SX,10*SY),10*R,3,,,,F,6
190 FOR C=0 TO U*R-5 STEP 10
200   Y0=R:X=-SX
210     FY=U*Y^3-C*Y*Y+U*(X*X-R*R)*Y-C*X*X
220     DY=FY/(3*U*Y*Y-2*C*Y+U*(X*X-R*R))
230     ADY=ABS(DY)
240     Y1=Y-DY
250     LOCATE 0,0:PRINT USING "count:## Ψ=####    x=####.##    y=####.##  |dy|=#.##^^^^";I,C,X,Y,ADY
260     IF ADY>.0001 THEN Y=Y1:GOTO 210
270   Y=R:X1=X
280   FOR I=0 TO N
290     FY=U*Y^3-C*Y*Y+U*(X*X-R*R)*Y-C*X*X
300     DY=FY/(3*U*Y*Y-2*C*Y+U*(X*X-R*R))
310     Y=Y-DY
320     ADY=ABS(DY)
330     LOCATE 0,0:PRINT USING "count:## Ψ=####    x=####.##    y=####.##  |dy|=#.##^^^^";I,C,X,Y,ADY
340     IF ADY>.0001 THEN 290
350     LINE (10*(X1+SX),10*(-Y1+SY))-(10*(X+SX),10*(-Y+SY)),5
360     LINE (10*(X1+SX),10*(Y1+SY))-(10*(X+SX),10*(Y+SY)),5
370     LOCATE 0,2:PRINT USING "Ψ=####    x=####.##    y=####.##";C,X,Y
380     X1=X:Y1=Y:X=X+DX:Y=R
390   NEXT I
400 NEXT C
401 PRINT "終わりますか?"
402 z$=inkey$
403 if z$="" then 402
404 if z$<>"n" then end else 140

図を見れば流線は円の最も上と下で間隔が最も狭くなっている。その狭さは見た目の通り、一様流れの場所の1/2である。
流線の間隔は流れの速さを表す。その関係はどんな関係なのかを明らかにすれば、円柱周りの流れの特徴が明白になる。

まず、一般論として任意の定常流れを考える。
流れの関数ψ(x,y)の値の一定である点(x,y)を結んでできた線は流線である。時間変化しない流線は流れの軌跡と一致する。
今、上(y)方向に流れの関数ψを微分する。
dψ/dy
これは右(x)方向の流れの速さuと同じで、
u=dψ/dy...(1)
ここで dψ/dy の意味を考える。
dψ/dy は、一点(x,y)から上(y方向)に向って、ψの値がどれだけ増えるかという増加量のことである。
図で dψ/dy はまさに流線の間隔の狭さのことである。
つまり(1)より、

いかなる流れの場合でも、時間変化しなければ、
右方向の速さ u は、流線の間隔の狭さに比例する。

これは天気図の等圧線の間隔が狭い程、その場所の風速は速くなるのと等価である。

ここまで分かれば、円柱周りの流れの特徴を明らかにできる。

結論…図より分かること

円柱周りの流れについて

�@ 円柱の上下の地点(円柱近くの、流れが一様流と平行になる部分)では、流れの速さは、もとの２倍になる。(理由…流線の間隔はもとの 1/2 だから)…これは「狭くなると流れが速くなる」という連続の式をも満たしている。

�A 円柱の影響は上下方向に遠くにまで及ぶ。

円柱周りの流線を描く[flow03.bas]からの発展です。プログラムリスト内に大まかな思考の背景を書いています。まずはリスト及び実行結果をじっくりご覧下さい。

なお、[flow03.bas]の中身は描画と刻み値DXに改良を加えて踏襲しています。
★流線形物体表面圧力分布も参照してみて下さい。後ほど考察をいれ、その時に補足します。

プログラムソースリスト
[flow031.bas]
100 '一様流中に置かれた静止円柱周りの流れ
110 '流れの関数Ψを使って流線を描く
120 'Ψ＝constで x,c:given  の時、ニュートン法を使ってy を得る
125 '複素速度ポテンシャルを f とすると
130 'ｆ＝Ｕ(ｚ＋ｒ^2／ｚ)　（ｚ＝ｘ＋iｙ　ｆ＝Φ＋iΨ）
135 '続いて円柱表面の速度ベクトルを求め、求めた速度から表面の圧力分布を計算する。
140 CLEAR
150 SX=300:SY=120:R=5:U=4
160 N=24:DX=8*R/N
170 n=16:Pi=3.14159265:dQ=Pi/n
180 dim x(n+n),y(n+n),u(n+n),v(n+n),p(n+n)
190 CLS 3
200 gosub *shp
210 gosub *ln
220 gosub *velo_s
230 gosub *prs_s
240 goto *ce
250 '==================
260 *ce
270 PRINT "終わりますか?"
280 z$=inkey$
290 if z$="" then 280
300 if z$="n" then 190
310 end
320 '======================
330 *shp:'物体形状
340 CIRCLE (SX,SY),10*R,3,,,,F,6
350 return
360 *ln:'流線
370 FOR C=0 TO 16 STEP 4
380   Y=R:X=-4*R
390     FY=U*Y^3-C*Y*Y+U*(X*X-R*R)*Y-C*X*X
400     DY=FY/(3*U*Y*Y-2*C*Y+U*(X*X-R*R))
410     ADY=ABS(DY)
420     Y1=Y-DY
430     LOCATE 0,0:PRINT USING "count:## Ψ=####    x=####.##    y=####.##  |dy|=#.##^^^^";I,C,X,Y,ADY
440     IF ADY>.0001 THEN Y=Y1:GOTO 390
450   Y=R:X1=-4*R
460   FOR I=0 TO N
470     FY=U*Y^3-C*Y*Y+U*(X*X-R*R)*Y-C*X*X
480     DY=FY/(3*U*Y*Y-2*C*Y+U*(X*X-R*R))
490     Y=Y-DY
500     ADY=ABS(DY)
510     LOCATE 0,0:PRINT USING "count:## Ψ=####    x=####.##    y=####.##  |dy|=#.##^^^^";I,C,X,Y,ADY
520     IF ADY>.0001 THEN 470
530     LINE (10*X1+SX,-10*Y1+SY)-(10*X+SX,-10*Y+SY),1
540     LINE (10*X1+SX,10*Y1+SY)-(10*X+SX,10*Y+SY),1
550     X1=X:Y1=Y:X=X+DX:Y=R
560   NEXT I
570 NEXT C
580 return
590 *velo_s:'表面速度計算
595 'z = r*exp(iQ) = r*(cos(Q)+i sin(Q)) で z が円周上なら、r=R なので、
600 'u-iv(表面)= df/dz(円上)= U(1-R^2/z^2)(円上) 
= U(1-exp(-i2Q))
= U(1-cos(2Q) + i sin(2Q)) ... 0<=Q<=2π
610 line (10,100)-(10+10*U,100),5:'一様流速度ベクトル表示
620 Q=0
630 for i=0 to n+n
640  u(i)=U*(1-cos(2*Q)) : v(i)=-U*sin(2*Q)
650 x(i)=R*cos(Q):y(i)=R*sin(Q)
660 line(SX+10*x(i),SY-10*y(i))-(SX+10*(x(i)+u(i)),SY-10*(y(i)+v(i))),5
670 Q=Q+dQ
680 next i
690 return
700 *prs_s:'圧力計算
720 for i=0 to n+n
730 p(i)=1-(u(i)*u(i)+v(i)*v(i))/U/U: '無次元化された圧力
760 next i
765 '円周に沿って圧力分布表示
770 for i=1 to n+n
780 line(SX+10*x(i-1),SY-20*p(i-1))-(SX+10*x(i),SY-20*p(i)),3
790 next i
795 '横軸に角度(rad)を取って、圧力分布をグラフ化
800 x=0:dx=R*dQ
810 for i=1 to n
820 line(SX-10*R+10*x,SY+20*R-20*p(i-1))-(SX-10*R+10*(x+dx),SY+20*R-20*p(i)),3
830 x=x+dx
840 next i
845 '軸線表示
850 line(SX-10*R,SY+20*R)-(SX-10*R+10*x,SY+20*R)
860 line(SX-10*R+5*x,SY+20*R+5)-(SX-10*R+5*x,SY+20*R)
870 line(SX-10*R,SY+20*R+5)-(SX-10*R,SY+20*R)
880 line(SX-10*R+10*x,SY+20*R+5)-(SX-10*R+10*x,SY+20*R)
890 line(SX-10*R,SY+20*R+60)-(SX-10*R,SY+20*R-20)
900 line(SX-10*R,SY+20*R+20)-(SX-10*R+5,SY+20*R+20)
910 line(SX-10*R,SY+20*R+40)-(SX-10*R+5,SY+20*R+40)
920 line(SX-10*R,SY+20*R+60)-(SX-10*R+5,SY+20*R+60)
930 line(SX-10*R,SY+20*R-20)-(SX-10*R+5,SY+20*R-20)
940 return

圧力とは絶対値ではなく、ある基準位置からの
圧力差を示す相対値です。
中でも大気の圧力は空気の比重量(→※1) によって変化し、その比重量は気圧と気温によって左右されます。
ここでは1気圧=1013[hPa]を基準にして、
現実に計測の比較的簡単な海抜0mでの気圧と気温を入力し、高さが変われば、気圧がどう変化するのかを計算してグラフ化します。

■空気の比重量と温度、圧力の関係
任意の位置での気圧をp、気温を絶対温度t[K]、比重量をγとして、空気を乾燥空気に見立てれば、気体の特性式より、

p=γRt　

式の R はガス定数で
R=29.27[kgf・m/kgf/K] です。

よって、pをPo=1[気圧]にすれば、その時の温度t=To、比重量γ=δとの関係式は、

Po=δ・R・To　…(1)



■大気の圧力と比重量の関係_1_
海抜0mでの気圧をP、比重量をξとすると、これは次式で表現できます。

P/ξⁿ = Po/δⁿ ( nは定数 )…(2) 

nはpの実測値より、n=1.23　と近似できます。
 
■大気の温度と圧力の関係
海抜0mでの気圧をP、温度をT[K]とすると、(1),(2)式から、

Po=δ・R・T・(Po/P)^(n-1)/n …(3)

■大気の圧力と比重量の関係_2_

(2),(3)式より、

ξ=δ・(P/Po)^1/n …(4)

(3),(4)式の連立により、δが消え、ξは既知のPo(=1気圧),P,T で表されます。

■静止流体の圧力
流体が重力だけを受けて静止している時には、圧力 p は、微分方程式

dp/dz=-γ …(5)
(zは鉛直上向きの座標、
γは比重量で変数)

が成り立ちます。

■大気の圧力と高さの関係
上の(2),(5)式から、

p=P・(1-(n-1)・ξ・z/n/Po)^n/(n-1) …(6)

上述の通り、(3),(4)式によって、ξは既知のPo(=1気圧),P,T で表されるので、高さzの圧力pは(6)式によって算出できます。


※1 「比重量」は工学単位系であり、物体の単位体積にかかる重力のことなので、重力を考えないSI単位系の「密度」とは意味は異なるが、地上付近においては、場所によらず重力加速度は不変だと考えることが出来るので、その場合２つは一方の単位系から他方を見たとき値は等しいと考えて良い。
この関係はちょうど工学単位の「重量」とSI単位の「質量」の関係に等しい。質量1kg(SI単位)の重量は地上付近で1kgf(工学単位)だから。

[press.bas]
90 '高度による気圧の変化
100 ON ERROR GOTO *JPE
102 CLEAR
103 defint C,I,r,g,b,Z,h
104 IX=SysMetric(SM_CXSCREEN)
105 IY=SysMetric(SM_CYSCREEN)
108 WndPos OwnerWnd%, 0, 0, IX, IY
120 CLS 3
125 LINE (430,10)-(610,470),6,B
130 P0=1013
140 N=1.23: U=N/(N-1)
142 r=&h70
143 g=&h20
144 b=&hc0
150 CLS 1
151 Z=0
152 r=&h90+(r+16) mod &h90
153 g=&h90+(g-16) mod &h90
154 b=&h90+(b+16) mod &h90
158 C=RGB(r,g,b)
190 INPUT "海抜 0mの気温は [℃]";T1
200 INPUT "海抜 0mの圧力p1は [hPa] ([enter]キーだけなら、P1= 1気圧= 1013 hPa)";P1
205 IF P1=0 THEN P1=1013
210 T0=(T1+273)*(P0/P1)^(1/U):'[K]
220 W0=10332/29.27/T0/P0*P1: 'W0..基準(P0,T0)での空気の比重量 [kgf/m3]
230 W1=W0*(P1/P0)^(1/N)
240 PRINT  USING "高度と気圧の関係 ( 海抜0mでの空気比重  #.###×0.001 g/cm3)";W1
250 Pago=P1/P0
260 P1=P1/P0*10332
270 COLOR 6: PRINT "高さ";
280 COLOR 4: PRINT "[m]";
290 COLOR 5: PRINT "  気圧";
300 COLOR 4: PRINT "[気圧(atm)]";
310 COLOR 5: PRINT "  気圧";
320 COLOR 4: PRINT "[hPa]"
330 COLOR 7
340  PRINT USING "#####    ##.###         ####";Z;P;P*P0
360 FOR Z=1 TO 180
370  P=P1*(1-W1/P1/U*Z*50)^U/10332
380  PRINT USING "#####    ##.###         ####";50*Z;P;P*P0
382  LN=LN+1
390  LINE (430+Z,640-620*Pago)-(431+Z,640-620*P),C,Bf
400  IF LN=24 THEN LN=0:GOSUB *STP
405  Pago=P
410 NEXT Z
430 PRINT "もう一度 ? ([y],[ent] / else key)"
440 Z$=INKEY$
450 IF Z$="y" OR Z$="Y" or asc(Z$)=&h0d THEN 150
460 IF Z$="" THEN 440
480 END
490 *STP
500 IF INKEY$="" THEN *STP
510 CLS 1
570 COLOR 6: PRINT "高さ";
580 COLOR 4: PRINT "[m]";
590 COLOR 5: PRINT "  気圧";
600 COLOR 4: PRINT "[気圧(atm)]";
610 COLOR 5: PRINT "  気圧";
620 COLOR 4: PRINT "[hPa]"
630 COLOR 7
640 RETURN
――――――――――

リストの210-240行を関数を変えれば、任意の物体を表現できる。
【例】以下のようにすると、上図のように連なる山が表れる。

210    X(I)=-4+8*I/M
220    FOR J=1 TO N
230      Y(J)=-3+8*J/N
240      Z(I,J)=2^(-X(I)*X(I)-Y(J)*Y(J))*2+3^(-(X(I)-1)^2-(Y(J)-2)^2)+1.5^(-(X(I)+2)^2-(Y(J)-2)^2)*3

[ 3dg_btn.bas ]
100 '３次元図−ボタン形
110 WndPos OwnerWnd%,  0, 0, 640, 450
120 CLEAR
130 CT=220
140 M=17:N=16
150 WX=40:WY=100
160 PI=3.14159
190 DIM X(M),Y(N),Z(M,N),X1(M,N),Y1(M,N),Y1X(M*N),Y1N(M*N),X1X(M*N),X1N(M*N)
:'XY平面は水平面、Z軸は垂線つまり高さ方向
195　'200-250行でボタン形を関数表現している
200 FOR I=1 TO M
210    X(I)=-20+40*I/M
220    FOR J=1 TO N
230      Y(J)=-20+40*J/N
240      ZS=100-X(I)*X(I)-Y(J)*Y(J):  IF ZS>0 THEN Z(I,J)=(SQR(ZS)-(ZS-110)/11)/2
250 NEXT J,I
260 *ST:  CLS 3
270 INPUT "俯角 (deg) ";Q0
280 INPUT "回転角 (deg) ";P0
290 Q=PI/180*Q0:P=PI/180*P0:'Qは俯角　PはＺ軸回りに回転した視線の角度
300 '以下3つのLINE関数で3軸を描く
310 LINE (CT,CT)-(CT,CT-WY*COS(Q))
320 LINE (CT,CT)-(CT-WY*SIN(P),CT+WY*COS(P)*SIN(Q))
330 LINE (CT,CT)-(CT+WY*COS(P),CT+WY*SIN(P)*SIN(Q))
340 X1X=0:X1N=0
:'以下の行列 X1() と Y1() 各要素は空間の点(X,Y,Z) を描画用平面上の点(X1,Y1) に変換したもの
350 FOR I=1 TO M
360    I1=M+1-I
370    FOR J=1 TO N
380      J1=N+1-J
390      X1(I1,J1)=8*(-X(I1)*SIN(P)+Y(J1)*COS(P))
400      Y1(I1,J1)=-8*((Z(I1,J1)*COS(Q)-Y(J1)*SIN(P)*SIN(Q)-X(I1)*COS(P)*SIN(Q)))
410 NEXT J,I
430 FOR J=1 TO N-1
440      J1=N+1-J
450      LINE (CT+X1(M,J1-1),CT+Y1(M,J1-1))-(CT+X1(M,J1),CT+Y1(M,J1)),3
460 NEXT J
470 FOR J=1 TO N
480      Y1X(J)=Y1(M,J):  X1X(J)=X1(M,J)
490      Y1N(J)=Y1(M,J):  X1N(J)=X1(M,J)
500 NEXT J
510 FOR J=2 TO M
520      Y1X(J+N)=Y1(J,N):  X1X(J+N)=X1(J,N)
530      Y1N(J+N)=Y1(J,N):  X1N(J+N)=X1(J,N)
540 NEXT J
550 FOR I=1 TO M-1
560   I1=M+1-I
570   FOR J=1 TO N-1
580      J1=N+1-J
590      GOSUB *NOVIEW
600 NEXT J:NEXT I
610 PRINT "俯角回転角変更してもう一度描画するか ? (y/n) "
620 Z$=INKEY$:  IF Z$="" THEN 620
630 IF Z$="y" OR Z$="Y" THEN *ST
640 END
650 *NOVIEW :'陰線処理
660 DIS=WX:
670 FOR JS=1 TO M+N
680    DISX=ABS(X1(I1,J1)-X1X(JS))
690    DISN=ABS(X1(I1,J1)-X1N(JS))
700     IF DISX<DIS THEN JD=JS:  DIS=DISX
710     IF DISN<DIS THEN JD=JS:  DIS=DISN
720 NEXT JS
730     IF Y1(I1,J1)>Y1X(JD) THEN Y1X(JD)=Y1(I1,J1):  X1X(JD)=X1(I1,J1):              GOTO *DRW
740     IF Y1(I1,J1)<Y1N(JD) THEN Y1N(JD)=Y1(I1,J1):  X1N(JD)=X1(I1,J1):              GOTO *DRW
750 *RET:  RETURN
760 *DRW:
  LINE (CT+X1(I1,J1-1),CT+Y1(I1,J1-1))-(CT+X1(I1,J1),CT+Y1(I1,J1)),3
770 GOTO *RET
――――――――――――

ライフゲームとは、シミュレーションの古典です。
ケンブリッジ大学の数学教授コンウェイ(John Horton Conway)の考案した数学パズルの一種で、生命体モデルの単純な活動を単純な環境下で模擬実験するものです。
その内容は簡単なので生命体モデルの1個体を碁石1個とし、環境を碁盤として行うことも出来ます。この時、碁盤の升目内に碁石を置くことにします。
生命体モデルの活動内容は誕生、生存、死滅の3つだけです。

初めに適当な形に碁石を配置します。
一マスごとにマス目全てに渡って一通り、3つのうちどの活動を判定を終えた後で石を置いたり取ったりの処理を実行して、世代が1つ交代したということになります。

3つの活動に分岐する条件は、着目している1つのマス目内の碁石の有無と、そのマス目の周りに隣接する碁石の数に依存し、以下の通りです。

誕生 … 着目したマス目に碁石が無く、その周りに隣接する碁石が3個あるなら、着目したマス目に碁石を置く。

生存 … 着目したマス目に碁石があり、その周りに隣接する碁石が2個、又は3個あれば、着目したマス目にある碁石をそのまま置く。

死滅 … 着目したマス目に碁石があり、その周りに隣接する碁石が4個以上、又は1個以下なら、着目したマス目にある碁石を取り除く。

以下に簡単なプログラムのソースを載せます。生命体モデルはピンク色の●です。●の初期配置は固定しているので、ユーザが自由に配置したい場合は、
�@ソース内のデータを直接変更するか
�A実行時にマウスを使って配置の指定が出来るようにマウス入出力関連の命令をソースに追加する

のどちらかを取って下さい。�@の場合はデータ変更だけなので簡単ですが、手間かかりますので何度も変更するには�Aがいいです。

�Aのマウス関連については、スプライン補間曲線のソースを参考にしてみて下さい。それは不親切だとか、一から考えるのが面倒だ、何か教えろといった意見がありましたら、メールで一言「頼む」とだけでも送って頂ければ、�Aのソースをこのページにて増補いたします。

[lifegame.bas] (58行)
100 'Life Game
110 WNDPOS OWNERWND%,  20, 20, 600, 400
120 CLEAR:  SCREEN 3:CLS 3
125 N=15
130 DIM A(N+1,N+1)
140  A(5,5)=1: A(5,9)=1: A(6,9)=1: A(7,8)=1
145  A(4,6)=1: A(4,7)=1: A(4,8)=1: A(4,9)=1
150  A(3,5)=1: A(2,5)=1
160 P=1
165 DX=580\(N+1): DY=380\(N+1)
170 LOCATE 10,1: PRINT USING "第###  世代";P
180 GOSUB *GRAPH
185 LOCATE 11,2: PRINT USING "###  個体";Q
190 *LP0
200 FOR I=1 TO N
210 FOR J=1 TO N
220     IF (A(I,J) AND 1)=1 THEN GOSUB *BE ELSE GOSUB *BENOT
230 NEXT J: NEXT I
250 FOR I=1 TO N
260 FOR J=1 TO N
270     A(I,J)=A(I,J)\2
280 NEXT J: NEXT I
282 P=P+1
283 LOCATE 10,1: PRINT USING "第###  世代";P
285 GOSUB *GRAPH
288 LOCATE 11,2: PRINT USING "###  個体";Q
289 PRINT "終わりますか？　(y/n) "
290 Z$=INKEY$
292 IF Z$="Y" OR Z$="y" THEN end
295 IF Z$="" THEN 290
296 LOCATE 0,3: PRINT "                     " 
299 GOTO *LP0
300 END
310 *BE: S=0
320 FOR K=I-1 TO I+1
330 FOR L=J-1 TO J+1
340     IF K=I AND L=J THEN 360
350     S=S+(A(K,L) AND 1)
360 NEXT L: NEXT K
370 IF S=3 OR S=2 THEN A(I,J)=A(I,J)+2
380 RETURN
390 *BENOT: S=0
400 FOR K=I-1 TO I+1
410 FOR L=J-1 TO J+1
420     IF K=I AND L=J THEN 440
430     S=S+(A(K,L) AND 1)
440 NEXT L: NEXT K
450 IF S=3 THEN A(I,J)=A(I,J)+2
460 RETURN
470 *GRAPH
472 CLS 2: Q=0
480 FOR I=1 TO N
490 FOR J=1 TO N
500     IF (A(I,J) AND 1)=0 THEN 520
505     Q=Q+1
510     CIRCLE (DX*(J+1),380-DY*(I+1)),DY/3,3,,,,F
520 NEXT J: NEXT I
530 RETURN
――――――――――――

本プログラムは Vector アップロード済みです。
フリーソフトです。

【ホイン法】
■過渡応答とは
電気回路で電源電圧がoffからonになる(onからoffも当てはまる)と、流れる電流は、ごく短時間のうちに、時間によらず安定した定常状態になる。
onにした瞬間から、定常状態になるまでの初期の短時間を過渡期と名付ける。
過渡期での回路中の素子に流れる電流と、かかる電圧
の変化の様子を過渡応答と言う。

■数値解法
微分方程式を解くのに、数値解法は極限:dt→0
を取らない。時間刻みdt に実数δtを与えたまま計算するので、時間的に連続せず、離散している。

　誤差が出るのに、なぜ数値解法か。
微分方程式を解いて連続した関数で表現する解析解は誤差が無い、正確な値。
しかし、解けるのは数学の教科書に載るような特殊な形だけで、実際には解けないのが普通。解けないと、方程式をただ眺めるだけで、そこから先に進めない。
そこで、たとえ誤差があっても何か値が欲しい、という時に数値解法を使う事になる。
連続関数が得られないなら、その関数を表す曲線グラフに近い折れ線グラフを引こう、というのが数値解法の発想である。
ある１つの値を得られたら、決められた手順に従って次の値との差を計算し２つ目の値を得、それを使って３つ目の値との差を求めて値を得、とどんどん進めていく。
微分方程式を解いた訳ではない為、得られた値に誤差が含まれるのは数値解法の宿命である。ここで出来得る事は誤差を極力小さくする事だけである。
誤差の大小は解法の選択、刻み値と初期値の取り方に大きく左右される。
解法の種類には、基本的なニュートン法、オイラー法、応用的物理計算にはルンゲクッタ法、弾性体を扱う有限要素法、境界要素法、流体力学では差分法等がある。
ここではオイラー法から発展し、誤差を抑え、収束性を向上させたホイン法を使う。

■直列LCR直流回路
　キルヒホッフ第１法則より、直列回路の任意の位置において
電流の値は同じ i[A](i は時間の関数 i(t))である。
ここで時間 t[sec] は電源onの瞬間を0とする。
電源電圧を E[V]、
コンデンサにかかる電圧を V[V](=V(t))、
電気容量を C[F:ファラド]とすると
帯電する電荷 Q[C:クーロン]は、

Q=∫idt=C・V
なので
i=C・(dV/dt)
変形して

dV/dt =i/C

また、コイルに生じる 自己誘導起電力は、
自己インダクタンスを L とすると、

-L・(di/dt)
　
マイナスなのは電流の変化を妨げるようにはたらく為である。
これを用いると
キルヒホッフ第２法則 電圧降下の和=起電力の和 は、

R・i+V=E-L・(di/dt)
となる。変形して、

di/dt =(E-V-R・i)/L

以上、i、v について2つの微分方程式

dV/dt =i/C
di/dt =(E-V-R・i)/L

が得られた。
ちなみにこれは解く事が出来て、i、vの関数をそれぞれ求められる。あえてここでは数値解法を使い、関数は数値解の正しさを裏付ける為のものとする。

■ホイン法導出
オイラー法は式中の dtをδtにして関数 i、v を離散化したものである。
　時間進行mステップ(= m番目の値 =m・δt[sec]後の値)を
添え字 _m
を付けて表し、
h =δt
f(i) =i/C
g(v,i) =(E-V-R・i)/L
とすると、２つの式はオイラー法ではそれぞれ以下の形になる。

v_n+1=v_n+h*f(i_n)
i_n+1= i_n+h*g(v_n , i_n)

ホイン法はこれを発展させて誤差を小さくし、収束性を向上させたものであり、
以下の２段階の手続きとなる。
　 

(a)第１段階
v0_n+1=v_n+h*f(i_n)
i0_n+1= i_n+h*g(v_n , i_n)

とし、
 
(b)第２段階
K1=f(i_n)
K2=f(i0_n+1)
とおくと、

v_n+1=v_n+(K1+K2)*h/2

また、
L1=g(v_n , i_n)
L2=g(v0_n+1 , i0_n+1)
とおくと、

i_n+1=i_n+(L1+L2)*h/2

(a)からここまでを繰り返す。

[ lcr-h.bas ] (80行)
100 'LCR CIRCUIT - Heun method
110 'LCR　回路　−　ホイン法による
120 CLEAR
130 key on
135 on key goto *en
140 DEFDBL A-I,K-R,T-Z
141 DEFINT J,M
142 *init
145 CLS 3
150 C=0.0000001:R=1000:L=.5:E=2
151 M=0
152 locate 10,0: print using"LCR直列回路、コンデンサに流れる電流と、かかる電圧　（電源 ## V）";E
153 locate 10,1: print using"L:##.##H  R:##.##kΩ  C:##.##μF 　";L;R/1000;C*10^6
154 color 6: locate 46,1: print "O.K.(y/n)";
155 color 4: print"?": color 7
157 z$=inkey$: if z$="" then 157
158 if z$="n" or z$="N" then gosub *nwdata
159 locate 60,1: print"強制終了は[Q]キー"
160 H=1/10000:  T=2
170 S=5
180 N=100
182 *strt
184 cls 2
190 line (0,300)-(800,300): line(10,0)-(10,600)
200 for J=0 to 79
210  line (10*J,200)-step(5,0): line(10*J,100)-step(5,0)
220 next J
230 I=0: I0=I:  I00=I:    V=0:  V0=V:  V00=V
240 A=0: JN=N*10^M
250 for J=1 to JN 
265 Z$=inkey$:  if Z$="q" or Z$="Q" then *rpt  
257     A=A+H
260     GOSUB *CAL
263    IF ABS(I-I0)+ABS(V-V0)>10^3 THEN *retry
266    LOCATE 20,2:PRINT USING "time ##.###### sec";A
267    color 3: LOCATE 20,3: print "電流i";:color 7
268    PRINT USING "##.######^^^^ A";I
269    color 5: locate 42,3: print "電圧v";: color 7
270 print using"#.######^^^^ V";V
271  if (J mod 10^M)>0 then *jp1
272 J0=J\10^M
273   LINE ((J0-1)*S+10,300-I00*100000)-(J0*S+10,300-I*100000),3
274   LINE ((J0-1)*S+10,300-V00*100)-(J0*S+10,300-V*100),5
275 I00=I: V00=V
277 *jp1
280  I0=I:  V0=V
285 NEXT J
290 *rpt
291 locate 40,4: color 6: print "繰り返すかな　(y/n)";
292 color 4: print"?"
293 color 7
295 *stp
296 z$=inkey$
297 if z$="y" or z$="Y" then *init
298 if z$="n" or z$="N" then *en
299 goto *stp
300 *en
310 END
320 *CAL
330   VD=V0+H*I0/C
340   ID=I0+H*(E-V0-I0*R)/L
342   I1=I:   V1=V:   GOSUB *FUNC:  K1=F:  L1=G
344   I1=ID:  V1=VD:  GOSUB *FUNC:  K2=F:  L2=G
345       V=V0+H/2*(K1+K2)
346       I=I0+H/2*(L1+L2)
350 RETURN
360  *FUNC
370        F=I1/C
375        G=(E-V1-I1*R)/L
380 RETURN
390 *nwdata
400 print "  ";:input "L ( H )";L: if L=0 then L=1e-10
410 print "  ";:input "R (kΩ)";R
420 print "  ";:input "C (μF)";C: if C=0 then C=1e-4 
430 R=R*1000: C=C/10^6: return
440 *retry
445 M=M+1
447 if M>4 then *rpt
450 H=H/10
460 goto *strt



――――――――――――――

【内　容】

　マウスクリックで確定した数個の節点を通る滑らかな曲線を描くプログラムです。


【補間多項式】

　未確定の線を描くには、点と点の間の曲線を、適当な関数で表現する必要があります。これを補間と呼びます。
　補間に使う関数は、計算の容易な多項式を用います。
補間多項式は数種ありますが、ここでは曲線の振動を抑える性質を持つスプライン関数を使います。
　スプライン関数とは、隣り合う２つの節点の間で作られる区間毎に比較的低次の多項式で近似し、区間毎に多項式の係数を特定する関数です。
　スプライン関数に１次多項式を使うなら、各節点が頂と底に来るような折れ線になります。２次以上なら曲線となります。しかし２次多項式では、自由度が不足して、曲線が滑らかにはなりません。
　そこで３次スプライン関数での補間が最適となります。

n 個の節点
　( x_i , y_i )  , i=1,2,...,n
部分区間
I_i : { x_i-1 <= x <= x_i } , i=2,...,n
スプライン関数   S(x)

とする。

  S(x) の２階微分 S"(x) は、x の１次式となり、区間 I_iの右端で定数

S"(x_i) = M_i , i=1,2,...,n

となる。今、

h_i = x_i - x_i-1 , i=2,...,n

とおくと、比例配分の関係から、区間 I_i において、

S"(x) = ( M_i-1*(x_i-x) + M_i*( x-x_i-1 ) )/h_i   ... (1)

という式が得られる。ここで * は積を意味する。( BASICの表記法 )
  (1)式を２回、連続積分する。この時現れる積分定数は、
節点の y 座標が

S(x_i) = y_i

であることから得られ、関数は以下のように表すことが出来る。

S(x) = M_i-1*(x_i-x)³/6/h_i + M_i*(x-x_i-1)³/6/hi　+ (y_i-1 - M_i-1*h_i²/6)*(x_i-x)/h_i　+ (y_i - M_i*h_i²/6)*(x-x_i-1)/h_i
... (2)

　S(x) を決定するには、(2) 式中の M_i (i=1,2,...,n) が得られればよい。
　つまり M_i を求めることが目標であることが明らかになった。以下、それを行う。

　まず、 各区間 I_i における境界条件を考える。境界条件とは、今の場合、
「S(x) は節点で滑らかにつながらなければならない。」ということである。

　「滑らか」とは、１つの節点に左と右からそれぞれ近づいたとき、「点の左右での傾きが等しい」ことである。
　S(x) を微分した S'(x) は、位置 x での曲線 S(x) の傾きを表す。
　そこで節点 x_i に左と右の両側からそれぞれ近づけて (2) 式を微分する。その値は等しくならなければならず、式で表すと、

h_i*M_i-1/6 + (h_i+h_i+1)*M_i/3 + h_i+1*M_i+1/6
=(y_i+1-y_i)/h_i+1 - (y_i-y_i-1 )/h_i       i=2,3,...,n-1

となる。両辺を６倍して簡単にすると、

b_i*M_i-1 + 2M_i + c_i*M_i+1 = d_i       i=2,3,...,n-1... (3)

という式になる。なお、式中に新出の b_i, c_i, d_i は

e_i=(y_i-y_i-1)/h_i      i=2,3,...,n
c_i=h_i+1/(h_i+h_i+1)       i=2,3,...,n-1
b_i=1-c_i       i=2,3,...,n-1
d_i=6(e_i+1-e_i)/(h_i+h_i+1)       i=2,3,...,n-1

.... (a)
と定義している。

こうして n 個の未知数 M_i に対して、n-2 個の方程式が得られた。
　後２つの方程式が得られれば、未知数と方程式の個数が一致して、未知数の解が求められる。
　２つの方程式は、x₁ と x_n での境界条件を

2M₁ + c₁*M₂ = d₁
b_n*M_n-1 + 2M_n = d_n ... (4)

とすることで与えられる。

ＬＵ分解とは、下三角行列 Ｌ を

 L₁₁  0   ...    0
 L₂₁ L₂₂ 0 ... 0
 :              :
 :                0
 L_n1 L_n2 ...  L_nn


上三角行列 Ｕ を

 U₁₁ U₁₂ ...U_1n
  0   U₂₂ .. U_2n
  :            :
  0  ...  0    U_nn

と定義し、任意の行列 Ｇ を

Ｇ＝ＬＵ

というように、ＬとＵの積の形にすることである。

　Ｇ が三重対角行列 Ａ の時、ＬＵ分解すると、
Ａ＝ＬＵ の Ｌ は

 a₁  0  ...    0
 b₂ a₂ 0 .. . 0
 0  b₃ a₃ 0..0
 :               0
 0 ...  0  b_n  a_n

Ｌ内の b_i は (a) 式で定義されている。

Ｕ は

 1 z₁  0  ...   0
 0 1   z₂ 0.. .0
 :            :
 0  ... 0  1  z_n-1
 0    ...  0    1 


　残る a_i と z_i については、次のように求める。

a₁=2 ...　(7.1)
次式を i=1,2,...,n-1 で逐次計算する。
z_i = c_i/a_i ...　(7.2)
a_i+1 = 2-b_i+1*z_i ...　(7.3)

　こうして得られた

Ａ=ＬＵ

を用いて (6) 式　ＡM=D の M を求める。

【連立方程式の解法】
  ＡM=D を解くのに、 新たに

ＬG=D  ,  ＵM=G

という G を導入する。G は縦１列の行列で、第i 要素を g_i とおく。

g₁=d₁/a₁ ...　(8.1)
次式を i=1,2,...,n-1 で逐次計算する。
g_i+1 = (d_i+1-b_i+1*g_i)/a_i+1  ...　(8.2)

  (7.1)と(8.1) は添字が全て 1 であり、
　(7.2)、(7.3)、(8.2)  は繰り返しの初期値と終値が同一で、しかも繰り返し処理の内部に登場する添字は、i と i+1 だけであるので、
(7.1)〜(8.2)  はまとめることができて、以下のようになる。

a₁=2 , g₁=d₁/a₁
 次式を i=1,2,...,n-1 で逐次計算する。
z_i = c_i/a_i
a_i+1 = 2-b_i+1*z_i
g_i+1 = (d_i+1-b_i+1*g_i)/a_i+1
  ...(b)

(b) の３つの式中の b、c、d は既に (a) 式によって得られている。
だから (b) からは、新たに g_i、z_i が決まる。
　こうして得られた G 及び Z を使って、次のように解 M が求められる。

M_n = g_n
M_j = g_j-z_j*x_j+1
    j=n-1,n-2,...,2,1

......(c)

　この M を (2) 式に入れれば、求めていたスプライン関数 S(x) が得られる。


■　まとめ

　以上、スプライン関数導出に必要な式は

(2),(5),(a),(b),(c)

に集約でき、この手続きをプログラムで組めば良い。

[ Spline.bas ]  (　50行　)
100 '３次スプライン補間曲線
　　  マウスクリックでデータ点を取る
110 WndPos OwnerWnd%,  0, 0, 600, 400
120 *STRT
130 CLEAR: SCREEN 3: CLS 3
150 ON MOUSE(2) GOSUB *PT
170 MOUSE(2) ON
180 INPUT"　節点の個数は";N
200 DIM X(N),Y(N),A(N),B(N),C(N),D(N),H(N),M(N),E(N),G(N),Z(N),X0(10),Y0(10)
210 LOCATE 2,1:PRINT USING"###  ###";MouseState(0),MouseState(1):IF I<N THEN 210
220  M1=(Y(2)-Y(1))/(X(2)-X(1))-.1:  MN=(Y(N)-Y(N-1))/(X(N)-X(N-1))+. 1:  '              m1...x1での傾き　　mn...xnでの傾き  最後の+α無いと 300行 で D1=DN=0 で ov
230 MOUSE(2) OFF
240 FOR  I=1 TO N:  X(I)=X(I)/16: Y(I)=Y(I)/16: NEXT  I:' 値を縮小しないと誤差が大きくなる。そのままでは範囲が狭隘なので計算中は値を 1/16 に圧縮している
250 FOR I=2 TO N:  H(I)=X(I)-X(I-1):  E(I)=(Y(I)-Y(I-1))/H(I): NEXT I
260 FOR I=2 TO N-1:  C(I)=H(I+1)/(H(I+1)+H(I))
270                     B(I)=1-C(I):  D(I)=6*(E(I+1)-E(I))/(H(I+1)+H(I))
280 NEXT I
285 'ここから 390行迄は三重対角行列のLU分解を使った連立１次方程式の解法
290 C(1)=0:  B(N)=0:  D(1)=0:  D(N)=0
300  'C(1)=1:  B(N)=1:  D(1)=6/H(2)*((Y(2)-Y(1))/H(2)-M1):   D(N)=6/H(N)*(MN-(Y(N)-Y(N-1))/H(N)): '両端つまりx1,xnでの傾きが既知 の場合はこちら
305 A(1)=2:  G(1)=D(1)/2
310 FOR I=1 TO N-1
320     Z(I)=C(I)/A(I)
330     A(I+1)=2-B(I+1)*Z(I)
340     G(I+1)=(D(I+1)-B(I+1)*G(I))/A(I+1)
350 NEXT I:     M(N)=G(N)
360 FOR I=1 TO N-1
370     K=N-I
380     M(K)=G(K)-Z(K)*M(K+1)
390 NEXT I
400 FOR I=2 TO N
410     PSET(16*X(I-1),16*Y(I-1)),2
420     FOR J=1 TO 10
430             X0(J)=X(I-1)+J*(X(I)-X(I-1))/10
440             Y0(J)=M(I-1)*(X(I)-X0(J))^3/6/H(I)+M(I)*(X0(J)-X(I-1))^3/6/H(I)+(Y(I-1)-M(I-1)*H(I)*H(I)/6)*(X(I)-X0(J))/H(I)+(Y(I)-M(I)*H(I)*H(I)/6)*(X0(J)-X(I-1))/H(I)
450     NEXT J
460     FOR J=1 TO 10
470             LINE-(16*X0(J),16*Y0(J)),2
480     NEXT J
490 NEXT I
500  PRINT "終わりますか？　(Y/N) "
502 Z$=INKEY$
504 IF Z$="Y" OR Z$="y" THEN end
510 IF Z$="" THEN 502
515 ERASE X,Y,A,B,C,D,H,M,E,G,Z,X0,Y0
518 GOTO *STRT
520 *PT
530 I=I+1
540 LOCATE 0,I+1
550 PRINT SPC(10);
555 PRINT USING"x##  y##";I;I;
560 X(I)=MouseState(0)
570 Y(I)=MouseState(1)
580 PRINT SPC(4);X(I);Y(I)
590 CIRCLE (X(I),Y(I)),8,2,,,,F
600 RETURN



――――――――――――

◆　実行ファイルの入手について

　実行ファイルをダウンロードすると、直ちに実行できます。
　入手先は、下記リンクの Vector です。

◆　実行すると
　
　的玉は黄一つで、直接当てるよう、マウスで向きを指定し、数字キーで力加減を入力します。

  詳細は、入手ファイル展開後の txtファイルを読んで下さい。

◆　リスト中の変数について

　下記の変数は、設定値を変えることで容易に変更が出来ます。
　�@　WID (290行)…台の穴の直径。
　�A　D (290行)…球の直径。
　　　　当然、WID>D でないといけません。
　　　　球径は大きい方が、やり易いです(当たり前)。

◆  注意
　もしマウスが表示されているにも関わらず、クリックが機能していないなら、
　180行　HD=1
として下さい。

[ grgame21.bas ] (188行)

100 'get@,put@を使った
110 *ST0:CLEAR
120 B=5:PI=3.14159:'  球数B は５個。今使用は４、５番（妙だが)
130　DIM　A(550),V(B),X(B),Y(B),X0(B),Y0(B),DX(B),DY(B),DV(B),H(B),G(B),M(B),N(B),CQ(B),SQ(B),B0X(B),B1X(B),B0Y(B),B1Y(B)
170 WndPos OwnerWnd%,  0, 0, 640, 450
180 HD=0
190 ON MOUSE(4+HD) GOSUB *M5: MOUSE(4+HD) OFF
220 RANDOMIZE
230 DVX=600:DVY=350: VX0=30:VY0=30
250 VX1=VX0+DVX:VY1=VY0+DVY
290 WID=28:D=22:R=SQR(WID*WID/4+4*WID+16-100):'WIDは穴径、Dは球径
300 *ST:CLS 3:S=0:N=0:M=0
305 LINE (VX0,VX0)-(VX1,VX1),3,BF
320 WX0=VX0+WID:WY0=VY0+WID:WX1=VX1-WID:WY1=VY1-WID
325 LINE (WX0,WY0)-(WX1,WY1),4,BF
330 CIRCLE (WX1-10,WY1-10),WID/2+8,4,,,,F
340 CIRCLE (WX0+10,WY1-10),WID/2+8,4,,,,F
350 CIRCLE (WX1-10,WY0+10),WID/2+8,4,,,,F
360 CIRCLE (WX0+10,WY0+10),WID/2+8,4,,,,F
370 CIRCLE (WX0,WY0),WID/2,0,,,,F
380 CIRCLE (WX1,WY0),WID/2,0,,,,F
390 CIRCLE (WX0,WY1),WID/2,0,,,,F
400 CIRCLE (WX1,WY1),WID/2,0,,,,F
410 CIRCLE ((WX0+WX1)/2,WY0),WID/2+8,4,,,,F
420 CIRCLE ((WX0+WX1)/2,WY1),WID/2+8,4,,,,F
430 CIRCLE ((WX0+WX1)/2,WY0-10),WID/2,0,,,,F
440 CIRCLE ((WX0+WX1)/2,WY1+10),WID/2,0,,,,F
450  K=5:' K:手球の番号
460  L=4:' l:的球の番号
470  GET@(WX0+20,WY0+20)-STEP(D,D),A(300)
480  FOR I=4 TO 5
490   X0(I)=WX0+(DVX-2*WID-D)*RND(1)
500   DV(I)=-.05
510  Y0(I)=WY0+(DVY-2*WID-D)*RND(1)
520   X(I)=X0(I):Y(I)=Y0(I)
530   CIRCLE(X0(I)+D/2,Y0(I)+D/2),D/2,I+10,,,,F
535 CIRCLE(X0(I)+D/2-D/16,Y0(I)+D/2-D/16),D*5/12,I+2,,,,F
540   GET@ (X0(I),Y0(I))-STEP(D,D),A(I*100-400)
550  NEXT I
590 *LP1
600 M5=0: CLS 1: S=0
640 MOUSE(4+HD) ON
650 *LPM1
660 LOCATE 0,0:PRINT "どの方向？"
670 IF M5<1 THEN *LPM1
680 MOUSE(4+HD) OFF
690 INPUT "力は (0-50)";V(K)
700 PUT@ (LX,LY),A(300),PSET
730 *LP00
740  IF V(K)=0 THEN VK=1:GOTO *AFCL
750  P=K
760  GOSUB *KL0
770  IF G(P)+H(P)=1 THEN GOSUB *HOLE
780  GOSUB *ML2
790  GOSUB *KL1
800  DS=SQR((X0(K)-LX)^2+(Y0(K)-LY)^2)
810 SV=SQR((X0(K)+DX(K)-LX)^2+(Y0(K)+DY(K)-LY)^2)
820  IF CL0=0 OR SV<DS THEN *LP00:'未衝突で今後も衝突しないなら*LP00へ
830  IF CL=1 THEN *AFCL
840  CL=CL+1
850  GOSUB *CLASH
860  *AFCL
870  IF V(L)=0 THEN VL=1:GOTO *S
880  P=L:AF=AF+1
890  GOSUB *KL0
900  IF G(P)+H(P)>=1 THEN GOSUB *HOLE:'台縁に当たった状態なら、穴に落ちてるか
910  GOSUB *ML2
920  GOSUB *KL1
930 *S
940 IF VK+VL<2 THEN *LP00:'どちらか動いていれば、以上の手続きを繰り返す
950 P=4
960 *EN
970  PRINT "もう一度やるか　(Y/N) "
980 Z$=INKEY$
990 IF Z$="Y" OR Z$="y" THEN *ST0
1000 IF Z$<>"N" AND Z$<>"n" THEN 980
1010 ON ERROR GOTO 0
1020 END
1030 *SR1
1040 T$=""
1050 FOR I=3 TO 5
1060  N(I)=0:M(I)=0
1070 NEXT I
1080 RETURN
1090 *M5:'初期状態で玉に当たるかどうかの判定
1100 M5=M5+1
1110 MX0=MouseState(0):'マウスはオリジナルスクリーン座標なので
1120 MY0=MOUSESTATE(1)
1130 MX=MX0-D/2
1140 MY=MY0-D/2
1150 LINE(X0(K)+D/2,Y0(K)+D/2)-(MX0,MY0)
1160  A=MY-Y0(K)
1170  G=MX-X0(K)
1180  CQ(K)=G/SQR(G^2+A^2)
1190  SQ(K)=A/SQR(G^2+A^2)
1200 I=4:L=4

1210  H=ABS(A*(X0(I)-X0(K))-G*(Y0(I)-Y0(K)))/SQR(A^2+G^2):' 手玉基点の軌道直線と的玉基点との距離
1220 VCD=(X0(L)-X0(K))*G+(Y0(L)-Y0(K))*A
1230 IF H>D OR VCD<0 THEN *N1:' [OR] 以降は内積を使ったCOSの判定
1240 'まずは衝突地点の近似計算、真正値は連立方程式を解かなければいけないので
1250 CL0=CL0+1
1260 DS0=(X0(L)-X0(K))^2+(Y0(L)-Y0(K))^2
1270 DDK=SQR(DS0-D*D):'現在の手玉から、衝突時の手玉の位置迄の距離
1280 DS0=SQR(DS0)
1290 B0=ATN(DDK/D)
1300 Q1=ATN((Y0(L)-Y0(K))/(X0(L)-X0(K))):'atn は角度のみで、±180゜は不明なので
1310 IF X0(K)>X0(L) THEN Q1=Q1+PI
1320 DB=.02
1330  B1=Q1+PI-B0:'B1は右向き水平軸からの角度
1340 H1=D
1350  WHILE B1<=Q1+PI+B0
1360   CBD=D*COS(B1)
1370   SBD=D*SIN(B1)
1380   LX0=X0(L)+D*COS(B1)
1390   LY0=Y0(L)+D*SIN(B1)
1400  CIRCLE (LX0+D/2,LY0+D/2),.2,L+2
1410   H0=ABS(SQ(K)*(LX0-X0(K))-CQ(K)*(LY0-Y0(K)))
1430  IF H0<H1 THEN H1=H0:LX=LX0:LY=LY0
1440   B1=B1+DB
1450 WEND
1460 *N1
1470 PUT@ (LX,LY),A(100),PSET
1480 RETURN
1490 *CLASH:'----- 手玉と的玉の衝突
1500 P=K
1520 CB=(X0(L)-LX)/D
1530 SB=(Y0(L)-LY)/D
1540 CP=CB*CQ(K)+SB*SQ(K)
1550 SP=SB*CQ(K)-CB*SQ(K)
1560 E=CP
1570 V(L)=V(K)*(1+E)*CP/2
1580 V1=V(K)*SQR(1+(E-3)*(E+1)*CP*CP/4)
1590 CS=(V(L)-V(K)*CP)/V1
1600 SS=V(K)/V1*SP
1610 CQ(K)=-CQ(K)*(CP*CS-SP*SS)+SQ(K)*(SP*CS+CP*SS):'三角関数の合成で表現
1620 SQ(K)=-SQ(K)*(CP*CS-SP*SS)-CQ(K)*(SP*CS+CP*SS)
1630 V(K)=V1
1640 X(K)=LX+V(K)*CQ(K)
1650 Y(K)=LY+V(K)*SQ(K)
1670 CQ(L)=CB
1680 SQ(L)=SB
1690 RETURN
1700 *ML2:'----- 玉の運動表示 -------
1710  PUT@ (X(P),Y(P)),A(100*P-400),PSET
1720  PUT@ (X0(P),Y0(P)),A(300),PSET
1730  PUT@ (X(P),Y(P)),A(100*P-400),PSET
1740 RETURN
1750 *KL0:'----- 動後の玉の座標 
1760  B0X(P)=X0(P)+DX(P)-WX0
1770  B1X(P)=X0(P)+DX(P)-WX1+D
1780  B0Y(P)=Y0(P)+DY(P)-WY0
1790  B1Y(P)=Y0(P)+DY(P)-WY1+D
1800  G(P)=-(B0X(P)<0)-(B1X(P)>0):'  G,H はクッションに当たる直前で１、
1810   H(P)=-(B0Y(P)<0)-(B1Y(P)>0):'                          離れれば０
1820  N(P)=N(P)+G(P):'Nは縦縁に当たった回数
1830  M(P)=M(P)+H(P):'V は速さ（スカラー量、向き無し）
1840  DX(P)=V(P)*CQ(P)*(-1)^N(P)
1850  DY(P)=V(P)*SQ(P)*(-1)^M(P)
1860  X(P)=(X0(P)+DX(P))*(G(P) XOR 1)-(WX0-B0X(P))*(B0X(P)<0)-(WX1-D-B1X(P))*(B1X(P)>0)
1870  Y(P)=(Y0(P)+DY(P))*(H(P) XOR 1)-(WY0-B0Y(P))*(B0Y(P)<0)-(WY1-D-B1Y(P))*(B1Y(P)>0)
1880 RETURN
1890 *KL1:' ---速さの変化、前位置の退避 
1900  V(P)=(V(P)+DV(P))*.5^(G(P)+H(P)):' クッションの反発係数　0.5
1910  IF V(P)<0 THEN V(P)=0
1920  X0(P)=X(P):Y0(P)=Y(P)
1930 RETURN

 

1940 *HOLE:'---穴に落ちたかどうか判定する---
1950 HOL=14+R
1960 DSX0=X0(P)+D/5*4-WX0:DSX1=WX1-X0(P)-D/5:HOL2=(WX0+WX1)/2-WID/2-4
1970 DSY0=Y0(P)+D/5*4-WY0:DSY1=WY1-Y0(P)-D/5
1980 IF DSX0<HOL*H(P) OR DSX1<HOL*H(P) OR DSY0<HOL*G(P) OR DSY1<HOL*G(P) OR (H(P)*(X0(P)+D/4)>HOL2 AND (X0(P)+3/4*D)<(HOL2+WID+8)*H(P)) THEN GOTO *DROP
1990 RETURN
2000 *DROP
2010 PUT@ (X0(P),Y0(P)),A(300),PSET
2020 PUT@ (X(P),Y(P)),A(300),PSET
2030 PRINT "落ちた"
2040 V(P)=0
2050 GOTO *LP00



――――――――――

◆　実行ファイルの入手について

実行ファイルをダウンロードすれば、直ちに実行できます。
入手先は、下記リンクの Vector です。
フリーソフトでアップロードしています。

◆実行すると

　色パネルのかき混ぜが自動で進行します。かき混ぜを止めるにはキーをどれか押して下さい。
　パネルの操作は全てキー入力です。皆さんからマウス操作の要求が聞こえたら追加するつもりです。
　パネルの周囲に書かれた文字と同じキーを押すと、三つのパネルがその方向に一つずつずれます。枠外にはみ出すパネルは、反対側の空いた所に回り込みます。
　押したキーの文字が画面に現れます。復習する時に使って下さい。
　手数のかかる場合でも、操作回数は、ほぼ７、８回が最短のようです。

[ puzzle33.bas ] (77行)

100 '３×３マス　数字並べ替えパズル
110 CLEAR:SCREEN 3:CLS 3
120 ShowWnd OwnerWnd%, SW_MAXIMIZE
130 DIM A(4,4)
140 RANDOMIZE
150 PRINT
160 PRINT "          　に　ならべてね"
170 FOR I=1 TO 3
180  FOR J=1 TO 3
190   A(I,J)=J
200 NEXT J,I
210 GOSUB *PR
220  R=INT(6*RND(1))+7
230  GOSUB *PROC
240  R$=HEX$(6*P+Q+3-6*U)+R$
250 N=1
260 *LP0
270  R=INT(12*RND(1))+1
280  GOSUB *PROC
290  R$=HEX$(6*P+Q+3-6*U)+R$
300  N=N+1
320 LOCATE 32,7:PRINT "1    2    3"
330 LOCATE 27,9:PRINT "7                   a"
340 LOCATE 27,12:PRINT "8                   b"
350 LOCATE 27,15:PRINT "9                   c"
360 LOCATE 32,17:PRINT "4    5    6"
368 LOCATE 55,3:PRINT "| stop on any key"
375 gosub *PR
377 if inkey$="" then 220
378 LOCATE 55,4:PRINT "| stop on [Q] key"
380 *LP1:Z$=INKEY$
390 IF Z$="q" OR Z$="Q" THEN LOCATE 5,9:PRINT R$:goto *En
400 IF Z$="" THEN *LP1
401 Z=ASC(Z$)
402 IF Z<&H31 OR Z>&H63 THEN *LP1
405 PRINT USING "& &";Z$;
430 R=VAL("&h"+Z$)
440 GOSUB *PROC
450 GOTO *LP1
460 *PR
461  FOR I=1 TO 3
462   FOR J=1 TO 3
463    LINE(10+10*J,25+10*I)-(18+10*J,33+10*I),6,BF,2^(J-1)
464   NEXT J
465  NEXT I
470  FOR I=1 TO 3
480   FOR J=1 TO 3
490    LINE(190+45*J,95+45*I)-(228+45*J,133+45*I),6,BF,2^(A(I,J)-1)
500   NEXT J
520  NEXT I
530 RETURN
540 *CI
550 DI=1-2*U
560 FOR I=1+2*U TO 3-2*U STEP DI
570   A(I-DI,V)=A(I,V)
580 NEXT I
590 A(3-2*U,V)=A(4*U,V)
600 RETURN
610 *CJ
620 DJ=1-2*U
630 FOR J=1+2*U TO 3-2*U STEP DJ
640   A(V,J-DJ)=A(V,J)
650 NEXT J
660 A(V,3-2*U)=A(V,4*U)
670 RETURN
680 *PROC
690  P=R\7:Q=(R MOD 7)+P
700  U=Q\4:V=(Q MOD 4)+U
710  IF P=0 THEN GOSUB *CI ELSE GOSUB *CJ
720 GOSUB *PR
730 RETURN
740 *En
742 LOCATE 55,5:PRINT "| end on  any key"
750 if inkey$="" then 750
760 end