function[]=Ass2_Question1_RK3_W5_Box()

clear all
close all

Res_min=3;      % Coarsest resolution
Res_max=10;     % Finest resolution
TF=4;           % Final time

L_Vector=Res_max-Res_min+1;     % Length of the vectors

% Defining the vectors that will be useful in the log-log-graph
Error_Inf=zeros(1,L_Vector);
Error_1=zeros(1,L_Vector);
Error_2=zeros(1,L_Vector);
Res=zeros(1,L_Vector);

c=0;
tic         % start of tic-toc

for k=Res_min:Res_max
    
    c=c+1;
    
    % Building the grid
    dx=2*pi/(2^k);
    X=0:dx:2*pi;
    N=length(X);
    
    dt=dx/3;    % CFL condition is satisfied
    M=ceil(TF/dt);
    
    % Initial u
    u=0*X;
    for i=1:N
        if X(i)<pi/4
            u(i)=1;
        elseif pi/4<X(i) && X(i)<pi/2
            u(i)=2;
        else
            u(i)=1;
        end
    end
    u=u';
    
    % RK3-TVD - WENO 5
    for t=dt:dt:M*dt
        space_u=WENO5(u,dx,N);
        psi_star=u-dt*space_u';
        space_u=WENO5(psi_star,dx,N);
        psi_star_star=(3/4)*u+(1/4)*psi_star-(1/4)*dt*space_u';
        space_u=WENO5(psi_star_star,dx,N);
        v=(1/3)*u+(2/3)*psi_star_star-2/3*dt*space_u';
        u=v;
    end
    
    % True solution at time T=TF
    u_true=0*X;
    for i=1:N
        if X(i)<pi/4+t
            u_true(i)=1;
        elseif pi/4+t<X(i) && X(i)<pi/2+t
            u_true(i)=2;
        else
            u_true(i)=1;
        end
    end
    
        % Storing the resolution and errors
        Res(c)=dx;
        Error_Inf(c)=max(abs(u-u_true'))
        Error_1(c)=dx*sum(abs(u-u_true'))
        Error_2(c)=(dx*sum(abs(u-u_true').^2))^(1/2)
end

%%% Plotting %%%
loglog(Res,Error_Inf,'+b',Res,Error_1,'+g',Res,Error_2,'+k')
hold on
% p_Inf=polyfit(log(Res),log(Error_Inf),1)
% loglog(Res,exp(p_Inf(2))*Res.^(p_Inf(1)),'-b')
% hold on
p_1=polyfit(log(Res(end-4:end)),log(Error_1(end-4:end)),1)
loglog(Res,exp(p_1(2))*Res.^(p_1(1)),'-g')
hold on
p_2=polyfit(log(Res(end-4:end)),log(Error_2(end-4:end)),1)
loglog(Res,exp(p_2(2))*Res.^(p_2(1)),'-k')
hold on
toc         % End of tic-toc


%%% Spatial discretization using WENO 5 %%%
function space_u=WENO5(u,dx,N)

space_u=zeros(1,N);

for i=1:N-1
    if i==1
        v=(1/dx)*[u(N-1)-u(N-2),u(i)-u(N-1),u(i+1)-u(i),u(i+2)-u(i+1),u(i+3)-u(i+2)];
    else if i==2
            v=(1/dx)*[u(i-1)-u(N-1),u(i)-u(i-1),u(i+1)-u(i),u(i+2)-u(i+1),u(i+3)-u(i+2)];
        else if i==N-2
                v=(1/dx)*[u(i-1)-u(i-2),u(i)-u(i-1),u(i+1)-u(i),u(1)-u(i+1),u(2)-u(1)];
            else if i==N-1
                    v=(1/dx)*[u(i-1)-u(i-2),u(i)-u(i-1),u(1)-u(i),u(2)-u(1),u(3)-u(2)];
                else
                    v=(1/dx)*[u(i-1)-u(i-2),u(i)-u(i-1),u(i+1)-u(i),u(i+2)-u(i+1),u(i+3)-u(i+2)];
                end
            end
        end
    end
    eps=10^(-6);
    s_1=(13/12)*(v(1)-2*v(2)+v(3))^2+(1/4)*(v(1)-4*v(2)+3*v(3))^2;
    s_2=(13/12)*(v(2)-2*v(3)+v(4))^2+(1/4)*(v(2)-v(4))^2;
    s_3=(13/12)*(v(3)-2*v(4)+v(5))^2+(1/4)*(3*v(3)-4*v(4)+v(5))^2;
    a_1=(1/10)/(eps+s_1)^2;
    a_2=(6/10)/(eps+s_2)^2;
    a_3=(3/10)/(eps+s_3)^2;
    s_a=a_1+a_2+a_3;
    w=(1/s_a)*[a_1,a_2,a_3];
    space_u(i+1)=(1/6)*(w(1)*(2*v(1)-7*v(2)+11*v(3))+w(2)*(-v(2)+5*v(3)+2*v(4))+w(3)*(2*v(3)+5*v(4)-v(5)));
end
space_u(1)=space_u(N);      % Periodic BC
