%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Novembre 2000
%
% schemas numeriques pour le probleme de Riemann avec une 
% equation scalaire a flux non-convexe
% u,t + (-1/4*pi^2)*cos(2*pi*u),x = 0
%
% si u est petit, on retrouve Burgers par developpement limite.
% si (ul-ur) est grand on trouve une onde complexe faite de 
% plusieurs chocs et detentes
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 1000 ;
% dx = pas d'espace
dx = lg/nx ;
% etat gauche et droit
ul = 1.5 ;
ur = 0. ;
liminf = min(ul,ur) - 0.1 ;
limmax = max(ul,ur) + 0.1 ;
% condition cfl
a = 1./(2*pi) ;
cfl = dx/a ;
% dt = pas de temps
dt = 0.9*cfl ;
% tt = temps de simulation
tt=lg/(2*a) ;
nt = round(tt/dt) ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
% x = points du maillage
% u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
     u(i) = ul ;
  else
     u(i) = ur ;
  end
end

plot(x,u) ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift3('1D',u,+1) ;
um = wshift3('1D',u,-1) ;
v  = (up+um)/2 - (-1./(4*pi^2))*(0.5*dt/dx)*(cos(2*pi*up)-cos(2*pi*um))/2 ;
u = v ; 
if mod(n,5) == 0
plot(x,u,'*') ;
axis([0 lg liminf limmax]) ;
title('Lax-Friedrichs, cfl=0.9')
pause(0) ;
end
%
end

