function de=stagr
% Lagrangian ideal gas - Riemann Problem, staggered mesh+artificial viscosity
global T t dt J x x2 mn mc u un Bl Br Bc ur ul 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]
% v u je [rho,u,eps]
% C je CFL cislo
% resime rovnici u_t + f_m = 0
Bl = 0; Br = 1; 
gamma = 1.4;
% Sod
ul = [1 0 1]; ur = [0.125 0 0.1]; Bc = 0.5; T = 0.2;
%Bl = -1; Br = 2; T=0.6;
% 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;
%C = 0.5;
% 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=-2; Br=2; Bc=0; T=1;gamma = 5/3;
%C=0.5;

% piston
%ul = [1 0 2/3*1e-4]; ur = [1 0 2/3*1e-4]; Bc = 0.5; T = 0.5;gamma = 5/3;C=0.0;dt = 0.00025; 
%ul = [1 0 2/3]; ur = [1 0 2/3]; Bc = 0.5; T = 0.5;gamma = 5/3;C=0.0;dt = 0.00025; 
% other single shock
%ul = [4 1 92/63]; Bc=0.5; T=0.3;
%ur = [1 0 8/63];

J = 100;
solve;

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

t=0;
nstep = 1;
init;
du = u;
while t < T
  time_step;
  one_step(nstep);
  plot(x,un(:,1),'bx-');hold on; plot(x,un(:,2),'rx-');
  plot(x,p,'gx-'); %plot(x,q/maxq/2,'mx-'); plot(x,s,'cx-');
  hold off;
  drawnow;
%  pause;
  u = un;
  t = t + dt; nstep = nstep + 1;
end
nstep = nstep

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

tmass = ul(1)*(Bc - Bl) + ur(1)*(Br - Bc);
dm = tmass/(J);
dx = (Br - Bl)/J;
x = zeros(J+1,1); mn = zeros(J+1,1);
x2 = zeros(J,1); mc = zeros(J,1);
u = zeros(J,3);
un = zeros(J,3);q = zeros(J,1);
x(1) = Bl;
for j=2:J+1
% mesh uniform in m  
  if(x(j-1) < Bc)
    dx = dm/uul(1);
    u(j-1,:) = uul;
  else
    dx = dm/uur(1);
    u(j-1,:) = uur;
  end
% mesh uniform in x
  dx = (Br - Bl)/J;
  x(j) = x(j-1) + dx;
  x2(j-1) = x(j-1) + dx/2;
  mc(j-1) = dx*u(j-1,1);
end
u(J+1,2) = uur(2);
% hmotnost uzlu
mn(1) = mc(1)/2;
for j=2:J
  mn(j) = (mc(j) + mc(j-1))/2;
end
mn(J+1) = mc(J)/2;
% for time_step
u(J+1,1) = u(J,1);
u(J+1,3) = u(J,3);

function stag=one_step(nstep)
% jeden casovy krok
global T t dt J x x2 mn mc u un Bl Br Bc ur ul gamma p q C

for j=1:J
% tlak
  p(j) = u(j,1)*u(j,3)*(gamma - 1);
  du = (u(j+1,2) - u(j,2));
  rho = u(j,1);
  cs2 = gamma*p(j)/rho;
  gamma4 = (gamma + 1)/4;
% Kuropatenko form of AV
  q2(j) = rho*(gamma4*du - sqrt(gamma4^2*du^2 + cs2))*du;
  divu = (u(j+1,2) - u(j,2));
  if divu >= 0
      q(j) = 0;
  else 
      q(j) = q2(j);
  end
end
q(J+1) = q(J);

% rovnice pro rychlost
for j=2:J
  un(j,2) = u(j,2) - dt*(p(j) + q(j) - p(j-1) - q(j-1))/mn(j);
end
% okrajove podminky pro rychlost
un(1,2) = un(2,2);
%un(1,2) = -1;  % for piston
un(J+1,2) = un(J,2);

% pohyb site
for j=1:J+1
  x(j) = x(j) + dt*(u(j,2) + un(j,2))/2;
end

% hustota z pohybu site
for j=1:J
  un(j,1) = mc(j)/(x(j+1) - x(j));
end

% rovnice pro vnitrni energii
for j=1:J
  un(j,3) = u(j,3) - dt*(p(j) + q(j))*  ...
            (un(j+1,2) + u(j+1,2) - un(j,2) - u(j,2))/2/mc(j);
end

% for time_step
un(J+1,1) = un(J,1);
un(J+1,3) = un(J,3);

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

p = (gamma - 1)*u(:,3).*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
if not(C==0) 
  dt = C*min(mn./  ...
             (abs(u(:,2)) + sqrt(gamma*p.*u(:,1))))
end






