clear;clc
N=10000;v0=1; %粒子个数，初始速度
v(1:N)=v0; %定义所有粒子的初始速度为v0
deltv=3*v0/20;
for i=1:20
vbin(i)=(i-0.5)*deltv; %定义一个速度分布区间，用于绘图
end
subplot(3,1,1),hist(v,vbin);
xlabel('v')
ylabel('N')
title('初始速度分布')
for i=1:10*N %碰撞总次数为10*N
j=ceil(rand*N); %随机选取粒子j
k=ceil(rand*N); %随机选取粒子k
while k==j %如果粒子k和j相同，则重新选粒子k，直到k和j不是同一个粒子
k=ceil(rand*N);
end
phij=2*pi*rand;
phik=2*pi*rand;
vj=[v(j)*cos(phij);v(j)*sin(phij)]; %粒子j在x方向和y方向的初始速度
vk=[v(k)*cos(phik);v(k)*sin(phik)]; %粒子k在x方向和y方向的初始速度
vc=(vj+vk)/2; %粒子i和j组成的质心速度（粒子质量为1）
wj=vj-vc; %粒子j在质心系里的初始速度
wk=vk-vc; %粒子k在质心系里的初始速度
phic=rand*2*pi;
cosphic=cos(phic);
sinphic=sin(phic);
vj=[cosphic,-sinphic;sinphic,cosphic]*wj+vc; %实验室系下，粒子j碰撞之后的速度，包括vx分量和vy分量
vk=[cosphic,-sinphic;sinphic,cosphic]*wk+vc; %实验室系下，粒子k碰撞之后的速度，包括vx分量和vy分量
v(j)=sqrt(vj'*vj); %求出粒子j在实验室系下的速度（x方向速度和y方向速度合成）
v(k)=sqrt(vk'*vk); %求出粒子k在实验室系下的速度（x方向速度和y方向速度合成）
if i==N %输出碰撞N次之后的速度分布
subplot(3,1,2),hist(v,vbin);
xlabel('v')
ylabel('N')
title([num2str(N) '次碰撞之后，速度分布'])
end
end
subplot(3,1,3),hist(v,vbin);
xlabel('v')
ylabel('N')
title([num2str(10*N) '次碰撞之后，速度分布'])
%% 欢迎交流学习： QQ:1960009019 微信公众号：大仙一品堂
%绘制与理论值比较的曲线
nv=hist(v,vbin);
p=nv/N/deltv;
for i=1:20
f(i)=exp(-vbin(i)^2/v0^2)*2*vbin(i)/v0^2;
end
figure
plot(vbin,p,'o',vbin,f,'-')
xlabel('v')
ylabel('f(v)')
title('理论曲线比较')
