% script file euler_sys.m
% solve the system of Euler equations in primitive
% and conservative variables
% with Lax-Friedrich's scheme and compare the results
% last changes: 16.06.03 by Vita
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear all; 

ax = -100; bx = 100;
n = 81; 
x = linspace(ax,bx,n);
h = x(2)-x(1);

tau = 0.1;
Tmax = 20;

% initial data
rho = ones(n,1);
u   = zeros(n,1);
p   = zeros(n,1);

lambda = tau/h;

left = find(x<0);
right = find(x>=0);

% init.1
%rho(left) = 5;

% init.2
%p(right) = 5;

% init.3
rho(left) = 5; u(left) = 5;

gamma = 7/5; % for a diatomic gas
E = p/(gamma-1) + (rho.*u.^2)/2;

% initial data for conservative variables
q1 = rho;
q2 = rho.*u;
q3 = E;


t = 0;

figure('name','Incompressible 1D Euler system')

while t<Tmax
  t = t+tau;
  
  E = p/(gamma-1) + (rho.*u.^2)/2;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % UPDATE IN PRIMITIVE VARIABLES  
  % does not make sense, as LxF scheme by it's construction
  % approximates equations in conservation form; and equation
  % in the primitive variables does not have it!
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % UPDATE IN CONSERVATIVE VARIABLES
  q1_temp = [q1(1); q1; q1(n)];
  q2_temp = [q2(1); q2; q2(n)];  
  q3_temp = [q3(1); q3; q3(n)];  
  
  % compute the f
  f1 = q2_temp;
  f2 = q2_temp.^2./q1_temp + (gamma-1).*(q3_temp-q2_temp.^2./q1_temp/2);
  f3 = q2_temp./q1_temp.*(gamma*q3_temp - (gamma-1)/2*q2_temp.^2./q1_temp);
  
  q1_new = (q1_temp(3:n+2)+q1_temp(1:n))/2 - ...
	   lambda/2*(f1(3:n+2)-f1(1:n));
  
  q2_new = (q2_temp(3:n+2)+q2_temp(1:n))/2 - ...
	   lambda/2*(f2(3:n+2)-f2(1:n));
  
  q3_new = (q3_temp(3:n+2)+q3_temp(1:n))/2 - ...
	   lambda/2*(f3(3:n+2)-f3(1:n));  

  q1 = q1_new;
  q2 = q2_new;
  q3 = q3_new;
  
  % calculate the corresponding rho, u, p:
  rho = q1;
  u = q2./q1;
  p = (gamma-1)*(q3-q2.^2./q1/2);
  
  
  subplot(3,1,1);  plot(x,rho,'b'); title('\rho');
  legend(['t=',num2str(t)],0)
  subplot(3,1,2);  plot(x,u,'b'); title('u');
  subplot(3,1,3);  plot(x,p,'b'); title('p')
  xlabel('x')
  
  drawnow;

  %pause
  
end

figure(2)
hold on;
P=plot(x,rho,'g'); title('\rho');
set(P,'linewidth',3);

figure(3)
hold on
P=plot(x,u,'g'); title('u');
set(P,'linewidth',3);

figure(4)
hold on
P=plot(x,p,'g'); title('p')
set(P,'linewidth',3);