%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ASDstudy.m - Plots shaft torque(Ts)and input power(PT) vs. speed(nm) 
%               for three-phase induction motor operation under
%               V/Hz control.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clf;
% Equivalent circuit values for rated voltage & frequency
R1=0.4; R2p=0.5; X1=1.0; X2p=1.0; Rc=485; Xm=60;
fR=60; VLR=460; p=4;   % Rated voltage & frequency; poles
Ll1=X1/2/pi/fR; Ll2p=X2p/2/pi/fR; Lm=Xm/2/pi/fR; % Inductance values
nf=35; fmax=70;  % No. frequencies for analysis; maximum frequency 
f=linspace(fmax/nf,fmax,nf); f=[f fR];  % Frequency array
% Coefficients & exponent of F&W equation for torque in N-m
a0=0.15; a1=3.62e-6; kF=1.8;
% Coefficients & exponent of load equation for torque in N-m
b0=0; b1=1.11e-3; kL=1.5;
a=0.6;          % Eddy current pu portion of core losses 
npts=250;       % No. points for torque-speed curve
for k=1:nf+1    % Calculation of performance for nf frequencies
   ws=2/p*2*pi*f(k); w=2*pi*f(k); % Synch. speed & radian frequency
   % Line voltage at frequency f(k)
   if f(k)/fR <= 1; VL=f(k)/fR*VLR; else; VL=VLR; end     
   wm=linspace(0,ws-0.0001,npts); % Speed points for analysis
   for m=1:npts       % Calculation of Ts-nm for frequency f(k)
      s=(ws-wm(m))/ws;
      Rcf=(a*fR^2+(1-a)*fR)/(a*f(k)^2+(1-a)*f(k))*Rc;
      Z11=R1+j*w*Ll1; Z22=R2p/s+j*w*Ll2p; Zm=j*w*Lm*Rcf/(Rcf+j*w*Lm);
      Zin=Z11+Z22*Zm/(Z22+Zm);
      I1=VL/sqrt(3)/Zin; I2p=Zm*I1/(Zm+Z22);
      TTd=3*(abs(I2p))^2*R2p/s/ws; % Total developed torque
      Ts(m)=TTd-(a0+a1*(wm(m)*30/pi)^kF); % Shaft torque
      if Ts(m)<0; Ts(m)=0; end
      PT(m)=sqrt(3)*VL*abs(I1)*cos(angle(I1)); % Input power
      %X rotloss(m)=3*(abs(I2p))^2*R2p/1000; % Rotor losses
      %X eff(m)=Ts(m)*wm(m)/PT(m)*100;  % Efficiency
   end
   % Plot complete curves for rated frequency, but only portions
   % break down torque speed(wmmax) for other frequencies.
   if f(k)==fR
      figure(1); plot(wm*30/pi,Ts); grid on; hold on
      title('SHAFT TORQUE-SPEED for VARIABLE FREQUENCY');
      xlabel('SPEED(n_m) - rpm'); ylabel('SHAFT TORQUE(T_s) - N-m');
      fdev=f(2)-f(1);
      text(0.05,0.95,['Frequency increment: ',num2str(fdev),' Hz'],'sc')
      text(0.05,0.90,['Base curve(solidline): ',num2str(fR),' Hz'],'sc')
      TL=b0+b1*(30/pi*wm) .^kL; nm=wm*max(f)/f(k)*30/pi;
      plot(nm, TL,'-.'); % Superimpose load torque plot
      figure(2); plot(wm*30/pi,PT/1000); grid on; hold on
      title('INPUT POWER - SPEED for VARIABLE FREQUENCY');
      xlabel('SPEED(n_m) - rpm'); ylabel('INPUT POWER(P_T) - KW');
      %X figure(3); plot(wm*30/pi,eff); grid on; hold on;
      %X title('EFFICIENCY - SPEED for VARIABLE FREQUENCY');
      %X xlabel('SPEED(n_m) - rpm'); ylabel('EFFICIENCY - %');
      %X figure(4); plot(wm*30/pi,rotloss); grid on; hold on;
      %X title('ROTOR LOSSES - SPEED for VARIABLE FREQUENCY');
      %X xlabel('SPEED(n_m) - rpm'); ylabel('ROTOR LOSSES - KW');
   else
      smax=R2p/sqrt(R1^2+w^2*(Ll1+Ll2p)^2); wmmax=(1-smax)*ws;
      for n=1:npts; if wm(n)>=wmmax; break; end; end
      n1=fix(1.1*n); n2=fix(0.9*n); if n2<1; n2=1; end
      figure(1); plot(wm(n:npts)*30/pi,Ts(n:npts),'--'); hold on
      figure(2); plot(wm(n:npts)*30/pi,PT(n:npts)/1000,'--'); hold on
      %X figure(3); plot(wm(n1:npts)*30/pi,eff(n1:npts),'--'); hold on;
      %X figure(4); plot(wm(n2:npts)*30/pi,rotloss(n2:npts),'--'); hold on;
   end
end; hold off;