%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    HM610.m  -  B-H interpolation routine
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = HM610(Bx)
% B-H values that follow are valid for M-6, 29 ga, 10 deg 
B=[0 38.7 45.2 52 58 64.5 71 77 84 90 97 100 103 106 ...
      110 113 116 119 123 126 129 135 142 ...
      200]*1000; % Klines/sq_in
H=[0 0.38 0.42 0.46 0.53 0.6 0.68 0.78 0.93 1.19 1.88 ...
      2.73 5.25 10.5 19.4 32.7 50.5 76.8 111 172 323 ...
      1576 8e3 1e9];  % A_t/in
% Activate to plot B-H curve
% m=21; plot(H(1:m),B(1:m)); grid; % Linear plot
% m=24; semilogx(H(2:m),B(2:m));grid; % Semilog plot
n=length(B); k=0;
if Bx==0; k=-1; y=0; end
if Bx<0; k=-1; y=0; disp('WARNING - Bx < 0, Hx = 0 returned'); end
if Bx>B(n); y=H(n); k=-1; disp('CAUTION - Beyond B-H curve'); end
for i=1:n
if k==0 & (Bx-B(i))<=0; k=i; break; end
end
if k>0;
y=H(k-1)+(Bx-B(k-1))/(B(k)-B(k-1))*(H(k)-H(k-1));
else;
end