clear;

%Number of wavelengths in the Box;
M=1;
%Nondimensional (maybe) advective speed C
C=.5;  
%Number of Full-domain Periods to Run
NT=2;
%Wave Amplitude
A=1;

%Spatial step size and temporal step duration
dx=1;
%dt=1, so that CFL=C (Comment for CFL=1)
dt=1;
%Calculate and show CFL number
CFL=C.*dt./dx

%Max domain size in space, and evolution duration time
Imin=0;
Imax=100;
Xmin=-dx.*(Imax+4)/2;
Xmax=dx.*(Imax+4)/2;
Xdomain=(Xmax-Xmin)-dx.*4;

%Spatial domain position data
x=[Xmin:dx:Xmax];

%Max spatial extent and temporal duration
Nmin=0;
Nmax=Xdomain.*NT./abs(C);
Tmin=dt.*Nmin;
Tmax=dt.*Nmax;

%Open Figure
figure(1);

%Calculate later times t...
for n=Nmin:1:Nmax   
    %Define analytical solution, also used for initial value
    Ua=A.*sign(cos(2*pi*(x./(Xdomain./M))-2*pi*((n.*dt)./(Xdomain./(C.*M)))));
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        U=Ua;   %Solution Init.
        U1=U;   %Initialize temporary U(n+1) from U(n)
        %Initialize flux values for Low and High order fluxes
        FL=0.*Ua;  %Set to zeros
        FH=0.*Ua;
        theta=0.*Ua;
        F=0.*Ua;
    else        %Evolve solution at later timesteps
        for im=Imin-1:1:Imax
            i=im+3;  %Convert to MATLAB Index
            %Obtain low-order flux from upwind, high-order from LW 
            FL(i)=(C/2)*(U(i)+U(i+1))-(abs(C)/2)*(U(i+1)-U(i));
            FH(i)=(C/2)*(U(i)+U(i+1))-(C^2*dt/(2*dx))*(U(i+1)-U(i));  
            %Set ratio of slopes
            theta(i)=(C>0)*(U(i)-U(i-1))/(U(i+1)-U(i)+1e-15)+(C<0)*((U(i+2)-U(i+1))/(U(i+1)-U(i)+1e-15));
        end;  
        
        %Uncomment only one of the following phi assignments
        
        %Linear Methods e.g., LeVeque [2002]
        %First-order upwind
        %phi=0;
        %Lax-Wendroff
        %phi=1;
        %Beam-Warming
        phi=theta;
        %Fromm
        %phi=0.5*(1+theta);
        
        %Highres Flux-Limiters e.g., LeVeque [2002]
        %Min-mod limiter
        %phi=max(0,min(1,theta));
        %Superbee limiter
        %phi=max(0,max(min(1,2*theta),min(2,theta)));
        %Van Leer
        %phi=(theta+abs(theta))./(1+abs(theta));
        %Monotized Center
        %phi=max(0,min(min((1+theta)/2,2),2*theta));
        
        %Obtain total flux
        F=FL+phi.*(FH-FL);
        for im=Imin-1:1:Imax+1
            i=im+3;  %Convert to MATLAB Index
            %Obtain U(n+1) for Flux-Limited Method
            U1(i)=U(i)-(dt/dx)*(F(i)-F(i-1));          
        end;
        
        %Assign value at right boundary, such that domain is "periodic"
        %This means, physically, that whatever goes out on the right will 
        %come back in on the left.
        %Update boundaries for Lax 
        U=U1;
        U(Imin+2)=U1(Imax+2);
        U(Imin+1)=U1(Imax+1);
        U(Imax+3)=U1(Imin+3);
        U(Imax+4)=U1(Imin+4);
end;
    %Show Results    
    plot(x,Ua,'k:',x,U1,'k-');
    axis([-50 50 -2 2]);
    xlabel('Horizontal Position');
    ylabel('Amplitude');
    title(['Numerical (Solid) and Exact (Dotted) Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    pause(0.0005); 
    
end;
