// resoud le pb dU/dt +A dU/dx=0
// avec le schema upwind ordre 1
// ou U=(p,v) et A correspond a l'acoustique

// fonction pour l'initialisation des champs p et v

function [U]=init_acoustique(a,b,N,c,x)
   U=zeros(2,N+5);

  // y=exp(-x.*x*10);
   // U(1,3:N+3)=y;

   y=sign(x);
   U(1,3:N+3)=y;
   U(2,3:N+3)=0;

endfunction;


a=-1; b=1;
N=200;
dx=(b-a)/N;
x=[a:dx:b];
Tf=4;

rho=1;
c=0.5;
CFL=.9;

 A=[0 rho*c^2;1/rho 0];
 lambda=[-c;c];
 r1=[-rho*c;1];
 r2=[rho*c;1];
 L1=0.5/c/rho*[-1 rho*c];
 L2=0.5/c/rho*[1 rho*c];
 P=[r1 r2];
 invP=[L1;L2];
 Aplus=P*[0 0;0 c]*invP;
 Amoins=P*[-c 0;0 0]*invP;

dt=CFL*dx/abs(c);

F=zeros(2,N+5);
div=zeros(2,N+5);

U=zeros(2,N+5);
U=init_acoustique(a,b,N,c);

p_init=U(1,3:N+3);

v_init=U(2,3:N+3);


t=0;
while t<Tf
 
 

 //conditions limites periodiques
    U(:,1)=U(:,N+1);
    U(:,2)=U(:,N+2);
    U(:,N+4)=U(:,4);
    U(:,N+5)=U(:,5);

 
 //calcul des flux Fi+1/2

 for i=2:N+3,
   F(:,i)=Aplus*U(:,i)+Amoins*U(:,i+1);
 end

 for i=3:N+3,
   div(:,i)=dt/dx*(F(:,i)-F(:,i-1));
 end

 U=U-div;
 
  xset("window",0)
  drawlater()
  clf();
  plot2d(x,p_init,2);
  plot2d(x,U(1,3:N+3));
  drawnow()


 t=t+dt;

  
end;

xset("window",0)
clf();
plot2d(x,p_init,2);
plot2d(x,U(1,3:N+3));

xset("window",1)
clf();
plot2d(x,v_init,2);
plot2d(x,U(2,3:N+3));



