/////////////////////////////////////////////////////////////
// modif des fichiers de greg
//
// ondes homogeneisées
// 
// B.D dernière version en avril 2008
/////////////////////////////////////////////////////////////
//



///////////////////////////
// parametres numeriques //
//////////////////////////

xset("font",2,3) ; 
xset("thickness",2) ;
//
pi = %pi ;
//
// lg = longueur du domaine
lg = 1. ;

// nx = nombre de mailles
nx =200 ;
// dx = pas d'espace
dx = lg/nx ;
// a = vitesse de transport
a = 1. ;
cfl = dx ;
// dt = pas de temps
dt =0.25*cfl ;

nu=dt/dx;
// nt = nombre de pas de temps effectues
ttotal=0.1;
nt = 400 ;
//


// initialisation
//
x=zeros(1,nx) ;

u=zeros(1,nx) ;
u0=zeros(1,nx) ;
u1=zeros(1,nx) ;
unew=zeros(1,nx) ;
ua=zeros(1,nx) ;
v=zeros(1,nx) ;
vnew=zeros(1,nx) ;
vb=zeros(1,nx) ;

p=zeros(1,nx) ;
q=zeros(1,nx) ;
z=zeros(1,nx) ;
pa=zeros(1,nx) ;
qb=zeros(1,nx) ;
pnew=zeros(1,nx) ;
qnew=zeros(1,nx) ;
pp=zeros(1,nx) ;
pm=zeros(1,nx) ;
qp=zeros(1,nx) ;
qm=zeros(1,nx) ;
fp=zeros(1,nx) ;
fq=zeros(1,nx) ;

fu=zeros(1,nx) ;
fv=zeros(1,nx) ;

w=zeros(1,nx) ;

up=zeros(1,nx) ;
um=zeros(1,nx) ;
upa=zeros(1,nx) ;
uma=zeros(1,nx) ;

vp=zeros(1,nx) ;
vm=zeros(1,nx) ;
vpb=zeros(1,nx) ;
vmb=zeros(1,nx) ;

rho=zeros(1,nx);
un_rho=zeros(1,nx);
sigma=zeros(1,nx);

coeffa=zeros(1,nx) ;
coeffb=zeros(1,nx) ;




////////////////////////////////////////
// Definitions des rho et sigma ///////
// C'est la que vous pouvez changer //
// les valeurs de rho et sigma     //
//  ainsi que la taille des couches/
///////////////////////////////////
nlg1=1; 
nlg2=1;

acoeff1=2.; // vous pouvez essayer plusieurs valeurs de ce coefficient, par
           // acoeff1=1, 2, 4  
bcoeff1=1./acoeff1;

acoeff2=1./acoeff1;    
bcoeff2=1./acoeff2;

cc=1;


for k=1:nx/(nlg1+nlg2)
     	for i=1+(k-1)*(nlg1+nlg2):(k-1)*(nlg1+nlg2)+nlg1
        rho(i)= acoeff1;
        sigma(i)= bcoeff1;
	end

	for i=1+(k-1)*(nlg1+nlg2)+nlg1:(k-1)*(nlg1+nlg2)+nlg1+nlg2
        rho(i)= acoeff2;
        sigma(i)= bcoeff2;
	end
        
end

for i=1:nx
      un_rho(i)=1./rho(i);
end 

////////////////////////////////////////
// Initialisation, notations          //
// en accord avec le sujet            //
///////////////////////////////////////
// x = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  u0(i)= 1e5* ( x(i)-7*lg/10)^2 *( x(i)-8*lg/10)^2 ;

end


for i=1:nx
  if x(i)<7*lg/10
      u0(i)=0 ;
  end 
  if x(i)> 8*lg/10
      u0(i)=0 ;
  end
end







///////////////////////////////////////////////
// construction d'un u1 particulier           //
// u1= derivee de u0 fois                     //
//      racine carre du produit 
//               rho(i) *sigma(i)            //
//////////////////////////////////////////////



for i=2:nx
u1(i)=((u0(i)-u0(i-1))/dx)*sqrt(  sigma(i)/rho(i));
end 
u1(1)=((u0(2)-u0(1))/dx)*sqrt(  sigma(1)/rho(1));



///////////////////////////////////////////////
// construction d'un u1 particulier           //
// u1= derivee de u0 fois                     //
//      racine carre du produit              //
//    rho  homogeneise *sigma homogeneise  //
//////////////////////////////////////////////

a2=nlg1/( (nlg1+nlg2) * acoeff1)+ nlg2/( (nlg1+nlg2) * acoeff2)    ;
//a2=1./a2;

a3= nlg1/( (nlg1+nlg2) * 1./bcoeff1)+ nlg2/( (nlg1+nlg2) * 1./bcoeff2)    ;

a4=sqrt(  a2 *a3); 



for i=2:nx
u1(i)=((u0(i)-u0(i-1))/dx)*a4;
end
u1(1)= ((u0(2)-u0(1))/dx)*a4;

//////////////////////////////
// p et q ///////////////////
/////////////////////////////
for i=1:nx
p(i)=rho(i)*u1(i);
u(i)=u0(i);
end 

for i=2:nx
q(i)=(u0(i)-u0(i-1))/dx;
end
q(1)=(u0(2)-u0(1))/dx;


for  i=1:nx
      z(i)=p(i)+q(i);
end


// affichage de la donnee initiale
xbasc() ;
//tics=[4,10,4,10];
//plotframe([0,-1.1,lg,1.1],tics);
plot2d(x,u0);
xtitle ('ondes, u0, condition initiale ') ;
halt() ;

xbasc() ;
plot2d(x,u1);
xtitle ('ondes, u1, condition initiale ') ;
halt() ;

xbasc() ;
plot2d(x,p);
xtitle ('ondes, p, condition initiale ') ;
halt() ;

xbasc() ;
plot2d(x,q);
xtitle ('ondes, q, condition initiale ') ;
halt() ;

xbasc() ;
plot2d(x,z);
xtitle ('ondes, z, condition initiale ') ;
halt() ;


/////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs 
/////////////////////////////////////////////////
//
// w est la solution exacte
// si u est le vecteur des valeurs u(i) alors 
//   up est le vecteur des u(i+1) 
//   um est le vecteur des u(i-1)
// 

T=0;
for n=1:nt;

//for n=1:1
T=T+dt;


for i=1:nx
u(i)=u(i)+p(i)*dt/rho(i);
end 



   for i=1:nx
       pa(i)=p(i)*un_rho(i);
       qb(i)=q(i)*sigma(i);
end

//
pp(1:nx-1)=pa(2:nx);
pp(nx)=pa(1);

pm(2:nx)=pa(1:nx-1);
pm(1)=pa(nx);

qp(1:nx-1)=qb(2:nx);
qp(nx)=qb(1);
qm(2:nx)=qb(1:nx-1);
qm(1)=qb(nx);




for i=1:nx
fp(i)=0.5*( qb(i) + qp(i))- (0.5* cc)*(pa(i)-pp(i));

fq(i)=0.5*( pa(i) + pp(i))- (0.5*cc)*(qb(i)-qp(i));
end

for i=2:nx
pnew(i)=p(i)+ nu *(fp(i)-fp(i-1)); 
qnew(i)=q(i)+ nu *(fq(i)-fq(i-1)); 
end

pnew(1)=p(1)+ nu *(fp(1)-fp(nx)); 
qnew(1)=q(1)+ nu *(fq(1)-fq(nx)); 

p=pnew;
q=qnew;

for  i=1:nx
      z(i)=p(i)+q(i);
end

////////////////////////////////////////
// affichage de la solution approchee //
//     tous les 10 pas de temps      //
//////////////////////////////////////
if modulo(n,10) == 0


  xbasc() ;
//  plotframe([0,-1.1,lg,1.1],tics);
//  xset("thickness",2) ;
 // plot2d(x,ua);
 //  plot2d(x,up);
  //plot2d(x,um);
//  plot2d(x,unew);

  plot2d(x,u);

//plot2d(x,v)
xtitle ('ondes, u, temps intermediaire  ') ;
T
//halt();
//  xset("thickness",1) ;
//  xbasc() ;
 // plot2d(x,v)
//xtitle ('ondes, v ') ;
//halt();
//  xtitle('transport, Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;


//