/////////////////////////////////////////////////////////////
// Paola Goatin, sett. 2006
//
// schemas numeriques sur l'equation de transport lineaire 
// avec condition aux limites de Neumann
// u,t + a u,x = 0
// 
/////////////////////////////////////////////////////////////
//
funcprot(0)
xset("font",2,3) ; 
xset("thickness",2) ;

exec('soltr.sci'); 
exec('shiftp.sci');

//
pi = %pi ;
//
// lg = longueur du domaine
lg = 1. ;
// nx = nombre de mailles
nx = 100 ;
// dx = pas d'espace
dx = lg/nx ;
// a = vitesse de transport
a = 1. ;
cfl = dx/a ;
// dt = pas de temps
dt = 0.5*cfl ;
// nt = nombre de pas de temps effectues
nt = 100 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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/4
  u(i) = 1. ;
  else
  u(i) = 0. ;
  end  
end

w = u ;

saisie = 4;


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 1: Lax-Friedrichs avec cfl=0.5
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
if saisie == 1
//  dt = 0.9*cfl ;
//  nt = 400 ;

  clf()
  tics=[4,10,4,10];
  plotframe([0,-0.3,lg,1.3],tics);
  plot2d(x,u,1,"000")
  xtitle ('transport, Lax-Friedrichs, cfl=0.5' ,' ',' ') ;

  xpause(800000) ;
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // marche en temps : Lax-Friedrichs
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  for n=1:nt
  //
  up = shiftp('+1',u) ;
  um = shiftp('-1',u) ;
  v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
  u = v ; 
  if modulo(n,5) == 0
    for i=1:nx
        xi = floor((lg/4+a*n*dt)/dx)  ;
        if i <= xi
           w(i) = 1. ;
        else
           w(i) = 0. ;
        end
    end
    clf()
    plotframe([0,-0.3,lg,1.3],tics);
    xset("thickness",2) ;
    plot2d(x,u,[1,1],"100","Lax-Friedrichs")
    xset("thickness",1) ;
    plot2d(x,w,[2,4],"100","exact")
    xtitle('transport, Lax-Friedrichs, cfl=0.5',' ',' ');
    xpause(100000) ;
  end
  //
  end
end


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 2: Upwind avec cfl=0.5
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
if saisie == 2
//  dt = 0.9*cfl ;
//  nt = 400 ;

  clf()
  tics=[4,10,4,10];
  plotframe([0,-0.3,lg,1.3],tics);
  plot2d(x,u,1,"000")
  xtitle ('transport, Upwind, cfl=0.5' ,' ',' ') ;

  xpause(800000) ;
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // marche en temps : Upwind
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  for n=1:nt
  //
  um = shiftp('-1',u) ;
  v = u - (dt/dx)*a*(u-um) ;
  u = v ; 
  if modulo(n,5) == 0
    for i=1:nx
        xi = floor((lg/4+a*n*dt)/dx)  ;
        if i <= xi
           w(i) = 1. ;
        else
           w(i) = 0. ;
        end
    end
    clf()
    plotframe([0,-0.3,lg,1.3],tics);
    xset("thickness",2) ;
    plot2d(x,u,[1,1],"100","Upwind")
    xset("thickness",1) ;
    plot2d(x,w,[2,4],"100","exact")
    xtitle('transport, Upwind, cfl=0.9',' ',' ');
    xpause(100000) ;
  end
  //
  end
end

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 3: Lax-Wendroff avec cfl=0.5
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
if saisie == 3
//  dt = 0.9*cfl ;
//  nt = 400 ;

  clf()
  tics=[4,10,4,10];
  plotframe([0,-0.3,lg,1.3],tics);
  plot2d(x,u,1,"000")
  xtitle ('transport, Lax-Wendroff, cfl=0.5' ,' ',' ') ;

  xpause(800000) ;
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // marche en temps : Lax-Wendroff
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  for n=1:nt
  //
  up = shiftp('+1',u) ;
  um = shiftp('-1',u) ;
  v = u - (0.5*dt/dx)*a*(up-um) + (0.5*(dt/dx)^2)*(a^2)*(up-2*u+um) ;
  u = v ; 
  if modulo(n,5) == 0
    for i=1:nx
        xi = floor((lg/4+a*n*dt)/dx)  ;
        if i <= xi
           w(i) = 1. ;
        else
           w(i) = 0. ;
        end
    end
    clf()
    plotframe([0,-0.3,lg,1.3],tics);
    xset("thickness",2) ;
    plot2d(x,u,[1,1],"100","Lax-Wendroff")
    xset("thickness",1) ;
    plot2d(x,w,[2,4],"100","exact")
    xtitle('transport, Lax-Wendroff, cfl=0.5',' ',' ');
    xpause(100000) ;
  end
  //
  end
end


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 4: Beam-Warming avec cfl=0.5
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
if saisie == 4
//  dt = 0.9*cfl ;
//  nt = 400 ;

  clf()
  tics=[4,10,4,10];
  plotframe([0,-0.3,lg,1.3],tics);
  plot2d(x,u,1,"000")
  xtitle ('transport, Beam-Warming, cfl=0.5' ,' ',' ') ;

  xpause(800000) ;
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  // marche en temps : Beam-Warming
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //
  for n=1:nt
  //
  um = shiftp('-1',u) ;
  umm = shiftp('-1',um) ;
  v = u - (0.5*dt/dx)*a*(3*u-4*um+umm) + (0.5*(dt/dx)^2)*(a^2)*(u-2*um+umm) ;
  u = v ; 
  if modulo(n,5) == 0
    for i=1:nx
        xi = floor((lg/4+a*n*dt)/dx)  ;
        if i <= xi
           w(i) = 1. ;
        else
           w(i) = 0. ;
        end
    end
    clf()
    plotframe([0,-0.3,lg,1.3],tics);
    xset("thickness",2) ;
    plot2d(x,u,[1,1],"100","Beam-Warming")
    xset("thickness",1) ;
    plot2d(x,w,[2,4],"100","exact")
    xtitle('transport, Beam-Warming, cfl=0.5',' ',' ');
    xpause(100000) ;
  end
  //
  end
end


//end