% Created 04/2010, Peter Gilbert
%
% Uses correlations from Lemmon and Jacobsen, 2000 Journal of Physical and  
% Chemical Reference Data and Lemmon and Jacobsen 2006, International 
% Journal of Thermophysics
%
% Input: Temperature (deg R), pressure (psi) column vectors
% Output: Dynamic viscosity (Pa-sec), thermal conductivity (W/m-K),
%         and specific heat (J/kg-K) for air in column vectors
%
% Note: Vector lengths must be >= 3 for finite difference calculation
%       of d(rho)/dp
%

function [mu thermcon CP] = airproperties(T, P)

if (length(T) < 3)
    ylen = 1;
    T(2:3,1) = T(1); 
    P(2:3,1) = P(1);
end
    

T = (T-459.67-32)*5/9+273.15; % Convert R to K

gamma = 1.4;
RK = 287.058;   % J/kg-K
Mol = 28.9586;  % Molecular weight for air, g/mol
sig = 0.360;

P = P .* 101325 ./ 14.7;  % Pa
rho = P./(RK*T)/Mol/1000; % Convert to mol/dm^3

%%%%%%%%%%%%%%%%%%%%%%%%
% Coefficients of the Collision Integral Equation (NIST props Table II, III)
k = 0:4;
b = [0.431; -.4623; 0.08406; 0.005341; -.00331]';
Ni = [10.72; 1.122; 0.002019; -8.876; -0.02916]';
ti = [0.2; 0.05; 2.4; 0.6; 3.6]';
di = [1; 4; 9; 1; 8]';
li = [0; 0; 0; 1; 1]';
gami = li;
Tc = 132.6312;      % K
rhoc = 10.4477;     % mol/dm^3
pc = 3.78502;       % MPa
tau = Tc./T;
del = rho./rhoc;
xi0 = 0.11;         % nm
qd = 0.31;          % nm
capgam = 0.055;
TrefK = 265.262;    % K
nu = 0.063;
lowgam = 1.2415;
R0 = 1.01;
bolt = 1.380658e-23;

% Them. con. constants (NIST props Table IV)
Nik = [1.308; 1.405; -1.036; 8.743; 14.76; -16.62; 3.793; -6.142; -.3778];
tik = [0; -1.1; -.3; .1; 0; .5; 2.7; .3; 1.3];
dik = [0 0 0 1 2 3 7 7 11]';
lik = [0 0 0 0 0 2 2 2 2]';
gamik = [0 0 0 0 0 1 1 1 1]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Cp in J/mol-K
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = [3.490888032;
    2.395525583E-6;
    7.172111248E-9;
    -3.115413101E-13;
    0.223806688;
    0.791309509;
    0.212236768;
    0.197938904;
    3364.011;
    2242.45;
    11580.4];

u = N(9)./T;
v = N(10)./T;
w = N(11)./T;

R = 8.314510;   % J/mol-K
cp = zeros(length(T),1); % Initialize cp vector

for m = 1:length(T)
    cp(m,1) = N(1).*T(m,1).^0+N(2).*T(m,1).^1+N(3).*T(m,1).^2+N(4).*T(m,1).^3+...
        N(5)./T(m,1).^(1.5)+N(6).*u(m,1).^2.*exp(u(m,1))./(exp(u(m,1))-1).^2+...
        N(7).*v(m,1).^2.*exp(v(m,1))./(exp(v(m,1))-1).^2+...
        (2/3)*N(8).*w(m,1).^2.*exp(-w(m,1))./((2/3).*exp(-w(m,1))+1).^2;
end
CP = cp.*1000/28.9586*R;  % Convert to J/kg-K
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute Viscosity in microPa-sec and Thermal Conductivity in mW/m-K
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize vectors
eta0 = zeros(length(tau),1);
etar = zeros(length(tau),1);
lam0 = zeros(length(tau),1);
lamr = zeros(length(tau),1);
lamc = zeros(length(tau),1);
drhodp = zeros(length(tau),1);

for n = 2:length(P)-1
    drhodpf = ((rho(n+1)-rho(n)))/((P(n+1)-P(n))*10^(-6));
    drhodpb = ((rho(n)-rho(n-1)))/((P(n)-P(n-1))*10^(-6));
    drhodp(n,1) = (drhodpf+drhodpb)/2;
end
drhodp(1,1) = drhodp(2);
drhodp(n+1,1) = drhodp(n);
if (sum(isnan(drhodp)) ~= 0)
    drhodp = zeros(length(drhodp),1);
end

if (T<=TrefK)
    Tindex = find(T<=TrefK,1,'last');
else
    Tindex = 1;
end
chiref = pc*10^(-6)*rho(Tindex)/rhoc^2*drhodp(Tindex);


for n = 1:length(tau)
    Omeg = exp(sum(b.*(log(T(n)/103.3)).^k,2));
    eta0(n,1) = 0.0266958*sqrt(Mol*T(n))./(sig^2.*Omeg);
    etar(n,1) = sum(Ni.*tau(n).^ti.*del(n).^di.*exp(-gami.*del(n).^li),2);
end
mu = (eta0+etar).*10^(-6);  % Convert to Pa-sec

for n = 1:length(tau)
    xiparam = ((pc*10^(-6)*rho(n)/rhoc^2*drhodp(n))-chiref*TrefK/T(n))/capgam; 
    xi = xi0*(((pc*10^(-6)*rho(n)/rhoc^2*drhodp(n))-chiref*TrefK/T(n))/capgam)^(nu/lowgam);
    Omegt = 2/pi*((1-1/gamma)*atan(xi/qd)+(1/gamma)*(xi/qd));
    Omegt0 = 2/pi*(1-exp(-1./(qd/xi+(1/3)*(xi/qd)^2.*(1./del(n)).^2)));
    lam0(n,1) = Nik(1)*eta0(n)+Nik(2)*tau(n).^tik(2)+Nik(3)*tau(n).^tik(3);
    lamr(n,1) = sum(Nik(4:end).*tau(n).^tik(4:end).*del(n).^dik(4:end).*exp(-gamik(4:end).*del(n).^lik(4:end)),1);
    if (xiparam <= 0)
        lamc(n,1) = 0;
    else
        lamc(n,1) = rho(n).*cp(n)*bolt*R0*T(n)*(Omegt-Omegt0)/(6*pi*xi*mu(n));
    end
end
lam = lam0+lamr+lamc; % mW/m-K
thermcon = lam/1000;  % W/m-K

if (ylen == 1)
    mu = mu(1);
    thermcon = thermcon(1);
    CP = CP(1);
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%







