///////////////////////////////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2002, Janvier 2012
//
// 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
//
///////////////////////////////////////////////////////////////////////////////////////////
//
xset("font",2,3) ; 
xset("thickness",2) ;
//
pi = %pi ;
//
// 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

clf()
tics=[4,10,4,10];
plotframe([0,liminf,lg,limmax],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Friedrichs, cfl=0.9' ,' ',' ') ;

halt() ;
///////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
///////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftn('+1',u) ;
um = shiftn('-1',u) ;
v  = (up+um)/2 - (-1./(4*pi^2))*(0.5*dt/dx)*(cos(2*pi*up)-cos(2*pi*um))/2 ;
u = v ; 
if modulo(n,5) == 0
  clf() 
  plotframe([0,liminf,lg,limmax],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xtitle('Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end