%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  dwleak2.m - Analysis of double window magnetic circuit
%              with leakage flux. Calculates air gap fluxes 
%              for specified value of coil current & plots
%              Flux vs. mmf. (Type 2 problem)
%              
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
% Units of input dimensions - 'MKS' or 'FPS' (British)
dim='FPS';
if dim=='MKS'; k=1; elseif dim=='FPS'; k=39.37; else; end
w1=3.50; h1=4.00;                % Window 1 width & height
w2=2.50; h2=3.00;                % Window 2 width & height
d=3.00; del1=0.100; del2=0.075;  % Core depth & air gap length
l1=3.00; l2=3.00; l3=2.00; l4=1.00; l5=1.00;  % Core widths
I=7.5;                   % Coil current (analysis point) - A
N=1000; SF=0.95;         % Coil turns, stacking factor
lm1=h1+l2; lm2=w1+(l1+l3)/2;     % Mean length paths
lm3=(h1+l2+h2+l4)/2-del1; lm4=w2+(l3+l5)/2; lm5=h2+l4-del2;
w1=w1/k; w2=w2/k; h1=h1/k; h2=h2/k; d=d/k;
l=[ l1 l2 l3 l4 l5]/k; del1=del1/k; del2=del2/k; 
lm=[lm1 lm2 lm3 lm4 lm5]/k; % All dimensions now in MKS units
Pg1=4e-7*pi*(l(3)*d/del1+0.52*(l(3)+d)+0.308*del1);
Pg2=4e-7*pi*(l(5)*d/del2+0.52*(l(5)+d)+0.308*del2);
Pl1=4e-7*pi*w1*d/h1; Pl2=4e-7*pi*w2*d/h2;
% Generate flux vs. Fc & Fr arrays (center & right legs)
phimax=2.2*d*min([l(4) l(5)])*SF; npts=250;
phir=linspace(0,phimax,npts);
for i=1:npts
   Fr(i)=phir(i)/Pg2+hm27(phir(i)/l(5)/d/SF)*lm(5)+ ...
      hm27(phir(i)/l(4)/d/SF)*lm(4);
   phil1(i)=Fr(i)*Pl1; phirt(i)=phir(i)+phil1(i);
   Fr(i)=Fr(i)+hm27(phirt(i)/l(4)/d/SF)*lm(4);
end
phimax=2.2*d*l(3)*SF; phic=linspace(0,phimax,npts);
for i=1:npts
   Fc(i)=phic(i)/Pg1+hm27(phic(i)/l(3)/d/SF)*lm(3);
end
% Generate coil flux(phit) vs. MMF(FT) array
Fmax=min([max(Fc) max(Fr)]); F=linspace(0,Fmax,npts);
for i=1:npts
   phi(i)=interp1(Fc,phic,F(i))+interp1(Fr,phirt,F(i));
   FT(i)=F(i)+hm27(phi(i)/l(2)/d/SF)*lm(2);
   phil2(i)=FT(i)*Pl2; phit(i)=phi(i)+phil2(i);
   FT(i)=FT(i)+hm27(phit(i)/l(2)/d/SF)*lm(2)+ ...
      hm27(phit(i)/l(1)/d/SF)*lm(1);
end
if N*I > max(FT); disp('CURRENT TOO HIGH'); end 
% Specified point analysis
flxp=interp1(FT, phit, N*I);
F=N*I-hm27(flxp/l(2)/d/SF)*lm(2)-hm27(flxp/l(1)/d/SF)*lm(1);
flxl1=F*Pl1; F=F-hm27((flxp-flxl1)/l(2)/d/SF)*lm(2);
flxg1=interp1(Fc, phic, F); phiR=flxp-flxl1-flxg1;
F=F-hm27(phiR/l(4)/d/SF)*lm(4); flxl2=F*Pl2; flxg2=phiR-flxl2;
disp([blanks(4) 'DOUBLE WINDOW MAGNETIC CIRCUIT, With Leakage']);
disp(' '); disp([blanks(8) ' COIL CURRENT = ', num2str(I)]);
disp(' '); disp([blanks(8) '    COIL FLUX = ', num2str(flxp)]);
disp(' '); disp([blanks(7) 'AIR GAP 1 FLUX = ', num2str(flxg1)]);
disp(' '); disp([blanks(7) 'AIR GAP 2 FLUX = ', num2str(flxg2)]);
plot(FT,phit,FT,phil1,'--',FT,phil2,'-.',N*I,flxp,'o'); grid
title('Double window magnetic circuit with leakage');
xlabel('Coil mmf, A-t'); ylabel('Flux, Wb');
legend('Coil flux','Leakage 1','Leakage 2', 'Specified pt.', 0);