Help calculating the thrust coefficient

Based on the code provided

My thrust coefficient is half of what is supossed to be. Can anyone point out what my mistakes are on my equation. The code starts from %%verification of conservation of momentum. I know based on my TSR my thrust coefficient is supposed to be somewhere around .7

Also my professor helped me calculate the power coefficient but i do not understand what this lines of code are doing?

What is tke1 and tke2 doing?

i know i have to do the same thing for the i,j and i,k but how would i do it?

could you do it for i,j component and comment how what is that it is doing it so that i can understand and do it fo i,k.

tke1=(u(istart,j,k)^2.0+v(istart,j,k)^2.0+w(istart,j,k)^2.0)/2;
tke2=(u(iend,j,k)^2.0+v(iend,j,k)^2.0+w(iend,j,k)^2.0)/2;

power_1=power_1+((-u(istart,j,k)))*((tke1)+pr(istart,j,k))*da;
power_2=power_2+(u(iend,j,k))*((tke2)+pr(iend,j,k))*da;

D = 126; %diamter
R = D/2;
Re = 7.3*10^7; %reynolds number
Rho = 1.225; %density
Mu = 1.785*10^-5 %viscosity of air
velocity = Re*Mu/Rho/D %velocity profile 
TSR = 7.5 %tip speed ration
Omega = TSR*velocity/R %angular velocity 
Power_in = 1/2*Rho*pi*R^2*velocity^3
figure 
utmpx(:,:)=u(:,:,n3/2);
contourf(xx,yy,utmpx',100, 'linestyle','none')
daspect([1,1,1])

figure 
[~,ind]=min(abs(xx-3));
utmpz(:,:)=u(ind,:,:);
contourf(zz,yy,utmpz,80, 'linestyle','none')
daspect([1,1,1])

%flux  going through the faces
f1=0;
f2=0;
for j=1: n2-1
    for k=1:n3-1
        da=(yy(j+1)-yy(j))*(zz(k+1)-zz(k));
        f1=f1+u(1,j,k)*da;
        f2=f2+u(n1,j,k)*da;
    end 
end
f3=0;
f4=0;
for i=1:n1-1
    for j=1:n2-1
        da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));
        f3=f3+w(i,j,1)*da;
        f4=f4+w(i,j,n3)*da;
    end
end
f5=0;
f6=0;
for i=1:n1-1
    for k=1:n3-1
        da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));
        f5=f5+v(i,n2,k)*da;
        f6=f6+v(i,1,k)*da;
    end
end
total_flux=-f1+f2+f3-f4+f5-f6

%%verification of conservation of momentum
f1=0;
f2=0;
p1=0;
p2=0;

f1=0;
f2=0;
for j=1: n2-1
    for k=1:n3-1
        da=(yy(j+1)-yy(j))*(zz(k+1)-zz(k));
        f1=f1+u(1,j,k)^2*da;
        f2=f2+u(n1,j,k)^2*da;
        p1=p1+pr(1,j,k)*da;
        p2=p2+pr(n1,j,k)*da;
    end 
end
f3=0;
f4=0;
for i=1:n1-1
    for j=1:n2-1
        da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));
        f3=f3+u(i,j,1)*w(i,j,1)*da;
        f4=f4+u(i,j,n3)*w(i,j,n3)*da;
    end
end
f5=0;
f6=0;
for i=1:n1-1
    for k=1:n3-1
        da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));
        f5=f5+u(i,n2,k)*v(i,n2,k)*da;
        f6=f6+u(i,1,k)*v(i,1,k)*da;
    end
end
trust_coefficient =-f1+f2+f3-f4+f5-f6+p1-p2

%Thrust Force 
TF = .5*trust_coefficient*1.225*velocity^2*pi*63^2

% Power extracted by the wind turbine
power_1=0;
power_2=0;
istart=75
iend=250
for j=1:n2-1
    for k=1:n3-1
        da=(zz(k+1)-zz(k))*(yy(j+1)-yy(j));
        tke1=(u(istart,j,k)^2.0+v(istart,j,k)^2.0+w(istart,j,k)^2.0)/2;
        tke2=(u(iend,j,k)^2.0+v(iend,j,k)^2.0+w(iend,j,k)^2.0)/2;

        power_1=power_1+((-u(istart,j,k)))*((tke1)+pr(istart,j,k))*da;
        power_2=power_2+(u(iend,j,k))*((tke2)+pr(iend,j,k))*da;
   end
end
power_3=0;
power_4=0;
for i=1:n1-1
    for j=1:n2-1
        da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));
        power_3=power_3+((-w(i,j,1)))*(((u(i,j,1)^2.0)/2)+pr(i,j,1))*da;
        power_4=power_4+(w(i,j,n3))*(((u(i,j,n3)^2.0)/2)+pr(i,j,n3))*da;
    end
end
power_5=0;
power_6=0;
for i=1:n1-1
    for k=1:n3-1
        da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));
        power_5=power_5+((v(i,100,k)))*(((u(i,100,k)^2.0)/2)+pr(i,100,k))*da;
        power_6=power_6+(-v(i,1,k))*(((u(i,1,k)^2.0)/2)+pr(i,1,k))*da;
    end
end
WT_power_1=4/3.14*(power_1+power_2)*2
WT_power_2=power_3+power_4
WT_power_3=power_5+power_6
c_p = -(WT_power_1+WT_power_2+WT_power_3)

%For visualisationn purpose Z and Y are transposed 
figure 
hold on
[X,Y,Z]=meshgrid(xx,zz,yy);
U=permute(u,[3,1,2]);
W=permute(v,[3,1,2]);
V=permute(w,[3,1,2]);

%Particles lines from 
[sx,sy,sz]=meshgrid(0,linspace(1.25,1.75,5),linspace(.4,.8,5));
h=streamline(X,Y,Z,U,V,W,sx,sy,sz);
set(h,'color','green');
daspect([1,1,1])
axis([min(xx),max(xx),min(zz),max(zz),min(yy),max(yy)])
z0=1.5;y0=.7031;
theta=0:.1:2*pi;
zc=.5*cos(theta)+z0;
yc=.5*sin(theta)+y0;
fill3(ones(1,length(zc))*3,zc,yc,[.6,.6,.6])
grid on
view(3)
xlabel('x')
ylabel('y')
zlabel('z')
        

% for visualization Z and Y are transposed 
figure 
hold on
[X,Y,Z]=meshgrid(xx,zz,yy);
U=permute(u,[3,1,2]);
W=permute(v,[3,1,2]);
V=permute(w,[3,1,2]);
%particles lines released from rotor plane 
[sx,sy,sz]=meshgrid(3,linspace(1.25,1.75,5),linspace(.4,.8,5));
h=streamline(X,Y,Z,U,V,W,sx,sy,sz);
set(h,'color','red');
daspect([1,1,1])
axis([min(xx),max(xx),min(zz),max(zz),min(yy),max(yy)])
z0=1.5;y0=.7031;
theta=0:.1:2*pi;
zc=.5*cos(theta)+z0;
yc=.5*sin(theta)+y0;
fill3(ones(1,length(zc))*3,zc,yc,[.6,.6,.6])
grid on
view(3)
xlabel('x')
ylabel('y')
zlabel('z')
%Pressure plot
figure 
pr_vsec(:,:)=pr(:,:,n3/2);
contourf(xx,yy,pr_vsec',180,'linestyle',' none')
daspect([1,1,1])
xlabel('x')
ylabel('y')

power = 1/2*Rho*velocity^3*pi*R^2*c_p