clear;
%Number of Full-domain Periods to Run
NT=1;
%Max domain size in space, and evolution duration time
Imin=0;
Imax=200;
Jmin=0;
Jmax=200;
%Spatial step size and temporal step duration
dx=0.5;
dy=0.5;
%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
CX=1; 
CY=1;

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


%Uncomment to set dt for CFL=
dCFL=1; %Desired CFL number
maxvels=sqrt(CX^2+CY^2);
dt=dCFL.*min(dx,dy)./(maxvels);



Nmin=0;
Nmax=Ydomain.*NT./(abs(CY).*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=A.*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.
    else        %Evolve solution at later timesteps
        %Initialize temporary U(n+1) from U(n)
        U1L=UL;
        ULa=UL;
        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(*) for Lax Method (x-direction)
                ULa(i,j)=(1./2).*(UL(i+1,j)+UL(i-1,j))-(CX.*UL(i+1,j)-CX.*UL(i-1,j)).*dt./(2.*dx);
            end;
        end;
        %Update boundaries for Lax-Wendroff intermediate
        ULa(Imin+1,:)=ULa(Imax+2,:);
        ULa(Imax+3,:)=ULa(2,:);     
        ULa(:,Jmin+1)=ULa(:,Jmax+2);
        ULa(:,Jmax+3)=ULa(:,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-f Method (y-direction)
                U1L(i,j)=(1./2).*(ULa(i,j+1)+ULa(i,j-1))-(CY.*ULa(i,j+1)-CY.*ULa(i,j-1)).*dt./(2.*dy);              
            end;
        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 
        UL=U1L;
        UL(Imin+1,:)=U1L(Imax+2,:);
        UL(Imax+3,:)=U1L(2,:);
        UL(:,Jmin+1)=U1L(:,Jmax+2);
        UL(:,Jmax+3)=U1L(:,2);
    end;
    %Show Results    
    figure(10);

    imagesc(x,y,UL');
    title(['Lax Friedrichs Solutions @ t=',num2str(n.*dt), ', n=',num2str(n),'']);
    axis([0 100 0 100]);
    axis square;
    caxis([-A A]); colorbar('horiz'); 
    xlabel('Horizontal Position');
    ylabel('Vertical Position');

end;
