%(Matlab)
%The following code allows to solve the Heat equation on the Field-Road space \R^{2}_{+} starting from the initial data (V_0,V_0).
%May be improved:
	%BDF2 scheme for time
	%Reactions function for the Field and the Road

%****: changeable parameters

clear()
close()

warning('off','all')        %To remove all warnings

%TIME PARAMETERS
Tmax=1000.;        %Final time****
dt=0.1;         %Time step****
TT=0:dt:Tmax;   %Time mesh
%DurationMovie=3000;   %Duration of the movie (in seconds)****

%SPACE PARAMETERS
H=100.;          %Height****
m=30;           %Number of vertical subdivisions (line)****
L=5*H;          %Width
n=5*m;          %Number of horizontal subdivisions (column)
h=H./m;         %Space step
XX=linspace(-L./2,L/2,n); %Horizontal space mesh
YY=linspace(H,0,m);       %Vertical space mesh

%PARAMETERS OF THE ROAD-FIELD SPACE
d = 1.;   %Diffusion intensity in the Field****
D = 2.;  %Diffusion intensity on the Road****
mu = 1.;  %Migratory equilibrium between Road and Field****
nu = 1.;  %Migratory equilibrium between Road and Field****

%SOME USEFULL CONSTANTS...
beta=d*dt/((h)^2);
gamma=D*dt/((h)^2);

%INITIAL DATUM
u0=@(x)indicator(x,-10,10);
v0=@(x,y)indicator(x,-10,10)*indicator(y,-10,10);

U0=zeros(1,n);
for j=1:n
    U0(1,j)=u0(XX(j));
end

V0=zeros(m,n);
for j=1:n
for i=1:m
    V0(i,j)=v0(XX(j),YY(i));
end
end

maxU0=max(U0);          %Get the max of U0
maxV0=max(max(V0));     %Get the max of V0

%We combine the two initial data (U0,V0) in a vector W for implicit resolution
W0=combineUV(m,n,U0,V0);    %Combine U0 and V0 in W0
W=W0;                   %Initialize the vector solution W


%CREATE THE LINEAR SYSTEM
    A=zeros((m+1)*n,(m+1)*n);
    TEST=zeros(m+1,n);  %Check if the index are good... (Take n=3 at least...)

    %Set the 6 corners
        %UP_LEFT_CORNER for V (Neumann) (1)
            i=1;
            j=1;
            TEST(i,j)=1;
            k=mat2vect(n,i,j);
            A(k,mat2vect(n,i,j))=1+4*beta;
            A(k,mat2vect(n,i+1,j))=-2*beta;
            A(k,mat2vect(n,i,j+1))=-2*beta;
        %UP_RIGHT_CORNER for V (Neumann) (2)
            i=1;
            j=n;
            TEST(i,j)=2;
            k=mat2vect(n,i,j);
            A(k,mat2vect(n,i,j))=1+4*beta;
            A(k,mat2vect(n,i+1,j))=-2*beta;
            A(k,mat2vect(n,i,j-1))=-2*beta;
        %DOWN_LEFT_CORNER for V (NEUMANN+INTERFACE) (3)
            i=m;
            j=1;
            TEST(i,j)=3;
            k=mat2vect(n,i,j);
            A(k,mat2vect(n,i,j))=1+4*beta+2*h*beta*nu/d;
            A(k,mat2vect(n,i,j+1))=-2*beta;
            A(k,mat2vect(n,i-1,j))=-2*beta;
            A(k,mat2vect(n,i+1,j))=-2*h*beta*nu/d;
        %DOWN_RIGHT_CORNER for V (NEUMANN+INTERFACE) (4)
            i=m;
            j=n;
            TEST(i,j)=4;
            k=mat2vect(n,i,j);
            A(k,mat2vect(n,i,j))=1+4*beta+2*h*beta*nu/d;
            A(k,mat2vect(n,i,j-1))=-2*beta;
            A(k,mat2vect(n,i-1,j))=-2*beta;
            A(k,mat2vect(n,i+1,j))=-2*h*beta*nu/d;
        %DOWN_LEFT_CORNER for U (NEUMANN+INTERFACE) (5)
            i=m+1;
            j=1;
            TEST(i,j)=5;
            k=mat2vect(n,i,j);
            A(k,mat2vect(n,i,j))=1+2*gamma+mu*dt;
            A(k,mat2vect(n,i,j+1))=-2*gamma;
            A(k,mat2vect(n,i-1,j))=-nu*dt;
        %DOWN_RIGHT_CORNER for U (NEUMANN+INTERFACE) (6)
            i=m+1;
            j=n;
            TEST(i,j)=6;
            k=mat2vect(n,i,j);
            A(k,mat2vect(n,i,j))=1+2*gamma+mu*dt;
            A(k,mat2vect(n,i,j-1))=-2*gamma;
            A(k,mat2vect(n,i-1,j))=-nu*dt;
    %Set the 5 edges
        %UP_EDGE for V (NEUMANN) (7)
            i=1;
            for j=2:(n-1)
                TEST(i,j)=7;
                k=mat2vect(n,i,j);
                A(k,mat2vect(n,i,j))=1+4*beta;
                A(k,mat2vect(n,i,j+1))=-beta;
                A(k,mat2vect(n,i,j-1))=-beta;
                A(k,mat2vect(n,i+1,j))=-2*beta;
            end
        %LEFT_EDGE for V (NEUMANN) (8)
            j=1;
            for i=2:(m-1)
                TEST(i,j)=8;
                k=mat2vect(n,i,j);
                A(k,mat2vect(n,i,j))=1+4*beta;
                A(k,mat2vect(n,i+1,j))=-beta;
                A(k,mat2vect(n,i-1,j))=-beta;
                A(k,mat2vect(n,i,j+1))=-2*beta;
            end
        %RIGHT_EDGE for V (NEUMANN) (9)
            j=n;
            for i=2:(m-1)
                TEST(i,j)=9;
                k=mat2vect(n,i,j);
                A(k,mat2vect(n,i,j))=1+4*beta;
                A(k,mat2vect(n,i+1,j))=-beta;
                A(k,mat2vect(n,i-1,j))=-beta;
                A(k,mat2vect(n,i,j-1))=-2*beta;
            end
        %DOWN_EDGE for V (INTERFACE) (10)
            i=m;
            for j=2:(n-1)
                TEST(i,j)=10;
                k=mat2vect(n,i,j);
                A(k,mat2vect(n,i,j))=1+4*beta+2*h*beta*nu/d;
                A(k,mat2vect(n,i,j+1))=-beta;
                A(k,mat2vect(n,i,j-1))=-beta;
                A(k,mat2vect(n,i-1,j))=-2*beta;
                A(k,mat2vect(n,i+1,j))=-2*h*beta*mu/d;
            end
        %DOWN_EDGE for U (INTERFACE) (11)
            i=m+1;
            for j=2:(n-1)
                TEST(i,j)=11;
                k=mat2vect(n,i,j);
                A(k,mat2vect(n,i,j))=1+2*gamma+mu*dt;
                A(k,mat2vect(n,i,j+1))=-gamma;
                A(k,mat2vect(n,i,j-1))=-gamma;
                A(k,mat2vect(n,i-1,j))=-nu*dt;
            end
    %Set the interior for V (12)
            for i=2:(m-1)
            for j=2:(n-1)
                TEST(i,j)=12;
                k=mat2vect(n,i,j);
                A(k,mat2vect(n,i,j))=1+4*beta;
                A(k,mat2vect(n,i,j+1))=-beta;
                A(k,mat2vect(n,i,j-1))=-beta;
                A(k,mat2vect(n,i+1,j))=-beta;
                A(k,mat2vect(n,i-1,j))=-beta;
            end
        end
%LINEAR SYSTEM HAS BEEN CREATED!

B=inv(A);   %Invert the linear system!

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


TLOG=[];
UMAXLOG=[];
VMAXLOG=[];

%SETUP THE MOVIE!
writerObj = VideoWriter('Movie.avi');
%writerObj.FrameRate = floor(length(TT)/DurationMovie); %speed of the film !
writerObj.FrameRate = 2; %Frame rate of the film
open(writerObj)

%SOLVING LOOP!
for t=1:length(TT)
if mod(t,10)==1 | t<1000 
    %PLOT SOLUTION
        clf();          %Clean the plot
        maxU=round(max(U0)*1000)/1000;
        maxV=round(max(max(V0))*1000)/1000;
        TLOG=[TLOG TT(t)];
        UMAXLOG=[UMAXLOG max(U0)];
        VMAXLOG=[VMAXLOG max(max(V0))];
        sgt=sgtitle(append('Heat on the Field-Road space $\bf{R}^{2}_{+} \quad | \quad t=$ ',string(TT(t)),' / ',string(Tmax),'$\quad | \quad \left( \mu , \nu \right)=$ $($',string(mu),',',string(nu),'$)\quad | \quad \left( d , D \right)=$ $($',string(d),',',string(D),'$)\quad | \quad \left\Vert\left( u\left(t,\cdot\right) , v\left(t,\cdot,\cdot\right) \right)\right\Vert_\infty=$ $($',string(maxU),' , ',string(maxV),'$)$'),'interpreter','latex','FontSize',18);
        mycolormapU=[linspace(0.3,0.9*maxU/maxU0,256);linspace(0.3,0.9*maxU/maxU0,256);0.7*ones(1,256)]';
        mycolormapV=[linspace(0.3,0.9*maxV/maxV0,256);linspace(0.3,0.9*maxV/maxV0,256);0.7*ones(1,256)]';
        %Plot Road u (1)
            ax(1)=subplot(3,2,1);
            plot(XX,U0,'LineWidth',2.0)                             %Show u
            hold on
            plot([0],[max(U0)],'o','LineWidth',2.0,'Color',[0.1 0.8 0])                 %Show the maximum position
            title('Curvy-plot of the Road: $u\left(t;\cdot\right)$','interpreter','latex','FontSize',16);    %Title
            aRoad=gca();                            %Get current axis for Road
            aRoad.YLim=[0 1.1];                     %Bounds for Y
            colormap(aRoad,mycolormapU);            %Color of the surface
            xlabel('$x$','interpreter','latex','FontSize',16);
            ylabel('$u\left( t ; x\right)$','interpreter','latex','FontSize',16);
            aRoad.YLabel.Rotation=45;              %Rotate the label
            aRoad.YLabel.Position=[-290.7237    0.5500   -1.0000];    %Shift the label on left
        %Plot Field v (2)
            ax(2)=subplot(3,2,2);
            surf(XX,YY,V0)                          %Show v
            hold on
            [MaxVindexX,MaxVindexY]=findMAX(V0);
            MaxVpositionX=XX(1,MaxVindexX);
            MaxVpositionY=YY(1,MaxVindexY);
            scatter3([MaxVpositionX]',[MaxVpositionY]',max(max(V0))*ones(length(MaxVpositionX),1),[100],[1 0 0],'filled','o','LineWidth',2.0)                 %Show the maximum position (parameters:X,Y,Z,size,color)
            title('Curvy-plot of the Field: $v\left(t; \cdot, \cdot\right)$','interpreter','latex','FontSize',16);    %Title
            aField=gca();                           %Get current axis for Field
            aField.View=[-32.4949   25.9708];       %Rotation around the surface
            aField.ZLim=[0 1];                      %Bounds for Z
            colormap(ax(2),mycolormapV);            %Color of the surface
            xlabel('$x$','interpreter','latex','FontSize',16);
            %aField.XLabel.Position=aField.XLabel.Position-[-50 -10 0];    %Shift the label
            ylabel('$y$','interpreter','latex','FontSize',16);
            %aField.YLabel.Position=aField.YLabel.Position;    %Shift the label
            zlabel('$v\left( t ; x,y\right)$','interpreter','latex','FontSize',16);
            aField.ZLabel.Rotation=45;              %Rotate the label
            aField.ZLabel.Position=[-380.0733  104.4665    0.5000];    %Shift the label on left
        %Plot Road u (3)
            ax(3)=subplot(3,2,3);
            pcolor(XX,[1 2],[U0;U0])                %Show u in 2D
            title('Colored-plot of the Road: $u\left(t; \cdot\right)$','interpreter','latex','FontSize',16);    %Title
            shading flat                            %To not show the mesh
            aRoad2=gca();                           %Get current axis for Road
            aRoad2.YTickLabel=[];
            colormap(ax(3),mycolormapU);            %Color of the pcolor
            colorbar('northoutside')
            xlabel('$x$','interpreter','latex','FontSize',16);
        %Plot Field v (4)
            ax(4)=subplot(3,2,4);
            pcolor(XX,YY,V0)                        %Show v in 2D
            title('Colored-plot of the Field: $v\left(t; \cdot, \cdot\right)$','interpreter','latex','FontSize',16);    %Title
            shading flat                            %To not show the mesh
            aField2=gca();                          %Get current axis for Field
            aField2.YLabel.Rotation=0;             %Rotate the label
            %aField2.YLabel.Position=aField2.YLabel.Position-[5 0 0];    %Shift the label on left
            colormap(ax(4),mycolormapU);            %Color of the pcolor
            colorbar
            xlabel('$x$','interpreter','latex','FontSize',16);
            ylabel('$y$','interpreter','latex','FontSize',16);
        %Plot Max u(t,.) in LOG/LOG (5)
            ax(5)=subplot(3,2,5);
            loglog(TLOG,UMAXLOG,'o-','LineWidth',2.0,'Color',[0.1 0.8 0]);
            hold on
            XXX=[300 1000];
            YYY=[1/3 0.1];
            regressionU=loglog(XXX,YYY,'-','LineWidth',3.0,'Color',[0.7 0.5 0]);
            legend(regressionU,'$y=-x$','interpreter','latex','FontSize',16);
            title('$\Vert u\left( t;\cdot \right)\Vert_\infty$ in LOG/LOG scale','interpreter','latex','FontSize',16);    %Title
            aLogU=gca();                            %Get current axis for MaxU in log scale
            aLogU.XLim=[1 max(TT)];   %Bounds for X
            aLogU.YLim=[0.001 2];   %Bounds for Y
            xlabel('$t$','interpreter','latex','FontSize',16);
            ylabel('$\Vert u\left( t;\cdot \right)\Vert_\infty$','interpreter','latex','FontSize',16);
            aLogU.YLabel.Rotation=45;              %Rotate the label
            aLogU.YLabel.Position=[0.37    0.0447   -1.0000];    %Shift the label on left
            grid on;
        %Plot Max v(t,.) in LOG/LOG (5)
            ax(6)=subplot(3,2,6);
            loglog(TLOG,VMAXLOG,'o-','LineWidth',2.0,'Color',[1 0 0]);
            hold on
            XXX=[300 1000];
            YYY=[1/3 0.1];
            regressionU=loglog(XXX,YYY,'-','LineWidth',3.0,'Color',[0.7 0.5 0]);
            legend(regressionU,'$y=-x$','interpreter','latex','FontSize',16);
            title('$\Vert v\left( t;\cdot ,\cdot \right)\Vert_\infty$ in LOG/LOG scale','interpreter','latex','FontSize',16);    %Title
            aLogV=gca();                            %Get current axis for MaxU in log scale
            aLogV.XLim=[1 max(TT)];   %Bounds for X
            aLogV.YLim=[0.001 2];   %Bounds for Y
            xlabel('$t$','interpreter','latex','FontSize',16);
            ylabel('$\Vert v\left( t;\cdot ,\cdot \right)\Vert_\infty$','interpreter','latex','FontSize',16);
            aLogV.YLabel.Rotation=45;              %Rotate the label
            aLogV.YLabel.Position=[0.35    0.0447   -1.0000];    %Shift the label on left
            grid on;
        drawnow;    %to have a continuous render!
        frame=getframe(gcf);             %Save the frame for movie!
        writeVideo(writerObj, frame);    %Save the frame for movie!
        hold on;
end
    %SOLVE THE PROBLEM
        W=B*W0;
    %SHIFT THE SOLUTION AND EXTRACT U AND V
        W0=W;               %Shift the solution
        U0=extractU(m,n,W0);    %Extract U
        V0=extractV(m,n,W0);    %Extract V
end

close(writerObj);

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%FUNCTIONS TO SWITCH BETWEEN MATRIX AND VECTOR
function k=mat2vect(n,i,j)    %IN: position in the matrix ; OUT: position in the vector
    k=n*(i-1)+j;
end

function [i,j]=vect2mat(n,k)    %IN: position in the vector ; OUT: position in the matrix
    i=floor((k-1)./n)+1;
    j=mod(k-1,n);
end

%TO BUILD INITIAL DATA
function z=indicator(x,a,b)
    if a<x && x<b
    z=1;
    else
    z=0;
    end
end

%COMBINE THE TWO INITIAL DATA INTO A VECTOR
function W=combineUV(m,n,U,V)
    M=[V ; U];   %Matrix which combines V and U
    W=zeros((m+1)*n);   %We convert M0 in the 
    for i=1:(m+1)
    for j=1:n
        W(mat2vect(n,i,j),1)=M(i,j);
    end
    end
end

%EXTRACT U AND V FROM THE COMBINATED VECTOR W
function U=extractU(m,n,W)
    U=zeros(1,n);
    for j=1:n
        U(1,j)=W(mat2vect(n,m+1,j));
    end
end

function V=extractV(m,n,W)
    V=zeros(m,n);
    for i=1:m
    for j=1:n
        V(i,j)=W(mat2vect(n,i,j));
    end
    end
end

%FIND THE INDEX OF THE MAX(s) IN A MATRIX
function [i,j]=findMAX(M)
    maxval=max(max(M));
    [j,i]=find(M == maxval);
end