%(Matlab)
%The following code allows to model particules motion on the Field-Road space \R^{2}_+
%****:changeble parameters
close()
clear()
warning('off','all') %To remove all warnings
%PARAMETERS
Tmax=1e-4; %**** Maximal time
d=1.; %**** Field diffusion
D=100; %**** Road diffusion
dt=1e-6; %**** Time step (the more dt is small, the more sapce step is too)
mu=1; %**** (mu=0 => Individuals only go from Field to Road (Aspiring Road))
nu=19; %**** (mu=0 => Individuals only go from Road to Field (Aspiring Field))
hF=2*sqrt(d*dt); %Assess space step following diffusion_coef=dt/(4h) (note hF<
hF
action="field";
end
if state(j)=="field" && y(j)<=hF
action="interface";
end
if state(j)=="road"
action="road";
end
if action=="field"
ran=4*rand();
if ran<=1
x(j)=x(j)+hF;
end
if ran<=2 && ran>1
x(j)=x(j)-hF;
end
if ran<=3 && ran>2
y(j)=y(j)+hF;
end
if ran>3
y(j)=y(j)-hF;
end
end
if action=="interface"
ran=(mu+nu)*rand();
if ran<=nu
state(j)="road";
ran=2*rand();
if ran<=1
x(j)=x(j)+hR;
end
if ran>1
x(j)=x(j)-hR;
end
else
ran=4*rand();
if ran<=1
x(j)=x(j)+hF;
end
if ran<=2 && ran>1
x(j)=x(j)-hF;
end
if ran>2
y(j)=y(j)+hF;
end
end
end
if action=="road"
ran=(mu+nu)*rand();
if ran<=nu
ran=2*rand();
if ran<=1
x(j)=x(j)+hR;
end
if ran>1
x(j)=x(j)-hR;
end
else
state(j)="field";
ran=4*rand();
if ran<=1
x(j)=x(j)+hF;
end
if ran<=2 && ran>1
x(j)=x(j)-hF;
end
if ran>2
y(j)=y(j)+hF;
end
end
end
end
%PLOT(Histograms)
ax(2)=subplot(2,2,3);
if length(xRoad)==0 %if xRoad is empty...
histoR=histogram([-1000]);
XRoadmax=1;
else
histoR=histogram(xRoad');
XRoadmax=max(abs(max(xRoad)),abs(min(xRoad)));
histoR.BinWidth=(xmax-xmin)/widthCOL;
histoR.FaceColor=[1 0.2 0.2];
end
title('Particule repartition on the Road','interpreter','latex','FontSize',16);
a2=gca();
a2.XLim=[xmin xmax];
a2.YLim(1,2)=a2.YLim(1,2)+1;
a2.XTick=[xmin 0 xmax];
a2.YTick=[];
ax(3)=subplot(2,2,4);
histoF=histogram2(xField',yField');
histoF.BinWidth=[(xmax-xmin)/widthCOL ymax/widthCOL];
histoF.FaceColor=[0.8 1. 0.8];
title('Particule repartition in the Field','interpreter','latex','FontSize',16);
a3=gca();
a3.XLim=[xmin xmax];
a3.YLim=[0 ymax/2];
a3.XTick=[xmin 0 xmax];
a3.YTick=[0 ymax/2];
a3.ZTick=[];
a3.View=[ -199.6931 36.4313];
drawnow;
frame=getframe(gcf); %Save the frame for movie!
writeVideo(writerObj, frame);
end
close(writerObj);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FUNCTIONS TO SWITCH BETWEEN MATRIX AND VECTOR
function drawcircle(x,y,radiusx,radiusy,color)
tt=0:pi/12:2*pi;
xx=x+radiusx*cos(tt);
yy=y+radiusy*sin(tt);
patch(xx,yy,color)
end