%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   sprBfld.m - Program plots the normalized B-field wave
%               form for a SP synchronous machine and gives
%               the Fourier spectrum for assessment of
%               harmonic content. Total harmonic distortion
%               (THD) is displayed to screen.
%               The pole arc design is nonconcentric with 
%               rotor center.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear, clf;
Dr=15;            % Rotor dia (in)
delc=0.125;       % Air gap length @ pole center (in)
k=0.95;           % pu dia of nonconcentric pole arc
psi=0.6667;       % Pole arc/pole pitch
b=k*Dr/2;         % Pole arc radius
a=Dr/2-b-delc;    % Offset for pole arc center
thetap=psi*pi/2;  % Span of half pole arc

% Build permeance over pole arc
theta=linspace(0,pi/2,256); npts=length(theta);
for i=1:npts; 
   if theta(i) >= thetap; nedge=i-1; break; end
end
for i=1:nedge
   x=a*cos(theta(i))+sqrt(a^2/2*(1-cos(2*theta(i)))+b^2);
   del=Dr/2-x;
   P(i)=1/del;
end
for i=nedge+1:npts  % Parabolic decay beyond pole shoe edge
   c=P(nedge)/(theta(nedge)-pi/2)^2;
   P(i)=c*(theta(i)-pi/2)^2;
end
P=[ fliplr(P) P]; P=[ P -fliplr(P) ]; % Balance of waveform
Bf=abs(fft(P)); Bf=Bf/max(Bf);  % Form normalized FFT
phi=linspace(0,360,1024);
y=Bf(2)*sin(phi*pi/180);   % Form fundamental of Bf
subplot(2,1,1); plot(phi,P/max(P),phi,y,'--'); grid;
title('Normalized rotor B-field waveform');
xlabel('Rotor position, deg'); ylabel('Br, pu');
legend('Total B-field','Fundamental B-field',1);
subplot(2,1,2); plot([0:25],Bf(1:26)); grid;
title('Rotor B-field Fourier spectrum');
xlabel('Harmonic number'); ylabel('Magnitude');

% Total harmonic distortion
THD=0; for i=3:50; THD=THD+(Bf(i)/Bf(2))^2; end
THD=sqrt(THD)*100; disp(' ');
disp(['TOTAL HARMONIC DISTORTION(%) = ',num2str(THD)]);
disp(' ');