卢瑟福散射演示动画

卢瑟福散射演示动画 动力学方程：⎪⎪⎩
⎪⎪⎨⎧++=++=22222212222222212211y x y y x m e
Z Z dt y d y x x y x m e Z Z dt x d ，初始条件：⎪⎩⎪⎨⎧========;
|,0|;0|,|00000d y x dt dy v dt dx t t t t 用matlab 求解微分方程前，进行转变： ⎪⎪⎩⎪⎪⎨⎧++=++=22222222222211y x y y x C dt y d y x x y x C dt x d ；令
;42;31;2,1x y y x y dt dy x x x dt dx →→→→→→；则动力学方程写为：⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=+==+=--;34;4)42(3;12;2)42(12/3222/322x dt
dx x x x C dt dx x dt dx x x x C dt dx 初始条件为：;)0(4;10)0(2;
0)0(3;)0(10d x x x v x =-===
积分时间区段取：;25,00==tf t
figure('name','散射演示');%设置标题名字
pausetime=.002;%设置暂停时间
set(gca,'xlim',[-15 20],'ylim',[-18 18]); %设置图片显示坐标范围
set(gcf,'doublebuffer','on') %消除抖动
%axis equal
hold on
plot(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=0;tf=25;x0=[v0,-10,0,d]';
[t,x]=ode45('xdot',t0,tf,x0);
len=length(t);
for m=1:len
set(pp,'xdata',x(m,2),'ydata',x(m,4));%入射粒子的即时位置
plot(x(m,2),x(m,4),'o','markersize',1);
if m==len
plot(x(m,2),x(m,4),'ok','markersize',5);
end%把入射粒子最后停留的位置记录下来pause(pausetime); %暂停一会
drawnow
end
end
%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');
求解动力学方程组所用的函数。

xdot
function 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 入射粒子数量n
2 入射粒子初速度v0
3 入射粒子位置d
4 入射粒子动力学参数C
散射分布结果如图：取[n v0 d C]=[100 1 20*dr-10 1]
可以看出散射角大致分布区域，以及存在大角度散射。

甚至存在沿原路返回的入（出）射粒子。

对于两个靶核的散射：
增大两个靶核间距，使得入射粒子可以从中间穿过。

多靶体系（两个靶原子）
动力学方程：
⎪⎪⎩⎪⎪⎨⎧++++++-+--+=+++++-+-+=2
222202222202222222022222022)()(1)()(1)()(1)()(1d y x d y d y x m e Z Z d y x d y d y x m e
Z Z dt y d d y x x d y x m e Z Z d y x x d y x m e Z Z dt x d B A B A
初始条件不变，变换后的动力学方程组为：
⎪⎪⎪⎪⎩
⎪⎪⎪⎪⎨⎧=++++--+==+++-+=----;34);4())4(2(2)4())4(2(13;12;2))4(2(22))4(2(112/3222/3222/3222/322x dt dx d x d x x C d x d x x C dt dx x dt dx x d x x C x d x x C dt dx
两个靶原子的情况：
figure('name','2靶核散射演示');%设置标题名字
pausetime=.002;%设置暂停时间
set(gca,'xlim',[-15 20],'ylim',[-18 18]);
set(gcf,'doublebuffer','on') %消除抖动
%axis equal
hold on
plot(0,4,'color','k','marker','.','markersize',40);%画靶原子A
plot(0,-4,'color','k','marker','.','markersize',40);%画靶原子B
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=0;tf=25;x0=[v0,-10,0,d]';
[t,x]=ode45('xdo',t0,tf,x0);
len=length(t);
for m=1:len
set(pp,'xdata',x(m,2),'ydata',x(m,4));%入射粒子的即时位置
plot(x(m,2),x(m,4),'o','markersize',1);
if m==len
plot(x(m,2),x(m,4),'ok','markersize',5);
end%把入射粒子最后停留的位置记录下来
pause(pausetime); %暂停一会
drawnow
end
end
%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);
如有侵权请联系告知删除，感谢你们的配合！。