%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  coenFe.m - evaluates developed force Fd for a magnetic
%             circuit for a specified current and air gap
%             value. Calls magckt1.m for flux linkage -
%             current array.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
x=input(' Air gap length = ');
if x==0; x=1e-06; else; end; del=0.01*x;
I=input('Current for Fd = ');
% Generate pair of flux linkage-current plots
LamI=magckt1(x-del); m=length(LamI);
Lam1(1:m)=LamI(1,1:m); I1(1:m)=LamI(2,1:m);
LamI=magckt1(x+del);
Lam2(1:m)=LamI(1,1:m); I2(1:m)=LamI(2,1:m);
if (I>max(I1)) | (I>max(I2)); disp(' ');disp('INVALID CURRENT'); end
Wf1=0; Wf2=0;
for i=2:m; % Calculate force
   if (I1(i)>I) | (I2(i)>I)
      Wf1=Wf1+(interp1(I1,Lam1,I)+Lam1(i-1))/2*(I-I1(i-1));
      Wf2=Wf2+(interp1(I2,Lam2,I)+Lam2(i-1))/2*(I-I2(i-1));
      break
   else
      Wf1=Wf1+(Lam1(i)+Lam1(i-1))/2*(I1(i)-I1(i-1));
      Wf2=Wf2+(Lam2(i)+Lam2(i-1))/2*(I2(i)-I2(i-1));
   end
end
   Fd=(Wf1-Wf2)/2/del;
disp(' '); disp([blanks(8) 'COENERGY FORCE ANALYSIS']);
disp(' '); disp([blanks(12) ['Air gap = ',num2str(x)]]);
disp(' '); disp([blanks(12) ['Current = ',num2str(I)]]);
disp(' '); disp([blanks(12) ['     Fd = ',num2str(Fd)]]);
% Activate to check lambda-i plots
% plot(I1,Lam1,I2,Lam2); grid;
% xlabel('Coil current, A'); ylabel('Flux linkage, Wb-t');