%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  phasemmf.m - Determines air gap mmfs for a p-pole, 3-phase
%               stator winding with integer slots/pole/phase. 
%               Values are normalized by Ni. First set of plots 
%               are mmfs for one coil pair. Second set of plots
%               are mmfs for phase winding.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
skew=input('Skew in no. slots (0-1)   ');
pitch=input('Coil pitch in no. slots   ');
p=input('No. poles   ');
S1=input('No. stator slots   '); clf;
x=ones(1,1024)/2; y=zeros(1,1024); z=y;
% Generate quasi-trapezoidal mmf of single coil.
m=p/2*round(skew/S1/2*1024); n=p/2*round(pitch/S1*1024);
x(n:1024)=-x(n:1024);
for i=1:m; if m==0; break; end
x(i)=(i-1)/m/2; x(n+1-i)=x(i); x(n+i)=-x(i); x(1025-i)=-x(i);
end
y=-shift(x,-0.5);   % Shift for mmf of second coil.
% Add mmfs to give coil pair resultant mmf.
ang=linspace(0,360,1024); z=x+y;  
% Plot component & resultant mmfs.
subplot(2,1,1); plot(ang,x,'--',ang,y,'-.',ang,z);
axis([0 360 -1.25 1.25]); grid;
title('Normalized air gap mmf for coil pair');
xlabel('Angular position, elect. deg.'); ylabel('mmf/Ni, pu');
legend('Coil 1', 'Coil 4', 'Coil pair',0);
z=abs(fft(z)); z=z/max(z)*100; % Determine Fourier coefficients.
subplot(2,1,2); axis; plot([0:15],z(1:16)); grid
title('Fourier coefficients of coil pair mmf');
xlabel('Harmonic number'); ylabel('Harmonic magnitude, %');
pause; figure;
% Construct mmf resulting from all coil pairs of a phase.
Nspp=fix(S1/2/3); z=x+y; y=z; vertscale=Nspp*max(z)+0.25;
subplot(2,1,1); plot(ang,z,'--');
axis([0,360,-vertscale,vertscale]); grid; hold;
title('Normalized air gap mmf for phase winding');
xlabel('Angular position, elect. deg.'); ylabel('mmf/Ni, pu');
for i=2:Nspp; if Nspp==1; break; end
x=shift(z,-1/S1); if Nspp<=4; plot(ang,x,'--'); else; end
y=y+x; z=x; end
plot(ang,y,'b');
legend('Coil pair mmf', 'Phase mmf');
% Determine phase mmf Fourier coefficients.
y=abs(fft(y)); y=y/max(y)*100;
subplot(2,1,2);hold; plot([0:15],y(1:16)); grid;
title('Fourier coefficients of phase mmf');
xlabel('Harmonic number'); ylabel('Harmonic magnitude, %');