clear all;
%Number of Full-domain Periods to Run
NT=1;
%Max domain size in space, and evolution duration time
Imin=0;
Imax=100;
Jmin=0;
Jmax=100;
%Spatial step size and temporal step duration
dx=1;
dy=1;
%Max spatial extent and temporal duration
Xmin=dx.*(Imin-1);
Xmax=dx.*(Imax+1);
Xdomain=(Imax-Imin+1).*dx;
Ymin=dy.*(Jmin-1);
Ymax=dy.*(Jmax+1);
Ydomain=(Jmax-Jmin+1).*dy;

x=[Xmin:dx:Xmax];
y=[Ymin:dy:Ymax];
[X,Y]=meshgrid(x,y);

%Number of wavelengths in the Box;
M=3;
%Nondimensional (maybe) advective speed C
A=1; 
B=1;

%Wave Amplitude
C=1;
%Gaussian Envelope? (1 yes, 0 no)
G=1;

%Uncomment to set dt for CFL=
dCFL=1; %Desired CFL number
maxvels=sqrt(A^2+B^2);
dt=dCFL.*min(dx,dy)./(maxvels);
%Set Time Parameters
Nmin=0;
Nmax=Ydomain.*NT./(abs(B).*dt);
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=C.*cos(2*pi*(X./(Xdomain./M)));
    if G==1 
        Ua=Ua.*exp(-.5*((X-Xmax./2).^2)./(Xdomain./10).^2).*exp(-.5*((Y-Ymax./2).^2)./(Ydomain./10).^2);
    end;
    %Calculate numerical solution (or set to initial value)
    if n==Nmin  %Set initial value at initial timestep
        UL=Ua;      %Lax Init., Unsplit
        U1L=0.*Ua;
        ULS=Ua;      %Lax Init., Split
        ULSa=0.*Ua;
        U1LS=0.*Ua;
        ULW=Ua;     %Lax-Wendroff Init., Unsplit
        U1LW=0.*Ua;
        ULWS=Ua;   %Lax-Wendroff Init., Split
        ULWSa=0.*Ua;
        U1LWS=0.*Ua;

    else        %Evolve solution at later timesteps
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            for jm=Jmin:1:Jmax
                j=jm+2;  %Convert to MATLAB Index
                
                %Obtain U(n+1) for Lax Unsplit Method
                U1L(i,j)=(1./4).*(UL(i+1,j)+UL(i-1,j)+UL(i,j+1)+UL(i,j-1))-(A.*UL(i+1,j)-A.*UL(i-1,j)).*dt./(2.*dx)-(B.*UL(i,j+1)-B.*UL(i,j-1)).*dt./(2.*dy);
                
                %Obtain U(n+1) for Lax Split Method (x-direction)
                ULSa(i,j)=(1./2).*(ULS(i+1,j)+ULS(i-1,j))-(A.*ULS(i+1,j)-A.*ULS(i-1,j)).*dt./(2.*dx);
                
                %Obtain U(n+1) for Lax-Wendroff Split Method (x-direction)
                ULWSa(i,j)=ULWS(i,j)-((A*dt)/(2*dx)).*(ULWS(i+1,j)-ULWS(i-1,j))+(((A^2)*(dt^2))/(2*(dx^2))).*(ULWS(i-1,j)-2.*ULWS(i,j)+ULWS(i+1,j));                
                
                %Obtain U(n+1) for Lax-Wendroff Unsplit Method
                U1LW(i,j)=ULW(i,j)-((A*dt)/(2*dx)).*(ULW(i+1,j)-ULW(i-1,j))+(((A^2)*(dt^2))/(2*(dx^2))).*(ULW(i-1,j)-2.*ULW(i,j)+ULW(i+1,j));
                %Terms for y-direction
                U1LW(i,j)=U1LW(i,j)-((B*dt)/(2*dy)).*(ULW(i,j+1)-ULW(i,j-1))+(((B^2)*(dt^2))/(2*(dy^2))).*(ULW(i,j-1)-2.*ULW(i,j)+ULW(i,j+1));
                %Cross-terms from Taylor Series
                U1LW(i,j)=U1LW(i,j)+(((A*B+B*A)*(dt^2))/(8*(dx*dy))).*((ULW(i+1,j+1)-ULW(i-1,j+1))-(ULW(i+1,j-1)-ULW(i-1,j-1)));

            end;
        end;
        %Update boundaries for Lax intermediate
        ULSa(Imin+1,:)=ULSa(Imax+2,:);
        ULSa(Imax+3,:)=ULSa(2,:);     
        ULSa(:,Jmin+1)=ULSa(:,Jmax+2);
        ULSa(:,Jmax+3)=ULSa(:,2);
        %Update boundaries for Lax-Wendroff intermediate
        ULWSa(Imin+1,:)=ULWSa(Imax+2,:);
        ULWSa(Imax+3,:)=ULWSa(2,:);     
        ULWSa(:,Jmin+1)=ULWSa(:,Jmax+2);
        ULWSa(:,Jmax+3)=ULWSa(:,2);
        for im=Imin:1:Imax
            i=im+2;  %Convert to MATLAB Index
            for jm=Jmin:1:Jmax
                j=jm+2;  %Convert to MATLAB Index
                %Obtain U(n+1) for Lax Split Method
                U1LS(i,j)=(1./2).*(ULSa(i,j+1)+ULSa(i,j-1))-(B.*ULSa(i,j+1)-B.*ULSa(i,j-1)).*dt./(2.*dy);
                
                %Obtain U(n+1) for Lax-Wendroff Split Method (y-direction)
                U1LWS(i,j)=ULWSa(i,j)-((B*dt)/(2*dy)).*(ULWSa(i,j+1)-ULWSa(i,j-1))+(((B^2)*(dt^2))/(2*(dy^2))).*(ULWSa(i,j-1)-2.*ULWSa(i,j)+ULWSa(i,j+1));                
            end;
        end;
        %Initialize temporary U(n+1) from U(n)
        %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 (Unsplit)
        UL=U1L;
        UL(Imin+1,:)=U1L(Imax+2,:);
        UL(Imax+3,:)=U1L(2,:);
        UL(:,Jmin+1)=U1L(:,Jmax+2);
        UL(:,Jmax+3)=U1L(:,2);
        %Update boundaries for Lax (Split)
        ULS=U1LS;
        ULS(Imin+1,:)=U1LS(Imax+2,:);
        ULS(Imax+3,:)=U1LS(2,:);
        ULS(:,Jmin+1)=U1LS(:,Jmax+2);
        ULS(:,Jmax+3)=U1LS(:,2);
        %Update boundaries for Lax-Wendroff (Unsplit)
        ULW=U1LW;
        ULW(Imin+1,:)=U1LW(Imax+2,:);
        ULW(Imax+3,:)=U1LW(2,:);     
        ULW(:,Jmin+1)=U1LW(:,Jmax+2);
        ULW(:,Jmax+3)=U1LW(:,2);
        %Update boundaries for Lax-Wendroff (split)
        ULWS=U1LWS;
        ULWS(Imin+1,:)=U1LWS(Imax+2,:);
        ULWS(Imax+3,:)=U1LWS(2,:);     
        ULWS(:,Jmin+1)=U1LWS(:,Jmax+2);
        ULWS(:,Jmax+3)=U1LWS(:,2);
    end;
    %Show Results    
    %figure(10);

    figure(3);
    title(['Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    subplot(2,2,1);
    imagesc(x,y,UL');
    title('2D Unsplit Lax');
    axis([0 100 0 100]);
    axis square;
    caxis([-A A]); colorbar('horiz'); 
    xlabel('Horizontal Position');
    ylabel('Vertical Position');
    
    subplot(2,2,2);
    imagesc(x,y,ULS');
    title('2D Split Lax');
    axis([0 100 0 100]);
    axis square;
    caxis([-A A]); colorbar('horiz'); 
    xlabel('Horizontal Position');
    ylabel('Vertical Position');
    
    subplot(2,2,3);    
    imagesc(x,y,ULW');
    title('2D Unsplit Lax-Wendroff 1');
    axis([0 100 0 100]);
    axis square;
    caxis([-A A]); colorbar('horiz');
    xlabel('Horizontal Position');
    ylabel('Vertical Position');
    
    subplot(2,2,4);    
    imagesc(x,y,ULWS');
    title('2D Split Lax-Wendroff 1');
    axis([0 100 0 100]);
    axis square;
    caxis([-A A]); colorbar('horiz');
    xlabel('Horizontal Position');
    ylabel('Vertical Position');
    
    pause(0.001);
    
end;
