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
U0=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.*sqrt(2));
%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=U0.*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
        %Initialize main U(n), U(*), and temporary U(n+1)
        ULW=Ua;
        U1LW=Ua;
        ULWa=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-Wendroff Method 
                %Terms for x-direction
                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;
        %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-Wendroff
        ULW=U1LW;
        ULW(Imin+1,:)=U1LW(Imax+2,:);
        ULW(Imax+3,:)=U1LW(2,:);     
        ULW(:,Jmin+1)=U1LW(:,Jmax+2);
        ULW(:,Jmax+3)=U1LW(:,2);
    end;

    %Show Results    
    figure(1);
    title(['Unsplit Lax-Wendroff Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    imagesc(x,y,ULW');
    axis([0 100 0 100]);
    axis square;
    caxis([-A A]); colorbar('horiz');
    xlabel('Horizontal Position');
    ylabel('Vertical Position');
    
end;
