%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  dwleak1.m - Analysis of double window magnetic circuit
%              with leakage flux. Calculates coil current 
%              for specified value of air gap 2 flux & plots
%              Flux vs. mmf. (Type 1 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
phig2=0.0030;            % Air gap 2 flux (analysis point) - Wb
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;
if phimax < phig2; disp('SPECIFIED FLUX TOO HIGH');end
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
Fp=phig2/Pg2+hm27(phig2/l(5)/d/SF)*lm(5)+ ...
   hm27(phig2/l(4)/d/SF)*lm(4); phig2l=Fp*Pl1;
Fp=Fp+hm27((phig2l+phig2)/l(4)/d/SF)*lm(4);
if Fp > max(Fc); disp('SPECIFIED POINT NOT POSSIBLE');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
phig1=interp1(Fc, phic, Fp); % Specified point analysis
phig1l=(Fp+hm27((phig2+phig2l+phig1)/l(2)/d/SF)*lm(2))*Pl2;
phip=phig2+phig2l+phig1+phig1l; Fp=interp1(phit,FT,phip);I=Fp/N;
disp([blanks(4) 'DOUBLE WINDOW MAGNETIC CIRCUIT, With Leakage']);
disp(' '); disp([blanks(8) ' COIL CURRENT = ', num2str(I)]);
disp(' '); disp([blanks(8) '    COIL FLUX = ', num2str(phip)]);
disp(' '); disp([blanks(7) 'AIR GAP 1 FLUX = ', num2str(phig1)]);
disp(' '); disp([blanks(7) 'AIR GAP 2 FLUX = ', num2str(phig2)]);
plot(FT,phit,FT,phil1,'--',FT,phil2,'-.',Fp,phip,'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);