%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   capcurv.m - plots the reactive capability curves for a
%               RR, 3-phase alternator at two rating values
%               dependent on H2 pressure. Constant Xs assumed.
%               0.9 PF leading assumed limit for stability
%               margin.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clf;
VLR=24000; S=400e6;        % Rated line voltage
S=[400e6 350e6];           % Rated apparent power - [HiPr LoPr]
Xs=1.6; Ra=0.08;           % Syn. reactance, phase resistance
PFR=0.85;                  % Rated PF
for k=1:length(S)
  IR=S(k)/sqrt(3)/VLR;     % Rated current
  Ef=VLR/sqrt(3)+IR*exp(-j*acos(PFR))*(Ra+j*Xs); % Rated condition
  % Ia-limit region
  theta=linspace(-acos(0.9),acos(PFR),200);
  m=length(theta); P=[0]; Q=[0];
  for i=1:m
    P=[P sqrt(3)*VLR*IR*cos(theta(i))];
    Q=[Q sqrt(3)*VLR*IR*sin(theta(i))];
  end
  % If-limit region
  delta=linspace(angle(Ef),0,200); n=length(delta);
  for i=2:n
    Ia=(abs(Ef)*exp(j*delta(i))-VLR/sqrt(3))/(Ra+j*Xs);I(i)=abs(Ia);
    PC=sqrt(3)*VLR*conj(Ia);
    P=[P real(PC)]; Q=[Q imag(PC)];
  end
plot(P(m),Q(m),'o',P,Q); legend('Rated point'); hold on
end
title('Reactive capability curve'); grid
xlabel('Output average power, W'); 
ylabel('Output reactive power, VARs');