///////////// Equation d'advection avec différents schemas //////////////
/////////////////////////////////////////////////////////////////////////

/////////////////////////////parametres//////////////////////////////////

T = 0.1;                      // Longueur de l'intervalle en temps
Xmax = 1;                      // Longueur de l'intervalle en espace

J = 50;                        // Nombre de pas en espace

dx = Xmax/J;                   // Pas en espace

nbC = 0.3;                     // Nombre de Courant.

x=(1:J)'*dx - 0.5*dx; // Grille en x

deff('[y]=vitesse(x,xmax)', 'y=1. + 0.*x;') 
//deff('[y]=vitesse(x,xmax)', 'y=2. + sin(2*pi*x);') 

//deff('[y]=CI(x,xmax)', 'y=1.*(x>0.5*xmax)+0.*(x<=0.5*xmax);') 
//deff('[y]=CI(x,xmax)', 'y=(x-0.5).*(x>0.5*xmax)+(0.5-x).*(x<=0.5*xmax);') 
deff('[y]=CI(x,xmax)', 'y=sin(2*pi*x);') 

u = zeros(J,1); // Solution numerique
uG = zeros(J,1); 
uD = zeros(J,1); 

t=0;
a = vitesse(x,Xmax);
u=CI(x,Xmax);          // On initialise la solution. 

A = max(a)
dt = min(dx*nbC/A,T);            // Pas de temps.

//// Boucle en temps

while (t<T)
    dt = min(dt,T-t);
    t = t + dt;            // Mise a jour de la date. 

//    for j=2:J,             // Décentré à gauche. 
//        uG(j) = u(j-1);   
//    end;
//    uG(1) = u(J);   
 
//    for j=2:J,             // Décentré à droite. 
//        uG(j) = u(j);   
//    end;
//    uG(1) = u(1);   
  
    for j=2:J,             // Centré. 
        uG(j) = 0.5*(u(j-1) + u(j));   
    end;
    uG(1) = 0.5*(u(J) + u(1));
    
    for j=1:J-1,
        uD(j) =  uG(j+1)   
    end;
    uD(J) =  uG(1)   
    
    for j=1:J,
        u(j) = u(j) - a(j)*dt/dx*(uD(j) - uG(j));
    end;

//  xset("wwpc");
    clf();
    plot2d(x,u);
//  clf();
//  plot2d(x,u,[1],"111","u(x,t)",[0,0,1,1]);
//  show_pixmap()
    sleep(10);
//pause;
//  halt;  

end;
