%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  vcurves.m - calculates data & plots V-curves for synchronous 
%              machine in motor mode.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
VL=2400; KVA=10000; PF=0.8; PR=KVA*PF; % Rated values
% Assume Xsu=1pu, Xl=0.1pu, 1% Cu losses
Xsu=VL^2/(KVA*1000); Xl=0.1*Xsu; Ra=1e-5*VL^2/KVA;

% OC sat curve data using line voltage
Voc=[0 1700 1950 2200 2300 2400 2500 2600 2700 2800 2900 ...
   3000 3100 3300]';
Ifoc=[ 0 70 83 98 105 115 126 140 160 180 205 235 280 400]';
        
% Set plot for rated Ia & If
IaR=1000*KVA/VL/sqrt(3)*(PF+j*sin(acos(PF)));
Er=abs(VL/sqrt(3)-j*IaR*Xl-IaR*Ra);
Ifs=interp1(Voc/sqrt(3),Ifoc,Er); Ifg=Ifoc(2)/(Voc(2)/sqrt(3))*Er;
Xss=(Xsu-Xl)*Ifg/Ifs+Xl;  % Saturated Xs
Ef=abs(VL/sqrt(3)-j*IaR*Xss-IaR*Ra);
IfR=interp1(Voc/sqrt(3),Ifoc,Er)*Ef/Er;
axis([0, IfR, 0, abs(IaR)]);

% Set family of output power values
ncurv=9; Po=linspace(0, PR, ncurv)*1000;

% Set PF angle range
ang=linspace(-60, 60, 50)*pi/180; n=length(ang);
for i=1:ncurv
   for k=1:n
   Ia(k)=Po(i)/sqrt(3)/VL/cos(ang(k));
   I=Ia(k)*cos(ang(k))+j*Ia(k)*sin(ang(k));
   Er=abs(VL/sqrt(3)-j*I*Xl-I*Ra);
   Ifs=interp1(Voc/sqrt(3),Ifoc,Er); Ifg=Ifoc(2)/Voc(2)/sqrt(3)*Er;
   Xss=(Xsu-Xl)*Ifg/Ifs+Xl;
   Ef=abs(VL/sqrt(3)-j*I*Xss-I*Ra);
   If(k)=interp1(Voc/sqrt(3),Ifoc,Er)*Ef/Er;
   end
Iau=Po(i)/sqrt(3)/VL; % Unity PF point
Er=abs(VL/sqrt(3)-j*Iau*Xl-Iau*Ra);
Ifs=interp1(Voc/sqrt(3),Ifoc,Er); Ifg=Ifoc(2)/Voc(2)/sqrt(3)*Er;
Xss=(Xsu-Xl)*Ifg/Ifs+Xl;
Ef=abs(VL/sqrt(3)-j*Iau*Xss-Iau*Ra);
Ifu=interp1(Voc/sqrt(3),Ifoc,Er)*Ef/Er;
Ial=Po(i)/sqrt(3)/VL/0.8; % 0.8 lead PF point
I=Ial*(0.8+j*0.6);
Er=abs(VL/sqrt(3)-j*I*Xl-I*Ra);
Ifs=interp1(Voc/sqrt(3),Ifoc,Er); Ifg=Ifoc(2)/Voc(2)/sqrt(3)*Er;
Xss=(Xsu-Xl)*Ifg/Ifs+Xl;
Ef=abs(VL/sqrt(3)-j*I*Xss-I*Ra);
Ifl=interp1(Voc/sqrt(3),Ifoc,Er)*Ef/Er;
plot(If,Ia,Ifu,Iau,'o',Ifl,Ial,'*'); grid; hold on;
end
title('V-curves for synchronous motor');
xlabel('Field current, A'); ylabel('Stator current, A');
legend(['KW incr. = ',num2str((Po(2)-Po(1))/1000)], ...
  'Unity PF', '0.8 PF leading',4);
hold off;