sigma=input('Enter solidity ratio= ');
0.1
lamda=input('Enter tip-speed ratio= ');
3
kappa=input('Enter angle of curvature: ');
0
delta=input('Enter pitch angle: ');
-2
theta=thetak';
phai=atan(sin(theta)./(lamda.+cos(theta)));
alpha=phai.-delta;
ct=(1.1*2*pi).*alpha;
geo=0.25*inv(pi)*sigma*sec(kappa).*ct;
Tol=input('Allowed tolerance value ');
0.05
dev=[Tol.*ones(N,2)];
nitermax=input('Number of maximum iterations = ');
100
wn=input('Enter initial value of the horizontal perturbation component: ');
-0.5
w_x=wn.*ones(N,1);
w_y=((sin(theta).*(1.+w_x).-cos(theta).*tan(phai).* (1.+w_x).-lamda.*tan(phai))./(sin(theta).*tan(phai).+cos(theta)));
v=[w_x,w_y];
w_x=v(:,1);
w_y=v(:,2);
F=(cos(theta).+w_x.*cos(theta).+w_y.*sin(theta).+lamda).^2;
f=geo.*(X*F).-sqrt(w_x.^2.+w_y.^2);
dFx=cos(theta).^2.*(2.*w_x.+2).+ 2.*lamda.* cos(theta).+w_y.*sin(2.*theta);
dFy=2.*w_y.*sin(theta).^2.+sin(2.*theta).+w_x.*sin(2.*theta).+2.*lamda.*sin(theta);
dF=[dFx,dFy];
J=geo.*(X*dF).-ones(N,2);
if J==0
disp('Error: (Division by zero, try another initial point) ')
return;
elseif J~=0
wi=v.-(J./f);
errxi=abs(wi(:,1).-v(:,1))./abs(v(:,1));
erryi=abs(wi(:,2).-v(:,2))./abs(v(:,2));
err=[errxi,erryi];
endif;
i=0;
while (err>=dev & 0<=abs(wi)<1 & i<=nitermax)
i=i+1
v=wi;
wi=v.-(J./f);
w_x=wi(:,1);
w_y=wi(:,2);
F=(cos(theta).+w_x.*cos(theta).+w_y.*sin(theta).+lamda).^2;
f=geo.*(X*F).-sqrt(w_x.^2.+w_y.^2);
dFx=cos(theta).^2.*(2.*w_x.+2).+ 2*lamda.* cos(theta).+w_y.*sin(2.*theta);
dFy=2.*w_y.*sin(theta).^2.+sin(2.*theta).+w_x.*sin(2.*theta).+2*lamda.*sin(theta);
dF=[dFx,dFy];
J=geo.*(X*dF).-ones(N,2);
wi=v.-(J./f);
errxi=abs(wi(:,1).-v(:,1))./abs(v(:,1));
erryi=abs(wi(:,2).-v(:,2))./abs(v(:,2));
err=[errxi,erryi];
if err<=dev
fprintf('Iteration complete at i= #%d\n', i)
else
i=i+1;
endif;
endwhile;
end;
