% program for calculation of three-phase winding distribution and mmf

% EMK 2007 A.R. IEA/LTH avo.reinap@iea.lth.se

% input parameters
K_s=0.5;            % slotting factor       
K_m=0.6;            % relative angular width of the surface mounted magnet(s)
N_s=6;              % number of slots
N_p=2;              % number of poles

%N_p=10; N_s=42;     % Lecture L8 example







% -------------------------------------------------- winding calculations
N_ph=3;             % number of phases
q=N_s/(N_ph*N_p);   % slot per phase and per pole
emf = [0 2*pi*N_p/2/N_s*[1:N_s-1]]; % emf angle

% winding factor
kd=zeros(1,21); kp=zeros(1,21); 
for h=1:2:21                              % the harmonic components
    kd(h) = 1/q*sin(h*q*N_p/2*2*pi/N_s/2)/(sin(h*N_p/2*2*pi/N_s/2));    % distribution factor % eq: 2.26
    % kd(h) = 1/q*sin(h*pi/6)/(sin(h*pi/6/q));    % distribution factor
    % a rule of thumb is that the maximum coil span should be equal to the
    % next lowest integer after dividing the number of rotor poles to stator slots
    span = fix(N_s/N_p);				        % coil span one coil in the phase belt
    kp(h) = cos(h/2*pi*(1-span/(N_s/N_p)));     % pitch factor % eq: 2.14 + 2.13 T.A.Lipo
    % kp(h) = sin(h*span*pi/6/q);                 % pitch factor
    if K_s==0 ks(h)=1; else ks(h)=sin(h*q*N_p/2*2*pi/N_s*K_s) / (h*q*N_p/2*2*pi/N_s*K_s); end 
    if K_s==1 ks1=1; end
    kw(h) = kd(h)*kp(h)*ks(h)/h;                % the winding factor   
    Bg(h)=sin((h)*K_m*pi/2)/(h);    % 4/pi*
end

   
   
   
   
   
% -------------------------------------------------- geometric modelling
Rm=[0.3 0.5 0.58 0.6 0.8 1]; % main radii of machine p.u.
% vizualisation settings
PoS=[0.02 0.04 0.36 0.92; 0.4 0.67 0.28 0.32; 0.4 0.34 0.28 0.32; 0.4 0.01 0.28 0.32; 0.7 0.31 0.28 0.68; 0.7 0.01 0.28 0.28];
figure(1); clf; 
scrsz = get(0,'ScreenSize'); 
set(gcf,'Position',[10 100 scrsz(3)-20 scrsz(4)/2],'Color',[1 1 1],'PaperPositionMode','auto')

% winding visualization settings
if ismember(q,[1:20])   % a full pitch winding 
    wdis=[ -2*ones(1,q) 1*ones(1,q) -3*ones(1,q) 2*ones(1,q) -1*ones(1,q) 3*ones(1,q) ];
elseif q==0.0       % a concentrated winding
    wdis=[1 -2 2 -3 3 -1]; 
else
    wdis=[];
end
wdist=[];
for ki=1:N_p/2
    wdist=[wdis wdist];
end
ColD(5,:)=[1 0 0]; ColD(3,:)=ColD(5,:)*0.5;   % Phase A
ColD(6,:)=[0 1 0]; ColD(2,:)=ColD(6,:)*0.5;   % Phase B
ColD(7,:)=[0 0 1]; ColD(1,:)=ColD(7,:)*0.5;   % Phase C

% -------------------------------------------------- emf
wdist=[];
LinS=[':'; ':'; ':'; '-'; '-'; '-'; '-'; ];
subplot('Position',PoS(5,:)); hold on
for ki=N_s:-1:1
    acoil=floor(emf(ki)/(1*pi/N_ph));
    if ismember(acoil,[0:6:100]+0), coil_nr = 1; cnr=-2; end
    if ismember(acoil,[0:6:100]+1), coil_nr = 5; cnr=1; end
    if ismember(acoil,[0:6:100]+2), coil_nr = 3; cnr=-3; end
    if ismember(acoil,[0:6:100]+3), coil_nr = 4; cnr=2; end
    if ismember(acoil,[0:6:100]+4), coil_nr = 2; cnr=-1; end
    if ismember(acoil,[0:6:100]+5), coil_nr = 6; cnr=3; end
    wdist=[cnr wdist];
    if sign(cnr)==-1, ang=N_p/2*pi/N_s+emf(ki)+pi; rad=1.2; lw=2; else ang=N_p/2*pi/N_s+emf(ki); rad=1.1; lw=1; end
    h=polar(ang*[1 1],[0 1.3]); set(h,'Color',ColD(cnr+4,:),'LineWidth',lw,'LineStyle',LinS(cnr+4,:))
    h=text(rad*cos(ang),rad*sin(ang),num2str(ki));
    set(h,'BackgroundColor',[1 1 .5],'Edgecolor',[.7 .7 .7],'FontSize',6,'HorizontalAlignment','center');
end
axis([-1 1 -1 1]*1.25)
axis off, axis equal, axis tight

if or(ismember(q,[1:20]),q==0.0)
else
    sp=span; wdist(2,:)=[ wdist(sp+1:length(wdist)) wdist(1:sp)]*-1;
end

% -------------------------------------------------- geometry of motor parts
subplot('Position',PoS(1,:)); hold on
% stator core
rs=[]; os=[];
if K_s==1; 
    [xs,ys]=pol2cart([0:pi/24:2*pi 2*pi:-pi/24:0],[Rm(5)*ones(1,49) Rm(6)*ones(1,49)]);
elseif K_s==0; 
    [xs,ys]=pol2cart([0:pi/24:2*pi 2*pi:-pi/24:0],[Rm(4)*ones(1,49) Rm(6)*ones(1,49)]);
else
	for ki=1:N_s
        rs=[rs Rm([5*ones(1,7) 4*ones(1,13) 5*ones(1,7)])];
        os=[os pi/N_s+[0:K_s/6:K_s K_s:(1-K_s)/6:(2-K_s) (2-K_s):K_s/6:2]*pi/N_s+(ki-1)*2*pi/N_s];
	end
	rs=[rs Rm(6)*ones(1,49)];
	os=[os (pi/N_s+2*pi):-pi/24:pi/N_s];
	[xs,ys]=pol2cart(os,rs);
end
fill(xs,ys,[0.8 0.8 0.9]-0.1,'EdgeColor',[0.8 0.8 0.9]-0.1)
%slot numeration
for ki=1:N_s
    h=text(mean(Rm(5:6))*cos(pi/N_s+(ki-1)*2*pi/N_s),mean(Rm(5:6))*sin(pi/N_s+(ki-1)*2*pi/N_s),int2str(ki));
    set(h,'BackgroundColor',[1 1 .5],'Edgecolor',[.7 .7 .7],'FontSize',6,'HorizontalAlignment','center');
end
% rotor core
[xs,ys]=pol2cart([0:pi/12:2*pi 2*pi:-pi/12:0],[Rm(1)*ones(1,25) Rm(2)*ones(1,25)]);
fill(xs,ys,[0.8 0.8 0.9]-0.1,'EdgeColor',[0.8 0.8 0.9]-0.1)
% magnets
for ki=1:N_p
    [xs,ys]=pol2cart(pi/N_p-K_m*pi/N_p+[0:pi/12:2*pi 2*pi:-pi/12:0]*K_m/N_p+(ki-1)*2*pi/N_p, ...
        [Rm(2)*ones(1,25) Rm(3)*ones(1,25)]);
    fill(xs,ys,[0.8 0.8 0.9]-0.4,'EdgeColor',[0.8 0.8 0.9]-0.4)
end

% -------------------------------------------------- geometry of windings
for ki=1:N_s
	if ismember(q,[1:20])   % a full pitch winding 
        [xs,ys]=pol2cart(pi/N_s-K_s*pi/N_s+[0:pi/12:2*pi 2*pi:-pi/12:0]*K_s/N_s+(ki-1)*2*pi/N_s, ...
        [Rm(4)*ones(1,25) Rm(5)*ones(1,25)]);
        fill(xs,ys,ColD(wdist(ki)+4,:),'EdgeColor',ColD(wdist(ki)+4,:));
	elseif q==0.0       % a concentrated winding
        wdis=[1 -2 2 -3 3 -1]; 
	else
        [xs,ys]=pol2cart(pi/N_s-K_s*pi/N_s+[0:pi/12:2*pi 2*pi:-pi/12:0]*K_s/N_s+(ki-1)*2*pi/N_s, ...
        [mean(Rm(4:5))*ones(1,25) Rm(5)*ones(1,25)]);
        fill(xs,ys,ColD(wdist(1,ki)+4,:),'EdgeColor',ColD(wdist(1,ki)+4,:));
        [xs,ys]=pol2cart(pi/N_s-K_s*pi/N_s+[0:pi/12:2*pi 2*pi:-pi/12:0]*K_s/N_s+(ki-1)*2*pi/N_s, ...
        [mean(Rm(4:5))*ones(1,25) Rm(4)*ones(1,25)]);
        fill(xs,ys,ColD(wdist(2,ki)+4,:),'EdgeColor',ColD(wdist(2,ki)+4,:));
	end
end
axis equal; axis off

% -------------------------------------------------- space spectrum
Fs=kw; Ts=Fs.*Bg;
subplot('Position',PoS(4,:)); hold on
bar(0:21,[[0 Bg]; [0 Fs]; [0 Ts] ]',1.5)
set(gca,'XLim',[0 20]);
h=legend('B(\theta)','MMF(\theta)','T(\theta)',1); legend boxoff; set(h,'FontSize',9)
colormap([[0.8 0.8 0.9]-0.4; 0 0 1; 0.9 0.7 0.2])
axis off, axis tight


% -------------------------------------------------- MMF distribution
subplot('Position',PoS(3,:)); hold on
set(gca,'XColor',[0 1 0],'FontSize',8) 
mmf1=[-1]; mmf2=[-1]; pos=[];
MMF=[1 -0.5 -0.5]; %MMF=[0 -sqrt(3/4) sqrt(3/4)]; %K_s=1;
for ki=1:N_s %length(wdist)
    pos=[pos 0.5*[1-K_s 1 1+K_s]+(ki-1)];
    mmf1=[mmf1 mmf1(length(mmf1))+MMF(abs(wdist(1,ki)))*[1 1 1]*sign(wdist(1,ki))/q];
    if or(ismember(q,[1:20]),q==0.0)
	else
        mmf2=[mmf2 mmf2(length(mmf2))+MMF(abs(wdist(2,ki)))*[1 1 1]*sign(wdist(2,ki))/q];
	end
end
mmf1=[-1 mmf1(1:length(pos)) ]-mean(mmf1(1:length(pos))); 
pos=[ 0.5 0.5+K_s*0.5 1+pos(1:length(pos)-1)];
%mmf=mmf([length(pos) 1:length(pos)-1]);
plot([0 N_s]+0.5,[0 0],'k-'); for ki=1:4:N_s, plot(ki,0,'+k'); end
plot([0 0]+0.5,[-1 1],'k-'); for ki=-1:0.5:1, plot(0.5,ki,'+k'); end
if or(ismember(q,[1:20]),q==0.0)
    plot(pos,mmf1,'.-b')
else
    mmf2=[-1 mmf2(1:length(pos)-1) ]-mean(mmf2(1:length(pos)));
    plot(pos,0.5*mmf1,'Color',[0.8 0.8 1])
    plot(pos,0.5*mmf2,'Color',[0.8 0.8 1])
    plot(pos,0.5*mmf1+0.5*mmf2,'.-','Color',[0 0 1])
end
for ki=1:4:N_s
    h=text(ki,-0.2,int2str(ki));
    set(h,'BackgroundColor',[1 1 .5],'Edgecolor',[.7 .7 .7],'FontSize',6,'HorizontalAlignment','center');
end
plot([0 N_s]+0.5,[0 0],'k-')
axis off, axis tight

    

% -------------------------------------------------- information box
subplot('Position',PoS(6,:)); hold on
text(0,9,['Ns=' num2str(N_s,'%1.0f')]); text(2,9,['Np=' num2str(N_p,'%1.0f')]); 
text(4,9,['q=' num2str(q,'%1.2f')]); text(6,9,['span=' num2str(span,'%1.1f')])
text(0,7,['kd=' num2str(kd(1),'%1.3f')]); text(2,7,['kp=' num2str(kp(1),'%1.3f')]);
text(4,7,['ks=' num2str(ks(1),'%1.3f')]); text(6,7,['k1=' num2str(kw(1),'%1.3f')]);
text(0,4,['Does not consider properly concentrated'],'Color','r'); 
text(0,2,['windings where q<1'],'Color','r');
axis([0 9 0 10])
axis off





% -------------------------------------------------- phase winding layout
subplot('Position',PoS(2,:)); hold on
%figure(7); clf; hold on
winax=[]; winay=[];
for ki=1:N_s
    if abs(wdist(1,ki))==1
        winax=[winax ki+span*((sign(wdist(1,ki))-1)*-0.5+[0.45 0.45 -0.03 -0.03 0.5 1.03 1.03 0.55 0.55]*sign(wdist(1,ki)))];
        winay=[winay [-0.6 -0.5 -0.4 0.4 0.5 0.4 -0.4 -0.5 -0.7]];
        if sign(wdist(1,ki))==1, 
            if and(ki>0,ki<N_s+1), text(ki+span*(-0.03),0,'\downarrow','HorizontalAlignment','center'); end
            if and(ki>0,ki+span<N_s+1), text(ki+span*(1.03),0,'\uparrow','HorizontalAlignment','center'); end
        else
            if and(ki>0,ki<N_s+1), text(ki+span*(-0.03),0,'\uparrow','HorizontalAlignment','center'); end
            if and(ki>0,ki+span<N_s+1), text(ki+span*( 1.03),0,'\downarrow','HorizontalAlignment','center'); end
        end
    end
end
plot(winax,winay,'r:')
plot(winax([1 length(winax)]),winay([1 length(winax)]),'r.')
winax=[]; winay=[];
for ki=1:N_s
    if abs(wdist(1,ki))==1
        winax=[winax ki-N_s+span*((sign(wdist(1,ki))-1)*-0.5+[0.45 0.45 -0.02 -0.02 0.5 1.02 1.02 0.55 0.55]*sign(wdist(1,ki)))];
        winay=[winay [-0.6 -0.5 -0.4 0.4 0.5 0.4 -0.4 -0.5 -0.7]];
    end
end
plot(winax,winay,'r:')
plot(winax([1 length(winax)]),winay([1 length(winax)]),'r.')
for ki=1:4:N_s
    plot(ki,-0.0,'k+');
    h=text(ki,-0.1,int2str(ki));
    set(h,'BackgroundColor',[1 1 .5],'Edgecolor',[.7 .7 .7],'FontSize',6,'HorizontalAlignment','center');
end
axis off, axis([0.5 N_s+0.5 -0.8 0.6])
