//
//  Propagation d'une onde de pression dans un réseau de trois arteres saine
//  avec un modele fluide/structure 1D
//  ici bifurcation:      Iliaque(1)--- >  extfem(2),intfem(3)
//          --------
//  --------  ->
//    ->    --------
//  --------  ->
//          --------
//  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----------
//
clear;
//
function [F1,F2]=F(A,Q,A0,bet);
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,bet);
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,A0,bet,dt);
[F1,F2]=F(A,Q,A0,bet);
[B1,B2]=B(A,Q);
[HU11,HU12,HU21,HU22]=H(A,Q,bet);
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
// conditions aux bifurcation:
// A1nplus1L,Q1nplus1L en z = L1
// A2nplus10,Q2nplus10 en z = 0
// A3nplus10,Q3nplus10 en z = 0
// 3 équations physiques:
// Q1nplus1L = Q2nplus1L + Q3nplus1L
// 0.5*U1nplus1L^2 + bet1(sqrt(A1nplus1L) - sqrt(A10)) = 0.5*U2nplus10^2 + bet2(sqrt(A2nplus10) - sqrt(A20))
// 0.5*U1nplus1L^2 + bet1(sqrt(A1nplus1L) - sqrt(A10)) = 0.5*U3nplus10^2 + bet3(sqrt(A3nplus10) - sqrt(A30))
// 3 équations numériques:
// W1nplus1L_1 = U1nplus1L + 2*sqrt(2*bet1*sqrt(A1nplus1L)/rho))
// W2nplus10_2 = U2nplus10 + 2*sqrt(2*bet2*sqrt(A2nplus10)/rho))
// W2nplus10_3 = U3nplus10 + 2*sqrt(2*bet3*sqrt(A3nplus10)/rho))
function y=BifSys(X)
    
    R = zeros(length(X));

    R(1) = X(2)-X(4)-X(6);
    R(2) = 0.5*((X(2)/X(1))^2-(X(4)/X(3))^2)+bet(1)*(sqrt(X(1)) - sqrt(A10(n1+2))) - bet(2)*(sqrt(X(3)) - sqrt(A20(1)));
    R(3) = 0.5*((X(2)/X(1))^2-(X(6)/X(5))^2)+bet(1)*(sqrt(X(1)) - sqrt(A10(n1+2))) - bet(3)*(sqrt(X(5)) - sqrt(A30(1)));
    R(4) = W1nplus1L_1 - (X(2)/X(1)) - 2*sqrt(2*bet(1)*sqrt(X(1))/rho);
    R(5) = W2nplus10_2 - (X(4)/X(3)) + 2*sqrt(2*bet(2)*sqrt(X(3))/rho);
    R(6) = W2nplus10_3 - (X(6)/X(5)) + 2*sqrt(2*bet(3)*sqrt(X(5))/rho);
    
    y = real(R);
endfunction
/////////////////////////////////////////////////////////
//-------------------------------------------------------------
// ---------------Programme principal -------------------------
//-------------------------------------------------------------
///////////////////////////////////////////////////////////
//
timer();
/// caracteristiques du sang ////////////////////
//
rho=1.05; // densité du sang
Kr=0.75;  // coefficient de viscosité

//
///////////////////////////////////////////////////////////
// valeurs de c0 (en m/s): vitesse de l'onde de pouls
//
//cmin=[5;5;5];
//cmax=[8;8;8];
//c=cmin+rand(3,1).*(cmax-cmin);
c =[6.5;6.5;6.5];
//load('Pg.sod');
//c = Pg(1:3);
//
c=100*c;   // valeur en cm/s

// caractéristiques des 3 artères
L=[14.4;45.3;33.3]; // longeurs
//A03min=0.2;
//A03max=0.4;
//A0val=[0.5327934;0.343;A03min+rand()*(A03max-A03min)]; // section au repos
A0val=[0.5327934;0.343;0.3];
//A0val=[0.5327934;0.343;Pg(6)];

bet=2*rho*(c.^2)./sqrt(A0val); // beta déduit de A0 et c

///////////////////////////////////////////////////////////////////
//
// résistances en sortie (valeur entre 0 et 1)
//
//Rmin=[0.3;0.3];
//Rmax=[0.7;0.7];
//Rt=Rmin+rand(2,1).*(Rmax-Rmin);
Rt=[0.5;0.5];
//Rt=Pg(4:5);

///////////////////////////////////////////////////////////////////
//
// discrétisation spatiale et temporelle //////////////
//
N1=5;                      // nombre de points intérieurs dans 1'artère
N2=7;
N3=6;

inti=1;     // temps initial
intf=60;    // temps final 
dt=2E-3;    // pas de temps

disp('c:');disp(c/100);
disp('section A0');disp(A0val);
disp('beta:');disp(bet);
disp('resistance:');disp(Rt);

//
z1=linspace(0,L(1),N1+1);
z2=linspace(0,L(2),N2+1);
z3=linspace(0,L(3),N3+1);
 
n1=length(z1)-2;l1=z1($);
n2=length(z2)-2;l2=z2($);
n3=length(z3)-2;l3=z3($);

//
//n1P=max(find(z1<0.25*L(1)));
n1P=1;
n1M=max(find(z1<0.5*L(1)));
n1D=length(z1);
//n1D=max(find(z1<0.75*L(1)));
//

//
//n2P=max(find(z2<0.25*L(2)));
n2P=1;
n2M=max(find(z2<0.5*L(2)));
n2D=length(z2);
//n2D=max(find(z2<0.75*L(2)));

//
//n3P=max(find(z3<0.25*L(3)));
n3P=1;
n3M=max(find(z3<0.5*L(3)));
n3D=length(z3);
//n3D=max(find(z3<0.75*L(3)));
//
A1=zeros(n1+2,1);
Q1=zeros(n1+2,1);
A10=A0val(1)*ones(z1');
//
A2=zeros(n2+2,1);
Q2=zeros(n2+2,1);
A20=A0val(2)*ones(z2');
//
A3=zeros(n3+2,1);
Q3=zeros(n3+2,1);
A30=A0val(3)*ones(z3');
//
////// données expérimentales (section) /////////////
//
Data=read('artere-iliac.txt',-1,15); 
diam = Data(inti:intf,8)./10000;
Abord0 = interp1(Data(inti:intf,15),diam',(Data(inti,15):dt:Data(intf,15)),'linear');
Tfin=Data(intf,15)-Data(inti,15);
Abord0 = (Abord0/2).^2*%pi;

Nt=int(Tfin/dt+1E-12);
t=(0:dt:Tfin);

//
A1init=A10;
Q1init=zeros(n1+2,1);
//
A2init=A20;
Q2init=zeros(n2+2,1);
//
A3init=A30;
Q3init=zeros(n3+2,1);
//
[W1init0_1,W2init0_1]=W12(A1init(1),Q1init(1),bet(1));
[W1initL_1,W2initL_1]=W12(A1init(n1+2),Q1init(n1+2),bet(1));
W1nplus1L_1=W1initL_1;
W2nplus1L_1=W2initL_1;
//
[W1init0_2,W2init0_2]=W12(A2init(1),Q2init(1),bet(2));
[W1initL_2,W2initL_2]=W12(A2init(n2+2),Q2init(n2+2),bet(2));
W1nplus1L_2=W1initL_2;
W2nplus1L_2=W2initL_2;
//
[W1init0_3,W2init0_3]=W12(A3init(1),Q3init(1),bet(3));
[W1initL_3,W2initL_3]=W12(A3init(n3+2),Q3init(n3+2),bet(3));
W1nplus1L_3=W1initL_3;
W2nplus1L_3=W2initL_3;
//
//g1=W2init0*ones(t)+4*sqrt(2/rho)*sqrt(bet*sqrt(Abord0));
g2_1=W2initL_1*ones(t);
g2_2=W2initL_2*ones(t);
g2_3=W2initL_3*ones(t);
//
M1=Masse(z1);
M1bis=M1(2:n1+1,:);
M1inv=inv(M1(2:n1+1,2:n1+1));
G1=Ggradient(ones(z1));
G1bis=G1(2:n1+1,:);
//
M2=Masse(z2);
M2bis=M2(2:n2+1,:);
M2inv=inv(M2(2:n2+1,2:n2+1));
G2=Ggradient(ones(z2));
G2bis=G2(2:n2+1,:);
//
M3=Masse(z3);
M3bis=M3(2:n3+1,:);
M3inv=inv(M3(2:n3+1,2:n3+1));
G3=Ggradient(ones(z3));
G3bis=G3(2:n3+1,:);
//

A1=A1init;
Q1=Q1init;
A1new=A1init;
Q1new=Q1init;
p1new=A1new;
pr1P=[0];pr1M=[0];pr1D=[0];A1P=[A0val(1)];A1M=[A0val(1)];A1D=[A0val(1)];Q1P=[0];Q1M=[0];Q1D=[0];
A1tot=A1;Q1tot=Q1;

A2=A2init;
Q2=Q2init;
A2new=A2init;
Q2new=Q2init;
p2new=A2new;
pr2P=[0];pr2M=[0];pr2D=[0];A2P=[A0val(2)];A2M=[A0val(2)];A2D=[A0val(2)];Q2P=[0];Q2M=[0];Q2D=[0];
A2tot=A2;Q2tot=Q2;

A3=A3init;
Q3=Q3init;
A3new=A3init;
Q3new=Q3init;
p3new=A3new;
pr3P=[0];pr3M=[0];pr3D=[0];A3P=[A0val(3)];A3M=[A0val(3)];A3D=[A0val(3)];Q3P=[0];Q3M=[0];Q3D=[0];
A3tot=A3;Q3tot=Q3;


/////////////////////////////////////////////////////////
//----------------debut boucle principale-------------
/////////////////////////////////////////////////////////
//
for i=1:Nt
//
tpsnew=i*dt;
A1=A1new;
Q1=Q1new;
p1=p1new;

A2=A2new;
Q2=Q2new;
p2=p2new;

A3=A3new;
Q3=Q3new;
p3=p3new;
//

// traitement de la bifurcation
[LAMB1_1,LAMB2_1]=LAMB12(A1(n1+2),Q1(n1+2),bet(1));
[W1_1,W2_1]=W12(A1,Q1,bet(1));

W1nplus1L_1=interpln([z1;W1_1'],l1-LAMB1_1*dt);

[LAMB1_2,LAMB2_2]=LAMB12(A2(1),Q2(1),bet(2));
[W1_2,W2_2]=W12(A2,Q2,bet(2));

W2nplus10_2=interpln([z2;W2_2'],-LAMB2_2*dt);

[LAMB1_3,LAMB2_3]=LAMB12(A3(1),Q3(1),bet(3));
[W1_3,W2_3]=W12(A3,Q3,bet(3));

W2nplus10_3=interpln([z3;W2_3'],-LAMB2_3*dt);

X0 = [A1(n1+2),Q1(n1+2),A2(1),Q2(1),A3(1),Q3(1)];
[X] = fsolve(X0,BifSys);

A1nplus1L = X(1);
Q1nplus1L = X(2);
A2nplus10 = X(3);
Q2nplus10 = X(4);
A3nplus10 = X(5);
Q3nplus10 = X(6);


// artère 1 
[LAMB1_1,LAMB2_1]=LAMB12(A1(1),Q1(1),bet(1));
[W1_1,W2_1]=W12(A1,Q1,bet(1));

W2nplus10_1=interpln([z1;W2_1'],-LAMB2_1*dt);
W1nplus10_1=W2nplus10_1 + 4*sqrt(2/rho)*sqrt(bet(1)*sqrt(Abord0(i+1)));

[A1nplus10,Q1nplus10]=invW12(W1nplus10_1,W2nplus10_1,bet(1));

// artère 2
[LAMB1_2,LAMB2_2]=LAMB12(A2(n2+2),Q2(n2+2),bet(2));

W1nplus1L_2=interpln([z2;W1_2'],l2-LAMB1_2*dt);
W2nplus1L_2=g2_2(i+1) - Rt(1)*(W1nplus1L_2 - W1initL_2);

[A2nplus1L,Q2nplus1L]=invW12(W1nplus1L_2,W2nplus1L_2,bet(2));

// artère 3
[LAMB1_3,LAMB2_3]=LAMB12(A2(n2+2),Q2(n2+2),bet(3));

W1nplus1L_3=interpln([z3;W1_3'],l3-LAMB1_3*dt);
W2nplus1L_3=g2_3(i+1) - Rt(2)*(W1nplus1L_3 - W1initL_3);

[A3nplus1L,Q3nplus1L]=invW12(W1nplus1L_3,W2nplus1L_3,bet(3));

// résolution:

// artère 1
A1new(1)=A1nplus10;
Q1new(1)=Q1nplus10;
A1new(n1+2)=A1nplus1L;
Q1new(n1+2)=Q1nplus1L;
//
[F1,F2]=F(A1,Q1,A10,bet(1));
[FLW1,FLW2]=FLW(A1,Q1,A10,bet(1),dt);
[BLW1,BLW2]=BLW(A1,Q1,dt);
[BU11,BU12,BU21,BU22]=BU(A1,Q1);
[HU11,HU12,HU21,HU22]=H(A1,Q1,bet(1));
//
RHS1part1=M1bis*A1+dt*G1bis*FLW1+dt*M1bis*BLW1;
RHS2part1=M1bis*Q1+dt*G1bis*FLW2+dt*M1bis*BLW2;
RHS1part2=-dt^2/2*(Llaplacien(HU11,z1)*F1+Llaplacien(HU12,z1)*F2)-dt^2/2*(Ddiv(BU11)*F1+Ddiv(BU12)*F2);
RHS2part2=-dt^2/2*(Llaplacien(HU21,z1)*F1+Llaplacien(HU22,z1)*F2)-dt^2/2*(Ddiv(BU21)*F1+Ddiv(BU22)*F2);//
//
RHS1=RHS1part1+RHS1part2;
RHS2=RHS2part1+RHS2part2;
//
  RHS1(1) =RHS1(1) -A1new(1)*M1(2,1);
  RHS1(n1)=RHS1(n1)-A1new(n1+2)*M1(n1+1,n1+2);
  RHS2(1) =RHS2(1) -Q1new(1)*M1(2,1);
  RHS2(n1)=RHS2(n1)-Q1new(n1+2)*M1(n1+1,n1+2);
//
A1prime=M1inv*RHS1;
Q1prime=M1inv*RHS2;
//
A1new(2:n1+1)=A1prime;
Q1new(2:n1+1)=Q1prime;
p1new=bet(1).*(sqrt(A1new)-sqrt(A10));
A1P=[A1P,A1new(n1P)];Q1P=[Q1P,Q1new(n1P)];pr1P=[pr1P,p1new(n1P)];
A1M=[A1M,A1new(n1M)];Q1M=[Q1M,Q1new(n1M)];pr1M=[pr1M,p1new(n1M)];
A1D=[A1D,A1new(n1D)];Q1D=[Q1D,Q1new(n1D)];pr1D=[pr1D,p1new(n1D)];

// artère 2
A2new(1)=A2nplus10;
Q2new(1)=Q2nplus10;
A2new(n2+2)=A2nplus1L;
Q2new(n2+2)=Q2nplus1L;
//
[F1,F2]=F(A2,Q2,A20,bet(2));
[FLW1,FLW2]=FLW(A2,Q2,A20,bet(2),dt);
[BLW1,BLW2]=BLW(A2,Q2,dt);
[BU11,BU12,BU21,BU22]=BU(A2,Q2);
[HU11,HU12,HU21,HU22]=H(A2,Q2,bet(2));
//
RHS1part1=M2bis*A2+dt*G2bis*FLW1+dt*M2bis*BLW1;
RHS2part1=M2bis*Q2+dt*G2bis*FLW2+dt*M2bis*BLW2;
RHS1part2=-dt^2/2*(Llaplacien(HU11,z2)*F1+Llaplacien(HU12,z2)*F2)-dt^2/2*(Ddiv(BU11)*F1+Ddiv(BU12)*F2);
RHS2part2=-dt^2/2*(Llaplacien(HU21,z2)*F1+Llaplacien(HU22,z2)*F2)-dt^2/2*(Ddiv(BU21)*F1+Ddiv(BU22)*F2);//
//
RHS1=RHS1part1+RHS1part2;
RHS2=RHS2part1+RHS2part2;
//
  RHS1(1) =RHS1(1) -A2new(1)*M2(2,1);
  RHS1(n2)=RHS1(n2)-A2new(n2+2)*M2(n2+1,n2+2);
  RHS2(1) =RHS2(1) -Q2new(1)*M2(2,1);
  RHS2(n2)=RHS2(n2)-Q2new(n2+2)*M2(n2+1,n2+2);
//
A2prime=M2inv*RHS1;
Q2prime=M2inv*RHS2;
//
A2new(2:n2+1)=A2prime;
Q2new(2:n2+1)=Q2prime;
p2new=bet(2).*(sqrt(A2new)-sqrt(A20));
A2P=[A2P,A2new(n2P)];Q2P=[Q2P,Q2new(n2P)];pr2P=[pr2P,p2new(n2P)];
A2M=[A2M,A2new(n2M)];Q2M=[Q2M,Q2new(n2M)];pr2M=[pr2M,p2new(n2M)];
A2D=[A2D,A2new(n2D)];Q2D=[Q2D,Q2new(n2D)];pr2D=[pr2D,p2new(n2D)];

// artère 3
A3new(1)=A3nplus10;
Q3new(1)=Q3nplus10;
A3new(n3+2)=A3nplus1L;
Q3new(n3+2)=Q3nplus1L;
//
[F1,F2]=F(A3,Q3,A30,bet(3));
[FLW1,FLW2]=FLW(A3,Q3,A30,bet(3),dt);
[BLW1,BLW2]=BLW(A3,Q3,dt);
[BU11,BU12,BU21,BU22]=BU(A3,Q3);
[HU11,HU12,HU21,HU22]=H(A3,Q3,bet(3));
//
RHS1part1=M3bis*A3+dt*G3bis*FLW1+dt*M3bis*BLW1;
RHS2part1=M3bis*Q3+dt*G3bis*FLW2+dt*M3bis*BLW2;
RHS1part2=-dt^2/2*(Llaplacien(HU11,z3)*F1+Llaplacien(HU12,z3)*F2)-dt^2/2*(Ddiv(BU11)*F1+Ddiv(BU12)*F2);
RHS2part2=-dt^2/2*(Llaplacien(HU21,z3)*F1+Llaplacien(HU22,z3)*F2)-dt^2/2*(Ddiv(BU21)*F1+Ddiv(BU22)*F2);//
//
RHS1=RHS1part1+RHS1part2;
RHS2=RHS2part1+RHS2part2;
//
  RHS1(1) =RHS1(1) -A3new(1)*M3(2,1);
  RHS1(n3)=RHS1(n3)-A3new(n3+2)*M3(n3+1,n3+2);
  RHS2(1) =RHS2(1) -Q3new(1)*M3(2,1);
  RHS2(n3)=RHS2(n3)-Q3new(n3+2)*M3(n3+1,n3+2);
//
A3prime=M3inv*RHS1;
Q3prime=M3inv*RHS2;
//
A3new(2:n3+1)=A3prime;
Q3new(2:n3+1)=Q3prime;
p3new=bet(3).*(sqrt(A3new)-sqrt(A30));
A3P=[A3P,A3new(n3P)];Q3P=[Q3P,Q3new(n3P)];pr3P=[pr3P,p3new(n3P)];
A3M=[A3M,A3new(n3M)];Q3M=[Q3M,Q3new(n3M)];pr3M=[pr3M,p3new(n3M)];
A3D=[A3D,A3new(n3D)];Q3D=[Q3D,Q3new(n3D)];pr3D=[pr3D,p3new(n3D)];

end
///////////////////////////////////////////////////////////
//----------fin boucle principale ----------------------
/////////////////////////////////////////////////////////
//
/////////////////////  tracé final  //////////////////////
disp('temps de calcul');disp(timer());

//
////**************************************************
//// TRACE DES RESULTATS
////**************************************************
//
Expe=read('artere-exp2.txt',-1,4);    // fichier expérimental 4 colonnes:
// vitesse artere 1, section artere 2, vitesse artere 3, temps
//
icouleur=5;
//
figure(0)
clf();
//
// Artère 1
//*********
subplot(2,3,1)
plot2d(t,A1P',icouleur) 
plot(t,Abord0,'b--')
legend(['A_Num','A-exp'])
title('A1')
//
subplot(2,3,4)
plot2d(t,(Q1P./A1P)',icouleur)
plot(Expe(:,4),Expe(:,1),'b--')
title('V1')
//
// Artère 2
//*********
subplot(2,3,2)
plot2d(t,A2M',icouleur)
plot(Expe(:,4),Expe(:,2),'b--')
title('A2')
//
subplot(2,3,5)
plot2d(t,(Q2M./A2M)',icouleur)
plot(Expe(:,4),Expe(:,3),'b--')

title('V2') 
//
// Artère 3
//*********
subplot(2,3,3)
plot2d(t,A3M',icouleur) 
title('A3')
//
subplot(2,3,6)
plot2d(t,(Q3M./A3M)',icouleur)
title('V3')
