clear;
clear all;
clc;
%constants
W=2*pi; %width of mesh grid
L=4*pi; %length of mesh grid
t_final=5; %max time
dx=2*pi/100; % point seperation in x-axis
dy=dx; %point seperation in y-axis
dt=0.01; %point seperation in time-axis
alpha=0.1; %thermal diffusivity constant in x-direction
beta=alpha; %thermal diffusivity constant in y-direction
n=round((L/dx)+1); %number of points in x-asix
m=round((W/dy)+1); %number of points in y-axis
tt=(t_final/dt)+1  %number of points in z-axis
%%
x=zeros(n,1);   %x vector
for i=1:n
    x(i)=(i-1)*dx;
end
y=zeros(m,1);   %y vector
for j=1:m
    y(j)=(j-1)*dy;
end
t=zeros(tt,1);    %t vector
for k=1:tt
    t(k)=(k-1)*dt;
end
%%
T=zeros(n,m,tt); % realisations
for k=1;
    for j=1:m
        for i = 1:n
        T(i,j,k)=20*cos(x(i))*sin(y(j));
        end
    end
end


%Setting up velocity
v=zeros(n,tt);
for k=1:tt
for i=1:n
    v(i,k)=5*(1-(((x(i))-2*pi)^2/(4*pi^2)))*cos(pi*t(k))*sin(x(i));
end
end

as=zeros(n,tt);
an=zeros(n,tt);
for k =1:tt
    for i=1:n
        as(i,k)=(-v(i,k)*dt/(2*dy))-(beta*dt/(dy^2));
        an(i,k)=(v(i,k)*dt/(2*dy))-(beta*dt/(dy^2));
    end
end

aw=(-alpha*dt/(dx^2));
ap=1+(2*alpha*dt/(dx^2))+(2*beta*dt/(dy^2));
ae=-alpha*dt/(dx^2);

b=zeros(n*m,tt); %Setting up matrix b
k=1;
for j =1:m
    for i=1:n
        P=(j-1)*n+i;
        b(P,k)=20*cos(x(i))*sin(y(j));
    end
end

Temp(:,k)=b(:,k);

A=zeros(n*m);  %setting up matrix A

%matrix A
for k=2:tt
%T1
for i=2:n-1
    for j=2:m-1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P+1)=ae;
        A(P,P-1)=aw;
        A(P,P+n)=an(i,k);
        A(P,P-n)=as(i,k);
    end
end
%T2
for i=1
    for j=2:m-1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P+1)=ae+aw;
        A(P,P-n)=as(i,k);
        A(P,P+n)=an(i,k);
    end
end

%T3
for i=n
    for j=2:m-1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P-1)=ae+aw;
        A(P,P-n)=as(i,k);
        A(P,P+n)=an(i,k);
    end
end

%T4
for j=1
    for i=2:n-1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P+1)=ae;
        A(P,P+(m-2)*n)=as(i,k);
        A(P,P+n)=an(i,k);
        A(P,P-1)=aw;
    end
end

%T5
for j=m
    for i=2:n-1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P+1)=ae;
        A(P,P-n)=as(i,k);
        A(P,P-(m-2)*n)=an(i,k);
        A(P,P-1)=aw;
    end
end

%T6
for i=1
    for j=1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P+1)=aw+ae;
        A(P,P+n)=an(i,k);
        A(P,P+(m-2)*n)=as(i,k);
    end
end

%T7
for i=n
    for j=1
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P-1)=aw+ae;
        A(P,P+n)=an(i,k);
        A(P,P+(m-2)*n)=as(i,k);
    end
end

%T8
for i=1
    for j=m
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P+1)=aw+ae;
        A(P,P-n)=as(i,k);
        A(P,P-(m-2)*n)=an(i,k);
    end
end

%T9
for i=n
    for j=m
        P=(j-1)*n+i;
        A(P,P)=ap;
        A(P,P-1)=aw+ae;
        A(P,P-n)=as(i,k);
        A(P,P-(m-2)*n)=an(i,k);
    end
end
A=sparse(A);
Temp(:,k)=A\b(:,k-1);
b(:,k)=Temp(:,k);
end
%% Graphs

%t=0
figure(1)
    T1_init=reshape(b(:,t==0),[n,m]);
    T1=flipud(T1_init');
    contourf(x,y,T1);
    c=colorbar;
    c.Label.String = ['Temperature in ' char(176) 'C'];
    title('Temperature distribution at 0 seconds');
    xlabel('Length (m)');
    ylabel('Width (m)');

%t=0.1
figure(2)
    T2_init=reshape(b(:,t==0.1),[n,m]);
    T2=flipud(T2_init');
    contourf(x,y,T2);
    c=colorbar;
    c.Label.String = ['Temperature in ' char(176) 'C'];
    title('Temperature distribution at 0.1 seconds');
    xlabel('Length (m)');
    ylabel('Width (m)');

%t=0.5
figure(3)
    T3_init=reshape(b(:,t==0.5),[n,m]);
    T3=flipud(T3_init');
    contourf(x,y,T3);
    c=colorbar;
    c.Label.String = ['Temperature in ' char(176) 'C'];
    title('Temperature distribution at 0.5 seconds');
    xlabel('Length (m)');
    ylabel('Width (m)');

%t=1
figure(4)
    T4_init=reshape(b(:,t==1),[n,m]);
    T4=flipud(T4_init');
    contourf(x,y,T4);
    c=colorbar;
    c.Label.String = ['Temperature in ' char(176) 'C'];
    title('Temperature distribution at 1 second');
    xlabel('Length (m)');
    ylabel('Width (m)');

%t=2
figure(5)
    T5_init=reshape(b(:,t==2),[n,m]);
    T5=flipud(T5_init');
    contourf(x,y,T5);
    c=colorbar;
    c.Label.String = ['Temperature in ' char(176) 'C'];
    title('Temperature distribution at 2 seconds');
    xlabel('Length (m)');
    ylabel('Width (m)');

%t=5
figure(6)
    T6_init=reshape(b(:,t==5),[n,m]);
    T6=flipud(T6_init');
    contourf(x,y,T6);
    c=colorbar;
    c.Label.String = ['Temperature in ' char(176) 'C'];
    title('Temperature distribution at =5 seconds');
    xlabel('Length (m)');
    ylabel('Width (m)');