% Example 13.4 Convection Equation
% Using 2nd Order Upwind Method

clear all
close all
clc

% Physical parameters
u = 0.10;
L = 1;
Lmax = 5;
tmax = 10;
T = 200;
x0 = 1;

% Mesh
dx = 0.05;
dt = 0.05;
c = u*dt/dx;

imax = Lmax/dx+1;
nmax = floor(tmax/dt)+1;
f = zeros(nmax,imax); % allocate solution matrix
x = [0:dx:Lmax];
t = [0:dt:tmax];

% initial condition
f(1,:) = initialwave(x0,L,T,x);


% Analytical Solution
xplot = [0:0.001:Lmax];
fexact(1,:) = initialwave(x0,L,T,xplot);
for n = 2:nmax-1
    fexact(n,:) = initialwave(x0+u*(n-1)*dt,L,T,xplot);
end

% 2nd Order Upwind Method
for n = 1:nmax-1
    for i = 3:imax-1
        f(n+1,i) = f(n,i) - c.*(f(n,i)-f(n,i-1)) - (c*(1-c)/2).*(f(n,i)-2.*f(n,i-1)+f(n,i-2));
    end
end


figure
ylim([-20,100])
xlabel('x')
ylabel('T ^oC')
title(['Convection Equation - 2nd Order Upwind Method (c = ',num2str(c),')'])
hold on
plot(xplot,fexact(1,:),'-','LineWidth',2)
plot(xplot,fexact(nmax-1,:),'-','LineWidth',2)
for j = 1:3:nmax-1
    if j ==1
        plot(x,f(j,:),'o-k')
    else
        plot(x,f(j,:),'ok')
    end
    %pause(dt)
    %saveas(gcf,[num2str(j),'.png'])
end
plot(x,f(nmax-1,:),'o-k')

figure
ylim([-20,100])
xlabel('x')
ylabel('T ^oC')
title(['Convection Equation - 2nd Order Upwind Method (c = ',num2str(c),')'])
%hold on
%plot(xplot,fexact(1,:),'-','LineWidth',2)
%plot(xplot,fexact(nmax-1,:),'-','LineWidth',2)
for j = 1:nmax-1
    clf
    plot(x,f(j,:),'ok')
    ylim([-20,100])
    xlim([0,5])
    xlabel('x')
    ylabel('T ^oC')
    title(['Convection Equation - 2nd Order Upwind Method (c = ',num2str(c),')'])
    pause(dt)
    %saveas(gcf,[num2str(j),'.png'])
end

figure
ylim([0,100])
xlabel('x')
ylabel('T ^oC')
title(['Convection Equation - Exact Solution'])
hold on
for j = 1:2:nmax-1
    plot(xplot,fexact(j,:),'-')
    %pause(dt)
    %saveas(gcf,[num2str(j),'.png'])
end
