%% Lax-Wendroff Two-Step Method 
% EP711 Example Code -- J.B. Snively

clear all;
close all;
clc;

%% Set Domain (Not all used, yet)

Xmin=-100000;
Xmax=100000;
Ymin=0;
Ymax=100000;
Xdomain=Xmax-Xmin;
Ydomain=Ymax-Ymin;

dx=1000;
dy=1000;

x_c=[Xmin-3*dx/2:dx:Xmax+3*dx/2];

y_c=[Ymin-3*dy/2:dy:Ymax+3*dy/2];

[X,Y]=meshgrid(x_c,y_c);

[J,I]=size(X);

F=cat(3,0.*X,0.*X,0.*X,0.*X);    %Fluxes for x-split
G=F;                        %Fluxes for y-split
Q=F;                        %Solution Variables
S=F;

%% Set Problem

%Initial Condition Parameters...

%Physical Parameters for Euler or Acoustics Equations
global g R P0 rho0 gamma C; 
global wind;
g=9.8;
R=287;
P0=1E5;
rho0=1.2;
gamma=1.4;
C=sqrt(gamma*P0/rho0);
H=P0/(rho0*g);
P0=P0.*exp(-Y./H);
rho0=rho0.*exp(-Y./H);
g=(P0(2:end,1)-P0(1:end-1,1))./(-0.5*dy*(rho0(2:end,1)+rho0(1:end-1,1)));
g=repmat(g,1,size(X,2)-1);
wind=100.*exp(-(Y-60000).^2./(2*5000^2));

%Time Parameters...

%CFL Constraints
dCFL=0.5; %Desired max CFL number (in either direction!)
dt=dCFL.*min(dx,dy)./(C);

%Number of Full-Domain Periods to Evolve (based on speed)
Nmin=0;
Nmax=500;
Tmin=dt.*Nmin;
Tmax=dt.*Nmax;

%Open Figure
fig1=figure(1);

%Wave Scale (box relative)
Mlength=3;
%Wave Amplitude
Amp=10000;
%Gaussian Envelope
Sigma_X=Xdomain/10;
Sigma_Y=Ydomain/10;

%Set Initial Condition (Gaussian Pressure Pulse)
%Pa=Amp.*exp(-.5*((X-Xmax./2).^2)./(Sigma_X).^2).*exp(-.5*((Y-Ymax./2).^2)./(Sigma_Y).^2);
Pa=Amp.*(((X).^2+(Y-10000).^2)<3000^2);
%Pa=0.*X;
Q(:,:,1)=rho0.*X.^0;
Q(:,:,2)=rho0.*wind;
Q(:,:,3)=0;
Q(:,:,4)=(Pa+P0.*X.^0)/(gamma-1)+0.5*rho0.*wind.^2;
%Set BCs for first time
Q=bc(Q);

% Set ranges for evaluations
iD=2:I-1;
jD=2:J-1;


%% Evolve Over Time

tic; % Measure Time

for n=Nmin:1:Nmax

    
    % Dimensional Splitting:
    
    
    %x-split
    F=Fflux(Q);
    Qs(jD,iD,:)=0.5*(Q(jD,iD,:)+Q(jD,iD+1,:))-(dt/(2*dx))*(F(jD,iD+1,:)-F(jD,iD,:));
    F=Fflux(Qs);
    Q(jD,iD,:) = Q(jD,iD,:)-(dt/dx)*(F(jD,iD,:)-F(jD,iD-1,:));
    %apply BCs
    Q=bc(Q);
    
    %y-split
    G=Gflux(Q);
    S=Source(0.5*(Q(jD,iD,:)+Q(jD+1,iD,:)),g(jD,iD));
    Qs(jD,iD,:)=0.5*(Q(jD,iD,:)+Q(jD+1,iD,:))-(dt/(2*dy))*(G(jD+1,iD,:)-G(jD,iD,:))+(dt/2)*S(:,:,:);
    G=Gflux(Qs);
    S=Source(Qs,g);
    Q(jD,iD,:) =Q(jD,iD,:)-(dt/dy)*(G(jD,iD,:)-G(jD-1,iD,:))+dt*0.5*(S(jD,iD,:)+S(jD-1,iD,:));
    %apply BCs
    Q=bc(Q);

        
    %store results
    Q_save(:,:,:,n+1)=Q;
    
end

toc; % Report Time

for n=Nmin:1:Nmax
    %plot results: Pressure
    subplot(2,1,1);
    KE=0.5*(Q_save(3:end-2,3:end-2,2,n+1).^2+Q_save(3:end-2,3:end-2,3,n+1).^2)./Q_save(3:end-2,3:end-2,1,n+1);
    PR=(Q_save(3:end-2,3:end-2,4,n+1)-KE).*(gamma-1)-P0(3:end-2,3:end-2);
    imagesc(x_c(3:end-2),y_c(3:end-2),PR); 
    axis xy; 
    colorbar; 
    xlabel('x'); ylabel('y'); title('Pressure Perturbation (Pa)');
    %plot results: Velocity
    subplot(2,1,2);
    imagesc(x_c(3:end-2),y_c(3:end-2),Q_save(3:end-2,3:end-2,2,n+1)./Q_save(3:end-2,3:end-2,1,n+1)); 
    axis xy; 
    colorbar; 
    xlabel('x'); ylabel('y'); title('Horizontal Velocity (m/s)');
    pause(0.001); 
end;


%% Functions

%Flux and Source Term Definitions


function F = Fflux(Q)
    global g R P0 rho0 gamma C; 
    P = (gamma-1)*(Q(:,:,4)-(0.5*(Q(:,:,2).^2+Q(:,:,3).^2)./Q(:,:,1)));
    F(:,:,1) = Q(:,:,2);
    F(:,:,2) = Q(:,:,2).*Q(:,:,2)./Q(:,:,1)+P;
    F(:,:,3) = Q(:,:,2).*Q(:,:,3)./Q(:,:,1);
    F(:,:,4) = (Q(:,:,4)+P).*Q(:,:,2)./Q(:,:,1);
end

function G = Gflux(Q)
    global g R P0 rho0 gamma C; 
    P = (gamma-1)*(Q(:,:,4)-(0.5*(Q(:,:,2).^2+Q(:,:,3).^2)./Q(:,:,1)));
    G(:,:,1) = Q(:,:,3);
    G(:,:,2) = Q(:,:,2).*Q(:,:,3)./Q(:,:,1);
    G(:,:,3) = Q(:,:,3).*Q(:,:,3)./Q(:,:,1)+P;
    G(:,:,4) = (Q(:,:,4)+P).*Q(:,:,3)./Q(:,:,1);
end

function S = Source(Q,g)
    S(:,:,1) = 0.*Q(:,:,1);
    S(:,:,2) = 0.*Q(:,:,1);
    S(:,:,3) = -Q(:,:,1).*g;
    S(:,:,4) = -Q(:,:,3).*g;
end

%Boundary Conditions

function Q = bc(Q)
    global g R P0 rho0 gamma C; 
    global wind;
    % Closed Bottom
    Q(1:2,:,1) = rho0(1:2,:)+(Q(3,:,1)-rho0(3,:)).*(rho0(1:2,:)./rho0(3,:)).^(0.5);
    Q(1:2,:,2) = rho0(1:2,:).*wind(1:2,:)+(Q(3,:,2)-rho0(3,:).*wind(3,:)).*(rho0(1:2,:)./rho0(3,:)).^(0.5);
    Q(1:2,:,3) = -Q(3,:,3).*(rho0(1:2,:)./rho0(3,:)).^(0.5);
    Q(1:2,:,4) = P0(1:2,:)./(gamma-1)+(Q(3,:,4)-P0(3,:)./(gamma-1)-0.5*rho0(3,:).*wind(3,:).^2).*(rho0(1:2,:)./rho0(3,:)).^(0.5)+0.5*rho0(1:2,:).*wind(1:2,:).^2;
    % Open Top
    Q(end-1:end,:,1) = rho0(end-1:end,:)+(Q(end-2,:,1)-rho0(end-2,:)).*(rho0(end-1:end,:)./rho0(end-2,:)).^(0.5);
    Q(end-1:end,:,2) = rho0(end-1:end,:).*wind(end-1:end,:)+(Q(end-2,:,2)-rho0(end-2,:).*wind(end-2,:)).*(rho0(end-1:end,:)./rho0(end-2,:)).^(0.5);
    Q(end-1:end,:,3) = Q(end-2,:,3).*(rho0(end-1:end,:)./rho0(end-2,:)).^(0.5);
    Q(end-1:end,:,4) = P0(end-1:end,:)./(gamma-1)+(Q(end-2,:,4)-P0(end-2,:)./(gamma-1)-0.5*rho0(end-2,:).*wind(end-2,:).^2).*(rho0(end-1:end,:)./rho0(end-2,:)).^(0.5)+0.5*rho0(end-1:end,:).*wind(end-1:end,:).^2;
    % Periodic Sides
    Q(:,2,:) = Q(:,end-2,:);
    Q(:,1,:) = Q(:,end-3,:);
    Q(:,end,:)  = Q(:,4,:);
    Q(:,end-1,:)= Q(:,3,:);   
    
end
