%陳金剛3100104568工程力學1001?修正版
clear?all；
N=50;%球數，可變?
rad=rand（N1）*2+4;%球半徑，[4，6]?
pos=[rand（N1）*90+5?rand（N1）*190+5];%初始位置：左半邊區域?
vel=rand（N2）*36-18;%各球初始速度==可變
color=rand（N3）;%各球顏色，隨機產生
mass=5*rad.^2;%各球質量，與半徑的平方成正比
Left=0;??%左邊界，下同
Right=200;
Up=200;
Down=0;
figure;
axis?manual;
axis?equal?square;?%固定坐標
axis（[0?200?0?200]）;
hold?on;
%===============畫不重疊的N個小球====
for?i=2:N
????index=0;
????while?index==0
?????????index2=0;
?????????for?j=1:i-1
????????????pp（j:）=[pos（i1）-pos（j1）pos（i2）-pos（j2）];
????????????rr（j）=rad（i）+rad（j）;?
????????????ppp（j）=norm（pp（j:））;
????????????if?ppp（j）?????????????????%修正一開始若球相切時，后面判斷可能誤認為相撞
????????????index2=1;break;
????????????end
?????????end
????????if?~index2
????????????index=1;
?????????else??pos（i:）=[rand（）*90+5?rand（）*190+5];
????????end
?????end
end
%====================================================================
couleft=0;%左壁面撞擊的總mv，以下類似
couright=0;
couup=0;
coudown=0;
dt=0;%最小時間
qq=0;%用于循環次數控制
while?qq<1000??
?%=========以下各球兩兩碰撞最小時間計算
?????k=1;
????for?i=1:N
????????for?j=i:N
??????????dis=[pos（j1）-pos（i1）pos（j2）-pos（i2）];
??????????vv=[vel（j1）-vel（i1）vel（j2）-vel（i2）];
??????????dist=norm（dis）;
??????????rr=rad（i）+rad（j）;
??????????cosAlpha?=abs（?sqrt（1-（rr/dist）^2））;
??????????cosTheta=（dot（disvv）/norm（vv）/dist）;
????????????if?cosTheta>=cosAlpha?&&?cosTheta<1?
????????????????dd=dist*cosTheta-sqrt（rr^2-（dist*sqrt（1-cosTheta^2））^2）;??
????????????????time（k）=dd/norm（vv）;
????????????????????k=k+1;
????????????end????
???????end
????end
?????tball=min（time）;
?%================各球碰墻最小時間計算
??for?i=1:N
???if?vel（i1）>0
??????tx（i）=（Right-pos（i1）+rad（i））/vel（i1）;%易錯點，少了半徑
???else?
??????tx（i）=（Left-pos（i1）+rad（i））/vel（i1）;
????end
???if?vel（i2）?>0
??????ty（i）=（Up-pos（i2）+rad（i））/vel（i2）;
???else?
??????ty（i）=（Down-pos（i2）+rad（i））/vel（i2）;
???end
???end
???twall=min（txty）;
???%====
???tBallWall=min（tballtwall）;
??dt=dt+tBallWall;
??%====與墻相撞改變速度
for?i=1:N
???if?pos（i1）-rad（i）<=Left?
???????pos（i1）=rad（i）;%修正邊
?????couleft=couleft+mass（i）*abs（vel（i1））;
???????vel（i1）=-vel（i1）;
?????elseif?pos（i1）+rad（i）>=Right?
?????????pos（i1）=200-rad（i）;%修正邊
????????couright=couright+mass（i）*abs（vel（i1））;
???????vel（i1）=-vel（i1）;
???????end
???if??pos（i2）-rad（i）<=Down?
???????pos（i2）=rad（i）;
?????coudown=coudown+mass（i）*abs（vel（i2））;
???????vel（i2）=-vel（i2）;
???elseif?pos（i2）+rad（i）>=Up?
???????pos（i2）=200-rad（i）;
???????couup=couup+mass（i）*abs（vel（i2））;
???????vel（i2）=-vel（i2）;
???end
end
%==========兩球碰撞改變速度
for?i=1:N
??for?j=i+1:N
??????twoBall=[pos（i1）-pos（j1）pos（i2）-pos（j2）];
??????D=norm（twoBall）;
??if?D-rad（i）-rad（j）<=0;
???????if?D??????????DRt=（rad（i）+rad（j）-D）/norm（[vel（i1）-vel（j1）?vel（i2）-vel（j2）]）;
??????????pos（i1）=pos（i1）-DRt*vel（i1）;
??????????pos（i2）=pos（i2）-DRt*vel（i2）;
????????

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件??????422291??2012-09-06?13:45??氣體擴散模擬實驗報告.pdf
?????文件????????4157??2012-09-07?01:24??a.m