//
//  Modele stratigraphique de remplissage des bassins par les sediments 
//
//
//  d_t h(x,t) + div(grad(psi(b(x,t)))) = 0 on (0,Lx)x(0,tf)
//
//  h(x,0) = hinit(x)   on (0,Lx)
//
//  grad(psi).n = g0 for x=0,  grad(psi).n = g1 for x = Lx  
//
//  b(x,t) = hmer(t) - h(x,t)
//
//
//
//  Sorties: * courbes h(x,t) fonction de x aux temps tk 
//
//   
//          * hs(x,t) = min_(t<=q<=tf) h(x,q) : surface du bassin aux temps t 
//                                              (tient compte des erosions)
//


clear

//donnees 
Lx=2; // longueur du bassin en km 


//Discretisation en espace 

n=10;
//n=100;

dx=Lx/n;

//Discretisation en temps 
tf = 1.5; // Temps final 

//ndt = 150
ndt = 1500


dt = tf/ndt; // Pas de temps 


nvisu = int(ndt/50) // plot tous les nvisu pas de temps 

// coefficients de diffusion 
Km = 1;
Kc = 10;

// flux aux bords gauche et droit 
g0 = -20;
g1 = 0;

function s = f_hmer(t)
    s = 25 + 5*cos(12*t) ;
endfunction

function s = f_psi(u)
    if (u<0) then
        s = Kc*u;
    else 
        s = Km*u;
    end
endfunction

function s = f_psip(u)
    if (u<0) then
        s = Kc;
    else 
        s = Km;
    end
endfunction


function s = f_hinit(x)
    s = 25*exp(-8*x/Lx) + 10;
endfunction


function F = f(h,t)

for i=1:n
 b(i) = f_hmer(t)-h(i);   
 psi(i) = f_psi(b(i));
end
for i=2:n-1 
 F(i) = (2*psi(i)-psi(i+1)-psi(i-1))/dx**2  
end
F(1) = ( (psi(1)-psi(2))/dx - g0 )/dx
F(n) = ( (psi(n)-psi(n-1))/dx - g1 )/dx

endfunction


function A = Df(h,t)

//A=sparse([],[],[n,n]);

A = zeros(n,n);

for i=1:n
 b(i) = f_hmer(t)-h(i);   
 psi(i) = f_psi(b(i));
 dpsi(i) = - f_psip(b(i))
end
for i=2:n-1 
 A(i,i) =   2*dpsi(i)/dx**2  
 A(i,i+1) = -dpsi(i+1)/dx**2  
 A(i,i-1) = -dpsi(i-1)/dx**2   
end
A(1,1) = dpsi(1)/dx**2
A(1,2) = -dpsi(2)/dx**2
A(n,n) = dpsi(n)/dx**2
A(n,n-1) = -dpsi(n-1)/dx**2

endfunction



//////////////////////////////////////////////
//
//    schemas 
//
//////////////////////////////////////////////
//
// fonction phi1(h,t,dt) schéma d'euler explicite 
function F = phi1(h,t,dt)

// A COMPLETER 

endfunction


// fonction phi2(h,t,dt) schéma du point milieu 
function F = phi2(h,t,dt)

// A COMPLETER

endfunction


// fonction phi3(h,t,dt) schéma de Heun 
function F = phi3(h,t,dt)

// A COMPLETER

endfunction



// schema Euler implicite 
//
///////////////////////////////////////////////////////////
//
//  Etant donnes h, t, dt, 
//
// on cherche x tel que G(x) = 0 avec G(x) = x-h - dt*f(x,t+dt)
//
//  en utilisant l'algorithme de Newton 
//
/////////////////////////////////////////////////////////
//
function F = phi4(h,t,dt)
eps = 1.0E-6;
kmax = 100;


// Algorithme de Newton pour la solution de G(x) = 0: 

x = h; // initialisation avec x = h 
r = x - h - dt*f(x,t+dt);
nr = norm(r);
nr0 = nr;
k = 0;

while (nr/nr0>eps)&(k<kmax)&(nr>eps)
    k = k +1;

// ALGORITHME DE NEWTON A COMPLETER 

end

if (k==kmax) then
    'newton non converge',x,nr
    abort
end

    F = (x-h)/dt;
    
endfunction



/////////////////////////////////////////////////////

for i=1:n
  X(i) = (i-1/2)*dx;
end  

/////////////////////////////////


//initialisation de h 
for i=1:n
 h(i)= f_hinit(X(i));    
end

hs(:,1) = h; // stockage de h init 

scf(2);
title('h')
plot (X,h,'r')  


//Boucle en temps  ------------------------------------------
t = 0;
for nt=1:ndt



/////////////////////////


// BOUCLE EN TEMPS A COMPLETER : mise a jour de h au pas de temps nt (au temps t=nt*dt)


////////////////////////

t = nt*dt 

//////// 
 
hs(:,nt) = h; // on stocke la solution au temps n 
     


 //plot h 
 scf(2);
 title('h') 
 if (int(nt/nvisu)*nvisu==nt)
  plot (X,h,'-r')  
 end

end // fin boucle en temps  ------------------------------------------------




// on enleve les parties erodees 
for i=1:n
for m=2:ndt
  j = ndt-m+1;
  hs(i,j) = min(hs(i,j),hs(i,j+1));  
end
end

scf(4);
title('hs')
for nt=1:ndt
if (int(nt/nvisu)*nvisu==nt)    
plot(X,hs(:,nt),'-b')
end
end





