% chaos2

% **************************************************************
% **************************************************************
% ***  aa)  Determination of d_n
% **************************************************************
% **************************************************************

tic

ww0=reshape(eye(4,4),1,16);

iiper1=1;
iihom1=iiper1+Nper1 ;
iihet12=iihom1+Nhom1;
iiper2=iihet12+Nhet12;
iihom2=iiper2+Nper2 ;
iihet21=iihom2+Nhom2;

% ***  almost orthonormal basis V_n and the error delta_4: 
% ***  || V^T*V - I ||  <  delta4
% ***  the matrices \tilde V_n
V = ortho(yy,mue.mid);
VV=sum(reshape(intval(V(:,[1:16 5:16 9:16 13:16])) ...
	.*V(:,[1:16 1:4 1:4 1:4 5:8 5:8 9:12]), NN, 4, 10),2);
VV=midrad(squeeze(VV.mid),squeeze(VV.rad));
VV(:,1:4)=VV(:,1:4)-1;
delta4=intval(max(sup(sqrt(sum(VV(:,[1:10 5:10]).^2,2)))));
clear VV;
radi=delta4/sqrt(1-delta4);
Vtild=midrad(V,radi.sup);

% ******************  Test: representation of the integral manifold Sigma
% ******************        as a function
ll0=5e-4; 
Uyy=midrad(yy,2*sup(sqrt(1+delta4)*ll0));
[fyy,grad] = funcgrad(Uyy,mue); 
Dy1=sum(grad.*V(:,[1:4]),2); 
Dy2=sum(grad.*V(:,[5:8]),2); 
Dy3=sum(grad.*V(:,[9:12]),2); 
Dy4=sum(grad.*V(:,[13:16]),2); 
dC = C0-jacint(intval(yy),mue);
llint=max(abs((abs(dC)+ll0*(abs(Dy1)+intval(abs(Dy3))+abs(Dy4)))./Dy2));
% ******************  Test: llint < ll0
if llint >= ll0, 
  text='******************  llint >= ll0  ******************'
  return
end

% ***  interval inclusion for x_n
ll00=1e-13;
Uyy=yy+infsup(-ll00,ll00)*V(:,5:8);
[fyy,grad]=funcgrad(Uyy,mue);
Dy2=sum(grad.*V(:,5:8),2);
ll=dC./Dy2;
llmax=max(abs(ll));
% ******************  Test: llmax < ll00
if llmax >= ll00, 
  text='******************  llmax >= ll00  ******************'
  return
end
xx=yy+(ll*ones(1,4)).*V(:,5:8);

% ***  interval inclusion xx1 for phi(h_n,xx)
XX = xtaylor(xx,hh,mue,ord-1);
intxx = xtstep(XX,infsup(0,1),ord-1);
rrx=intval(1e-8);
intxx=midrad(intxx.mid,sup(intxx.rad+rrx));
intXX = xtaylor(intxx,hh,mue,ord);
rrmax=sup(sqrt(intval(max(sup(sum(intXX(:,:,ord+1).^2,2))))));
% ******************  Test: rrmax < rrx
if rrmax >= rrx, 
  text='******************  rrmax >= rrx  ******************'
  return
end
XX(:,:,ord+1)=intXX(:,:,ord+1);
clear intXX;
xx1 = xtstep(XX,intval(1),ord);

% ***  neighbourhood Uxx1 (= C_n of paper) containing phi(tn,xx)
z1=sum((xx1-yy(index,:)).*Vtild(index,1:4),2);
z2=sum((xx1-yy(index,:)).*Vtild(index,5:8),2);
z3=sum((xx1-yy(index,:)).*Vtild(index,9:12),2);
z4=sum((xx1-yy(index,:)).*Vtild(index,13:16),2);
Uxx1=midrad(xx1.mid,sup(xx1.rad+sqrt(2*(1+delta4))*abs(z1)*ones(1,4)));
fUxx1=funcgrad(Uxx1,mue);
fq1=sum(fUxx1.*Vtild(index,1:4),2);
fq2=sum(fUxx1.*Vtild(index,5:8),2);
fq3=sum(fUxx1.*Vtild(index,9:12),2);
fq4=sum(fUxx1.*Vtild(index,13:16),2);
% ******************  Test: Uxx1 contains phi(tn,xx)
minfq1=min(inf(fq1));
maxfq234=sqrt(intval(max(sup(fq2.^2+fq3.^2+fq4.^2))));
if minfq1 < maxfq234 , 
  text='******************  phi(tn,xx) not in C_n  ******************'
  return
end

% ***  time step tau_n 
tau=-z1./fq1;
isneg=tau.inf.*tau.sup;
% ******************  Test: 0 in [tau]
if max(isneg)>=0,
  text='******************  tau.inf.*tau.sup > 0  ******************'
  return
end

% ***  estimate for |P_n(0)|: epsn
zeta1=[z3+fq3.*tau z4+fq4.*tau];
epsn=sup(sqrt(sum(zeta1.^2,2)));

% **************************************************************
% **************************************************************
% ***  bb)  Determination of d_n
% **************************************************************
% **************************************************************

% ***  neighbourhood Uxx (= \hat U_n of paper)
Uxx=midrad(yy,1.0001*rho.sup);
[fUxx,gradUxx] = funcgrad(Uxx,mue);
dGq2=sum(gradUxx.*V(:,5:8),2);
dGq3=sum(gradUxx.*V(:,9:12),2);
dGq4=sum(gradUxx.*V(:,13:16),2);
sprime=-[dGq3./dGq2 dGq4./dGq2];
Lambda=sqrt(intval(max(sup(sum(sprime.^2,2)))));
llmax=max(abs(ll));
% ******************  Test: Uxx contains V_n S_n(Omega_n)
if (1+delta4)*(rho^2+(llmax+Lambda*rho)^2) > (1.0001*rho)^2,
  text='******************  V_n S_n(Omega_n) not in Uxx  ******************'
  return
end
VnSnprim=sprime(:,[1 1 1 1 2 2 2 2]).*V(:,[5:8 5:8])+V(:,9:16);

% ***  interval inclusion Uxx1 for phi(h_n,Uxx)
[UXX,UWW] = taylor(Uxx,intval(ones(NN,1)*ww0),hh,mue,ord-1);
[intUxx,intUww] = tstep(UXX,UWW,infsup(0,1),ord-1);
rrx=intval(1e-8);
rrw=intval(1e-6);
intUxx=midrad(intUxx.mid,sup(intUxx.rad+rrx));
intUww=midrad(intUww.mid,sup(intUww.rad+rrw));
[intUXX,intUWW] = taylor(intUxx,intUww,hh,mue,ord);
rxmax=sup(sqrt(intval(max(sup(sum(intUXX(:,:,ord+1).^2,2))))));
rwmax=sup(sqrt(intval(max(sup(sum(intUWW(:,:,ord+1).^2,2))))));
% ******************  Test: rxmax < rrx & rwmax < rrw
if (rxmax >= rrx) | (rwmax >= rrw), 
  text='******************  rxmax >= rrx or rwmax >= rrw ******************'
  return
end
UXX(:,:,ord+1)=intUXX(:,:,ord+1);
UWW(:,:,ord+1)=intUWW(:,:,ord+1);
clear intUXX;
clear intUWW;
[Uxx1,Uww1] = tstep(UXX,UWW,intval(1),ord);


% ***  neighbourhood UUxx1 (\hat C_n of paper) containing phi(tn,Uxx)
z1=sum((Uxx1-yy(index,:)).*Vtild(index,1:4),2);
z2=sum((Uxx1-yy(index,:)).*Vtild(index,5:8),2);
z3=sum((Uxx1-yy(index,:)).*Vtild(index,9:12),2);
z4=sum((Uxx1-yy(index,:)).*Vtild(index,13:16),2);
UUxx1=midrad(Uxx1.mid,sup(Uxx1.rad+sqrt(2*(1+delta4))*abs(z1)*ones(1,4)));
fUUxx1=funcgrad(UUxx1,mue);
Ufq1=sum(fUUxx1.*Vtild(index,1:4),2);
Ufq2=sum(fUUxx1.*Vtild(index,5:8),2);
Ufq3=sum(fUUxx1.*Vtild(index,9:12),2);
Ufq4=sum(fUUxx1.*Vtild(index,13:16),2);
Uminfq1=min(inf(Ufq1));
Umaxfq234=sqrt(intval(max(sup(Ufq2.^2+Ufq3.^2+Ufq4.^2))));
% ******************  Test: UUxx1 contains phi(tn,Uxx)
if Uminfq1 < Umaxfq234 , 
  text='******************  phi(tn,Uxx) not in \hat C_n  ******************'
  return
end

% ***  time step tau_n 
tau=-z1./Ufq1;
isneg=tau.inf.*tau.sup;
% ******************  Test: 0 in [tau]
if max(isneg)>=0,
  text='******************  tau.inf.*tau.sup > 0  ******************'
  return
end
  
% ***  zetaquer = \bar zeta: P_n(zeta) in zetaquer
dz2=tau.*Ufq2;
dz3=tau.*Ufq3;
dz4=tau.*Ufq4;
Uzet2=z2+dz2;
Uzet3=z3+dz3;
Uzet4=z4+dz4;
zetaquer=sup(sqrt(intval(max(sup(Uzet2.^2 + Uzet3.^2 + Uzet4.^2)))));
% ******************  Test: P_n(zeta) wohldefiniert
if zetaquer > ll0 , 
  text='******************  zetaquer > ll0  ******************'
  return
end

% ***  U1xx1: neighbourhood containing phi(t,Uxx) for t in [[hn,tn]]
U1xx1=yy(index,:) + (Uzet2*ones(1,4)).*V(index,5:8) + (Uzet3*ones(1,4)).*V(index,9:12) + (Uzet4*ones(1,4)).*V(index,13:16);

% ***  U1ww1: monodromy matrix M_n(Omega_n)
d0ww1=1e-7*infsup(-1,1)+infsup(0,1.02)*(tau*ones(1,16)).*fprime(UUxx1,Uww1,mue);
Xihat=Uww1+d0ww1;
d1ww1=(tau*ones(1,16)).*fprime(UUxx1,Xihat,mue);
% ******************  Test: in(d1ww1,d0ww1)
if sum(sum(~in(d1ww1,d0ww1))) > 0, 
  text='******************  NOT in(d1ww1,d0ww1) ******************'
  return
end
U1ww1=Uww1+d1ww1;

% ***  U1xx1: set containing phi(tn(x),x) for all x in \hat U_n
fU1xx1 = funcgrad(U1xx1,mue);
fu1nn=fU1xx1./(sum(Vtild(index,1:4).*fU1xx1,2)*ones(1,4));

% ***  AA: midpoint of [DP_n(Omega_n)]
mm1=reshape(reshape(1:16,4,4)',1,16);
mm1=[mm1 mm1];
mm2=[1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 5 6 7 8 5 6 7 8 5 6 7 8 5 6 7 8];
wwVnSn=reshape(sum(reshape(U1ww1(:,mm1).*VnSnprim(:,mm2),NN,4,8),2),NN,8);
q1wwVnSn=reshape(sum(reshape(Vtild(index,[1 2 3 4 1 2 3 4]).*wwVnSn,NN,4,2),2) ,NN,2);
% fUaxx1 = funcgrad(UUxx1,mue);
% fu1nn=fUaxx1./(sum(Vtild(index,1:4).*fUaxx1,2)*ones(1,4));
% fu1nn=fUUxx1./(sum(Vtild(index,1:4).*fUUxx1,2)*ones(1,4));
MatwwVnSn=wwVnSn-fu1nn(:,[1 2 3 4 1 2 3 4]).*q1wwVnSn(:,[1 1 1 1 2 2 2 2]);
AAint=reshape(sum(reshape(Vtild(index,[9:16 9:16]).*MatwwVnSn(:,[1 2 3 4 1 2 3 4 5 6 7 8 5 6 7 8]),NN,4,4),2),NN,4);
AA=AAint.mid;
Delta=max(sup(sqrt(sum(intval(AAint.rad).^2,2))));

% ***  Triangularisation of the matrices A_n
[Qper1,Rper1]= triangdb(AA(iiper1:iiper1+Nper1-1,:),Nper1);
[Qper2,Rper2]= triangdb(AA(iiper2:iiper2+Nper2-1,:),Nper2);
[Qhom1,Rhom1]= triang0(AA(iihom1:iihom1+Nhom1-1,:),Qper1.mid(1,:),Nhom1);
[Qhom2,Rhom2]= triang0(AA(iihom2:iihom2+Nhom2-1,:),Qper2.mid(1,:),Nhom2);
[Qhet12,Rhet12]= triang0(AA(iihet12:iihet12+Nhet12-1,:),Qper1.mid(1,:),Nhet12);
[Qhet21,Rhet21]= triang0(AA(iihet21:iihet21+Nhet21-1,:),Qper2.mid(1,:),Nhet21);
Qn=[Qper1; Qhom1; Qhet12; Qper2; Qhom2; Qhet21];
Rn=[Rper1; Rhom1; Rhet12; Rper2; Rhom2; Rhet21];

% ***  delta3
delta3=intval(max(abs(Qn(:,1).^2+Qn(:,2).^2-1)));
% ******************  Test: delta3 < 1
if (delta3 >= 1), 
  text='******************  delta3 >= 1 ******************'
  return
end

% ***  delta2
[Kinv,Kinvd]= Kinvers(Rn,epsn,Nper1,Nhom1,Nhet12,Nper2,Nhom2,Nhet21, ...
	iiper1,iihom1,iihet12,iiper2,iihom2,iihet21);
Sn=[AA(:,1).*Qn(:,1)+ AA(:,3).*Qn(:,2)- Qn(index,1).*Rn(:,1)- ...
	 Qn(index,3).*Rn(:,2), AA(:,2).*Qn(:,1)+ AA(:,4).*Qn(:,2)- ...
	 Qn(index,2).*Rn(:,1)- Qn(index,4).*Rn(:,2), AA(:,1).*Qn(:,3)+ ...
	 AA(:,3).*Qn(:,4)- Qn(index,1).*Rn(:,3)- Qn(index,3).*Rn(:,4), ...
	 AA(:,2).*Qn(:,3)+ AA(:,4).*Qn(:,4)- Qn(index,2).*Rn(:,3)- ...
	 Qn(index,4).*Rn(:,4)];
Snorm=sqrt(intval(max(sup(sum(Sn.^2,2)))));
delta2=Kinv*Snorm/sqrt(1-delta3);
% ******************  Test: delta2 < 1
if (delta2 >= 1), 
  text='******************  delta2 >= 1 ******************'
  return
end

% ***  delta0, delta1
delta0=Kinvd*sqrt((1+delta3)/(1-delta3))/(1-delta2);
delta1=Delta*Kinv*sqrt((1+delta3)/(1-delta3))/(1-delta2);
% ******************  Test: delta1 < 1
if (delta1 >= 1), 
  text='******************  delta1 >= 1 ******************'
  return
end

% ***  determine r
rr=delta0/(1-delta1);
% ******************  Test: rr < rho
if (rr >= rho), 
  text='******************  rr >= rho ******************'
  return
end

'********* test ok **********'
toc


return


