clear;
%Boundary Condition bc=0 for open, bc=1 for reflecting, bc=2 for periodic
bc=1;
%Number of wavelengths in the Box;
M=4;
%Nondimensional (maybe) advective speed C
C=1;  
%Number of Wave Periods to Run
NT=2;
%Wave Amplitude
A=1;
%Gaussian Envelope? (1 yes, 0 no)
G=1;
%V Envelope Multiplier (i.e., -1 for backward, 1 for forward)
B=1;
%Phase shift of second wave
phishift=0;

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

%Max domain size in space, and evolution duration time
Imin=0;
Imax=100;
Xmin=-(50+1); %Add one for BCs
Xmax=50+1;
Xdomain=(Imax-Imin+1).*dx;
Nmin=0;
Nmax=Xdomain.*NT./C;
Tmin=dt.*Nmin;
Tmax=dt.*Nmax;

%Domain X-axis Reference
x=[Xmin:dx:Xmax];

%Open Figure
figure(1);

%Calculate later times t...
for n=Nmin:1:Nmax   

    
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        %Define analytical solution, also used for initial value
        Ua=A.*cos(2*pi*(x./(Xdomain./M))-2*pi*((n.*dt)./(Xdomain./(C.*M))));
        Va=A.*cos(2*pi*(x./(Xdomain./M))-2*pi*((n.*dt)./(Xdomain./(C.*M)))-phishift);
        if G==1 
            Ua=Ua.*exp(-.5*((x./2).^2)./(Xdomain./20).^2);
            Va=B.*Va.*exp(-.5*((x./2).^2)./(Xdomain./20).^2);
        end;
        
        UL=Ua;  %Lax Init.
        VL=Va;
        ULW=Ua; %Lax-Wendroff Init.
        VLW=Va;
        
    else        %Evolve solution at later timesteps
        
        %Initialize temporary U(n+1) from U(n)
        U1L=UL;
        U1LWa=ULW;
        V1L=VL;
        V1LWa=VLW;
        
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1) for Lax Method
            U1L(i)=(1./2).*(UL(i+1)+UL(i-1))-(C.*VL(i+1)-C.*VL(i-1)).*dt./(2.*dx);
            V1L(i)=(1./2).*(VL(i+1)+VL(i-1))-(C.*UL(i+1)-C.*UL(i-1)).*dt./(2.*dx);  
        end;
        for im=Imin-1:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1/2) for Lax-Wendroff 2-Step Method
            U1LWa(i)=(1./2).*(ULW(i+1)+ULW(i))-(C.*VLW(i+1)-C.*VLW(i)).*dt./(2.*dx);       
            V1LWa(i)=(1./2).*(VLW(i+1)+VLW(i))-(C.*ULW(i+1)-C.*ULW(i)).*dt./(2.*dx);     
        end;

        %Perform second Lax-Wendroff step:
        %Update results for n+1/2 (Lax-Wendroff)
        ULWa=U1LWa; 
        VLWa=V1LWa; 
        %Initialize Array
        U1LW=ULWa;
        V1LW=VLWa;
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            %Obtain U(n+1) for Lax-Wendroff 2-Step Method
            U1LW(i)=ULW(i)-(C.*(VLWa(i)-VLWa(i-1))).*(dt./dx);   
            V1LW(i)=VLW(i)-(C.*(ULWa(i)-ULWa(i-1))).*(dt./dx);            
        end;
        %Set solution to updated values
        UL=U1L;
        VL=V1L;
        ULW=U1LW;
        VLW=V1LW;
        %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 
        UL(Imin+1)=U1L(Imin+2+(Imax-Imin)*(bc==2));
        UL(Imax+3)=U1L(Imax+2+(-Imax)*(bc==2));
        VL(Imin+1)=(1-2*(bc==1)).*V1L(Imin+2+(Imax-Imin)*(bc==2));
        VL(Imax+3)=(1-2*(bc==1)).*V1L(Imax+2+(-Imax)*(bc==2));
        %Update boundaries for Lax-Wendroff
        ULW(Imin+1)=ULW(Imin+2+(Imax-Imin)*(bc==2));
        ULW(Imax+3)=ULW(Imax+2+(-Imax)*(bc==2));
        VLW(Imin+1)=(1-2*(bc==1)).*VLW(Imin+2+(Imax-Imin)*(bc==2));
        VLW(Imax+3)=(1-2*(bc==1)).*VLW(Imax+2+(-Imax)*(bc==2));
    end;
    
    %Show Results    
    plot(x,Ua,'k:',x,ULW,'k--',x,UL,'k-',x,Va,'b:',x,VLW,'b--',x,VL,'b-');
    axis([-50 50 -2.*A 2.*A]);
    xlabel('Horizontal Position');
    ylabel('Amplitude');
    title(['Lax (Solid), Lax-Wendroff(Dashed) and Initial Condition (Dotted) @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    pause(0.005); 
    
end;
