%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  soldes.m - solenoid design program that uses the data
%             provided by soldata.m & calls magckt2.m for
%             calculation of flux linkages-current at each
%             specified value of air gap position.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
soldata,
I=input('   Current for Force Study = ');
% Approximate flux density checks
flx=Bmax*pi*d^2/4;  % Plunger for specified Bmax
A=[A1 A2 A3 A4 A5 A6 A7 A8 A9 A10]; mA=length(A);
figure(1);
for i=1:mA; Bck(i)=flx/A(i); end
plot([1:mA],Bck,'p',[0 mA],[Bmax Bmax]); grid;
ylabel('Flux density, T'); xlabel('Member no.');
X1=(Bmax+min(Bck))/2; X2=0.9*X1;
text(1,X1,'Consider dimension adjustment if flux density');
text(1,X2,'of any member varies significantly from Bmax.');
% Weight, size & current density estimates
J=N*I/((db-dc)/2*lc)/39.37^2/0.65;  % 65% copper fill
WTFe=pi*df^2*le+pi*(df^2-db^2)/4*lc+pi*d^2/4*1.1*lc;
WTFe=WTFe*0.283*39.37^3;
WTCu=0.75*pi*(db^2-dc^2)/4*lc*0.322*39.37^3;
Ac=0.65*(db-dc)/2*lc/N*39.37^2;
Rc=1.025e-06*N*pi*(db+dc)*39.37/2/Ac;
Lc=4e-07*pi*N^2*pi/4*d^2/(0.125/39.37);
disp(' '); disp( '  WEIGHT, SIZE & COIL CHARACTERISTIC ESTIMATES');
disp(' '); disp(['       Core Weight(lbs) = ' num2str(WTFe)]);
disp(' '); disp(['       Coil Weight(lbs) = ' num2str(WTCu)]);
disp(' '); disp(['           Frame OD(in) = ' num2str(df*39.37)]);
disp(' '); disp(['       Frame length(in) = ' num2str(lf*39.37)]);
disp(' '); disp(['  Conductor area(sq in) = ' num2str(Ac)]);
disp(' '); disp(['  Current Density(Apsi) = ' num2str(J)]);
disp(' '); disp(['  Resistance(Ohm,150 C) = ' num2str(Rc)]);
disp(' '); disp(['     Inductance(H,1/8") = ' num2str(Lc)]);
disp(' ');
% Generate a set of flux-current plots for air gap lengths
% specified by x. Enter dimensions in inches.
x=[1e-12 0.1 0.2 0.3 0.4 0.5 1.0]/39.37; m=length(x);
figure(2);
for i=1:m
   LamI=magckt2(x(i)); n=length(LamI);
   plot(LamI(2,1:n),LamI(1,1:n)/N); grid on; hold on
   xlabel('Coil current, A'); ylabel('Flux, Wb');
   title('Cylindrical plunger solenoid');
end
% Force study for specified current I
for i=1:m
   % Generate pair of flux linkage - current plots
   del=x(i)/100; LamI=magckt2(x(i)-del); n=length(LamI);
   Lam1(1:n)=LamI(1,1:n); I1(1:n)=LamI(2,1:n);
   LamI=magckt2(x(i)+del);
   Lam2(1:n)=LamI(1,1:n); I2(1:n)=LamI(2,1:n);
   if (I>max(I1)) | (I>max(I2));
      disp(' ');disp('INVALID CURRENT');
   end
   Wf1=0; Wf2=0;
   for j=2:n; % Calculate force
      if (I1(j)>I) | (I2(j)>I)
         Wf1=Wf1+(interp1(I1,Lam1,I)+Lam1(j-1))/2*(I-I1(j-1));
         Wf2=Wf2+(interp1(I2,Lam2,I)+Lam2(j-1))/2*(I-I2(j-1));
         break
      else
         Wf1=Wf1+(Lam1(j)+Lam1(j-1))/2*(I1(j)-I1(j-1));
         Wf2=Wf2+(Lam2(j)+Lam2(j-1))/2*(I2(j)-I2(j-1));
      end
   end
   Fd(i)=(Wf1-Wf2)/2/del;
end
figure(3); plot(x*39.37,Fd*0.2248); grid
title('Force vs. displacement');
xlabel('Displacement, in'); ylabel('Force, lb');