%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  dwckt1.m - Analysis of double window magnetic circuit.
%             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);
% Generate flux vs. Fc & Fr arrays (center & right legs)
phimax=2.5*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)+ ...
      2*hm27(phir(i)/l(4)/d/SF)*lm(4);
end
phimax=2.5*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=interp1(phir, Fr, phig2);
if Fp>max(Fc); disp('SPECIFIED POINT NOT POSSIBLE');end
% Generate coil flux(phi) 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, phir, F(i));
   FT(i)=F(i)+2*hm27(phi(i)/l(2)/d/SF)*lm(2)+ ...
      hm27(phi(i)/l(1)/d/SF)*lm(1);
end
% Specified point analysis
phig1=interp1(Fc, phic, Fp); phip=phig1+phig2;
Fp=interp1(phi, FT, phip);I=Fp/N;
disp([blanks(4) 'DOUBLE WINDOW MAGNETIC CIRCUIT, No 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,phi,Fr,phir,'--',Fc,phic,'-.',Fp,phip,'o'); grid
title('Double window magnetic circuit');
xlabel('Coil mmf, A-t'); ylabel('Flux, Wb');
legend('Coil flux','Right flux','Center flux', 'Specified pt.', 4);