%(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<<hR)
hR=2*sqrt(D*dt);

T=0:dt:Tmax;        %Time mesh

sizedots=0.008;     %**** Size of the dots when ploting
numberdots=300;    %**** Size of the population

time_multiplicator=1000;   %**** Resacling in time if needed...

widthCOL=50;        %**** Width of the columns in the histograms

xmin=-1;            %**** Size of the viewed domain
xmax=1;             %****
ymin=-0.15;         %****
ymax=0.25;          %****

Vshift=-2*ymin/3;   %To shift particules down (reps. up) if they are on the Road (resp. Field)

%INITIALIZE VECTORS OF POSITION
    %Dirac distribution
        x=zeros(1,numberdots);
        y=0.01*ones(1,numberdots);
    %Ponctual distribution
        %x=[-0.6 -0.6 0.0 0.6];
        %y=[+0.2 +0.01 +0.01 +0.01];
        
state=[];   %State of each particule (field or road)
for j=1:length(x)
    state=[state "field"];
end

%Permute randomly the colors of individuals! (esthetics)
colors0=jet(length(x));
colors=colors0;
J=randperm(length(x));
for j=1:length(x)
    k=J(j);
    colors(j,:)=colors0(k,:);
end

%SETUP THE MOVIE!
writerObj = VideoWriter('Movie.avi');
writerObj.FrameRate = 24;       %Frame rate of the film
open(writerObj)

fig=figure();                   %To get the figure full-screened
fig.WindowState='maximized';

action=[];  %Vector of size numberdots with values field/road/interface to know what kind of movement to do

%MOVING LOOP!
for i=1:length(T)
    clf();
    text_time=sprintf('%.2f', time_multiplicator*T(i));       %To force one decimal showing
    sgt=sgtitle(append('Particule motions on the $\bf{R}^{2}_{+}$-Field-Road space$\quad | \quad t=$ ',text_time,' / ',string(time_multiplicator*Tmax),'$\quad | \quad \left( \mu , \nu \right)=$ $($',string(mu),',',string(nu),'$)\quad | \quad \left( d , D \right)=$ $($',string(d),',',string(D),'$)$'),'interpreter','latex','FontSize',18);
    ax(1)=subplot(2,2,[1,2]);
    area([xmin xmax],[ymax,ymax],'FaceColor',[0.8 1. 0.8],'LineWidth',0.001,'EdgeColor','none');
    hold on;
    plot([xmin xmax],[-Vshift -Vshift],'red','LineWidth',2);
    hold on;
    a1=gca;
    a1.Box='off';
    a1.XLim=[xmin xmax];
    a1.YLim=[ymin ymax]; 
    a1.YTick=linspace(0,ymax,3);
    ratio1=2*a1.Position(3)/a1.Position(4);    

    xField=[];
    xRoad=[];
    yField=[];
    yRoad=[];
    
    for j=1:length(x)
        t=T(i);
        
	%PLOT!
        ax(1)=subplot(2,2,[1,2]);
        if state(j)=="field"
            drawcircle(x(j),y(j),sizedots,ratio1*((ymax-ymin)/(xmax-xmin))*sizedots,colors(j,:))
            xField=[xField x(j)];
            yField=[yField y(j)];
        else    %ie. state is road
            drawcircle(x(j),y(j)-Vshift,sizedots,ratio1*((ymax-ymin)/(xmax-xmin))*sizedots,colors(j,:))
            xRoad=[xRoad x(j)];
            yRoad=[yRoad y(j)];
        end
        a1.XLim=[xmin xmax]; 
        a1.YLim=[ymin ymax];
        
	%MOVE!
        if state(j)=="field" && y(j)>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