ti0_TN=0.16438;

tim_TN=fzero(@wf_TN_1_1,[0.0 ti0_TN]);

function [gy]=wf_TN_1_1(x)

ti0_TN=0.16438;

D_TN=6.6E-6;

tim_TN=x;

gy1=ab_TN_1(tim_TN,ti0_TN,@TN_f1);

gy=(gy1-(D_TN/2))/(D_TN/2);

end

function =ab_TN_1(down,up,f)                     %integral subroutine, where double(*f)(double) is the integral function, down and up is the limit of integal

xi=[0.0,-0.90618,-0.538469,0.0,0.538469,0.90618];

ai=[0.0,0.236927,0.478629,0.568889,0.478629,0.236927];

tol=1e-3;

mm=1;

c1=abs(down)+abs(up);

r1=0;

while(1)

s=0;

h1=(up-down)/mm;

for i=1:mm

aa=down+h1*(i-1);

bb=aa+h1;

for j=1:6

x_TN=0.5*((bb-aa)*xi(j)+aa+bb);

s=s+ai(j)*f(x_TN);                  %double(*f)(double) is the integral function

end

end

s=0.5*s*h1;

if(tol-abs(s-r1)/(abs(s)+1)>=0)

break;

end

r1=s;

mm=mm+1;

if((c1+h1==c1)||(mm<=100))

break;

end

end

end

function [f1_TN]=TN_f1(ti_TN)

tim_TN=0.032077;

K11=10e-12;

K33=15e-12;

r1_TN=(K33-K11)/K33;

f1_TN=sqrt(1-r1_TN*cos(ti_TN)*cos(ti_TN))/sqrt((cos(ti_TN)*cos(ti_TN)-cos(tim_TN)*cos(tim_TN))*g_TN(ti_TN));

end

function [g5_TN]=g_TN(ti_TN)

tw_d=5.8677e5;

DZ=1.0908e-6;

tim_TN=0.032077;

PI=3.1415926;

K22=6e-12;

K33=15e-12;

E0=8.85e-12;

ER1=28.0;

ER2=8.0;

MAX_TW_ANG=(250.0*PI/180.0);

D_TN=6.6e-6;

P0=(2*PI/MAX_TW_ANG)*D_TN;

q0=2*PI/P0;

r2_TN=(K33-K22)/K33;

a_TN= (ER1-ER2)/ER1;

g1_TN=cos(tim_TN)*cos(tim_TN)*(1-r2_TN*cos(tim_TN)*cos(tim_TN))*(1-r2_TN*(cos(ti_TN)*cos(ti_TN)+cos(tim_TN)*cos(tim_TN)))*tw_d/(cos(ti_TN)*cos(ti_TN)*(1-r2_TN*cos(ti_TN)*cos(ti_TN)));

g2_TN=2*q0*K22*cos(tim_TN)*cos(tim_TN)*(1-r2_TN*cos(tim_TN)*cos(tim_TN))*tw_d/(K33*cos(ti_TN)*cos(ti_TN)*(1-r2_TN*cos(ti_TN)*cos(ti_TN)));

g3_TN=q0*q0*K22*K22*(cos(ti_TN)*cos(ti_TN)-cos(tim_TN)*cos(tim_TN))/(K33*K33*cos(ti_TN)*cos(ti_TN)*(1-r2_TN*cos(ti_TN)*cos(ti_TN)));

g4_TN=DZ*DZ*a_TN/(ER1*E0*K33*(1-a_TN*cos(ti_TN)*cos(ti_TN))*(1-a_TN*cos(tim_TN)*cos(tim_TN)));

g5_TN=g1_TN-g2_TN-g3_TN+g4_TN;

end