%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  rdrotmmf.m - Program plots the resulting normalized B-field
%               wave over two pole pitches for a round 
%               rotor synchronous machine and gives the Fourier
%               spectrum for assessment of harmonic content. 
%               Effect of rotor slot width is included. Total 
%               harmonic distortion (THD)is calculated and
%               displayed to screen.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clf;
p=2;     % No. poles
S2=20;   % No. rotor slots
bf=2;    % Rotor slot width (in.)
Dr=41;   % Rotor diameter (in.)
taur=pi*Dr/p; Wf=0.25*taur; nr=S2/2/p;
gam2=(pi*Dr-p*Wf)/2/Dr/nr; lamr=Dr/p*gam2;
nlamr=fix(gam2/2/pi*1024);nrmp=fix(bf/pi/Dr*1024);
if nrmp/2~=fix(nrmp/2); nrmp=nrmp-1; end; % Even integer
mag=0; k=1;
for i=1:256  % Build first quarter of B-field wave
   if (i==(2*k-1)/2*nlamr)&(k<nr+1); k=k+1; mag=mag+1/nr;
   else; end
   Ff(i)=mag;
end
k=1;  % Add slot ramps
for i=1:256
   if (i==(2*k-1)/2*nlamr)&(k<nr+1); k=k+1;
      for j=i-nrmp/2+1:i+nrmp/2-1
         Ff(j)=Ff(i-nrmp/2)+1/nr/nrmp*(j-i+nrmp/2);
      end
   else; end
end
Ff=[Ff fliplr(Ff)]; Ff=[Ff -Ff];  % Form balance of waveform
theta=linspace(0,360,1024);
x=abs(fft(Ff));        % Determine Fourier spectrum
y=x(2)/512*sin(theta*pi/180);   % Form fundamental of B-field
subplot(2,1,1); plot(theta,Ff,theta,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],x(1:26)/512); grid;
title('Rotor B-field Fourier spectrum');
xlabel('Harmonic Number'); ylabel('Magnitude');

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