%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  cleak2.m - Analysis of C-shaped magnetic circuit with
%             leakage flux. Calculates air gap flux for 
%             specified value of coil current. Plots Flux 
%             vs. mmf. (Type 2 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
I=20.0;              % Coil current for analysis point - A
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;
% Flux vs. mmf curve
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
% Specified point analysis
phig=interp1(F, phip, N*I);
phill=interp1(F, phil, N*I);
disp([blanks(4) 'C-SHAPED MAGNETIC CIRCUIT, With Leakage']);
disp(' '); disp([blanks(8) 'AIR GAP FLUX = ', num2str(phig)]);
disp(' '); disp([blanks(8) 'LEAKAGE FLUX = ', num2str(phill)]);
disp(' '); disp([blanks(8) 'COIL CURRENT = ', num2str(I)]);
plot(F,phip,F,phil,'--',N*I,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);