clear;

%Number of wavelengths in the Box;
M=1;
%Number of Full-domain Periods to Run
NT=2;
%Wave Amplitude
ampl=1;

%Pressure
P0=1e5;
%Density
rho0=1;
%Gamma
gamma=1.4;

%Speed of sound
C=sqrt(gamma*P0/rho0);

%Constant Coefficient A-matrix
A=[[0,gamma*P0];[1/rho0,0]];
eigVec=[[-rho0*C,rho0*C];[1,1]];
eigVal=[[-C,0];[0,C]];

%The following aren't necessary for LeVeque's method! 
%Calculate them in case an alternate low-order flux is specified...
%|A|
absA=eigVec*abs(eigVal)*(eigVec^-1);
%A+
Aplus=eigVec*max(0,eigVal)*(eigVec^-1);
%A-
Aminus=eigVec*min(0,eigVal)*(eigVec^-1);  

%Spatial step size and temporal step duration
dx=.2833;
%Target CFL
CFL=0.5;
%dt=1, so that CFL=C (Comment for CFL=1)
dt=CFL*dx./C;
%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(dt*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=ampl.*sign(cos(2*pi*(x./(Xdomain./M))));
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        U=[[Ua];[zeros(size(Ua))]];%[Ua/(rho0*C)]];   %Solution Init.
        U1=U;   %Initialize temporary U(n+1) from U(n)
        %Initialize flux values for Low and High order fluxes
        FL=0.*U;  %Set to zeros
        FH=0.*U;
        theta=0.*U;
        alpha=theta;
        phi=theta;
        apdq=theta;
        amdq=theta;
        F=0.*U;
    else        %Evolve solution at later timesteps
        for im=Imin-1:1:Imax+1
            i=im+3;  %Convert to MATLAB Index
            deltaP=U(1,i)-U(1,i-1);
            deltaU=U(2,i)-U(2,i-1);
            %Calculate eigencoefficients from Riemann Problem
            alpha(1,i)=(-deltaP+rho0*C*deltaU)/(2*rho0*C);
            alpha(2,i)=(deltaP+rho0*C*deltaU)/(2*rho0*C);
        end;  
        for im=Imin-1:1:Imax+1
            i=im+3;  %Convert to MATLAB Index
            %Calculate slope ratios
            %Negative wave 1
            theta(1,i)=alpha(1,i+1)/(alpha(1,i)+1e-15);
            %Positive wave 2
            theta(2,i)=alpha(2,i-1)/(alpha(2,i)+1e-15);
        end;  

        %-----------------------------------------------------------------
        %Uncomment only one of the following phi assignments
        %
        %Linear Methods e.g., LeVeque [2002]
        %First-order upwind
        %phi=zeros(size(theta));
        %Lax-Wendroff
        %phi=ones(size(theta));
        %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));
        %-----------------------------------------------------------------
        
        %Apply limiters to characteristic coefficients
        alphalimit=alpha.*phi;    
        
        %Obtain flux differences / second order corrections
        for im=Imin-1:1:Imax+1
            i=im+3;  %Convert to MATLAB Index
            %Obtain Fluctuations for First-order Godunov Method
            %via Characteristics e.g. LeVeque
            %Negative wave contribution to pressure
            amdq(1,i)=eigVal(1,1)*(alpha(1,i)*eigVec(1,1));
            %Positive wave contribution to pressure
            apdq(1,i)=eigVal(2,2)*(alpha(2,i)*eigVec(1,2));
            %Negative wave contribution to velocity
            amdq(2,i)=eigVal(1,1)*(alpha(1,i)*eigVec(2,1));
            %Positive wave contribution to velocity
            apdq(2,i)=eigVal(2,2)*(alpha(2,i)*eigVec(2,2));
            %Via the A-matrices (This is equivalent)
            %apdq(:,i)=Aplus*(U(:,i)-U(:,i-1));
            %amdq(:,i)=Aminus*(U(:,i)-U(:,i-1));
            %Obtain High-order Fluxes e.g. LeVeque
            FH(1,i)=(1/2)*(abs(eigVal(1,1))*(1-(dt/dx)*abs(eigVal(1,1)))*alphalimit(1,i)*eigVec(1,1)+(abs(eigVal(2,2)))*(1-(dt/dx)*abs(eigVal(2,2)))*alphalimit(2,i)*eigVec(1,2));
            FH(2,i)=(1/2)*(abs(eigVal(1,1))*(1-(dt/dx)*abs(eigVal(1,1)))*alphalimit(1,i)*eigVec(2,1)+(abs(eigVal(2,2)))*(1-(dt/dx)*abs(eigVal(2,2)))*alphalimit(2,i)*eigVec(2,2));
        end;
        
        %Obtain total flux
        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)*(apdq(:,i)+amdq(:,i+1))-(dt/dx)*(FH(:,i+1)-FH(:,i));          
        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,U(1,:),'k-',x,U(2,:)*rho0*C/2,'k--');
    axis([-50*dx 50*dx -2 2]);
    xlabel('Horizontal Position');
    ylabel('Amplitude');
    title(['P (Solid), V*(Z_0/2) (Dashed), and Initial (Dotted) Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    pause(0.0005); 
    
end;