% The INTLAB code to verify the numerical estimates in the paper:
% Parabolic and near-parabolic renormalizations for local degree three
% arXiv: 1510.00043v5, 2024.

% First step: set the precise of numbers: long

>> intvalinit('DisplayInfsup')

>> c6=infsup(6,6); 
mu=11-4*sqrt(c6)
cp_P=1-2*sqrt(c6)/3
cv_P=-16/(3*(8*sqrt(c6)+3))
v1P=9*sqrt(c6)-47/2-sqrt(996*sqrt(c6)-2439)/2
v2P=9*sqrt(c6)-47/2+sqrt(996*sqrt(c6)-2439)/2


>> cp=(1+4*sqrt(c6)+2*sqrt(2*(9+sqrt(c6))))/5
cpp=(1+4*sqrt(c6)-2*sqrt(2*(9+sqrt(c6))))/5
omega=(8*sqrt(c6)-3)/25+i*(6*sqrt(c6)+4)/25
cv=3*(8*sqrt(c6)+3)/4
v1=(2+v1P)/(-v1P)+i*2*sqrt(-1-v1P)/(-v1P)
v2=(2+v2P)/(-v2P)+i*2*sqrt(-1-v2P)/(-v2P)


>> a0=infsup(-0.06,-0.06);
r0=infsup(1.07,1.07);
epsilon1=infsup(0.128,0.128);
epsilon2=infsup(0.007,0.007);
epsilon3=infsup(0.014,0.014);


>> eta=infsup(3,3);
epsilon1^6*(2-epsilon1)^4/((1+epsilon1)*(2+epsilon1)^8)
cv*exp(-2*pi*eta)
(2-epsilon2)^6*epsilon2^4/((1+epsilon2)*(1+epsilon2)^8)
(2+epsilon3)^6*(1+epsilon3)^4/((1-epsilon3)*epsilon3^4*(1-epsilon3)^4)
cv*exp(2*pi*eta)
(real(omega)+0.06)^2+imag(omega)^2-(r0-epsilon3)^2


>> a4=infsup(-0.22+0.69*i,-0.22+0.69*i);
epsilon4=abs(a4-omega);
a5= infsup(0.78+0.21*i,0.78+0.21*i);
epsilon5=abs(a5-omega);
a6=omega;
epsilon6=infsup(0.41,0.41);
epsilon7=infsup(0.82,0.82);
r1=infsup(1.2,1.2);


>> addpath('C:\Users\yangfei\Desktop\functions'); %add function file


>> x60=infsup(0.54,0.54);
x61=real(omega);
x62p=infsup(1.07,1.07);
x62m=infsup(1.067,1.067);
x63=real(omega)+epsilon6;

(x60-real(omega))^2+(y4p(x60)-imag(omega))^2-epsilon6^2
(x60-real(a4))^2+(y6p(x60)-imag(a4))^2-epsilon4^2
(x62p-real(omega))^2+(y5p(x62p)-imag(omega))^2-epsilon6^2
(x62m-real(a5))^2+(y6m(x62m)-imag(a5))^2-epsilon5^2


>> y4p(x60)-23*(x60+1)/26
x601=infsup(0.6,0.6);
Xi4p(x60,x601)
Xi4p(x601,x61)


>> s61=imag(a5+omega)/real(a5-omega);
real(a5)+epsilon5/sqrt(1+s61^2)
imag(a5)-imag(omega)*(real(a5)-1)/(real(omega)-1)
x611=infsup(0.99,0.99);

xi55pmax=(sqrt((real(a5-omega))^2+(imag(a5+omega))^2)+epsilon5)^4;
xi51p(x61)*xi52p(x611)/( xi53p(x611)*xi54p(x611)*xi55pmax)

x612=infsup(1.05,1.05);
x613=infsup(1.06,1.06);

HatXi5p(x611,x612)
HatXi5p(x612,x613)
HatXi5p(x613,x62p)


>> s62=imag(omega)/(real(omega)+1);
real(omega)+epsilon6/sqrt(1+s62^2)
s63=imag(omega)/(real(omega)-1);
real(omega)-epsilon6/sqrt(1+s63^2)
s64=imag(omega)/real(omega);
real(omega)+epsilon6/sqrt(1+s64^2)

x621=real(omega);
x622=0.99;
x623=1.07;
x63=real(omega)+epsilon6;

xi61p(x60)*xi62p(x621)/(xi63p(x621)*epsilon6^4*xi65p(x621))
xi61p(x621)*xi62p(x622)/((1+epsilon6)*epsilon6^4*xi65p(x621))
xi61p(x622)*xi62p(x623)/(xi63p(x622)*epsilon6^4*xi65p(x622))
xi61p(x623)*xi62p(x623)/(xi63p(x622)*epsilon6^4*xi65p(x622))

((x63+1)^2+(imag(omega))^2)^3*((x63-1)^2+(imag(omega))^2)^2/(xi63p(x623)*epsilon6^4*xi65p(x623))

x631=1.069;
x632=1.072;

Xi6m(x631,x62m)
Xi6m(x632,x631)
xi61m(x632)*xi62m(x632)/((x63^2+(imag(omega))^2)^(1/2)*epsilon6^4*((x63-real(omega))^2+(2*imag(omega))^2)^2)


>> x70=-1-epsilon7;
x71=-1.095;
x72=real(a4)-epsilon4
x73=-0.77;

(x71-real(a4))^2+(y7p(x71)-imag(a4))^2-epsilon4^2
(x71+1)^2+(y4p(x71))^2-epsilon7^2
y4m(x73)

s71=-imag(omega)/(real(omega)+1);
-1-epsilon7/sqrt(1+s71^2)

x701=-1.3;
(x70+1)^6*(x70-1)^4/(xi73p(x701)*xi74p(x701)*((x70-real(omega))^2+(imag(omega))^2)^2)
(x70+1)^6*(x70-1)^4/(xi73p(x701)*xi74p(x701)*xi75p(x701))
Xi7p(x71,x701)

s72=imag(a4)/(real(a4)-1);
real(a4)-epsilon4/sqrt(1+s72^2)

xi41p(x71)*xi42p(x71)/((x72^2+(imag(a4))^2)^(1/2)*xi44p(x71)*((x72-real(omega))^2+(imag(a4)+imag(omega))^2)^2)

s73=imag(a4)/(real(a4)+1);
real(a4)-epsilon4/sqrt(1+s73^2)

s74=imag(omega-a4)/real(omega-a4);
real(a4)-epsilon4/sqrt(1+s74^2)

x721=-1;

((x72+1)^2+(imag(a4))^2)^3*((x72-1)^2+(imag(a4))^2)^2/(xi43m(x721)* ((x72-real(omega))^2+(imag(a4)-imag(omega))^2)^2*xi45m(x721))

((x72+1)^2+(imag(a4))^2)^3*((x72-1)^2+(imag(a4))^2)^2/(xi43m(x721)*xi44m(x721)*xi45m(x721))
Xi4m(x73,x721)
xi41m(x73)*xi42m(x721)/(xi43m(x73)*xi44m(x73)*xi45m(x73))


>> x4=a0+r0*r1;
x1=-1.1;
x2=0.54;
x3=1.03;

(x1+1)^2+(yy(x1))^2-epsilon7^2

(x1-real(a4))^2+(yy(x1)-imag(a4))^2-epsilon4^2
(x2-real(a4))^2+(yy(x2)-imag(a4))^2-epsilon4^2

(x2-real(omega))^2+(yy(x2)-imag(omega))^2-epsilon6^2
(x3-real(omega))^2+(yy(x3)-imag(omega))^2-epsilon6^2

(x3-real(a5))^2+(yy(x3)-imag(a5))^2-epsilon5^2
(x4-real(a5))^2+(imag(a5))^2-epsilon5^2

x1p=-0.5;
x2p=0.3;

(x1p-real(a4))^2+(imag(a4))^2-epsilon4^2
(x2p-real(a4))^2+(imag(a4))^2-epsilon4^2
(x2p-real(a5))^2+(imag(a5))^2-epsilon5^2


>> b0=2*(13+32*sqrt(c6))/25
b1=(2029+256*sqrt(c6))/125

a11=2*(617+688*sqrt(c6)); 
a01=25*(119+16*sqrt(c6));
a12=3889250+837000*sqrt(c6); 
a02=2755539+487396*sqrt(c6);
a13=31356325+8965425*sqrt(c6); 
a03=66811702+23697378*sqrt(c6);
a14=102142212+38104768*sqrt(c6); 
a04=240990025+94826600*sqrt(c6);


>> Q2max(11)
64*(617+688*sqrt(c6))/3125
LogDQmax(5.6)


>> c00=infsup(0.06,0.06);
c01max=infsup(2.14,2.14);


>> c11=infsup(11,11);
log(c11)+0.06*pi/3
14*sqrt(c6)/25+0.45


>> alpha1=atan(imag(a5)/(real(a5)-a0))
alpha2m=infsup(0.54,0.54);
alpha2p=infsup(0.55,0.55);
alpha3=atan(imag(omega)/(real(omega)-a0))
alpha4=atan((imag(omega)+epsilon6)/(real(omega)-a0))

h0=14*sqrt(c6)/25+0.45;
alpha5=pi-atan((h0+1)/(a0+1))


>> r5(0)-r2(alpha1)
(a0+r0*r5(alpha2p)*cos(alpha2p)-real(omega))^2 +(r0*r5(alpha2p)*sin(alpha2p)-imag(omega))^2-epsilon6^2
r5(alpha2p)-r2(alpha2p)

r6(alpha4)-r2(alpha4)
(a0+r0*r6(alpha2m)*cos(alpha2m)-real(a5))^2 +(r0*r6(alpha2m)*sin(alpha2m)-imag(a5))^2-epsilon5^2
r6(alpha2m)-r2(alpha3)


>> tt4=ones(1,24);
for k = 2:10
  tt4(k)=1+0.01*k;
end

for k = 11:15
  tt4(k)=1.1+0.02*(k-10);
end

for k = 16:23
  tt4(k)=1.2+0.05*(k-15);
end

r7(alpha4)-r2(tt4(2))
for k = 2:22
  r7(tt4(k))-r2(tt4(k+1))
end
r7(tt4(23))-r2(alpha5)


>> theta51=infsup(2.38,2.38);
theta52=infsup(2.6,2.6);
r8(alpha5)-r2(theta51)
r8(theta51)-r2(theta52)


>> c00+c01max+varphi1max(11)
(cv-11)*sin(pi/6)


>> b0-c00-b1/(2*6.1)
betamax(6.1)


>> c00+c01max+phi1max(cv-2.35)
2.35/exp(-LogDvarphimax(cv-2.35-3.5))
2.4*exp(LogDQmax(cv-2.4))

2.75+Q(cv)-25.5
infsup(1,1)*(25.5-22)*sin(7*pi/20)
b1/20+Q2max(20)


>> (phi1max(125)+b1/122+Q2max(122))/5
LogDFmax(5*25)


>> theta1=infsup(3*pi/20,3*pi/20);
theta2=infsup(pi/4,pi/4);
u1theta1=infsup(8.5,8.5);
u2theta1=infsup(6.1,6.1);
u3=infsup(22,22)*cos(theta1);
u4=infsup(17.3,17.3);
u1theta2=infsup(9,9);
u2theta2=infsup(6.6,6.6);

u0theta1=u1theta1/cos(theta1)
u0theta2=u1theta2/cos(theta2)

u2theta1+c00*cos(theta1)+c01max+varphi1max(u2theta1)
u2theta2+c00*cos(theta2)+c01max+varphi1max(u2theta2)

u4+c00*cos(theta1)+c01max+varphi1max(u4)
u3

ArgDeltaFmax(u2theta1,theta1)
ArgDeltaFmax(u2theta1,-theta1)
-ArgDeltaFmin(u2theta1,theta1)
-ArgDeltaFmin(u2theta1,-theta1)

ArgDeltaFmax(u2theta2,theta2)
ArgDeltaFmax(u2theta2,-theta2)
-ArgDeltaFmin(u2theta2,theta2)
-ArgDeltaFmin(u2theta2,-theta2)

pi/2-theta2


>> r4=infsup(0.34,0.34);
u5=u3-u1theta1;

abs(b0-c00+b1*exp(-theta1*i)/(2*u4)-2*u5*r4^2*exp(theta1*i)/(1-r4^2))+betamax(u4)-2*u5*r4/(1-r4^2)

-ArgDeltaFmin(u4,theta1)+LogDFmax(u4)/2-log(1-r4^2)/2
pi/infsup(5,5)

-ArgDeltaFmax(u4,theta1)-LogDFmax(u4)/2+log(1-r4^2)/2
-3*pi/infsup(20,20)

exp(LogDFmax(u4)/2)/(AbsDeltaFmin(u4,theta1)*sqrt(1-r4^2))
sqrt(1-r4^2)/(AbsDeltaFmax(u4,theta1)*exp(LogDFmax(u4)/2))


>> infsup(1,1)*tan(1.245)
infsup(1,1)*sqrt(1+8^2)/0.083

u6=infsup(10.7,10.7);
b1/u6+Q2max(u6)
(22-b0)*cos(theta1)-u6
LogDQmax(u6)
LogDvarphimax(u6)


>> t0=6.5*sqrt(2)-cp;

tt0=ones(1,12);
for k = 1:8
  tt0(k)=0.5*k;
end

for k = 9:10
  tt0(k)=0.8+0.4*k;
end

for k = 11:12
  tt0(k)=2.8+0.2*k;
end


-3*pi/infsup(4,4)
vartheta(0,tt0(1))
for k = 1:11
  vartheta(tt0(k),tt0(k+1))
end

pi/infsup(2,2)
vartheta(tt0(1),0)
for k = 1:11
  vartheta(tt0(k+1),tt0(k))
end


>> atan(t0/cp)
pi/infsup(4,4)
atan((t0-imag(omega))/(cp-real(omega)))
atan((t0+imag(omega))/(cp-real(omega)))


>> (6.5+b0)*cos(pi/4)+2^4*a11/(5^5*(6.5-1)^2)*cos(pi/4)+Q3max(6.5)
cv*cos(pi/4)


>> cp+c00-c01max-varphi1max(cp)
6.5+c00*cos(pi/4)-c01max-varphi1max(6.5)


>> xx1=ones(1,16);
for k = 1:16
  xx1(k)=0.04*k;
end

-5*pi/infsup(4,4)

theta21(xx1(1))+(theta22(xx1(1))+4*pi)-pi/2-(theta24(0)+4*pi)-(theta25(0)+4*pi)

for k = 1:15
  Theta2(xx1(k),xx1(k+1))-4*pi
end

theta21(real(omega))+theta22(real(omega))-theta23(xx1(16))-theta24(xx1(16))-theta25(xx1(16))-4*pi

theta21(0)+theta22(0)-theta23(xx1(1))-theta24(xx1(1))-theta25(xx1(1))-4*pi

for k = 1:15
  Theta2(xx1(k+1),xx1(k))-4*pi
end

theta21(0.64)+(theta22(0.64)+4*pi)-theta23(real(omega))-2*pi-2*pi
-3*pi/infsup(4,4)


>> Theta2(-1,-0.975)-5*pi
Theta2(-0.975,-0.95)-5*pi
Theta2(-0.95,-0.925)-5*pi
Theta2(-0.925,-0.9)-5*pi

6*pi/2+theta22(-1)-theta23(-0.975)-theta24(-0.975)-theta25(-0.975)-5*pi
Theta2(-0.95,-0.975)-5*pi
Theta2(-0.925,-0.95)-5*pi
Theta2(-0.9,-0.925)-5*pi


>> Xi2(-0.975,-1)*cos(6*pi/2+theta22(-1)-theta23(-0.975)-theta24(-0.975)-theta25(-0.975)-5*pi)
h2(-0.975,-0.95)
h2(-0.95,-0.925)
h2(-0.925,-0.9)

-Xi2(-0.975,-1)-(1.7-4.2*sqrt(2))
-Xi2(-0.95,-0.975)-(1.7-4.2*sqrt(2))
-Xi2(-0.925,-0.95)-(1.7-4.2*sqrt(2))
-Xi2(-0.9,-0.925) -(1.7-4.2*sqrt(2))


>> xx2=ones(1,22);
for k = 4:22
  xx2(k)=-1.1+0.05*k;
end

for k = 4:21
  Theta2(xx2(k),xx2(k+1))-5*pi
end

for k = 4:20
  Theta2(xx2(k+1),xx2(k))-5*pi
end

theta21kmax=theta21(-0.05)+(theta22(-0.05)+4*pi)-pi/2-(theta24(0)+4*pi)-(theta25(0)+4*pi)


>> for k = 4:20
  h2(xx2(k),xx2(k+1))
end

Xi2(0,-0.05)*cos(theta21kmax)

for k = 4:20
  vv2(xx2(k),xx2(k+1))
end

Xi2(0,-0.05)*sin(theta21kmax)-(Xi2(0,-0.05)*cos(theta21kmax)-4.2*sqrt(2))


>> xx3=ones(1,40);
for k = 1:20
  xx3(k)=-1-0.04*k;
end

for k = 21:29
  xx3(k)=2.2-0.2*k;
end

for k = 30:40
  xx3(k)=8-0.4*k;
end

theta31kmin=6*pi/2+(theta32(-1)+4*pi)-(theta33(-1.04)+pi)-(theta34(-1.04)+4*pi)-(theta35(-1.04)+4*pi)

for k = 1:7
  Theta3(xx3(k+1),xx3(k))
end

theta31kmax=Theta3(-1,-1.04)

for k = 1:7
  Theta3(xx3(k),xx3(k+1))
end


>> 4.2*sqrt(infsup(2,2))-1.7
Xi3(-1,-1.04)

for k = 1:7
  Xi3(xx3(k),xx3(k+1))
end


>> Xi3(-1,-1.04)*cos(Theta3(-1,-1.04))

for k = 1:7
  Xi3(xx3(k),xx3(k+1))*cos(Theta3(xx3(k),xx3(k+1)))
end

for k = 8:28
  Xi3(xx3(k),xx3(k+1))
end


>> for k = 29:39
  Theta3(xx3(k+1),xx3(k))
end

for k = 29:39
  Theta3(xx3(k),xx3(k+1))
end


>> for k = 29:39
  h3(xx3(k),xx3(k+1))
end

for k = 29:39
  vv3(xx3(k),xx3(k+1))
end


>> 2*Q2max(8)-8-1.6+b0


>> x42m=(2*real(a4)*imag(a4)*imag(omega)+(real(a4)^2-imag(a4)^2)*real(omega))/(real(a4)^2+imag(a4)^2)
y4m(x42m)
x40m=-0.22-epsilon4
x41m=-1;

s41=-((1+real(a4))*imag(a4)+epsilon4*sqrt((1+real(a4))^2+(imag(a4))^2-epsilon4^2))/(epsilon4^2-(1+real(a4))^2)
x41tilde=(real(a4)+imag(a4)*s41-s41^2)/(1+s41^2)

s42=-((1-real(a4))*imag(a4)+epsilon4*sqrt((1-real(a4))^2+(imag(a4))^2-epsilon4^2))/((1-real(a4))^2- epsilon4^2);
x42tilde=(real(a4)+imag(a4)*s42+s42^2)/(1+s42^2)


>> x401m=-1.096;
theta41m(x401m)+theta42m(x401m)-atan(imag(a4)/x40m)-4*atan((imag(a4)-imag(omega))/(x40m-real(omega)))-4*atan((imag(a4)+imag(omega))/(x40m-real(omega)))+pi+2*pi

6*atan(imag(a4)/(x40m+1))+ 4*atan(imag(a4)/(x40m-1))-theta43m(x401m)-theta44m(x401m)-theta45m(x401m)+pi+2*pi


6*atan(s41)+theta42m(-1)-theta43m(x401m)-theta44m(x401m)- theta45m(x401m)+pi+2*pi
6*(-pi/2)+theta42m(x401m)-theta43m(-1)-theta44m(-1)- theta45m(-1)+pi+2*pi

x411m=-1.04;
Theta41mm(x41tilde,x411m)
Theta41mm(x411m,x41tilde)
Theta41mm(x411m,-1)
Theta41mm(-1,x411m)


>> x40p=-0.22-epsilon4;
x41p=-1;
x42p=0;
x43p=real(omega);

xx40=ones(1,14);
for k = 1:6
  xx40(k)=-1.11+0.005*k;
end

for k = 7:14
  xx40(k)=-1.14+0.01*k;
end

6*atan(imag(a4)/(x40p+1))+ 4*atan(imag(a4)/(x40p-1))-theta43p(xx40(1))-theta44p(xx40(1))-theta45p(xx40(1))+pi+2*pi

for k = 1:13
  Theta4p(xx40(k+1),xx40(k))+pi+2*pi
end

theta41p(xx40(1))+theta42p(xx40(1))- atan(imag(a4)/ x40p)- 4*atan((imag(a4)-imag(omega))/(x40p-real(omega)))-4*atan((imag(a4)+imag(omega))/(x40p-real(omega)))+pi+2*pi

for k = 1:12
  Theta4p(xx40(k),xx40(k+1))+pi+2*pi
end

-3*pi+theta42p(xx40(14))-theta43p(xx40(13))-theta44p(xx40(13))-theta45p(xx40(13))+pi+2*pi


>> xx41=ones(1,50);
for k = 1:20
  xx41(k)=-1+0.01*k;
end

for k = 21:40
  xx41(k)=-1.2+0.02*k;
end

for k = 41:50
  xx41(k)=-2+0.04*k;
end

3*pi+theta42p(-1)-theta43p(xx41(1))-theta44p(xx41(1))-theta45p(xx41(1))-5*pi+2*pi

for k = 1:48
  Theta4p(xx41(k+1),xx41(k))-5*pi+2*pi
end

theta41p(-0.04)+theta42p(-0.04)+pi/2-theta44p(0)-theta45p(0)-5*pi+2*pi

Theta4p(-1,-0.99)-5*pi+2*pi
for k = 1:49
  Theta4p(xx41(k),xx41(k+1))-5*pi+2*pi
end


>> xx42=ones(1,26);
for k = 1:11
  xx42(k)=0.05*k;
end

xx42(12)=0.575;

Theta4p(xx42(1),0)-4*pi+2*pi

for k = 1:11
  Theta4p(xx42(k+1),xx42(k))-4*pi+2*pi
end

Theta4p(x42tilde,xx42(12))-4*pi+2*pi

theta41p(xx42(1))+theta42p(xx42(1))-pi/2-theta44p(0)-theta45p(0)-4*pi+2*pi

for k = 1:11
  Theta4p(xx42(k),xx42(k+1))-4*pi+2*pi
end

Theta4p(xx42(12), x42tilde)-4*pi+2*pi


>> for k = 14:20
  xx42(k)=0.625+0.005*(k-14);
end

for k = 21:25
  xx42(k)=0.659+0.001*(k-21);
end

xx42(26)=0.6635;

Theta4phat(xx42(14),x42tilde)
for k = 14:25
  Theta4phat(xx42(k+1),xx42(k))
end

Theta4phat(x42tilde,xx42(14))
for k = 14:25
  Theta4phat(xx42(k),xx42(k+1))
end

theta41p(0.6635)+(theta42p(real(omega))+4*pi)-theta43p(real(omega))-4*(atan((imag(omega)-0.69)/(real(omega)+0.22))+pi/2)-4*pi/2+2*pi

Theta4phat(0.6635,real(omega))


>> x51m=real(a5)+sqrt(1+(real(a5))^2-2*(real(a5)*real(omega)+imag(a5)*imag(omega)))
x52m=0.78+epsilon5

>> infsup(3,3)*pi/4

xx50=ones(1,12);
for k = 1:11
  xx50(k)=0.66+0.03*k;
end

xx50(12)=1;

for k = 1:10
  Theta5p(xx50(k),xx50(k+1))+4*pi
end

theta51p(1)+4*(-pi/2)-theta53p(xx50(11))-theta54p(xx50(11))-theta55p(xx50(11))+4*pi

for k = 1:11
  Theta5p(xx50(k+1),xx50(k))+4*pi
end


>> Theta5p(xx50(1),real(omega))+4*pi
theta51p(xx50(1))+theta52p(xx50(1))-theta53p(real(omega))-4*atan((0.78-real(omega))/(imag(omega)-0.21))-4*pi/2+4*pi

infsup(3,3)*pi/4


>> xx51p=ones(1,23);
xx51p(1)=1.04;
xx51p(2)=1.08;

Theta5p(1,xx51p(1))
Theta5p(xx51p(1),xx51p(2))

theta51p(1)+4*pi/2-theta53p(xx51p(1))-theta54p(xx51p(1))-theta55p(xx51p(1))

Theta5p(xx51p(2),xx51p(1))


>> s51=imag(a5)/(real(a5)+1);
xtilde51=real(a5)+epsilon5/sqrt(1+s51^2)

s52=imag(a5)/real(a5);
xtilde52=real(a5)+epsilon5/sqrt(1+s52^2)

s53=imag(a5+omega)/real(a5-omega);
xtilde53=real(a5)+epsilon5/sqrt(1+s53^2)


>> for k = 3:10
  xx51p(k)=1.1+0.02*(k-3);
end

for k = 11:16
  xx51p(k)=1.25+0.01*(k-11);
end

xx51p(17)=1.305;

for k = 2:16
  Theta5p(xx51p(k),xx51p(k+1))
end

Theta5p(xx51p(17),xtilde52)

for k = 2:16
  Theta5p(xx51p(k+1),xx51p(k))
end

Theta5p(xtilde52,xx51p(17))

for k = 2:16
  h5p1(xx51p(k),xx51p(k+1))
end

h5p1(xx51p(17),xtilde52)

for k = 2:16
  vv5p1(xx51p(k),xx51p(k+1))
end

vv5p1(xx51p(17),xtilde52)


>> xx51p(19)=1.317;
xx51p(20)=1.323;

Theta5p(xtilde52,xx51p(19))
Theta5p(xx51p(19),xx51p(20))
Theta5p(xx51p(20),xtilde51)

Theta5p(xx51p(19),xtilde52)
Theta5p(xx51p(20),xx51p(19))
Theta5p(xtilde51,xx51p(20))

h5p2(xtilde52,xx51p(19))
h5p2(xx51p(19),xx51p(20))
h5p2(xx51p(20),xtilde51)

vv5p2(xtilde52,xx51p(19))
vv5p2(xx51p(19),xx51p(20))
vv5p2(xx51p(20),xtilde51)


>> xx51p(22)=1.329;

Theta5p(xtilde51,xx51p(22))
Theta5p(xx51p(22),xtilde51)
h5p3(xtilde51,xx51p(22))
vv5p3(xtilde51,xx51p(22))

Theta5p1(xx51p(22),0.78+epsilon5)
Theta5p2(0.78+epsilon5,xx51p(22))

h5p33(xx51p(22),0.78+epsilon5)
vv5p33(xx51p(22),0.78+epsilon5)


>> Q(x51m)
0.25*0.7/DQ(x51m)

x511m=1.293;
x511m-x51m
y5m(x511m)

infsup(1,1)*0.0175/sqrt(2)

s54=imag(omega)/real(x51m-omega);
(1.2-real(a5))^2+(s54*(1.2-x51m)-imag(a5))^2-epsilon5^2


>> xx51m=ones(1,11);
xx51m(1)=1.293;
xx51m(2)=1.297;
xx51m(3)=1.3;
xx51m(4)=1.305;
xx51m(5)=1.31;
xx51m(6)=1.315;
xx51m(7)=1.32;
xx51m(8)=1.325;
xx51m(9)=1.327;
xx51m(10)=1.329;

for k = 1:9
  Theta5m(xx51m(k+1),xx51m(k))
end

Theta5m1110=theta51m(xx51m(10))+theta52m(xx51m(10))-atan(imag(a5)/(0.78+epsilon5))-4*atan(imag(a5-omega)/(0.78+epsilon5-real(omega)))-4*atan(imag(a5+omega)/(0.78+epsilon5-real(omega)))

for k = 1:9
  Theta5m(xx51m(k),xx51m(k+1))
end

6*atan(imag(a5)/(0.78+epsilon5+1))+4*atan(imag(a5)/(0.78+epsilon5-1))-theta53m(xx51m(10))-theta54m(xx51m(10))-theta55m(xx51m(10))


>> for k = 1:9
  Xi5mtilde(xx51m(k),xx51m(k+1))
end

Xi5mtilde1011=((0.78+epsilon5+1)^2+(imag(a5))^2)^3*((0.78+epsilon5-1)^2+(imag(a5))^2)^2/(xi53m(xx51m(10))*((0.78+epsilon5-real(omega))^2+(imag(a5-omega))^2)^2*xi55m(xx51m(10)))

>> for k = 1:9
  Xi5mtilde(xx51m(k),xx51m(k+1))*cos(Theta5m(xx51m(k+1),xx51m(k)))
end

Xi5mtilde1011*cos(Theta5m1110)