1、卢瑟福散射演示动画动力学方程： ，初始条件：22212 22212 yxmeZdtyxtx ;|,0|0dyxtvdt用 matlab 求解微分方程前，进行转变：；令2222221yxCdtyxtx；则动力学方程写为：;42;31,xyxydt ;34;4)2(;1;2)42(2/3/3xdtxCxdtx初始条件为： ;)0(4;1)(230dxxv积分时间区段取： 25,tffigure(name,散射演示);%设置标题名字pausetime=.002;%设置暂 停时间set(gca,xlim,-15 20,ylim,-18 18); %设置图片显示坐标范围set(gcf,doublebuf
2、fer,on) %消除抖动%axis equalhold onplot(0,0,color,k,marker,.,markersize,40);%画靶原子line(-10 -10,-10 10,LineStyle,-); %入射粒子初始位置范围u=0:0.01*pi:2*pi;X=14*cos(u);Y=14*sin(u);plot(X,Y,-k); %环形探测器len=0;for n=1:100 %入射粒子数量v0=1; %入射粒子的初速度dr=rand;d=20*dr-10; %入射粒子的位置pp=plot(0,d,r,marker,o,markersize,5);%入射粒子的初始位置t0
3、=0;tf=25;x0=v0,-10,0,d;t,x=ode45(xdot,t0,tf,x0);len=length(t);for m=1:lenset(pp,xdata,x(m,2),ydata,x(m,4);%入射粒子的即时位置plot(x(m,2),x(m,4),o,markersize,1);if m=lenplot(x(m,2),x(m,4),ok,markersize,5);end %把入射粒子最后停留的位置 记录 下来pause(pausetime); %暂停一会drawnowendend%subplot(1,3,1);%plot(t,x(:,1),:b,t,x(:,2),-r)
4、;%subplot(1,3,2);%plot(t,x(:,3),:b,t,x(:,4),-r);%subplot(1,3,3);%plot(x(:,2),x(:,4),-b);求解动力学方程组所用的函数。xdotfunction xd=xdot(t,x)C=1;xd=zeros(4,1);xd(1)=C*(x(2)2+(x(4)2)(-1.5)*x(2);xd(2)=x(1);xd(3)=C*(x(2)2+(x(4)2)(-1.5)*x(4);xd(4)=x(3);可调整参数为：1 入射粒子数量 n2 入射粒子初速度 v03 入射粒子位置 d4 入射粒子动力学参数 C散射分布结果如图：取n v
5、0 d C=100 1 20*dr-10 1可以看出散射角大致分布区域，以及存在大角度散射。甚至存在沿原路返回的入（出）射粒子。对于两个靶核的散射：增大两个靶核间距，使得入射粒子可以从中间穿过。多靶体系（两个靶原子）动力学方程： 222022202 222022202 )()(1)()(1 dyxdyxmeZdyxdyxmeZdty xxtx BA初始条件不变，变换后的动力学方程组为： ;34 );4()4(2)4()4(21; ;2)4(2)4(21 2/32/3 /32/3xdt dxdxCdxdCxdt xdxxd两个靶原子的情况：figure(name,2靶核散射演示 );%设置标题名
6、字pausetime=.002;%设置暂 停时间set(gca,xlim,-15 20,ylim,-18 18);set(gcf,doublebuffer,on) %消除抖动%axis equalhold onplot(0,4,color,k,marker,.,markersize,40);%画靶原子Aplot(0,-4,color,k,marker,.,markersize,40);%画靶原子Bline(-10 -10,-10 10,LineStyle,-);u=0:0.01*pi:2*pi;X=14*cos(u);Y=14*sin(u);plot(X,Y,-k);len=0;for n=1
7、:100v0=1; %入射粒子的初速度dr=rand;d=20*dr-10; %入射粒子的位置pp=plot(0,d,r,marker,o,markersize,5);%入射粒子的初始位置t0=0;tf=25;x0=v0,-10,0,d;t,x=ode45(xdo,t0,tf,x0);len=length(t);for m=1:lenset(pp,xdata,x(m,2),ydata,x(m,4);%入射粒子的即时位置plot(x(m,2),x(m,4),o,markersize,1);if m=lenplot(x(m,2),x(m,4),ok,markersize,5);end %把入射粒子
8、最后停留的位置 记录 下来pause(pausetime); %暂停一会drawnowendend%subplot(1,3,1);%plot(t,x(:,1),:b,t,x(:,2),-r);%subplot(1,3,2);%plot(t,x(:,3),:b,t,x(:,4),-r);%subplot(1,3,3);%plot(x(:,2),x(:,4),-b);用于求解微分方程组的函数：function xdd=xdo(t,x)C1=1;C2=1;dd=4;xdd=zeros(4,1);xdd(1)=C1*(x(2)2+(x(4)-dd)2)(-1.5)*x(2)+C2*(x(2)2+(x(4)+dd)2)(-1.5)*x(2);xdd(2)=x(1);xdd(3)=C1*(x(2)2+(x(4)-dd)2)(-1.5)*(x(4)-dd)+C2*(x(2)2+(x(4)+dd)2)(-1.5)*(x(4)+dd);xdd(4)=x(3);
