% Solving the Euler Equations
%    d      d
%    -- U = -- F(U)
%    dt     dx
%
% using Lax-Friedrichs 
% 
C1='Solving the Euler Equations';
C2='    d      d               ';
C3='    -- U = -- F(U)         ';
C4='    dt     dx              ';
C5='                           ';
C6=' using Lax-Friedrichs      ';
C=[C5;C1;C2;C3;C4;C5;C6;C5];
disp(C);
clear all
clf

M=input('Give number of grid-points (M): ');
dx=20/(M+1);

x=[0:dx:20];
disp('The integration is from 0 to 2, 2/dt should be an integer')
dt=input('Give time-step, dt: ');
N=ceil(2/dt);

gamma=1.4;

MM=M+2;
U=zeros(3,MM);
U_LEFT=zeros(3,1);
U_RIGHT=zeros(3,1);

U(1,1:MM/2)=1.0;
U(2,1:MM/2)=0.0;
U(3,1:MM/2)=1.0/(gamma-1.0);

U(1,MM/2+1:MM)=0.125;
U(2,MM/2+1:MM)=0.0;
U(3,MM/2+1:MM)=0.1/(gamma-1.0);

U_LEFT(1)=1.0;
U_LEFT(2)=0.0;
U_LEFT(3)=1.0/(gamma-1.0);

U_RIGHT(1)=0.125;
U_RIGHT(2)=0.0;
U_RIGHT(3)=0.1/(gamma-1.0);

for n=1:N
  UNY(:,2:MM-1)=0.5*(U(:,3:MM)+U(:,1:MM-2))-...
                 dt/(2*dx)*(F(U(:,3:MM),gamma)-...
                    F(U(:,1:MM-2),gamma));
  UNY(:,1)=0.5*(U(:,2)+U_LEFT)-...
                 dt/(2*dx)*(F(U(:,2),gamma)-F(U_LEFT,gamma));
  UNY(:,MM)=0.5*(U_RIGHT+U(:,MM-1))-...
                  dt/(2*dx)*(F(U_RIGHT,gamma)-...
                        F(U(:,MM-1),gamma));
  U=UNY;
  t(n)=n*dt;
  rho(n,:)=U(1,:);
  u(n,:)=U(2,:)./U(1,:);
end
tf=N*dt;
TF=num2str(tf);
[T X]=meshgrid(t,x);
figure(1)
surf(X,T,rho')
view(2)
shading interp
drawnow
xlabel('x')
ylabel('t')
title('\rho')

figure(2)
surf(X,T,u')
view(2)
shading interp
drawnow
xlabel('x')
ylabel('t')
title('u')

figure(3)
plot(x,rho(N,:))
xlabel('x')
title(['\rho at t=',TF])

figure(4)
plot(x,u(N,:))
xlabel('x')
title(['u at t=',TF])