% Function to interpolate data from thermal property tables 
function [varargout] = prop(varargin)
  J=length(varargin);
  if ischar(varargin{J})             % with units
    [mat,T,Units]=varargin{J-2:J};
    switch Units
      case 'C'                       % Celsius
        T=T+273.15;
      case 'F'                       % Fahrenheit
        T=(T+459.67)*5/9;
      case 'K'                       % Kelvin
      case 'R'                       % Rankine
        T=T*5/9;
      otherwise
        warning = 'prop: Unknown temperature units'
    end
    nout = J-3;
  else                               %  Units=='K'
    [mat,T]=varargin{J-1:J};
    nout = J-2;
  end

  for n=1:nout
    switch varargin{n}               % flag
      case 'k'                       % thermal conductivity
        val = get_prop(2,mat,T);
      case 'Dm'                      % mass diffusivity
        val = get_prop(2,mat,T);
      case 'Psat'                    % saturated pressure
        val = get_prop(2,mat,T);
      case 'den'                     % density
        val = get_prop(3,mat,T);
      case 'rho'                     % density
        val = get_prop(3,mat,T);
      case 'hlg'                     % heat of vaporization
        val = get_prop(3,mat,T);
      case 'Cp'                      % specific heat at constant pressure
        val = get_prop(4,mat,T);
      case 'Cv'                      % specific heat at constant volume
        val = get_prop(5,mat,T);
      case 'mu'                      % dynamic viscosity
        val = get_prop(6,mat,T);
      case 'expan'                   % thermal expansion coefficient
        val = get_prop(7,mat,T);
      case 'alpha'                   % thermal diffusivity
        val = get_prop(2,mat,T)./get_prop(3,mat,T)./get_prop(4,mat,T);
      case 'gam'                     % specific heat ratio
        val = get_prop(4,mat,T)./get_prop(5,mat,T);
      case 'nu'                      % kinematic viscosity
        val = get_prop(6,mat,T)./get_prop(3,mat,T);
      case 'Pr'                      % Prandtl number
        val = get_prop(6,mat,T).*get_prop(4,mat,T)./get_prop(2,mat,T);
      case 'Sc'                      % Schmitt number
        val = get_prop(6,mat,T).*get_prop(4,mat,T)./get_prop(2,mat,T);
      case 'R'                       % gas constant
        val = get_prop(4,mat,T)-get_prop(5,mat,T);
      otherwise
        flag = varargin{n}
        error('Unknown flag');
    end
    varargout{n} = val;
  end
end

function val = get_prop(i,mat,T)
  if (mat(1,1)==T)
    val = mat(1,i);
  elseif (T>mat(end,1))
    warning = 'mat: T > Tmax'
    val = mat(end,i);
  elseif (T<mat(1,1))
    warning = 'mat: T < Tmin'
    val = mat(1,i);
  else
    val = interp1(mat(:,1),mat(:,i),T);
  end
end
