function Wave_Equation
% Wave_Equation
%   Q_t + [F(Q)]_x = 0
%   Q(0,t) = Q(2,t);
%   Q(x,0) = g(x);  i.e. g(x) = sin(pi*x)
%   x = [-1,1]
%   Final Time = 2 (Problem = 0), 0.45 (Problem = 1)
% ====================================================
% Input
Problem          = 1      ;  % Problem_Type : 0-Wave ; 1-Burgers
Amplitude        = 1      ;

x0               = 0      ;
x1               = 2      ;
N                = 81     ;

Final_Time       = 0.45   ;
CFL              = 0.45   ;

Movie_Resolution = 0.1   ;
% ====================================================
% WENO Setup
WENO_Order       = 5     ;
WENO_Flag        = 1     ;
WENO_Epsilon     = 1e-12 ;
WENO_Power       = 2     ;
% ====================================================
% Index Setup
r  = (WENO_Order + 1)/2  ;
N0 = 1      	         ;
N2 = N0 + r              ;
N1 = N2 - 1              ;
N3 = N2 + (N-1) - 1      ;
N4 = N3 + 1              ;
N5 = N3 + r              ;
% ====================================================
[ x , dx ] =   Domain_Setup ( N      ) ;

Q =  Initial_Condition ( x, 0 ) ;
Q = Boundary_Condition ( Q    ) ;

Step = 0 ; Time = 0 ; dt = Time_Step ( Q ) ;

while ( Time < Final_Time )
    dt   = Time_Step ( Q )   ;
    Step = Step + 1          ;
    Time = Time + dt         ;
    
    if Time > Final_Time, dt = dt - (Time - Final_Time); Time = Final_Time; end
    
    Q = Runge_Kutta ( Q )    ;
    
    Output( x, Q )           ;
end

if Problem == 0
    Q_E = Exact_Solution( x, Time )   ;
    Err = norm(Q-Q_E, inf); disp(Err) ;
end

fprintf('Done ... \n');

% ================================================================
    function [x, dx] = Domain_Setup ( N )
        
        dx = ( x1 - x0 ) / (N-1)  ;
        x  = x0 + ((N0:N5)'-N2) * dx + 1/2 * dx ;
        
    end
% ================================================================
    function Q = Runge_Kutta ( Q )
        
        %  Stage 1
        t = Time                        ;
        
        D_F = D_Fluxes( Q )             ;
        Q1  = Q + dt * D_F              ;
        
        Q1  = Boundary_Condition ( Q1 )	;
        
        % Stage 2
        t = Time - dt/2                 ;
        
        D_F = D_Fluxes( Q1 )            ;
        Q1  = (3.0*Q + Q1 + dt*D_F)/4.0 ;
        
        Q1  = Boundary_Condition ( Q1 ) ;
        
        % Stage 3
        t = Time + dt                   ;
        
        D_F = D_Fluxes( Q1 )            ;
        Q   = (Q + 2*Q1 + 2*dt*D_F)/3.0 ;
        
        Q   = Boundary_Condition ( Q  )	;
    end
% ================================================================
    function dt = Time_Step ( Q )
        if Problem == 0, dt = CFL*dx                    ; end
        if Problem == 1, dt = CFL*dx/(max(abs(Q(:,1)))) ; end
    end
% ================================================================
    function D_F = D_Fluxes ( Q )
        D_F = zeros(N5-N0+1,1);
        
        if Problem == 0, F =     Q    ; end
        if Problem == 1, F = 1/2*Q.^2 ; end
        
        [Flux_P, Flux_M] = Flux_Spliting(F, Q);

        F_Half = WENO_Flux ( Flux_P, Flux_M );
        
        D_F(N2:N3) = -(F_Half(N2:N3) - F_Half(N2-1:N3-1))/dx;
        
    end
% ================================================================
    function F_Half = WENO_Flux ( Flux_P, Flux_M )
        F_P = zeros(N5-N0+1,1); F_M = zeros(N5-N0+1,1);
        
        i0 = -2; i1 = 3;
        for i = N2-1:N3
            F_P(i) = WENO5_P( Flux_P(i+i0:i+i1) );
            F_M(i) = WENO5_M( Flux_M(i+i0:i+i1) );
        end
        F_Half = F_P + F_M;
 
    end
% ================================================================
    function [Flux_P, Flux_M] = Flux_Spliting(F, Q)
        if Problem == 0, Lambda =     1       ; end
        if Problem == 1, Lambda = max(abs(Q)) ; end
        
        Flux_P = 1/2 * (F + Lambda * Q);
        Flux_M = 1/2 * (F - Lambda * Q);
    end
% ================================================================
    function fh = WENO5_P( f )
        r    = (WENO_Order + 1)/2   ;
        ep   = WENO_Epsilon         ;
         p   = WENO_Power           ;
        flag = WENO_Flag            ;
        
        d  = [1/10, 6/10, 3/10];
        C1 = 1/4; C2 = 13/12;
        
        f  = f(:);
        
        P    = zeros(r,1);
        beta = zeros(1,r);
        
        P(1) = [ 2 -7 11]/6 * f(r-2:r  );
        P(2) = [-1  5  2]/6 * f(r-1:r+1);
        P(3) = [ 2  5 -1]/6 * f(r  :r+2);
        
        beta(1) = C1 * ([ 1 -4  3] * f(r-2:r  )).^2 + C2 * ([1 -2 1] * f(r-2:r  )).^2;
        beta(2) = C1 * ([-1  0  1] * f(r-1:r+1)).^2 + C2 * ([1 -2 1] * f(r-1:r+1)).^2;
        beta(3) = C1 * ([-3  4 -1] * f(r  :r+2)).^2 + C2 * ([1 -2 1] * f(r  :r+2)).^2;
        
        if flag == 1
            tau = abs(beta(3) - beta(1));
            alpha = d .* (1 + (tau ./ (ep+beta)).^p);
        else
            alpha = d ./ (ep+beta).^p;
        end
        
        omega = alpha / sum(alpha);
        fh    = omega * P;
    end
% ================================================================
    function fh = WENO5_M( f )
        f = f(end:-1:1);
        fh = WENO5_P( f );
        %{
        r    = (WENO_Order + 1)/2   ;
        ep   = WENO_Epsilon         ;
         p   = WENO_Power           ;
        flag = WENO_Flag            ;

        d  = [1/10, 6/10, 3/10];
        C1 = 1/4; C2 = 13/12;
        
        f  = f(:);
        
        P    = zeros(r,1);
        beta = zeros(1,r);
        
        P(1) = [11 -7  2]/6 * f(r+1:r+3);
        P(2) = [ 2  5 -1]/6 * f(r  :r+2);
        P(3) = [-1  5  2]/6 * f(r-1:r+1);
        
        beta(1) = C1 * ([ 3 -4  1] * f(r+1:r+3)).^2 + C2 * ([1 -2 1] * f(r+1:r+3)).^2;
        beta(2) = C1 * ([ 1  0 -1] * f(r  :r+2)).^2 + C2 * ([1 -2 1] * f(r  :r+2)).^2;
        beta(3) = C1 * ([-1  4 -3] * f(r-1:r+1)).^2 + C2 * ([1 -2 1] * f(r-1:r+1)).^2;
        
        if flag == 1
            tau = abs(beta(3) - beta(1));
            alpha = d .* (1 + (tau ./ (ep+beta)).^p);
        else
            alpha = d ./ (ep+beta).^p;
        end
        
        omega = alpha / sum(alpha);
        fh    = omega * P;
        %}
    end
% ================================================================
    function f = Exact_Solution( x, Time )
        f =  Amplitude*sin(pi*(x-Time)) ;
    end
% ================================================================
    function f = Initial_Condition( x, Time )
        f = Exact_Solution (x, Time) ;
    end
% ================================================================
    function f = Boundary_Condition( f )
        L1 = N1 - N0; L2 = N5 - N4 ;
        f(N0:N1) = f(N3-L1:N3) ;
        f(N4:N5) = f(N2:N2+L2) ;
    end
% ================================================================
    function Output( x, Q )
        fprintf('Plot ... \n') ; plot(x, Q, 'o-'); grid on
        title(['t = ', num2str(Time, '%1.2f')]);
        axis([x0, x1, -(Amplitude+0.5), Amplitude+0.5]);
        pause(Movie_Resolution)
    end
% ================================================================
%     function Stability_Check
%         if dt < 1.0e-16, error('=========  STOP ! Unstable  =============='); end
%     end
% ================================================================
end