%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Te_wm.m - calculates total developed torque ( 3Td ) 
%           using [6.52]. Also determines VTh, RTh, & XTh.
%           Assumes V1 on reference. Requires equivalent
%           circuit parameter values.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
p=4; f=60;                              % Poles, frequency
R1=0.2268; R2pr=0.12528; X1=0.51252;    % Equiv. ckt. values
X2pr=0.76878;Rc=185.1; Xm=9.8554; 
V1=230/sqrt(3);                         % Phase voltage
RcXm=Rc*j*Xm/(Rc+j*Xm);
VTh=RcXm*V1/(R1+j*X1+RcXm); ang=angle(VTh)*180/pi; VTh=abs(VTh);
ZTh=RcXm*(R1+j*X1)/(RcXm+R1+j*X1); RTh=real(ZTh); XTh=imag(ZTh);
Thevenin_Values=[ '       VTh' '         Theta'... 
 '        RTh' '          XTh' ]
format short e
[ VTh ang RTh XTh ]
pause;
npts=99; s=linspace(0.00001,1,npts); s=fliplr(s); ws=2/p*2*pi*f;
R2pr0=R2pr; X2pr0=X2pr; smax=R2pr/sqrt(RTh^2+(XTh+X2pr)^2);
for i=1:npts
% Activate indented code below to frequency-adjust R2pr & X2pr.
%    if i ==1; disp('R2pr & X2pr empirically adjusted for fr'); end  
%    if s(i)>=smax
%    R2pr=(0.5+0.5*sqrt(s(i)/smax))*R2pr0;
%    X2pr=(0.4+0.6*sqrt(smax/s(i)))*X2pr0;
%    else; R2pr=R2pr0; X2pr=X2pr0;
%    end
TTd(i)=3*VTh^2*R2pr/s(i)/ws/((RTh+R2pr/s(i))^2+(XTh+X2pr)^2);
nm(i)=(1-s(i))*ws*30/pi;
end
plot(nm,TTd);grid;
title('Speed-torque curve');
xlabel('Shaft speed, rpm'); ylabel('Torque (3T_d), N-m');