HSPTV!掲示板

HSPを使わせていただいています。
授業で多粒子の衝突（反発係数が１）を説明したいと考えていますが、
壁との衝突まで出来ましたが、粒子同士の衝突が出来ません。
アルゴリズム等をご教示お願いいたします。
（obaqでやってみましたが、反発係数１が出来ませんでしたので諦めました。）
よろしくお願いいたします。

#include "hgimg3.as"
//
wx=1366:wy=768
bgscr 0,wx,wy,0,0,0
hgini
//
s=16 //円テクスチャの作成
buffer 1,s,s,0:color:boxf
color 255,255,255
circle 0,0,s,s,0
settex:tex1=stat
//
ddim px,100 //座標用
ddim py,100
ddim pvx,100 //初速度用
ddim pvy,100
//
wallx1=0
wallx2=wx
wally1=0
wally2=wy
//
v1=1.0
//

////////初期座標////
num=0
wire=0
for i,0,10 //横の並び
for j,0,10 //縦の並び
ddx=100.0
ddy=100.0
x0=50.0
y0=50.0
x=x0+ddx*i
if wire=0:y=y0+ddy*j
if wire=1:y=y0+ddy*(1.0*j+0.5)

px.num=x //初期座標
py.num=y

th=0.02*rnd(315)
pvx.num=1.0*v1*sin(th) //初速度
pvy.num=1.0*v1*cos(th)
num++
next
if wire=0:wire=1:goto *a100
if wire=1:wire=0:goto *a100
*a100
next
sum=num //粒子の総数

//////// main ////////////////////////////////////////////////
*main
hgdraw

for j,0,sum //座標を計算
px.j=px.j+pvx.j
py.j=py.j+pvy.j
next
//球の衝突の処理
//？？？？？？？？


for j,0,sum
//左右壁の衝突の処理
if px.j<=wallx1 or px.j>=wallx2 :pvx.j=-pvx.j
//上下壁の衝突の処理
if py.j<=wally1 or py.j>=wally2 :pvy.j=-pvy.j
next

//描画
for j,0,sum
gmode 2,s,s:pos px.j,py.j
hgrotate tex1,0,0,0,s,s
next
stick k,0
if k&128 : goto *owari ; [ESC]で終了

hgsync 2
goto *main
////////////////////////////////
*owari
end

お疲れ様です。こんなのでよかったらどうぞ。

#define s	16	//粒子直径
#const double r		0.5*s	//粒子半径
r2 = r*r	//半径^2
#define T_mainLoop	2	//メインループ周期[ms]
#define cReb	1.0	//反発係数


#include "hgimg3.as"
//
wx=1366:wy=768
bgscr 0,wx,wy,0,0,0
hgini
//
//円テクスチャの作成
buffer 1,s,s,0:color:boxf
color 255,255,255
circle 0,0,s,s,0
settex:tex1=stat
//
ddim px,100 //座標用
ddim py,100
ddim pvx,100 //(初)速度用
ddim pvy,100
//
wallx1=0
wallx2=wx
wally1=0
wally2=wy
//
v1=1.0
//

////////初期座標////
num=0
wire=0
for i,0,10 //横の並び
	for j,0,10 //縦の並び
		ddx=100.0
		ddy=100.0
		x0=50.0
		y0=50.0
		x=x0+ddx*i
		if wire=0:y=y0+ddy*j
		if wire=1:y=y0+ddy*(1.0*j+0.5)
		
		px.num=x //初期座標
		py.num=y
		
		th=0.02*rnd(315)
		pvx.num=1.0*v1*sin(th) //初速度
		pvy.num=1.0*v1*cos(th)
		num++
	next
	if wire=0:wire=1:goto *a100
	if wire=1:wire=0:goto *a100
	*a100
next
sum=num //粒子の総数

//////// main ////////////////////////////////////////////////
*main
hgdraw

for j,0,sum //座標を計算
	px.j=px.j+pvx.j
	py.j=py.j+pvy.j
next

/* 衝突処理 */
	repeat sum
		i=cnt
		repeat sum-i-1
			j=i+cnt+1
			/* 衝突判定 */
				dist = sqrt(powf(px(i)-px(j),2)+powf(py(i)-py(j),2))
				if (dist <= s) {	//衝突
					dup xi,px(i) : dup yi,py(i)
					dup vxi,pvx(i) : dup vyi,pvy(i)
					dup xj,px(j) : dup yj,py(j)
					dup vxj,pvx(j) : dup vyj,pvy(j)
					_X = xi-xj : _Y = yi-yj
					_Vx = vxi-vxj : _Vy = vyi-vyj
					if (absf(_Vx*_Vx+_Vy*_Vy) == 0.0) : continue	//あまりにも小さい場合は断念する
					t1 = (_X*_Vx+_Y*_Vy + sqrt(powf(_X*_Vx+_Y*_Vy,2)-(_Vx*_Vx+_Vy*_Vy)*(_X*_X+_Y*_Y-4.0*r2)))/(_Vx*_Vx+_Vy*_Vy)	//今からt1(秒? あなたの世界設定による。)前に衝突した
					/* 両者の速度を衝突面の法線,接線成分に分解する */
						inx = -_X/dist : iny = -_Y/dist	//法線方向単位ベクトル (i番ボール基準)
						itx = -iny : ity = inx	//接線方向〃
						vi_n = inx*vxi + iny*vyi	//i番ボール速度の法線成分
						vi_t = itx*vxi + ity*vyi	//〃接線成分
						vj_n = inx*vxj + iny*vyj	//j番ボール速度の法線成分
						vj_t = itx*vxj + ity*vyj	//〃接線成分
					/* 両者を接触位置に戻す */
						xi -= vxi*t1 : yi -= vyi*t1
						xj -= vxj*t1 : yj -= vyj*t1
					/* 速度を更新 */
						vi_n_new = 0.5*((1.0-cReb)*vi_n + (1.0+cReb)*vj_n)
						vj_n_new = 0.5*((1.0-cReb)*vj_n + (1.0+cReb)*vi_n)
						vxi = vi_n_new*inx + vi_t*itx : vyi = vi_n_new*iny + vi_t*ity
						vxj = vj_n_new*inx + vj_t*itx : vyj = vj_n_new*iny + vj_t*ity
					/* 移動 */
						xi += vxi*t1 : yi += vyi*t1
						xj += vxj*t1 : yj += vyj*t1
				}
		loop
	loop


for j,0,sum
	//左右壁の衝突の処理
	if px.j<=wallx1 or px.j>=wallx2 :pvx.j=-pvx.j
	//上下壁の衝突の処理
	if py.j<=wally1 or py.j>=wally2 :pvy.j=-pvy.j
next
	
//描画
for j,0,sum
	gmode 2,s,s:pos px.j,py.j
	hgrotate tex1,0,0,0,s,s
next
stick k,0
if k&128 : goto *owari ; [ESC]で終了

hgsync T_mainLoop
goto *main
////////////////////////////////
*owari
end

教えていただき、誠にありがとうございました。
お陰様で、理科の授業「物質の三態」の構成粒子の熱運動のイメージを
上手く見せる事が出来ました。感謝いたします。

#include "hgimg3.as"
//#epack "ball-sb80.bmp"
#define s 32 //16 //粒子直径
#const double r 0.5*s //粒子半径
r2 = r*r //半径^2
#define cReb 1.0 //反発係数
//
wx=1366:wy=768
bgscr 0,wx,wy,0,0,0
hgini
//
ddim px,100 //座標用
ddim py,100
ddim pvx,100 //(初)速度用
ddim pvy,100
//
wallx1=0 //壁
wallx2=wx
wally1=0
wally2=wy
v1=1.0　　//速度用
//球の読み込み
//円テクスチャの作成
texs=s
buffer 1,texs,texs,0:boxf
color 255,255,255
circle 0,0,texs,texs,1
settex:tex1=stat
//texload "ball-sb80.bmp":texs=80:tex1=stat //実際は画像を使っています。

//メニュー用
color 255,255,255
font "ＭＳ ゴシック",16
texmake 512,16:mest1 = stat:texmes "F1:固体　 　F2:液体　 　F3:気体　　↑↓:温度",mest1,0,0
/////////////////////////////////
*reset
menu=1
rc=0.0　//背景　赤
ac=1.0 //速度の倍率
///////////////////粒子の初期座標///////////////////////////////
num=0
wire=0
for i,0,16//13 //横の並び
for j,0,9//7　　//縦の並び
ddx=80.0
x0=80.0
y0=50.0
x=x0+ddx*i
if wire=0:y=y0+ddx*j
if wire=1:y=y0+ddx*(1.0*j+0.5)
px.num=x //初期座標
py.num=y
num++
next
if wire=0:wire=1:goto *a100
if wire=1:wire=0:goto *a100
*a100
next
sum=num //粒子の総数
///////////////// main ////////////////////////////////////////////////
*main
hgdraw
color 0,255,255
//hgrect 200,100,0,200,100
hgline 0+4,0+4,wx-3,0+4
hgline 0+4,wy-4,wx-3,wy-4
hgline 3,4,3,wy-4
hgline wx-3,4,wx-3,wy-4
/////////////背景の色
gc=0//green
bc=0//blue
if rc<0:rc=0.0
if rc>255:rc=255.0

col="$" +strf("%02x",int(rc)) + strf("%02x",gc) + strf("%02x",bc)
clscolor $+col
//文字テクスチャの描画
color 255,0,0:gmode 2,512,16 　
pos 300,16:hgrotate mest1

////////////////////////固体の振動
if menu=1 {
for j,0,sum
th=0.02*rnd(315)
pvx.j=1.0*ac*v1*sin(th)
pvy.j=1.0*ac*v1*cos(th)
x=px.j+pvx.j
y=py.j+pvy.j
gmode 2,texs,texs:pos x,y :hgrotate tex1,0,0,0,s,s
next
}
////////////////////////液体
if menu=2 {
for j,0,sum
th=0.02*rnd(315)
pvx.j=0.3*ac*v1*sin(th)
pvy.j=0.3*ac*v1*cos(th)
px.j=px.j+pvx.j
py.j=py.j+pvy.j
gmode 2,texs,texs:pos px.j,py.j :hgrotate tex1,0,0,0,s,s
next
}
//////////////////////////気体
if menu=3 {
for j,0,sum //座標
px.j=px.j+5.0*ac*pvx.j
py.j=py.j+5.0*ac*pvy.j

//左右壁の衝突の処理
if px.j<=wallx1+s :px.j=1.0*wallx1+s :pvx.j=-pvx.j
if px.j>=wallx2-s :px.j=1.0*wallx2-s :pvx.j=-pvx.j
//上下壁の衝突の処理
if py.j<=wally1+s :py.j=1.0*wally1+s :pvy.j=-pvy.j
if py.j>=wally2-s :py.j=1.0*wally2-s :pvy.j=-pvy.j
next

/* 衝突処理 */
repeat sum
i=cnt
repeat sum-i-1
j=i+cnt+1
/* 衝突判定 */
dist = sqrt(powf(px(i)-px(j),2)+powf(py(i)-py(j),2))
if (dist <= s) { //衝突
dup xi,px(i) : dup yi,py(i)
dup vxi,pvx(i) : dup vyi,pvy(i)
dup xj,px(j) : dup yj,py(j)
dup vxj,pvx(j) : dup vyj,pvy(j)
_X = xi-xj : _Y = yi-yj
_Vx = vxi-vxj : _Vy = vyi-vyj
if (absf(_Vx*_Vx+_Vy*_Vy) == 0.0) : continue //あまりにも小さい場合は断念する
t1 = (_X*_Vx+_Y*_Vy + sqrt(powf(_X*_Vx+_Y*_Vy,2)-(_Vx*_Vx+_Vy*_Vy)*(_X*_X+_Y*_Y-4.0*r2)))/(_Vx*_Vx+_Vy*_Vy) //今からt1(秒? あなたの世界設定による。)前に衝突した
/* 両者の速度を衝突面の法線,接線成分に分解する */
inx = -_X/dist : iny = -_Y/dist //法線方向単位ベクトル (i番ボール基準)
itx = -iny : ity = inx //接線方向〃
vi_n = inx*vxi + iny*vyi //i番ボール速度の法線成分
vi_t = itx*vxi + ity*vyi //〃接線成分
vj_n = inx*vxj + iny*vyj //j番ボール速度の法線成分
vj_t = itx*vxj + ity*vyj //〃接線成分
/* 両者を接触位置に戻す */
xi -= vxi*t1 : yi -= vyi*t1
xj -= vxj*t1 : yj -= vyj*t1
/* 速度を更新 */
vi_n_new = 0.5*((1.0-cReb)*vi_n + (1.0+cReb)*vj_n)
vj_n_new = 0.5*((1.0-cReb)*vj_n + (1.0+cReb)*vi_n)
vxi = vi_n_new*inx + vi_t*itx : vyi = vi_n_new*iny + vi_t*ity
vxj = vj_n_new*inx + vj_t*itx : vyj = vj_n_new*iny + vj_t*ity
/* 移動 */
xi += vxi*t1 : yi += vyi*t1
xj += vxj*t1 : yj += vyj*t1
}
loop
loop

//描画
for j,0,sum
gmode 2,texs,texs:pos px.j,py.j
hgrotate tex1,0,0,0,s,s
next

}
/////
hgsync 2
//////////////////////////////////////////////////////////////////////
stick k,0
if k&128 : goto *owari ; [ESC]で終了

if k&16 {
repeat
stick k,0
if k&128 : goto *owari ; [ESC]で終了
if k&16 : break ; スペースでぬける
loop
}
// ↑
drc=0.4
getkey k,38:if k&1{
if menu=1:ac=1.0*ac*1.008:rc=rc+drc:goto *main
if menu=2:ac=1.0*ac*1.009:rc=rc+drc:goto *main
if menu=3:ac=1.0*ac*1.02:rc=rc+drc:goto *main
}
//　↓
getkey k,40:if k&1{
if menu=1:ac=1.0*ac/1.008:rc=rc-drc:goto *main
if menu=2:ac=1.0*ac/1.009:rc=rc-drc:goto *main
if menu=3:ac=1.0*ac/1.02:rc=rc-drc:goto *main
}
//
getkey k,112:if k&1:menu=1:goto *reset;F1
getkey k,113:if k&1:menu=2:goto *main;F2
getkey k,114:if k&1:menu=3:goto *main ;F3

/////////////////////////////////////////////////////
goto *main
*owari
end

解決済みのチェックを忘れていました。