%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  imthd.m - Calculates percent total harmonic distortion
%            and percent 3rd & 5th harmonic content of the
%            input current to a single-phase transformer
%            based on harmonic profile determined by hyst.m.
%            Assumes that magnetizing current is 3% of full-
%            load current and that a 2% of full-load core loss
%            fundamental component of current exists in phase
%            with terminal voltage. Must have executed hyst.m 
%            in past so that Imag.mat exists with stored Im 
%            harmonic profile for desired percentage flux.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; load Imag.mat Im; PF=0.80;  % Neg PF indicates leading
npts=50; puload=linspace(0.01,1,npts); % Per unit load
% Add core loss component & fundamental of Im to load current to
% form total fundamental per unit input current.
for i=1:npts
   puI(i)=abs(0.02+puload(i)*exp(-j*PF/abs(PF)*acos(abs(PF)))- ...
      j*Im(2)/max(Im)*0.03);
end
for i=1:npts
   z=0;
   for j=4:2:32   % Odds thru 31st harmonic
      z=z+(Im(j)/max(Im)*0.03/puI(i))^2;
   end
   thd(i)=sqrt(z)*100;               % Total harmonic distortion
   hd3(i)=Im(4)/max(Im)*0.03/puI(i)*100; % Percentage 3rd harmonic
   hd5(i)=Im(6)/max(Im)*0.03/puI(i)*100; % Percentage 5th harmonic
end
load=puload*100;
plot(load,thd,load,hd3,'--',load,hd5,'-.'); grid
title('Harmonic content of transformer input current');
ylabel('Percentage'); xlabel('Percent load current');
legend('Total harmonic distortion','Percent 3rd harmonic', ...
   'Percent 5th harmonic',1);
text(65,0.7*max(thd),['PF = ',num2str(PF)]);