clear all
close all
lbox=50;
vs=0.05;
dt=1.;
r=1.;
rr=2;
%random fluctuation
eta = (2.*pi).*.1;
nbird=300;
nmolecules=600;
birdl=[1:nbird];
onebd=ones(1,nbird)';
moll=[1:nmolecules];
onem=ones(1,nmolecules)';
figure
axis([0 lbox 0 lbox])
axis('square')
hold on
xb=rand(nbird,1).*lbox;  %first possition
yb=rand(nbird,1).*lbox;    %first possition
ang=pi.*2.*rand(nbird,1);  %first oriantation
xm=rand(nmolecules,1).*lbox;  %first possition
ym=rand(nmolecules,1).*lbox;    %first possition
vxb = vs.*cos(ang); %first velo
vyb = vs.*sin(ang);
vxm = 3.*rand(nmolecules,1);
vym = 3.*rand(nmolecules,1);
for nsteps=1:500;
xb = xb + vxb.*dt;
yb = yb + vyb.*dt;
xm = xm + vxm.*dt;
ym = ym + vym.*dt;
for mol1 = 1:nmolecules;
%periodic boundary condition
if(xm(mol1)<0);xm(mol1)=xm(mol1)+lbox; end
if (ym(mol1)<0);ym(mol1)=ym(mol1)+lbox;end
if (xm(mol1)>lbox);xm(mol1)=xm(mol1)-lbox;end
if (ym(mol1)>lbox);ym(mol1)=ym(mol1)-lbox;end
%sep(1:nbird)=(onebd.*xb(bird1)-xb(1:nbird)).^2+ (onebd.*yb(bird1)-   yb(1:nbird)).^2;
%nearang=ang(sep<r);
%mang(bird1)=mean(nearang);
end
for bird1 = 1:nbird;
%periodic boundary condition
if(xb(bird1)<0);xb(bird1)=xb(bird1)+lbox; end
if (yb(bird1)<0);yb(bird1)=yb(bird1)+lbox;end
if (xb(bird1)>lbox);xb(bird1)=xb(bird1)-lbox;end
if (yb(bird1)>lbox);yb(bird1)=yb(bird1)-lbox;end
sep(1:nbird)=(onebd.*xb(bird1)-xb(1:nbird)).^2+ (onebd.*yb(bird1)-yb(1:nbird)).^2;
sepm(1:nmolecules)=(onem.*xb(bird1)-xm(1:nmolecules)).^2+     (onem.*yb(bird1)-ym(1:nmolecules)).^2;
nearang=ang(sep<r);
nearm=sepm<rr^2
xmeanm=dot(nearm,xm)/(sum(nearm));
ymeanm=dot(nearm,ym(1:nmolecules))/(sum(nearm));
mang2(bird1)= atan(ymeanm/xmeanm);
mang1(bird1)=mean(nearang);
mang = (mang1(bird1)+ mang2(bird1) )/2;
end
ang=mang'+eta.*(rand(nbird,1)-0.5);
vxb = vs.*cos(ang);
vyb = vs.*sin(ang);
vxm = 3.*(rand(nmolecules,1)-0.5);
vym =3.*(rand(nmolecules,1)-0.5);
cla
set(gcf,'doublebuffer','on')
plot(xb,yb,'.b')
drawnow
end