%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  cleak1.m - Analysis of C-shaped magnetic circuit with
%             leakage flux.Calculates coil current for 
%             value of air gap flux & plots Flux vs. mmf.
%             (Type 1 problem)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
% Units of input dimensions - 'MKS' or 'FPS' (British)
dim='MKS';
if dim=='MKS'; k=1; elseif dim=='FPS'; k=39.37; else; end
w=0.15; h=0.10;      % Window width & height
d=0.05; del=0.002;   % Core depth & air gap length
l1=0.05; l2=0.05; l3=0.05; l4=0.05; % Core widths
phig=0.0040;         % Air gap flux for analysis point - Wb
N=500; SF=0.95;      % Coil turns, stacking factor
w=w/k; h=h/k; d=d/k; del=del/k;
l=[ l1 l2 l3 l4]/k;  % All dimensions now in MKS units
lm=[h+(l(2)+l(4))/2 w+(l(1)+l(3))/2 h+(l(2)+l(4))/2-del ...
      w+(l(1)+l(3))/2];
Pg=4e-7*pi*(l(3)*d/del+0.52*(l(3)+d)+0.308*del);
Pl=4e-7*pi*w*d/h;
% Specified point analysis
F1=phig/Pg;
F1=F1+hm27(phig/l(3)/d/SF)*lm(3)+hm27(phig/l(2)/d/SF) ...
   *lm(2)/2+hm27(phig/l(4)/d/SF)*lm(4)/2;
phil=F1*Pl; phigl=phig+phil;     % Leakage & coil fluxes
F1=F1+hm27(phigl/l(1)/d/SF)*lm(1)+hm27(phigl/l(2)/d/SF) ...
   *lm(2)/2+hm27(phigl/l(4)/d/SF)*lm(4)/2;
I1=F1/N;
disp([blanks(4) 'C-SHAPED MAGNETIC CIRCUIT, With Leakage']);
disp(' '); disp([blanks(8) 'COIL CURRENT = ', num2str(I1)]);
disp(' '); disp([blanks(8) 'AIR GAP FLUX = ', num2str(phig)]);
disp(' '); disp([blanks(8) 'LEAKAGE FLUX = ', num2str(phil)]);
% Flux vs. mmf plot
phimax=2.0*d*min(l)*SF;
phip=linspace(0,phimax,100); npts=length(phip);
for j=1:npts
   F(j)=phip(j)/Pg;
   F(j)=F(j)+hm27(phip(j)/l(3)/d/SF)*lm(3)+hm27(phip(j)/ ...
      l(2)/d/SF)*lm(2)/2+hm27(phip(j)/l(4)/d/SF)*lm(4)/2;
   phil(j)=F(j)*Pl; phigl=phip(j)+phil(j); 
   F(j)=F(j)+hm27(phigl/l(1)/d/SF)*lm(1)+hm27(phigl/l(2)/d/SF) ...
      *lm(2)/2+hm27(phigl/l(4)/d/SF)*lm(4)/2;
   
end
plot(F,phip,F,phil,'--',F1,phig,'o'); grid
title('C-shaped magnetic circuit with leakage');
xlabel('Coil mmf, A-t'); ylabel('Flux, Wb');
legend('Air gap flux', 'Leakage flux', 'Specified pt.', 0);