/////////////////////////////////////////
//
////******VERSION  N arteres  *************
//        
/////////////////////////////////////////
// 
//  ici bifurcation:      Iliaque(1)--- >  extfem(2),intfem(3)
//
//
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//
clear all
//
function [F1,F2]=F(A,Q,A0,bet);
F1=Q;
F2=alpha*(Q.^2)./A+1/3/rho*bet.*(A.^1.5-A0.^1.5);
endfunction
///////////////////////////////////////////////////
function [B1,B2]=B(A,Q,A0,dbet);
B1=zeros(A);
B2=-Kr*Q./A+1/rho*dbet.*(sqrt(A0).*(A-A0)-2/3*(A.^1.5-A0.^1.5));
endfunction
///////////////////////////////////////////////////
function [HU11,HU12,HU21,HU22]=H(A,Q,bet);
HU11=zeros(A);
HU12=ones(A);
HU21=-alpha*(Q.^2)./(A^2)+1/2/rho*bet.*sqrt(A);
HU22=2*alpha*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 [BU1,BU2]=GU(A,Q,A0,dbet)
        BU1=zeros(A);
        BU2=-Kr*Q./A+A/rho*dbet.*(sqrt(A0)-2/3*sqrt(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 [BU11,BU12,BU21,BU22]=BU(A,Q,A0,dbet);
BU11=zeros(A);
BU12=zeros(A);
BU21=Kr*Q./(A.^2)+1/rho*dbet.*(sqrt(A0)-sqrt(A));
BU22=-Kr./A;
endfunction
///=////////////////////////////////////////////////
function [FLW1,FLW2]=FLW(A,Q,A0,bet,dbet);
[F1,F2]=F(A,Q,A0,bet);
[B1,B2]=B(A,Q,A0,dbet);
[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,A0,dbet);
[B1,B2]=B(A,Q,A0,dbet);
[BU11,BU12,BU21,BU22]=BU(A,Q,A0,dbet);
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
//
/////////////////////////////////////////////////
function [nP,nM,nD]=npt(z)                      // position des points P,M et D
  n=length(z);L=z(n);
  nP=1;
  nM=max(find(z<0.5*L));
  nD=n;
endfunction
///////////////////////////////////////////////////
//
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//
/// Systeme de 6 équations a resoudre pour les conditions aux bords au niveau de la bifurcation
//
function y = BifurcSystEq(X)
  
  Anplus1L_L = X(1) ; Qnplus1L_L = X(2) ; 
  Anplus10_R1 = X(3) ; Qnplus10_R1 = X(4) ;
  Anplus10_R2 = X(5) ; Qnplus10_R2 = X(6) ;
  
  u_L = Qnplus1L_L/Anplus1L_L; u_R1 = Qnplus10_R1/Anplus10_R1; u_R2 = Qnplus10_R2/Anplus10_R2; 
  
  R = zeros(X);
  
  R(1) = Qnplus1L_L - Qnplus10_R1 - Qnplus10_R2;
  
  R(2) = (1/2)*(u_L^2 - u_R1^2) + (bet(Conect(j,1)).entries(N(Conect(j,1))+2)/rho)*(sqrt(Anplus1L_L)-sqrt(A0(Conect(j,1)).entries(N(Conect(j,1))+2)))  - (bet(Conect(j,2)).entries(1)/rho)*(sqrt(Anplus10_R1)-sqrt(A0(Conect(j,2)).entries(1))) ;
  
  R(3) = (1/2)*(u_L^2 - u_R2^2) + (bet(Conect(j,1)).entries(N(Conect(j,1))+2)/rho)*(sqrt(Anplus1L_L)-sqrt(A0(Conect(j,1)).entries(N(Conect(j,1))+2)))  - (bet(Conect(j,3)).entries(1)/rho)*(sqrt(Anplus10_R2)-sqrt(A0(Conect(j,3)).entries(1)));
  
  R(4) =  u_L + 2*sqrt( 2/(rho)*bet(Conect(j,1)).entries(N(Conect(j,1))+2).*sqrt(Anplus1L_L)) - W1nplus1L_L;
  
  R(5) =  u_R1 - 2*sqrt( 2/(rho)*bet(Conect(j,2)).entries(1).*sqrt(Anplus10_R1)) - W2nplus10_R1;
  
  R(6) =  u_R2 - 2*sqrt( 2/(rho)*bet(Conect(j,3)).entries(1).*sqrt(Anplus10_R2)) - W2nplus10_R2;
  
  y = real(R);
endfunction
//
//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------
//           PROGRAMME  PRINCIPAL
//--------------------------------------------------------------------//------------------
//--------------------------------------------------------------------//------------------
//
//
timer();
//
///////////////////////////////////////////////////////////////////
//
// connectivité
//
NbrArt=3;
Conect = [1 2 3];
Entree = 1;
Sorties = [2  3];
//
// ////////////////////////////////////////////////////////////////////
//------------ Paramètres ------------
/////////////////////////////////////////////////////////////////////
//
//   section d'entree expérimentale
//
Data=read('artere-iliac.txt',-1,16);    // fichier de données
//
A0min1=1E-8*%pi*(Data(1,8)/2)^2;           // section initiale A0
////
///////////////////////////////////////////////////////////////////
//
// paramètres physiques
//
  Kr=0.75;;         // Contient la Viscosité du sang
  rho=1.05;         
  alpha=1;          
//
///////////////////////////////////////////////////////////////////
//
// longueur des arteres (en cm)
//
  Long=[14.4;45.3;33.3];
//
///////////////////////////////////////////////////////////
// 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=100*c;   // valeur en cm/s
//
///////////////////////////////////////////////////////////////////
//
// 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);
//
///////////////////////////////////////////////////////////////////
//
// sections au repos (en cm^2)
//  
A03min=0.2;
A03max=0.4;
A0val=[0.5327934;0.343;A03min+rand()*(A03max-A03min)];
//
//
///////////////////////////////////////////////////////////////////
//
// valeur de beta déduite de C0 et A0
// 
 bet_val=2*rho*(c.^2)./sqrt(A0val)  
//
///////////////////////////////////////////////////////////////////
//
// nombre de points de discrétisation par artère
//  
  Ndis=[8;13;11];
  //
////////////////////////////////////////////////////////////////////
//
// discrétisation en temps 
//   
  inti=1;      // temps initial
  intf=130;    // temps final
  dt=1E-3      // pas de temps
// 
/////////////////////////////////////////////////////////////////////
  disp('c0:');disp(c/100);
  disp('section A0');disp(A0val);
  disp('beta:');disp(bet_val);
  disp('resistance:');disp(Rt);
//
///////////////////////////////////////////////////////////
//
  // ------- discrétisation spatiale ----//
  //****************************//
 
 z = cell(NbrArt,1);
  for i=1:NbrArt 
  z(i).entries = linspace(0,Long(i),Ndis(i)+1)';
  end
// 
  N = zeros(NbrArt,1); L = zeros(NbrArt,1);
  for i=1:NbrArt
    N(i)=length(z(i).entries)-2; L(i) = z(i).entries($);
  end
  PMD = zeros(NbrArt,3); 
  for i=1:NbrArt
    [nP,nM,nD]=npt(z(i).entries);
    PMD(i,:) = [nP,nM,nD];
  end
////////////////////////////////////////////////////////////
  //
  //----- calcul de beta ---------------------
  //
  bet=cell(NbrArt,1);
  dbet=cell(NbrArt,1);      // dérivée de beta (nulle ici)
  for i=1:NbrArt   
    bet(i).entries = bet_val(i).*ones(N(i)+2,1); 
    dbet(i).entries = zeros(N(i)+2,1);     
  end
  //
  //----- calcul de A0 ---------------------
  //
  A0 = cell(NbrArt,1);
  for i=1:NbrArt
    A0(i).entries=A0val(i)*ones(N(i)+2,1) ;
  end 
//
/////////////////////////////////////////////////
//
 // Profil de section expérimental à l'entrée
 //------------------------------------
 //
 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);
 Nt=int(Tfin/dt+1E-12)
 t=(0:dt:Tfin);
 STOP = ceil(Nt);
 Abord0 = (Abord0/2)^2*%pi; // section à l'entree
//
/////////////////////////////////////////////////
//  
// initialisation à t=0: A = A0 et Q=0
//
  A = cell(NbrArt,1);
  Q = cell(NbrArt,1);
  for i=1:NbrArt
    A(i).entries=zeros(N(i)+2,1);
    Q(i).entries=zeros(N(i)+2,1);
  end 
//
  A = A0;
  for i=1:NbrArt
    Q(i).entries=zeros(N(i)+2,1);
  end
/////////////////////////////////////////////////
  
 // conditions aux limites initiales pour W1 et W2
  
  for i=1:NbrArt
    [W1init0,W2init0]=W12(A(i).entries(1),Q(i).entries(1),bet(i).entries(1));
    [W1initL,W2initL]=W12(A(i).entries(N(i)+2),Q(i).entries(N(i)+2),bet(i).entries(N(i)+2));
    W(i,:) = [W1init0,W2init0, W1initL,W2initL];
  end
  
 /////////////////////////////////////////////////
  
 // conditions à l'entrée pour W1 et à la sortie pour W2

     for i=1:NbrArt
         g1(i).entries=W(i,2)*ones(t) + 4*sqrt(2/rho)*sqrt(bet(i).entries(1)*sqrt(Abord0));
         g2(i).entries=W(i,4)*ones(t);
       end

  
  //-------- matrices pour le calcul EF ------------------------

  
   for i=1:NbrArt
      M(i).entries=Masse(z(i).entries);
      Mbis(i).entries=M(i).entries(2:N(i)+1,:);
      Minv(i).entries=inv(M(i).entries(2:N(i)+1,2:N(i)+1));
      G = Ggradient(ones(z(i).entries));
      Gbis(i).entries = G(2:N(i)+1,:);
   end
 //-------- initialisation  ------------------------
  
  for i=1:NbrArt
    p(i).entries=bet(i).entries.*(sqrt(A(i).entries)-sqrt(A0(i).entries));

    
    pr_PMD(i).entries=zeros(Nt+1,3);
    A_PMD(i).entries=zeros(Nt+1,3);
    
    A_PMD(i).entries(1,:)=[A0(i).entries(1) A0(i).entries(1) A0(i).entries(1)];
    Q_PMD(i).entries=zeros(Nt+1,3);
    
  end
  //-------------------------------------------------------
  //----------------debut boucle principale-------------
  // -----------------------------------------------------
//
Anew=A;
Qnew=Q;
//
  for i=1:STOP  

 
    ////// Conditions aux bords au niveau des bifurcations //////
  
  
  
    for j=1:size(Conect,1)  //
        
      [LAMB1_0_R1,LAMB2_0_R1]=LAMB12(A(Conect(j,2)).entries(1),Q(Conect(j,2)).entries(1),bet(Conect(j,2)).entries(1));
      [W1_R1,W2_R1]=W12(A(Conect(j,2)).entries,Q(Conect(j,2)).entries,bet(Conect(j,2)).entries);
      W2nplus10_R1=interpln([z(Conect(j,2)).entries';W2_R1'],-LAMB2_0_R1*dt);
  
      [LAMB1_0_R2,LAMB2_0_R2]=LAMB12(A(Conect(j,3)).entries(1),Q(Conect(j,3)).entries(1),bet(Conect(j,3)).entries(1));
      [W1_R2,W2_R2]=W12(A(Conect(j,3)).entries,Q(Conect(j,3)).entries,bet(Conect(j,3)).entries);
      W2nplus10_R2=interpln([z(Conect(j,3)).entries';W2_R2'],-LAMB2_0_R2*dt);
  
      [LAMB1_L_L,LAMB2_L_L] = LAMB12(A(Conect(j,1)).entries(N(Conect(j,1))+2),Q(Conect(j,1)).entries(N(Conect(j,1))+2),bet(Conect(j,1)).entries(N(Conect(j,1))+2));
      [W1_L,W2_L]=W12(A(Conect(j,1)).entries,Q(Conect(j,1)).entries,bet(Conect(j,1)).entries);
      W1nplus1L_L=interpln([z(Conect(j,1)).entries';W1_L'],L(Conect(j,1))-LAMB1_L_L*dt); 
      
  
      X0 = [A(Conect(j,1)).entries(N(Conect(j,1))+2) Q(Conect(j,1)).entries(N(Conect(j,1))+2) A(Conect(j,2)).entries(1) Q(Conect(j,2)).entries(1) A(Conect(j,3)).entries(1) Q(Conect(j,3)).entries(1)];
  

      [Solution VS] = fsolve(X0,BifurcSystEq);

      
      Anew(Conect(j,1)).entries(N(Conect(j,1))+2)=Solution(1);
      Qnew(Conect(j,1)).entries(N(Conect(j,1))+2)=Solution(2);
      Anew(Conect(j,2)).entries(1)=Solution(3);
      Qnew(Conect(j,2)).entries(1)=Solution(4);
      Anew(Conect(j,3)).entries(1)=Solution(5);
      Qnew(Conect(j,3)).entries(1)=Solution(6);
      
  
    end
    
    ////// Conditions aux bords au niveau de l'entree et des sorties //////
    
    // Entree du réseau //
    
    [LAMB1_0_L,LAMB2_0_L]=LAMB12(A(Entree).entries(1),Q(Entree).entries(1),bet(Entree).entries(1));
    [W1_L,W2_L]=W12(A(Entree).entries,Q(Entree).entries,bet(Entree).entries);
    
    
    W2nplus10_L=interpln([z(Entree).entries';W2_L'],-LAMB2_0_L*dt);
    W1nplus10_L=W2nplus10_L + 4*sqrt(2/rho)*sqrt(bet(1).entries(1)*sqrt(Abord0(i+1)));
    
    

    [Anplus10_L,Qnplus10_L]=invW12(W1nplus10_L,W2nplus10_L,bet(Entree).entries(1));
    
    Anew(Entree).entries(1)=Anplus10_L;
    Qnew(Entree).entries(1)=Qnplus10_L;
    
    // Sorties du réseau //
    

        for j=1:length(Sorties)
      
      [LAMB1_L_S,LAMB2_L_S] = LAMB12(A(Sorties(j)).entries(N(Sorties(j))+2),Q(Sorties(j)).entries(N(Sorties(j))+2),bet(Sorties(j)).entries(N(Sorties(j))+2));
      [W1_S,W2_S]=W12(A(Sorties(j)).entries,Q(Sorties(j)).entries,bet(Sorties(j)).entries);
      W1nplus1L_S=interpln([z(Sorties(j)).entries';W1_S'],L(Sorties(j))-LAMB1_L_S*dt);
      

      [l1,l2]=Lvect(A(Sorties(j)).entries(N(Sorties(j))+2),Q(Sorties(j)).entries(N(Sorties(j))+2), bet(Sorties(j)).entries(N(Sorties(j))+2))
      W2nplus1L_S=g2(Sorties(j)).entries(i+1) - Rt(j)*(W1nplus1L_S - W(Sorties(j),3)); 
      [Anplus1L_S,Qnplus1L_S]=invW12(W1nplus1L_S,W2nplus1L_S,bet(Sorties(j)).entries(N(Sorties(j))+2));
    
      Anew(Sorties(j)).entries(N(Sorties(j))+2)=Anplus1L_S;
      Qnew(Sorties(j)).entries(N(Sorties(j))+2)=Qnplus1L_S;
      
    end
    
    ////// Résolution du Pb élément fini pour chaque artre (Taylor-Galerkin) //////
     
  
    for j=1:NbrArt

      [F1,F2]=F(A(j).entries,Q(j).entries,A0(j).entries,bet(j).entries);
      [FLW1,FLW2]=FLW(A(j).entries,Q(j).entries,A0(j).entries,bet(j).entries,dbet(j).entries);
      [BLW1,BLW2]=BLW(A(j).entries,Q(j).entries,A0(j).entries,dbet(j).entries);
      [BU11,BU12,BU21,BU22]=BU(A(j).entries,Q(j).entries,A0(j).entries,dbet(j).entries);
      [HU11,HU12,HU21,HU22]=H(A(j).entries,Q(j).entries,bet(j).entries);

      //
      RHS1part1=Mbis(j).entries*A(j).entries+dt*Gbis(j).entries*FLW1+dt*Mbis(j).entries*BLW1;
      RHS2part1=Mbis(j).entries*Q(j).entries+dt*Gbis(j).entries*FLW2+dt*Mbis(j).entries*BLW2;
      RHS1part2=-dt^2/2*(Llaplacien(HU11,z(j).entries)*F1+Llaplacien(HU12,z(j).entries)*F2)-dt^2/2*(Ddiv(BU11)*F1+Ddiv(BU12)*F2);
      RHS2part2=-dt^2/2*(Llaplacien(HU21,z(j).entries)*F1+Llaplacien(HU22,z(j).entries)*F2)-dt^2/2*(Ddiv(BU21)*F1+Ddiv(BU22)*F2);
      //
      RHS1=RHS1part1+RHS1part2;
      RHS2=RHS2part1+RHS2part2;
      //
      RHS1(1)=RHS1(1)-Anew(j).entries(1)*M(j).entries(2,1);
      RHS1(N(j))=RHS1(N(j))-Anew(j).entries(N(j)+2)*M(j).entries(N(j)+1,N(j)+2);
      RHS2(1)=RHS2(1)-Qnew(j).entries(1)*M(j).entries(2,1);
      RHS2(N(j))=RHS2(N(j))-Qnew(j).entries(N(j)+2)*M(j).entries(N(j)+1,N(j)+2);
      //
      Aprime=Minv(j).entries*RHS1;
      Qprime=Minv(j).entries*RHS2;
      //
      Anew(j).entries(2:N(j)+1)=Aprime;
      Qnew(j).entries(2:N(j)+1)=Qprime;
      
      pnew(j).entries=bet(j).entries.*(sqrt(Anew(j).entries)-sqrt(A0(j).entries));

      
      A_PMD(j).entries(i+1,:)=[Anew(j).entries(PMD(j,1)), Anew(j).entries(PMD(j,2)), Anew(j).entries(PMD(j,3))];
      Q_PMD(j).entries(i+1,:)=[Qnew(j).entries(PMD(j,1)), Qnew(j).entries(PMD(j,2)), Qnew(j).entries(PMD(j,3))];
      pr_PMD(j).entries(i+1,:)=[pnew(j).entries(PMD(j,1)), pnew(j).entries(PMD(j,2)), pnew(j).entries(PMD(j,3))];
      
    end    
    
    
    A=Anew;
    Q=Qnew;
    p=pnew;
    
  

end  
//////////////////////////////////////////////////////////////////
//
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 2, temps
//
icouleur=5;
//
figure(0)
clf();
//
// Artère 1
//*********
subplot(2,3,1)
plot2d(t,A_PMD(1).entries(:,1),icouleur) 
plot(t,Abord0,'b--')
legend(['A-exp','A_Num'])
title('A1')
//
subplot(2,3,4)
plot2d(t,Q_PMD(1).entries(:,1)./A_PMD(1).entries(:,2),icouleur)
plot(Expe(:,4),Expe(:,1),'black--')
title('V1')
//
// Artère 2
//*********
subplot(2,3,2)
 plot2d(t,A_PMD(2).entries(:,2),icouleur)
 plot(Expe(:,4),Expe(:,2),'black--')
 title('A2')
//
subplot(2,3,5)
 plot2d(t,Q_PMD(2).entries(:,2)./A_PMD(2).entries(:,2),icouleur)
  plot(Expe(:,4),Expe(:,3),'black--')
 title('V2') 
//
// Artère 3
//*********
subplot(2,3,3)
plot2d(t,A_PMD(3).entries(:,2),icouleur) 
title('A3')
//
subplot(2,3,6)
plot2d(t,Q_PMD(3).entries(:,2)./A_PMD(3).entries(:,2),icouleur)
title('V3')
