function artere=principalN()
global dt Kr rho artere Conect j N W1nplus1L_L W2nplus10_R1 W2nplus10_R2
% connectivité

NbrArt=3;
Conect = [1 2 3];
Entree = 1;
Sorties = [2  3];
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------ Paramètres ------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
%
%   section d'entree expérimentale
%
Data=load('artere-illiac.txt');    % fichier de données
%
tic;
A0min1=1E-8*pi*(Data(1,8)/2)^2;           % section initiale A0

  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];
cmin=[5;5;5];
cmax=[8;8;8];
%c=cmin+rand(3,1).*(cmax-cmin);
%
c=[ 6.5;6.5;6.5];

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);
Rt=[ 0.5500;
    0.4500
   ];

% sections au repos (en cm^2)
%  
A03min=0.2;
A03max=0.4;
A0val=[0.5327934;0.343; 0.2481];%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=[5;7;6];
% discrétisation en temps 
%   
  inti=1;      % temps initial
  intf=60;    % temps final
  dt=2E-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 ----%
  %****************************%
 
 
  for i=1:NbrArt 
  artere(i).z= linspace(0,Long(i),Ndis(i)+1)';
  end
% 
  N = zeros(NbrArt,1); L = zeros(NbrArt,1);
  for i=1:NbrArt
    N(i)=length(artere(i).z)-2; L(i) = artere(i).z(end);
  end
  PMD = zeros(NbrArt,3); 
  for i=1:NbrArt
    [nP,nM,nD]=npt(artere(i).z);
    PMD(i,:) = [nP,nM,nD];
  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %----- calcul de beta ---------------------
  %
      
  for i=1:NbrArt   
    artere(i).bet = bet_val(i).*ones(N(i)+2,1); 
    artere(i).dbet = zeros(N(i)+2,1);     
  end

  %----- calcul de A0 ---------------------
  %
  
  for i=1:NbrArt
   artere(i).A0=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=floor(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
 
  for i=1:NbrArt
    artere(i).A=zeros(N(i)+2,1);
    artere(i).Q=zeros(N(i)+2,1);
  end 
for i=1:NbrArt
     artere(i).A = artere(i).A0
  end

  for i=1:NbrArt
    artere(i).Q=zeros(N(i)+2,1);
  end
%%%%%%%%%%%%%%%%%%%%%%%% 
 % conditions aux limites initiales pour W1 et W2
  
  for i=1:NbrArt
    [W1init0,W2init0]=W12(artere(i).A(1),artere(i).Q(1),artere(i).bet(1));
    [W1initL,W2initL]=W12(artere(i).A(N(i)+2),artere(i).Q(N(i)+2),artere(i).bet(N(i)+2));
    W(i,:) = [W1init0,W2init0, W1initL,W2initL];
  end
  
 %%%%%%%%%%%%%%%%%%%%%%%%/
  
   % conditions à l'entrée pour W1 et à la sortie pour W2
[pp qq]=size(t);
     for i=1:NbrArt
         artere(i).g1=W(i,2)*ones(pp,qq) + 4*sqrt(2/rho)*sqrt(artere(i).bet(1)*sqrt(Abord0));
         artere(i).g2=W(i,4)*ones(pp,qq);
       end
  %-------- matrices pour le calcul EF ------------------------

  
   for i=1:NbrArt
      artere(i).M=Masse(artere(i).z);
      artere(i).Mbis=artere(i).M(2:N(i)+1,:);
      artere(i).Minv=inv(artere(i).M(2:N(i)+1,2:N(i)+1));
      [pp qq]=size(artere(i).z);
      G = Ggradient(ones(pp,qq));
      artere(i).Gbis = G(2:N(i)+1,:);
   end
   %-------- initialisation  ------------------------
  
  for i=1:NbrArt
    artere(i).p=artere(i).bet.*(sqrt(artere(i).A)-sqrt(artere(i).A0));

    
    artere(i).pr_PMD=zeros(Nt+1,3);
    artere(i).A_PMD=zeros(Nt+1,3);
    
    artere(i).A_PMD(1,:)=[artere(i).A0(1) artere(i).A0(1) artere(i).A0(1)];
    artere(i).Q_PMD=zeros(Nt+1,3);
    
  end
   %-------------------------------------------------------
  %----------------debut boucle principale-------------
  % -----------------------------------------------------
for k=1:NbrArt
artere(k).Anew=artere(k).A;
artere(k).Qnew=artere(k).Q;
end
  for i=1:STOP  
   %%% Conditions aux bords au niveau des bifurcations %%%
  
  
  
    for j=1:size(Conect,1)  %
      [LAMB1_0_R1,LAMB2_0_R1]=LAMB12(artere(Conect(j,2)).A(1),artere(Conect(j,2)).Q(1),artere(Conect(j,2)).bet(1));
      [W1_R1,W2_R1]=W12(artere(Conect(j,2)).A,artere(Conect(j,2)).Q,artere(Conect(j,2)).bet);
      W2nplus10_R1=interp1(artere(Conect(j,2)).z',W2_R1',-LAMB2_0_R1*dt);
  
      [LAMB1_0_R2,LAMB2_0_R2]=LAMB12(artere(Conect(j,3)).A(1),artere(Conect(j,3)).Q(1),artere(Conect(j,3)).bet(1));
      [W1_R2,W2_R2]=W12(artere(Conect(j,3)).A,artere(Conect(j,3)).Q,artere(Conect(j,3)).bet);
      W2nplus10_R2=interp1(artere(Conect(j,3)).z',W2_R2',-LAMB2_0_R2*dt);
  
      [LAMB1_L_L,LAMB2_L_L] = LAMB12(artere(Conect(j,1)).A(N(Conect(j,1))+2),artere(Conect(j,1)).Q(N(Conect(j,1))+2),artere(Conect(j,1)).bet(N(Conect(j,1))+2));
      [W1_L,W2_L]=W12(artere(Conect(j,1)).A,artere(Conect(j,1)).Q,artere(Conect(j,1)).bet);
      W1nplus1L_L=interp1(artere(Conect(j,1)).z',W1_L',L(Conect(j,1))-LAMB1_L_L*dt); 
    
  
      X0 = [artere(Conect(j,1)).A(N(Conect(j,1))+2) artere(Conect(j,1)).Q(N(Conect(j,1))+2) artere(Conect(j,2)).A(1) artere(Conect(j,2)).Q(1) artere(Conect(j,3)).A(1) artere(Conect(j,3)).Q(1)];
      

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

      
      artere(Conect(j,1)).Anew(N(Conect(j,1))+2)=Solution(1);
      artere(Conect(j,1)).Qnew(N(Conect(j,1))+2)=Solution(2);
      artere(Conect(j,2)).Anew(1)=Solution(3);
      artere(Conect(j,2)).Qnew(1)=Solution(4);
      artere(Conect(j,3)).Anew(1)=Solution(5);
      artere(Conect(j,3)).Qnew(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(artere(Entree).A(1),artere(Entree).Q(1),artere(Entree).bet(1));
    [W1_L,W2_L]=W12(artere(Entree).A,artere(Entree).Q,artere(Entree).bet);
    
    
    W2nplus10_L=interp1(artere(Entree).z',W2_L',-LAMB2_0_L*dt);
    W1nplus10_L=W2nplus10_L + 4*sqrt(2/rho)*sqrt(artere(1).bet(1)*sqrt(Abord0(i+1)));
    
    

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

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

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

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

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

      
      artere(jj).A_PMD(i+1,:)=[artere(jj).Anew(PMD(jj,1)), artere(jj).Anew(PMD(jj,2)), artere(jj).Anew(PMD(jj,3))];
      artere(jj).Q_PMD(i+1,:)=[artere(jj).Qnew(PMD(jj,1)), artere(jj).Qnew(PMD(jj,2)), artere(jj).Qnew(PMD(jj,3))];
      artere(jj).pr_PMD(i+1,:)=[artere(jj).pnew(PMD(jj,1)), artere(jj).pnew(PMD(jj,2)), artere(jj).pnew(PMD(jj,3))];
      
    end    
    
    for kk=1:NbrArt
    artere(kk).A=artere(kk).Anew;
    artere(kk).Q=artere(kk).Qnew;
    artere(kk).p=artere(kk).pnew;
    end
  

end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%**************************************************
% TRACE DES RESULTATS
disp 'temps de calcul'
disp(toc)
Expe=load('artere-exp2.txt');    % fichier expérimental 4 colonnes:
% vitesse artere 1, section artere 2, vitesse artere 2, temps


% Artère 1

subplot(2,3,1)
plot(t,artere(1).A_PMD(:,1)) 
hold on
plot(t,Abord0,'b--')
legend(['A-exp','A_Num'])
title('A1')


subplot(2,3,4)
plot(t,artere(1).Q_PMD(:,1)./artere(1).A_PMD(:,2))
hold on
plot(Expe(:,4),Expe(:,1),'black--')
title('V1')

% Artère 2

subplot(2,3,2)
 plot(t,artere(2).A_PMD(:,2))
 hold on
 plot(Expe(:,4),Expe(:,2),'black--')
 title('A2')

subplot(2,3,5)
 plot(t,artere(2).Q_PMD(:,2)./artere(2).A_PMD(:,2))
 hold on
  plot(Expe(:,4),Expe(:,3),'black--')
 title('V2') 
% Artère 3

subplot(2,3,3)
plot(t,artere(3).A_PMD(:,2)) 
title('A3')

subplot(2,3,6)
plot(t,artere(3).Q_PMD(:,2)./artere(3).A_PMD(:,2))
title('V3')

 
  
end

function [F1,F2]=F(A,Q,A0,bet)
global rho
F1=Q;
F2=(Q.^2)./A+1/3/rho*bet.*(A.^1.5-A0.^1.5);
end

function [B1,B2]=B(A,Q,A0,dbet);
global Kr rho
[aa,bb]=size(A);
B1=zeros(aa,bb);
B2=-Kr*Q./A+1/rho*dbet.*(sqrt(A0).*(A-A0)-2/3*(A.^1.5-A0.^1.5));
end

function [HU11,HU12,HU21,HU22]=H(A,Q,bet);
global rho
[aa,bb]=size(A);
HU11=zeros(aa,bb);
HU12=ones(aa,bb);

HU21=-(Q.^2)./(A.^2)+1/2/rho*bet.*sqrt(A);      % loooooooooooooooooool puissance
HU22=2*Q./A;
end

function [W1,W2]=W12(A,Q,bet);
global rho

W1=Q./A+2*sqrt(2/(rho)*bet.*sqrt(A));
W2=Q./A-2*sqrt(2/(rho)*bet.*sqrt(A));
end
function [A,Q]=invW12(W1,W2,bet);
global rho
A=(rho./bet).^2*(W1-W2).^4/4^5;
Q=A.*(W1+W2)/2;
end

function [LAMB1,LAMB2]=LAMB12(A,Q,bet);
global rho
LAMB1=Q./A+sqrt(1/(2*rho)*bet.*sqrt(A));
LAMB2=Q./A-sqrt(1/(2*rho)*bet.*sqrt(A));
end


%%%%%%%%%%%%%%%%%%%%%%%%%/
function [BU11,BU12,BU21,BU22]=BU(A,Q,A0,dbet);
global Kr rho
[aa,bb]=size(A);
BU11=zeros(aa,bb);
BU12=zeros(aa,bb);
BU21=Kr*Q./(A.^2)+1/rho*dbet.*(sqrt(A0)-sqrt(A));
BU22=-Kr./A;
end

%Left eigenvalue
function [l1,l2]=Lvect(A,Q,bet);     
[LAMB1,LAMB2]=LAMB12(A,Q,bet);
l1= [-LAMB2,1];
l2=[-LAMB1,1];
end
function [FLW1,FLW2]=FLW(A,Q,A0,bet,dbet);
global dt
[F1,F2]=F(A,Q,A0,bet);
[B1,B2]=B(A,Q,A0,dbet);
[HU11,HU12,HU21,HU22]=H(A,Q,bet);
FLW1=F1+dt/2*(HU11.*B1+HU12.*B2);
FLW2=F2+dt/2*(HU21.*B1+HU22.*B2);
end
function [BLW1,BLW2]=BLW(A,Q,A0,dbet);
global dt
[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);
end

%%-----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;
end

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;
end

function M=Ddiv(f)
n=length(f);
M=Ggradient(f)';
M=M(2:n-1,:);
end

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,:)
end
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;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function y = BifurcSystEq(X)
  global rho artere Conect j N W1nplus1L_L W2nplus10_R1 W2nplus10_R2
  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; 
  [pp qq]=size(X);
  R = zeros(pp,qq);
  
  R(1) = Qnplus1L_L - Qnplus10_R1 - Qnplus10_R2;
  
  R(2) = (1/2)*(u_L^2 - u_R1^2) + (artere(Conect(j,1)).bet(N(Conect(j,1))+2)/rho)*(sqrt(Anplus1L_L)-sqrt(artere(Conect(j,1)).A0(N(Conect(j,1))+2)))  - (artere(Conect(j,2)).bet(1)/rho)*(sqrt(Anplus10_R1)-sqrt(artere(Conect(j,2)).A0(1))) ;
  
  R(3) = (1/2)*(u_L^2 - u_R2^2) + (artere(Conect(j,1)).bet(N(Conect(j,1))+2)/rho)*(sqrt(Anplus1L_L)-sqrt(artere(Conect(j,1)).A0(N(Conect(j,1))+2)))  - (artere(Conect(j,3)).bet(1)/rho)*(sqrt(Anplus10_R2)-sqrt(artere(Conect(j,3)).A0(1)));
  
  R(4) =  u_L + 2*sqrt( 2/(rho)*artere(Conect(j,1)).bet(N(Conect(j,1))+2).*sqrt(Anplus1L_L)) - W1nplus1L_L;
  
  R(5) =  u_R1 - 2*sqrt( 2/(rho)*artere(Conect(j,2)).bet(1).*sqrt(Anplus10_R1)) - W2nplus10_R1;
  
  R(6) =  u_R2 - 2*sqrt( 2/(rho)*artere(Conect(j,3)).bet(1).*sqrt(Anplus10_R2)) - W2nplus10_R2;
  
  y = real(R);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


