碰撞系统的matlab程序,一个基于蒙特卡罗的碰撞程序，请高手指点

%M-file, valve1.m

d=input('diameter:');

l=input('the centraline length:');

k=input('hatch length of the valve (less than radius) to the radius ratio:');

b=2*l/d;  %b是具有单位半径的管长，即中心线与半径长度之比

if b<2

error('不符合实际');

end

tnt=input('the total number of trajectories:');  %给轨迹总数tnt赋值

np=0;nd=0;nr=0; %初始化计数器

for nt=1:tnt

ranf=rand(1);rm=sqrt(ranf); %半径均匀分布

ranf=rand(1);a=2*pi*ranf; %位置均匀分布

xmi=-b/2;ymi=rm*sin(a);zmi=rm*cos(a);

[al,am,an]=dcs();  %溢出是在x方向

[s1,s2]=quad(am^2+an^2,am*ymi+an*zmi,ymi^2+zmi^2-1); %碰到入口部分圆柱壁面的距离

s=s1;

if s<0

s=s2;  %有一个正根和一个负根，s要取正根

end

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s;

s0=(b/2-1+k)/al;x0=k-1;y0=ymi+s0*am;z0=zmi+s0*an;  %碰到阀芯的距离

if s0<=s  %碰到前边阀芯的条件

i=2;

xm=x0;ym=y0;zm=z0;

else if xm<=-ym %碰到入口部分

i=-1;

else   %进入出口部分

[s1,s2]=quad(al^2+an^2,al*xmi+an*zmi,xmi^2+zmi^2-1);

s=s1;

if s2>s

s=s2; %两个根中较大的一个是合适的

end

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s; %撞击点坐标

if ym

i=1; %反射将是从出口部分

else

i=0;np=np+1;nd=nd+1; %分子直接通过弯管

end

end

end

while i~=0  %当分子溢出管道时，i会被置零

nr=nr+1;xmi=xm;ymi=ym;zmi=zm; %从xmi,ymi,zmi反射

if i==-1 %起始点在入口部分

[a,c,al]=dcs(); %a在与起始点垂直的方向上

am=-a*ymi+c*zmi;an=-a*zmi-c*ymi; %转化到轴上

s=-2*(am*ymi+an*zmi)/(am^2+an^2);

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s; %撞击点坐标

s0=(-xmi-1+k)/al;x0=k-1;y0=ymi+s0*am;z0=zmi+s0*an;  %碰到阀芯的距离

if s0<=s  %碰到前边阀芯的条件

i=2;

xm=x0;ym=y0;zm=z0;

else if xm<=-ym %碰到入口部分

if xm<=-b/2

i=0; %分子通过入口离开

end

else  %碰撞肯定是在出口部分

[s1,s2]=quad(al^2+an^2,al*xmi+an*zmi,xmi^2+zmi^2-1);

s=s1;

if s2>s

s=s2; %两个根中较大的一个是合适的

end

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s;

if ym

i=1; %下一个反射将是从出口部分

else

np=np+1;i=0; %分子通过出口离开

end

end

end

else if i==1  %起始点在出口部分

[a,c,am]=dcs(); %a在与起始点垂直的方向上

al=-a*xmi-c*zmi;an=-a*zmi+c*xmi; %旋转到轴

s=-2*(al*xmi+an*zmi)/(al^2+an^2);

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s; %撞击点坐标

s0=k/al; x0=k-1; y0=ymi+s0*am; z0=zmi+s0*an;

[st1,st2]=quad(am^2+an^2,am*ymi+an*zmi,ymi^2+zmi^2-1);

if st1

st=st1;

st1=st2;

st2=st;  %st1是两个根里较大的一个

end

xmt=xmi+al*st2;ymt=ymi+am*st2;zmt=zmi+an*st2;

if isreal(st2)==1

if st2

if xmt>k-1

i=3;  %碰到阀芯的壁面部分

xm=xmt;ym=ymt;zm=zmt;

else if xmt>-1

if s0<=st1

if s0>=st2

i=2;  %碰到阀芯前边部分

xm=x0;ym=y0;zm=z0;

else

xm=xmi+al*st1;ym=ymi+am*st1;zm=zmi+an*st1;

if xm<=-b/2

i=0;

else

i=-1;

end

end

else

xm=xmi+al*st1;ym=ymi+am*st1;zm=zmi+an*st1;

if xm<=-b/2

i=0;

else

i=-1;

end

end

end

end

else

if ym>=b/2

i=0;np=np+1;

end

end

else

if ym>=b/2

i=0;np=np+1;

end

end            else if i==2  %起始点在阀芯前边部分

[al,am,an]=dcs();al=-al;

[s1,s2]=quad(am^2+an^2,am*ymi+an*zmi,ymi^2+zmi^2-1); %碰到入口部分圆柱壁面的距离

s=s1;

if s<0

s=s2; %取正的一个根

end

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s;

if xm<=-ym

if xm<=-b/2

i=0;

else

i=-1;

end

else

[s1,s2]=quad(al^2+an^2,al*xmi+an*zmi,xmi^2+zmi^2-1);

s=s1;

if s<0

s=s2;

end

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s;

if ym>=b/2

i=0;np=np+1;

else

i=1;

end

end

else %起始点在阀芯壁面上

[a,c,al]=dcs();

am=a*ymi-c*zmi;an=a*zmi+c*ymi;

[s1,s2]=quad(al^2+an^2,al*xmi+an*zmi,xmi^2+zmi^2-1);

s=s1;

if s<0

s=s2;

end

xm=xmi+al*s;ym=ymi+am*s;zm=zmi+an*s;

if ym>=b/2

i=0;np=np+1;

else

i=1;

end

end

end

end

end

end

w=np/tnt

个人感觉出问题的应该是红色的部分，输入参数后matlab就没有反应了，好像进入死循环了，也不知道是为什么，还请各位高手指点一下，谢谢啦