function T = Ktemp(e,Tref)
% Ttemp  Temperature of K-type thermocouple
%
% Synopsis:  T = Ktemp(e)
%            T = Ktemp(e,Tref)
%
% Input:  e    = emf of thermocouple junction, Volts:  -5.891e-3 <= e <= 54.866e-3
%         Tref = (optional) temperature of reference junction
%                Default:  Tref = 0 C (ice point)
%
% Output:  T = temperature, degrees C

if nargin<2
  T = TofKtype(e);  return;
elseif nargin==2
  eref = emfKtype(Tref);     %  emf of ice-point referenced T-couple at Tref
  T = TofKtype(e+eref);      %  correct for nonzero reference temperature
  return;
else
  error('only one or two inputs allowed');
end

% ==============================================================
function T = TofKtype(e)
% TofKtype  Temperature of K type thermocouple with ice-point compensation
%
% Synopsis:  T = TofKtype(e)
%
% Input:   e = (vector or scalar) emf of thermocouple junction in Volts
%              T-couple is assumed to have a reference junction at 0 C
%
% Output:  T = temperature, degrees C
em = e*1e3;             %  emf in millivolts, possibly a vector
T = zeros(size(em));    %  preallocate so T vector will be contiguous in memory
n = 1:length(em);       %  vector to monitor which em values have been converted
imd = find(0<=em & em<20.644);      %  vector of indices such that 0 <= e(i) <= 20.644
ihi = find(20.644<=em & em<54.866); %   20.644 <= e(i) <= 54.866
ilo = find(-5.891<=em & em<0);      %   -5.891 <= e(i) <= 0
if ~isempty(ihi)
  c = [-3.110810e-8  8.802193e-6  -9.650715e-4  5.464731e-2  -1.646031 ...
       48.30222  -131.8058];
  T(ihi) = polyval(c,em(ihi));
  n(ihi) = 0;   %  these i values are done, set elements in monitor vector to zero
end
if ~isempty(imd)
  c = [-1.052755e-8   1.057734e-6  -4.413030e-5  9.804036e-4  -1.228034e-2 ...
        8.315270e-2  -2.503131e-1   7.860106e-2  25.08355  0];
  T(imd) = polyval(c,em(imd));
  n(imd) = 0;   %  these i values are done, set elements in monitor vector to zero
end
if ~isempty(ilo)
  c = [-5.1920577e-4  -1.0450598e-2  -8.6632643e-2  -3.7342377e-1  -8.9773540e-1 ...
       -1.0833638     -1.1662878     25.173462  0];
  T(ilo) = polyval(c,em(ilo));
  n(ilo) = 0;   %  these i values are done, set elements in monitor vector to zero
end
% --- Check to see if all elements of em vector have been converted
if any(n)
  nn = n(find(n~=0));   %  select only the nonzero indices
  warning(sprintf('e = %g value outside of range -5.891 <= e <= 20.644 mV\n',em(nn)));
end