function z=fctz(L,N,a1,a2,delta,tauraff)
  n1=int((a1-delta)/L*(N-1))+1;
  n2=int((a1+delta)/L*(N-1))+1;
  n3=int((a2-delta)/L*(N-1))+1;
  n4=int((a2+delta)/L*(N-1))+1;
  z1=linspace(0,a1-delta,n1)
  z2=linspace(a1-delta,a1+delta,(n2-n1)*tauraff+1)
  z3=linspace(a1+delta,a2-delta,n3-n2+1)
  z4=linspace(a2-delta,a2+delta,(n4-n3)*tauraff+1)
  z5=linspace(a2+delta,L,N-n4+1)
  z=[z1,z2(2:$),z3(2:$),z4(2:$),z5(2:$)];
endfunction
//
function [nP,nM,nD]=npt(z)
  n=length(z);L=z(n);
  nP=max(find(z<0.25*L));
  nM=max(find(z<0.5*L));
  nD=max(find(z<0.75*L));
endfunction
//
function pinit=fctpinit(dt,Tfin,Timp,pamp);
Nt=int(Tfin/dt);
t=(0:dt:Tfin);
pinit=zeros(t);
Nimp=int(Timp/dt);
timp=(0:dt:Nimp*dt);
pinit(1:Nimp+1)=pamp*sin(%pi/Timp*timp);
endfunction
//
function bet=fctbet(bet1,bet2,delta,a1,a2,z);
  n=length(z);L=z(n);
  n1=max(find(z<a1-delta));
  n2=max(find(z<a1+delta));
  n3=max(find(z<a2-delta));
  n4=max(find(z<a2+delta));
  zbis=[a1-delta,a1,a1+delta];
  betbis=[bet1,(bet1+bet2)/2,bet2];
  d=splin(zbis,betbis,'clamped',[0,0])
  betbis1=interp(z(n1:n2), zbis, betbis, d)';
  zbis=[a2-delta,a2,a2+delta];
  betbis=[bet2,(bet1+bet2)/2,bet1];
  d=splin(zbis,betbis,'clamped',[0,0])
  betbis2=interp(z(n3:n4), zbis, betbis, d)';
  bet=[bet1*ones(1:n1-1)';betbis1;bet2*ones(n2+1:n3-1)';betbis2;bet1*ones(n4+1:n)'];
endfunction
//
function dbet=diffbet(bet,z);
  n=length(z);
  dbet=bet;
  for i=2:n-1
      dbet(i)=(bet(i+1)-bet(i))/(z(i+1)-z(i));
//    dbet(i)=(bet(i+1)-bet(i-1))/(z(i+1)-z(i-1));
  end
    dbet(1)=(bet(2)-bet(1))/(z(2)-z(1)); 
    dbet(n)=(bet(n)-bet(n-1))/(z(n)-z(n-1));
endfunction
//
//--------------------------------------------------------------
//
function [F1,F2]=F(A,Q);
F1=Q;
F2=alpha*(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+1/rho*dbet.*(sqrt(A0).*(A-A0)-2/3*(A.^1.5-A0.^1.5));
endfunction
//
function [HU11,HU12,HU21,HU22]=H(A,Q);
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 [BU11,BU12,BU21,BU22]=BU(A,Q);
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);
[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);
[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()
//
//L1=15;N1=80;delta=0.5;a1=5;a2=10;    // cas 1
L1=60;N1=240;delta=1;Lpro=14;a1=L1/2-Lpro/2;a2=L1/2+Lpro/2;  // cas 2
//
dt=1E-4       // cas sain
//dt=8E-6;   // cas patho
//Timp=0.0025;Tfin=0.035;
Timp=0.165;Tfin=0.25;
//
tauraff=1;betamp=1;     //cas sain
//tauraff=4;betamp=10;     //cas patho
z=fctz(L1,N1,a1,a2,delta,tauraff);
N=length(z)-2;L=z($);
[nP,nM,nD]=npt(z)
//
A=zeros(N+2,1);
Q=zeros(N+2,1);
//
E0=3E6;
hw=0.05;
R0=0.5;
A0val=%pi*R0^2;
rho=1;
nu=0.035;
pamp=20E3;
Kr=8*%pi*nu;
alpha=1;
//
bet1=4*sqrt(%pi)*hw*E0/(3*%pi*R0^2);
bet2=betamp*bet1;
//
bet=fctbet(bet1,bet2,delta,a1,a2,z);
dbet=diffbet(bet,z);
xset('window',0);
xbasc()
xsetech([0,0,1,0.5])
plot2d(z,bet)
plot2d(z,bet,-1)
xsetech([0,0.5,1,0.5])
plot2d(z,dbet)
plot2d(z,dbet,-1)
//
A0=A0val*ones(N+2,1);
//
pinit=fctpinit(dt,Tfin,Timp,pamp);
//
xset('window',1);
xbasc()
Nt=int(Tfin/dt);
t=(0:dt:Tfin);
plot2d(t,pinit)
//
xset('window',3)
xbasc()
xset('window',4)
xbasc()
//
Ainit=A0;
Qinit=zeros(N+2,1);
//
[W1init0,W2init0]=W12(Ainit(1),Qinit(1),bet(1));
[W1initL,W2initL]=W12(Ainit(N+2),Qinit(N+2),bet(1));
//
g1=W2init0*ones(t)+4*sqrt(2/rho)*sqrt(pinit+bet(1)*sqrt(A0(1)));
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;
p=bet.*(sqrt(A)-sqrt(A0));
prP=[0];prM=[0];prD=[0];AP=[A0val];AM=[A0val];AD=[A0val];QP=[0];QM=[0];QD=[0];
//
for i=1:Nt
//
tpsnew=i*dt;
A=Anew;
Q=Qnew;
p=pnew;
//
[LAMB1,LAMB2]=LAMB12(A(1),Q(1),bet(1));
[W1,W2]=W12(A,Q,bet);
W1nplus10=g1(i+1);
W2nplus10=interpln([z;W2'],-LAMB2*dt);
[Anplus10,Qnplu10]=invW12(W1nplus10,W2nplus10,bet(1));
//
[LAMB1,LAMB2]=LAMB12(A(N+2),Q(N+2),bet(N+2));
W1nplus1L=interpln([z;W1'],L-LAMB1*dt);
W2nplus1L=g2(i+1);
[Anplus1L,Qnplus1L]=invW12(W1nplus1L,W2nplus1L,bet(N+2));
//
Anew(1)=Anplus10;
Qnew(1)=Qnplu10;
Anew(N+2)=Anplus1L;
Qnew(N+2)=Qnplus1L;
//
[F1,F2]=F(A,Q);
[FLW1,FLW2]=FLW(A,Q);
[BLW1,BLW2]=BLW(A,Q);
[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)];
//
xset('window',2)
xbasc()
xsetech([0,0,1,0.5])
plot2d(z,Anew);
xsetech([0,0.5,1,0.5])
plot2d(z,Qnew);
//
xset('window',3)
plot2d([tpsnew-dt;tpsnew],[prP($-1),prM($-1),prD($-1);prP($),prM($),prD($)])
xset('window',4)
xsetech([0,0,1,0.5])
plot2d([tpsnew-dt;tpsnew],[AP($-1),AM($-1),AD($-1);AP($),AM($),AD($)])
xsetech([0,0.5,1,0.5])
plot2d([tpsnew-dt;tpsnew],[QP($-1),QM($-1),QD($-1);QP($),QM($),QD($)])
end
