%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    Hcs.m  -  B-H interpolation routine
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = Hx(Bx)
% B-H values that follow are valid for Cast Steel
Bcs=[0 0.04 0.1 0.2 0.25 0.3 0.5 0.73 1.15 1.38 1.51 1.63 1.8];
Hcs=[0 50 90 140 170 182 275 400 1000 2000 3000 6000 20000];
B=linspace(0,max(Bcs),200); H=spline(Bcs,Hcs,B);
% Extending beyond range that spline yields monotonic curve
B=[B 1.85 1.9]; H=[H 40000 150000]; m=length(B);
% Activate to plot B-H curve
% plot(H(1:m),B(1:m)); grid; pause; % Linear plot
% semilogx(H(2:m),B(2:m));grid; pause; % 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)+(Bx-B(n))/(pi*4e-07); % Use after complete saturation
   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