function de=pde
% Lagrangian ideal gas - Riemann Problem
global T t dt J x x2 m2 u un Bl Br Bc ur ul nlwlf gamma p C
% T je konecny cas
% dt je casovy  krok
% J je pocet bunek
% x(j) je prostorova souradnice j-teho bodu
% u(j) je reseni v bode x(j) na casove vrstve n
% un(j) je reseni v bode x(j) na casove vrstve n+1
% resime na intervalu (Bl,Br), Bc je bod ve kterem je pocatecni
%       nespojitost
% ur,ul jsou leva a prava poc. hodnota RP [rho,u,p]
% C je CFL cislo
% resime rovnici u_t + f_m = 0
Bl = 0; Br = 1; 
gamma = 1.4;
% test 1
ul = [1 0.75 1]; ur = [0.125 0 0.1]; Bc = 0.3; T = 0.2;
C = 1;
% test 2
%ul = [1 -2 0.4]; ur = [1 2 0.4]; Bc = 0.5; T = 0.3;
% test 3
%ul = [1 0 1000]; ur = [1 0 0.01]; Bc = 0.5; T = 0.012;
% test 4
%ul = [5.99924 19.5975 460.894]; Bl = -0.4; Br = 1.2;
%ur = [5.99924 -6.19633 46.095 ]; Bc = 0.4; T = 0.035;
%C = 0.5;
% Noh
%ul = [1 1 1.e-6]; ur = [1 -1 1.e-6]; Bl=-3; Br=3; Bc=0; T=2;gamma = 5/3;
%C=0.75;

J = 200;
nlwlf = 1;
solve;

function so=solve
% resi PDR u_t = a u_x
global T t dt J x x2 m2 u un Bl Br Bc ur ul nlwlf gamma p C

t=0;
nstep = 1;
init;
du = u;
while t < T
  time_step;
  one_step(nstep);
  plot(x2,ones(J,1)./un(:,1),'x-');
%  plot(x2,un(:,3));
%  axis([-B B 0 1]);
  drawnow;
%  pause;
  u = un;
  t = t + dt; nstep = nstep + 1;
end

function ii=init
% inicializace, pocatecni podminky
global T t dt J x x2 m2 u un Bl Br Bc ur ul nlwlf gamma p C
uul = [1/ul(1), ul(2), ul(3)./((gamma - 1)*ul(1)) + ul(2)^2/2];
uur = [1/ur(1), ur(2), ur(3)./((gamma - 1)*ur(1)) + ur(2)^2/2];

tmass = ul(1)*(Bc - Bl) + ur(1)*(Br - Bc);
dm = tmass/(J);
dx = (Br - Bl)/J;
x = zeros(J+1,1);
x2 = zeros(J,1); m2 = zeros(J,1);
u = zeros(J,3);
un = zeros(J,3);
x(1) = Bl;
for j=2:J+1
  if(x(j-1) < Bc)
    dx = dm*uul(1);
    u(j-1,:) = uul;
  else
    dx = dm*uur(1);
    u(j-1,:) = uur;
  end
  x(j) = x(j-1) + dx;
  x2(j-1) = x(j-1) + dx/2;
  m2(j-1) = dx/u(j-1,1);
end

function lf=one_step(nstep)
% jeden casovy krok Lax-Friendrichsovo schema
global T t dt J x x2 m2 u un Bl Br Bc ur ul nlwlf gamma p C

% predictor
% f is flux, u2 = U^{n+1/2}_i, us = u^n_i - interpolated velocity
f = zeros(J,3); u2 = zeros(J+1,3); us = zeros(J+1);
p = (gamma - 1)*(u(:,3) - u(:,2).*u(:,2)/2)./u(:,1);
f(:,1) = u(:,2);
f(:,2) = -p;
f(:,3) = -p.*u(:,2);

for j=2:J
  u2(j,:) = (m2(j-1)*u(j-1,:) + m2(j)*u(j,:))/(m2(j-1) + m2(j)) ...
	    + dt/(m2(j-1) + m2(j))*(f(j,:) - f(j-1,:));
% interpolate velocities on time level n
  us(j) = ((x2(j) - x(j))*u(j-1,2) + (x(j) - x2(j-1))*u(j,2))/ ...
	     (x2(j) - x2(j-1));
end
us(1) = u(1,2);
us(J+1) = u(J,2);
u2(1,:) = u2(2,:); u2(J+1,:) = u2(J,:);
% advance x_i
for j=1:J+1
  x(j) = x(j) + dt/2*(us(j) + u2(j,2));
%  x(j) = x(j) + dt*u2(j,2);
end
for j=1:J
  x2(j) = (x(j) + x(j+1))/2;
end

%plot(x,ones(J+1,1)./u2(:,1)); pause;
%plot(x,u2(:,3)); pause;
% corrector
% density from geometry
for j=1:J
  un(j,1) = (x(j+1) - x(j))/m2(j);
end
f = zeros(J+1,3); 
p = (gamma - 1)*(u2(:,3) - u2(:,2).*u2(:,2)/2)./u2(:,1);
f(:,1) = u2(:,2);
f(:,2) = -p;
f(:,3) = -p.*u2(:,2);

m1 = 2;  % defines if density from geometry or from pde
if(mod(nstep,nlwlf)==0 | t >= T)
  for m=m1:3
    for j=2:J-1
      un(j,m) = (u2(j+1,m) + u2(j,m))/2 ...
		+ dt/(2*m2(j))*(f(j+1,m) - f(j,m));
    end
  end
else
  for m=m1:3
    for j=2:J-1
      un(j,m) =  u(j,m) ...
		+ dt/m2(j)*(f(j+1,m) - f(j,m));
    end
  end
end
% okrajove podminky
un(1,:) = un(2,:);
un(J,:) = un(J-1,:);

function ti=time_step
global T t dt J x x2 m2 u un Bl Br Bc ur ul nlwlf gamma p C

p = (gamma - 1)*(u(:,3) - u(:,2).*u(:,2)/2)./u(:,1);
sp = size(p);
for i=1:sp(1)
  if(p(i,1)<1.e-6) 
    p(i,1) = 1.e-6;
  end
end

dt = C*min(m2./  ...
            ((abs(u(:,2)) + sqrt(gamma*p.*u(:,1))) ./ u(:,1)));






