%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Tphckt.m - solves for currents & voltages of a 3-phase load
%            supplied by a 3-phase source. Source & load can be 
%            unbalanced with or without neutral connection. Both
%            source & load are wye connected. a.b.c phase 
%            sequence is assumed. Line impedance is modeled.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; disp(' ');
disp(' THREE-PHASE CIRCUIT ANALYSIS'); disp(' ');
BalV=input('  Source balanced? ( 1=Yes, 2=No )  ');
if BalV==1
   magVan=input('  magVan = ');
   angVan=input('  angVan = ')*pi/180;
   magVbn=magVan; angVbn=angVan+2*pi/3;
   magVcn=magVan; angVcn=angVan-2*pi/3;
elseif BalV==2
   magVan=input('  magVan = ');
   angVan=input('  angVan = ')*pi/180;
   magVbn=input('  magVbn = ');
   angVbn=input('  angVbn = ')*pi/180;
   magVcn=input('  magVcn = ');
   angVcn=input('  angVcn = ')*pi/180;
else; disp(' INPUT ERROR on VOLTAGE BALANCE');
end
Van=magVan*exp(j*angVan); Vbn=magVbn*exp(j*angVbn);
Vcn=magVcn*exp(j*angVcn); ;disp(' ');
BalL=input('  Load balanced? ( 1=Yes, 2=No )  ');
if BalL==1
   magZa=input('  magZa = ');
   angZa=input('  angZa = ')*pi/180;
   magZb=magZa; angZb=angZa;
   magZc=magZa; angZc=angZa;
elseif BalL==2
   magZa=input('  magZa = ');
   angZa=input('  angZa = ')*pi/180;
   magZb=input('  magZb = ');
   angZb=input('  angZb = ')*pi/180;
   magZc=input('  magZc = ');
   angZc=input('  angZc = ')*pi/180;
else; disp(' INPUT ERROR on LOAD BALANCE');
end;
disp(' '); LZ=input(' Line impedance? ( 1=Yes, 2=No ) ');
if LZ==2; magZl=0; angZl=0;
else
   magZl=input('  magZl = ');
   angZl=input('  angZl = ')*pi/180;
end
Za=magZa*exp(j*angZa); Zb=magZb*exp(j*angZb);
Zc=magZc*exp(j*angZc); Zl=magZl*exp(j*angZl); disp(' ');
N=input(' Neutral connection? ( 1=Yes, 2=No ) '); disp(' ');
if N~=1 & N~=2; disp(' INPUT ERROR on NEUTRAL'); end
if N==1
   A=[Za+Zl 0 0; 0 Zb+Zl 0; 0 0 Zc+Zl];
   I=inv(A)*[Van Vbn Vcn].';
else
   A=[Za+Zl -Zb-Zl 0; 0 Zb+Zl -Zc-Zl; 1 1 1];
   I=inv(A)*[Van-Vbn Vbn-Vcn 0].';
end
VN=[Za 0 0; 0 Zb 0; 0 0 Zc]*I;
disp(' '); disp(' LOAD CURRENT in POLAR FORM');
disp([abs(I) angle(I)*180/pi]);
if N==1
   disp(' '); In=I(1)+I(2)+I(3); angIn=angle(In)*180/pi;
   disp([' In = ', num2str(abs(In)), '|_', ...
         num2str(angIn), ' Load->Source']);
end
disp(' '); disp(' LOAD VOLTAGE in POLAR FORM');
disp([abs(VN) angle(VN)*180/pi]);