%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Td_nmvf.m - Plots total developed torque(TTd) vs speed(nm)
%              for variable frequency induction motor operation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clf;
R1=0.2; R2=0.3; X1=0.754; X2=0.754; 
Rc=110; Xm=33.9;
p=4; fR=60; VR=460; nmLim=2100;  % nmLim is max allowable speed

% Thevenin resistance & inductance
Zm=j*Xm*Rc/(Rc+j*Xm);
ZTh=Zm*(R1+j*X1)/(R1+j*X1+Zm);
RTh=real(ZTh); LTh=imag(ZTh)/2/pi/fR;
L1=X1/2/pi/fR; L2=X2/2/pi/fR; Lm=Xm/2/pi/fR;

% Frequencies for analysis in addition to rated
nf=14; f=linspace(nmLim/nf,nmLim,nf)*p/120;

% Rated frequency calculations( Base curve )
npts=50; ws=2/p*2*pi*fR; wm=linspace(0,ws-0.01,npts);
w=2*pi*fR; VTh=abs(VR/sqrt(3)*Zm/(R1+j*X1+Zm));
for i=1:npts
   s=(ws-wm(i))/ws;
   TTd(i)=3*VTh^2*R2/s/ws/( (RTh+R2/s)^2+w^2*(LTh+L2)^2 );
end
plot(wm*30/pi,TTd);grid;
title('Induction motor torque-speed for variable frequency');
xlabel('Speed, rpm'); ylabel('Total developed torque, N-m');
fdev=f(2)-f(1);
text(0.65,0.95,['Frequency increment: ',num2str(fdev),' Hz'],'sc')
text(0.65,0.92,['Base curve(solidline): ',num2str(fR),' Hz'],'sc')
hold on;  % Hold base curve for overplotting

for k=1:nf  % Other than rated frequency calculations
ws=2/p*2*pi*f(k); wm=linspace(0,ws-0.01,npts);
w=2*pi*f(k);

% Set voltage-frequency control
if f(k)/fR <= 1; b=0.06*VR; VL=(VR-b)*f(k)/fR+b; else; VL=VR; end

% Empirical adjustment of core loss resistance
Rcf=0.5*Rc*(VR/VL*f(k)/fR)^2*(fR/f(k)+(fR/f(k))^2);
Zm=j*w*Lm*Rcf/(Rcf+j*w*Lm);
ZTh=Zm*(R1+j*w*L1)/(R1+j*w*L1+Zm);
RTh=real(ZTh); LTh=imag(ZTh)/2/pi/f(k);
VTh=abs(VL/sqrt(3)*Zm/(R1+j*w*L1+Zm));
for i=1:npts
   s=(ws-wm(i))/ws;
   TTd(i)=3*VTh^2*R2/s/ws/( (RTh+R2/s)^2+w^2*(LTh+L2)^2 );
end

% Determination of points above nmmax for plot
smax=R2/sqrt(R1^2+w^2*(LTh+L2)^2); wmmax=(1-smax)*ws;
for m=1:npts; if wm(m)>=wmmax; break; end; end
plot(wm(m:npts)*30/pi, TTd(m:npts),'--');
end
% Activate to superimpose load torque plot
%TL=50+0.004052847*wm .^2;
%plot(wm*30/pi, TL,'-.');
hold off;