%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&&&
%
%  TAperf.m  -  Calculates magnetization & open circuit sat-
%               uration curves for turboalternator. Then cal-
%               culates load performance.
%               Calls TAdata.m for input data.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clf;
TAdata,
npts=50; Fp=zeros(1,npts+1);  % npts=points for sat curve
Ns=S1/p/m; Cphi=Ns*Cs*p/a; lam1=pi*D/S1; del=(D-Dr)/2;
tau=pi*D/p; taur=pi*Dr/p; Wf=0.25*taur; nr=S2/2/p;
gam2=(pi*Dr-p*Wf)/2/Dr/nr; lamr=Dr/p*gam2;
% Carter coefficients
qty=lam1*(5*del+bs); ks=qty/(qty-bs^2);
qty=lamr*(5*del+bf); kr=(qty/(qty-bf^2)-1)*(1-Wf/taur)+1;
% Stator winding factor
gam=p*pi/S1; kd=sin(Ns*gam/2)/Ns/sin(gam/2);
rho=pitch*gam; kp=cos(pi/2-rho/2); kw=kd*kp;

% Magnetization curve calculations
PHI=1.10e08*VL/sqrt(3)/2.22/kw/Cphi/f;
phi=linspace(0,1.1*PHI,npts); phim=[phi';PHI];
for i=1:npts+1;    % Calculation loop
Fg=0.09972*p*phim(i)*ks*kr*del/D/ln/SF;  % Air gap mmf
t13=pi*(D+2*d1s/3)/S1-bs;                % Stator teeth mmf
Bts=p*phim(i)/S1/t13/ln/SF;Fts=HM690(Bts)*d1s;
t13=(Dr-4*d1f/3)*gam2/p-bf;              % Rotor teeth mmf
Ar13=pi*(Dr-4*d1f/3)/p-2*nr*bf;
Btr=phim(i)/Ar13/l; Ftr=H1010(Btr)*d1f;
ls=pi*(Do+D+2*d1s)/12/p;                 % Stator yoke mmf
Bys=0.5*phim(i)/(Do-D-2*d1s)/ln/SF;
Fys=(HM610(Bys)+HM610(sqrt(3)*Bys)+HM610(2*Bys))*ls;
lr=pi*(Dr-2*d1f)/12/p;       % Rotor core mmf
Brc=0.5*phim(i)/(Dr-2*d1f-2*hc)/(l+2*(1+nr*lamr));
Frc=(H1010(Brc)+H1010(sqrt(3)*Brc)+H1010(2*Brc))*lr;
Fp(1,i)=Fg+Fts+Ftr+Fys+Frc;
end

% Plot magnetization & open circuit saturation curves
flux=phim'/1000; Ifnl=Fp(npts+1)/nr/Csf;
plot(Fp(1,1:npts),flux(1,1:npts),Fp(1,npts+1),flux(1,npts+1),'o');
title(TTL); grid;     % 'o' indicates rated point
xlabel('MMF/pole, Amp-turns'); ylabel('FLUX/pole, kilolines');
text(0.7*max(Fp),0.15*max(flux),['No-load I_f(A) = ',num2str(Ifnl)]);
Voc(1:npts)=flux(1:npts)*sqrt(3)*2.22e-5*kw*Cphi*f;
Ifoc(1:npts)=Fp(1:npts)/nr/Csf;
figure(2); plot(Ifoc,Voc); title(TTL); grid;
xlabel('Field current, A'); ylabel('O-C line voltage, V');

% Resistance & reactance calculations
Zb=VL^2*1e-6/MVA;
alfa=asin((bs+0.5)/lam1); lalf=rho*tau/(2*pi*cos(alfa));
Ra=1.065e-6*(l+2+2*(lalf+d1s))*Cphi/a/sa;
sum=0; for i=1:nr; sum=sum+2*i; end
Rf=1.014e-6*nr*Csf*(2*p*(Wf+l)+pi*Dr*gam2*sum)/sf;
pupitch=p*pitch/S1;       % phase factor
if pupitch <= 0.5; Ks=pupitch;
elseif pupitch <= 2/3; Ks=-0.25+1.5*pupitch;
else; Ks=0.25+0.75*pupitch; end
if S2<S1; ks1=1; else; ks1=S1/S2; end
a1=(pi*D/S1-bs-bf)/2; if a1<0; a1=0; end
Xphi=19.15e-8*f*(kw*Cphi)^2*D*ln*SF/p^2/(ks*kr*del);
Xslt=2.0055e-7*m*f*Cphi^2*l*Ks/S1*(d1c/3/bs+d0s/bs);
Xe=0.3192e-7*f*(kw*Cphi)^2*(0.5+lalf)*Ks/p*(1+ ...
   log(3.54*D/S1/sqrt(bs*d1s)));
Xl=Xslt+Xe; Xsu=Xphi+Xl;
Xphipu=Xphi/Zb; Xlpu=Xl/Zb; Xspu=Xphipu+Xlpu;    % pu reactances

% Performance calculations
ang=-acos(PF);       % Lagging PF assumed
Pw=3e-11*pi*Dr*(l+2*(1+nr*lamr))*(2*pi*f*Dr/p)^3;   % F&W losses
Pf=2*4.62e-6*db^2*lb*(120*f/p)^2; PFW=Pf+Pw;
Wtt=(pi/4*((D+2*d1s)^2-D^2)-S1*bs*d1s)*ln*SF*0.283; % Teeth weight
Wys=pi/4*(Do^2-(D+2*d1s)^2)*ln*SF*0.283;            % Yoke weight
npts=120; Po=linspace(0, 1.2*MVA*PF, npts)*1e06; % Set output power 
for i=1:npts
   Ia(i)=Po(i)/sqrt(3)/VL/cos(ang);
   I=Ia(i)*cos(ang)+j*Ia(i)*sin(ang);
   Er=abs(VL/sqrt(3)+j*I*Xl+I*Ra);
   Ifs=interp1(Voc/sqrt(3),Ifoc,Er);
   Ifu=Ifoc(2)/(Voc(2)/sqrt(3))*Er;
   Xss=(Xsu-Xl)*Ifu/Ifs+Xl; Xs(i)=Xss;
   Ef=VL/sqrt(3)+j*I*Xss+I*Ra;
   If(i)=abs(Ef)/Er*Ifs;
   PFWT(i)=PFW; PCuT(i)=3*Ia(i)^2*Ra; PfT(i)=If(i)^2*Rf;
   phil=Er/(2.22e-8*kw*f*Cphi);
   t13=pi*(D+2*d1s/3)/S1-bs; Bts=p*phil/S1/t13/ln/SF;
   Bys=0.866*phil/(Do-D-2*d1s)/ln/SF;
   PcT(i)=(PM690(Bts)*Wtt+PM610(Bys)*Wys)*2.5;
   Pin=3*real(Ef*conj(I))+PFW+PcT(i)+PfT(i);
   Pout=3*real(VL/sqrt(3)*conj(I));
   eff(i)=Pout/Pin*100;
end

figure(3); subplot(211); plot(Po/1e6,eff); grid
title('Turboalternator performance');
xlabel('Output power, MW'); ylabel('Efficiency, %');
text(0.7*max(Po/1e6),0.15*max(eff),'No fan & pump losses');
subplot(212); plot(Po/1e6,If); grid
title('Field current requirement');
xlabel('Output power, MW'); ylabel('I_f, A');

figure(4); subplot(211); plot(Po/1e6,Ia); grid
title('Turboalternator performance');
xlabel('Output power, MW'); ylabel('Current, A');
figure(5); subplot(111);
plot(Po/1e6,PFWT,Po/1e6,PcT,'--',Po/1e6,PCuT,'-.', ...
   Po/1e6,PfT,':'); grid
title('Core, copper, field, and F&W losses');
xlabel('Output power, MW'); ylabel('Losses, W');
figure(6); plot(0,0,Po/1e6,Xs,'-.'); grid
title('Synchronous reactance');
xlabel('Output power, MW'); ylabel('X_s, Ohms');

% Display equivalent circuit parameters
disp('  '); disp('  ')
disp(['    Xlpu   ' '   Xphipu   ' ' Xspu   ']);
disp([ Xlpu Xphipu Xspu ]);
disp('  '); disp('  ')
disp(['    Rapu ' '     Rf (ohm) ']);
disp([ Ra/Zb Rf ]);

% Display rated point current densities
Deltas=Ia(100)/a/sa; Deltaf=If(100)/sf;
disp('  '); disp('  ')
disp(['  Deltas(Apsi) ' 'Deltaf(Apsi)     ']);
disp([ Deltas Deltaf ]);

% Display length over windings for stator & rotor
LosW=l+1+2*(d1s+lalf*sin(alfa)); LofW=l+2*(1+nr*lamr);
disp('  '); disp('  ')
disp(['  LosW (in) ' 'LofW (in)     ']);
disp([ LosW LofW ]);