
//clear; close;
//da=gda(); da.auto_clear = 'off';
//
// Donnes d entree

T=input("Duree totale en temps ")
a=0
b=1.
Va=1
// choix en espace

J=input("nombre de points en espace ")
delta_x=(b-a)/J
disp(delta_x,'delta x')
co_a=input("CFL ")
//disp (co_a,'co_a')
delta_t = abs(co_a)*delta_x/abs(Va)
disp(delta_t,'delta_t')
N=floor(T/delta_t)
disp(N,'nombre d iterations')

L=input("numero du mode de Fourier ")
// donnees initiales regulieres
D=1
function y=u0(x)
       x1=x-floor(x)
       y=sin(2*%pi*L*x1/D)
endfunction

// donnees initiales regulieres
//function y=u0(x)
//    x1=x-floor(x)
//    y=sin(2*%pi*x1)
//endfunction

// donnees chocs
//function y=u0(x0)
//    x1=x0-floor(x0)
// y=double(x1<0.25)+double(x1>0.75)
//endfunction
//

x=linspace(a,b,J+1)

Vc=zeros(J+1)
Vs=zeros(J+1)

W0=feval(x,u0)

Vc=W0

Nt=input("nombre d instantanés pour la visu ")
pasN=N/(Nt-1)
disp(pasN,"pas en temps des visus de sortie ")
// W_d=spzeros(J+1,Nt)
W_d=zeros(size(x,2),Nt)
count=0
pascount=0
nn=1

 W_d(:,nn)=W0(:)
// W_d=W0
for n=1:1:N
    count=count+1    
    if (count==floor(pasN)) then
        pascount=pascount+1
        count=0
        nn=nn+1
        xr=x-Va*n*delta_t
        W_d(:,nn)=feval(xr(:),u0)
       end
end
plot2d(x,W_d(:,1:Nt))



count=0
pascount=0
nn=1
V_d=zeros(size(x,2),Nt)
V_d(:,nn)=Vc(:)

for n=1:1:N
    count=count+1
// schema centre explicite : instable
 //  Vs=Vc+0.5*co_a*[Vc(J),Vc(1:J)] -0.5*co_a*  [Vc(2:J+1),Vc(2)]
//// schema decentre amont : stable 0 inf coa inf 1

  Vs=(1-co_a)*Vc+co_a*[Vc(J),Vc(1:J)]

// schema Lax Wendroff: stable 0 inf coa inf 1
//
//    Vs=(1-co_a*co_a)*Vc+0.5*co_a*(1+co_a)*[Vc(J),Vc(1:J)]-0.5*co_a*(1-co_a)*[Vc(2:J+1),Vc(2)]
//  
    if (count==floor(pasN)) then
        pascount=pascount+1
        disp(pascount*floor(pasN),'iteration')
        count=0
        nn=nn+1
        V_d(:,nn)=Vs(:)
      
    end
    Vc=Vs

end

plot2d(x,V_d(:,1:Nt))
scf
//for i=1:1:Nt
//subplot(Nt,1,Nt-i+1)
//plot(x',W_d(:,i),'b',x',V_d(:,i),'r')
//end
