//
//  Propagation d'une onde de pression dans une artere saine
//  avec un modele fluide/structure 1D
//
//  Attention: par rapport a certaines references
//  la loi de pression est p=beta*(A^0.5-A0^0.5)
//  
//
//-------------------------------------------------------------
// ---------------  Fonctions auxiliaires ---------------------
//-------------------------------------------------------------
//
//-----fonctions associees a l'EDP----------
//
function [F1,F2]=F(A,Q);
F1=Q;
F2=(Q.^2)./A+1/3/rho*bet.*(A.^1.5-A0.^1.5);
endfunction
//
function [B1,B2]=B(A,Q);
B1=zeros(A);
B2=-Kr*Q./A;
endfunction
//
function [HU11,HU12,HU21,HU22]=H(A,Q);
HU11=zeros(A);
HU12=ones(A);
HU21=-(Q.^2)./(A^2)+1/2/rho*bet.*sqrt(A);
HU22=2*Q./A;
endfunction
//
function [W1,W2]=W12(A,Q,bet);
W1=Q./A+2*sqrt(2/(rho)*bet.*sqrt(A));
W2=Q./A-2*sqrt(2/(rho)*bet.*sqrt(A));
endfunction
//
function [A,Q]=invW12(W1,W2,bet);
A=(rho./bet).^2*(W1-W2).^4/4^5;
Q=A.*(W1+W2)/2;
endfunction
//
function [LAMB1,LAMB2]=LAMB12(A,Q,bet);
LAMB1=Q./A+sqrt(1/(2*rho)*bet.*sqrt(A));
LAMB2=Q./A-sqrt(1/(2*rho)*bet.*sqrt(A));
endfunction
//
function [BU11,BU12,BU21,BU22]=BU(A,Q);
BU11=zeros(A);
BU12=zeros(A);
BU21=Kr*Q./(A.^2);
BU22=-Kr./A;
endfunction

// Left eigenvalue
//-----------------
function [l1,l2]=Lvect(A,Q,bet)     
[LAMB1,LAMB2]=LAMB12(A,Q,bet);
l1= [-LAMB2,1];
l2=[-LAMB1,1]
endfunction
//
function [FLW1,FLW2]=FLW(A,Q,dt);
[F1,F2]=F(A,Q);
[B1,B2]=B(A,Q);
[HU11,HU12,HU21,HU22]=H(A,Q);
FLW1=F1+dt/2*(HU11.*B1+HU12.*B2)
FLW2=F2+dt/2*(HU21.*B1+HU22.*B2)
endfunction
//
function [BLW1,BLW2]=BLW(A,Q,dt);
[B1,B2]=B(A,Q);
[BU11,BU12,BU21,BU22]=BU(A,Q);
BLW1=B1+dt/2*(BU11.*B1+BU12.*B2)
BLW2=B2+dt/2*(BU21.*B1+BU22.*B2)
endfunction
//
//
//-----matrices EF--------------------------
//
function M=Masse(z)
n=length(z)
M=zeros(n,n);
M(1,1)=(z(2)-z(1))/3;
M(1,2)=(z(2)-z(1))/6;
M(2,1)=(z(2)-z(1))/6;
for i=2:n-1
  M(i,i)=(z(i)-z(i-1))/3+(z(i+1)-z(i))/3;
  M(i,i+1)=(z(i+1)-z(i))/6;
  M(i+1,i)=(z(i+1)-z(i))/6;
end
M(n,n)=(z(n)-z(n-1))/3;
endfunction
//
function M=Ggradient(f)
n=length(f);
M=zeros(n,n);
M(1,1)=(-2*f(1)-f(2))/6;
M(1,2)=(-2*f(2)-f(1))/6;
M(2,1)=(2*f(1)+f(2))/6;
for i=2:n-1
  M(i,i)=(2*f(i-1)+f(i))/6-(2*f(i)+f(i+1))/6;
  M(i,i+1)=-(2*f(i+1)+f(i))/6;
  M(i+1,i)=(2*f(i)+f(i+1))/6
end
M(n,n)=(2*f(n-1)+f(n))/6;
endfunction
//
function M=Ddiv(f)
n=length(f);
M=Ggradient(f)';
M=M(2:n-1,:)
endfunction
//
function M=Llaplacien(f,z)
n=length(z);
M=zeros(n,n);
M(1,1)=(f(1)+f(2))/2/(z(2)-z(1));
M(1,2)=-(f(1)+f(2))/2/(z(2)-z(1));
M(2,1)=-(f(1)+f(2))/2/(z(2)-z(1));
for i=2:n-1
  M(i,i)=(f(i)+f(i-1))/2/(z(i)-z(i-1))+(f(i+1)+f(i))/2/(z(i+1)-z(i));
  M(i,i+1)=-(f(i+1)+f(i))/2/(z(i+1)-z(i));
  M(i+1,i)=-(f(i+1)+f(i))/2/(z(i+1)-z(i));
end
M(n,n)=(f(n)+f(n-1))/2/(z(n)-z(n-1));
//
M=M(2:n-1,:)
endfunction
/////////////////////////////////////////////////////////
//-------------------------------------------------------------
// ---------------Programme principal -------------------------
//-------------------------------------------------------------
///////////////////////////////////////////////////////////
//
clf()
clear all;
//
//
/// caracteristiques de l'artere  ////////////////////
//
L1=5.9;
rho=1.05;
A0val=0.5227943;
Kr=0.75;
bet=649E3;
//
// discrétisation spatiale et temporelle //////////////
//
N1=10;                      // nombre de points intérieurs
dt=8E-4;                    // pas de temps
//
z=linspace(0,L1,N1);  
N=length(z)-2;L=z($);
//
nP=max(find(z<0.25*L1));
nM=max(find(z<0.5*L1));
nD=max(find(z<0.75*L1));
//
A=zeros(N+2,1);
Q=zeros(N+2,1);
A0=A0val*ones(z');
//
////// données expérimentales (section) /////////////
//
Data=read('artere-illiac.txt',-1,15); 
diam = Data(16:126,8)./10000;
Tfin=Data(126,15)-Data(16,15);
texp=Data(16:126,15)-Data(16,15);
Abord0 = interp1(texp,diam',0:dt:Tfin,'linear');
Abord0 = (Abord0/2)^2*%pi;
//
//
//////////tracé de la donnée initiale (section)///
Nt=int(Tfin/dt+1E-12);
t=(0:dt:Tfin);
xset('window',0)
clf()
plot2d(t,Abord0)
xtitle('section en entrée t-->A(t)');
///////////////////////////////////////////////////
//
Ainit=A0;
Qinit=zeros(N+2,1);
//
[W1init0,W2init0]=W12(Ainit(1),Qinit(1),bet);
[W1initL,W2initL]=W12(Ainit(N+2),Qinit(N+2),bet);
W1nplus1L=W1initL;
W2nplus1L=W2initL;
//
g1=W2init0*ones(t)+4*sqrt(2/rho)*sqrt(bet*sqrt(Abord0));
g2=W2initL*ones(t);
//
M=Masse(z);
Mbis=M(2:N+1,:);
Minv=inv(M(2:N+1,2:N+1));
G=Ggradient(ones(z));
Gbis=G(2:N+1,:);
//
A=Ainit;
Q=Qinit;
Anew=Ainit;
Qnew=Qinit;
pnew=Anew;
prP=[0];prM=[0];prD=[0];AP=[A0val];AM=[A0val];AD=[A0val];QP=[0];QM=[0];QD=[0];
Atot=A;Qtot=Q;
/////////////////////////////////////////////////////////
//----------------debut boucle principale-------------
/////////////////////////////////////////////////////////
timer();
xset('window',1)
//
for i=1:Nt
//
tpsnew=i*dt;
A=Anew;
Q=Qnew;
p=pnew;
//
[LAMB1,LAMB2]=LAMB12(A(1),Q(1),bet);
[W1,W2]=W12(A,Q,bet);

W2nplus10=interpln([z;W2'],-LAMB2*dt);
W1nplus10=W2nplus10 + 4*sqrt(2/rho)*sqrt(bet*sqrt(Abord0(i+1)));


[Anplus10,Qnplu10]=invW12(W1nplus10,W2nplus10,bet);
//
[LAMB1,LAMB2]=LAMB12(A(N+2),Q(N+2),bet);
W1nplus1L=interpln([z;W1'],L-LAMB1*dt);
W2nplus1L=g2(i+1);
[Anplus1L,Qnplus1L]=invW12(W1nplus1L,W2nplus1L,bet);
//
Anew(1)=Anplus10;
Qnew(1)=Qnplu10;
Anew(N+2)=Anplus1L;
Qnew(N+2)=Qnplus1L;
//
[F1,F2]=F(A,Q);
[FLW1,FLW2]=FLW(A,Q,dt);
[BLW1,BLW2]=BLW(A,Q,dt);
[BU11,BU12,BU21,BU22]=BU(A,Q);
[HU11,HU12,HU21,HU22]=H(A,Q);
//
RHS1part1=Mbis*A+dt*Gbis*FLW1+dt*Mbis*BLW1;
RHS2part1=Mbis*Q+dt*Gbis*FLW2+dt*Mbis*BLW2;
RHS1part2=-dt^2/2*(Llaplacien(HU11,z)*F1+Llaplacien(HU12,z)*F2)-dt^2/2*(Ddiv(BU11)*F1+Ddiv(BU12)*F2);
RHS2part2=-dt^2/2*(Llaplacien(HU21,z)*F1+Llaplacien(HU22,z)*F2)-dt^2/2*(Ddiv(BU21)*F1+Ddiv(BU22)*F2);//
//
RHS1=RHS1part1+RHS1part2;
RHS2=RHS2part1+RHS2part2;
//
  RHS1(1)=RHS1(1)-Anew(1)*M(2,1);
  RHS1(N)=RHS1(N)-Anew(N+2)*M(N+1,N+2);
  RHS2(1)=RHS2(1)-Qnew(1)*M(2,1);
  RHS2(N)=RHS2(N)-Qnew(N+2)*M(N+1,N+2);
//
Aprime=Minv*RHS1;
Qprime=Minv*RHS2;
//
Anew(2:N+1)=Aprime;
Qnew(2:N+1)=Qprime;
pnew=bet.*(sqrt(Anew)-sqrt(A0));
AP=[AP,Anew(nP)];QP=[QP,Qnew(nP)];prP=[prP,pnew(nP)];
AM=[AM,Anew(nM)];QM=[QM,Qnew(nM)];prM=[prM,pnew(nM)];
AD=[AD,Anew(nD)];QD=[QD,Qnew(nD)];prD=[prD,pnew(nD)];
//
///////////////// tracé dynamique en fonction du temps /////////
drawlater
clf()
xsetech([0,0,1,0.5])
plot2d(z,Anew,rect=[0,A0val,L,A0val+0.1])
xtitle('fonction z--> A(t,z)')
xsetech([0,0.5,1,0.5])
plot2d(z,Qnew,rect=[0,-10,L,50])
xtitle('fonction z--> Q(t,z)')
drawnow
//
end
///////////////////////////////////////////////////////////
//----------fin boucle principale ----------------------
/////////////////////////////////////////////////////////
//
/////////////////////  tracé final  //////////////////////
disp('temps de calcul');disp(timer());
xset('window',2)
clf()
plot2d(t,[prP',prM',prD'],[1,2,3])
xtitle('fonctions t--> p(t,z) en 3 points P, M et D')
xset('window',3)
clf()
xsetech([0,0,1,0.5])
plot2d(t,[AP',AM',AD'],[1,2,3])
xtitle('fonctions t--> A(t,z) en 3 points P, M et D')
xsetech([0,0.5,1,0.5])
plot2d(t,[QP',QM',QD'],[1,2,3])
xtitle('fonctions t--> Q(t,z) en 3 points P, M et D')

//
