classdef Diffusivity
    %DIFFUSIVITY Summary of this class goes here
    %   Detailed explanation goes here

    properties
    end

    methods(Static)
        function [results] = getParams(specie)
            switch(specie)
                % epsK = epsilon/k, Kelvins
                % sigma, Angstroms
                % Values from Poling pp. B.1-B.2, unless otherwise stated

                case 'H2'
                    sigma = 2.827;
                    epsK = 59.7;
                case 'CO'
                    sigma =  3.690;
                    epsK = 91.7;
                case 'CO2'
                    sigma = 3.941;
                    epsK = 195.2;
                case 'H2S'
                    sigma = 3.623;
                    epsK = 301.1;
                case 'SO2'
                    sigma = 4.112;
                    epsK = 335.4;
                case 'N2'
                    sigma = 3.798;
                    epsK = 71.4;
                case 'H2O'
                    sigma = 2.641;
                    epsK = 809.1;
                case 'COS'
                    sigma = 4.130;
                    epsK = 336.0;
                case 'CS2'
                    sigma = 4.483;
                    epsK = 467;
                case 'NH3'
                    sigma = 2.900;
                    epsK = 558.3;
                case 'S6'
                    % Values from Meyer
                    Tc = 1313; % K
                    Pc = (179.7+200)/2; % atm
                    % Estimations from Welty, p409-410
                    sigma = 2.44 * (Tc/Pc)^(1/3);
                    epsK = 0.77*Tc;
            end

            results = [sigma, epsK];
        end

        function dba = D_BA(specieA, specieB, T, P)
            M_A = molarMass.getMolarMass(specieA);
            M_B = molarMass.getMolarMass(specieB);

            B = 10^-4*(10.85 - 2.50*(1/M_A+1/M_B)^0.5);

            [sigmaA, epsK_A] = Diffusivity.getParams(specieA);
            [sigmaB, epsK_B] = Diffusivity.getParams(specieB);
            kappa = 1.38064852e-23; % Boltzman constant,  m^2.kg/(s^2.K)
            epsilonA = epsK_A * kappa;
            epsilonB = epsK_B * kappa;

            epsilon_AB = (epsilonA*epsilonB)^0.5;
            Omega = kappa/epsilon_AB * T;

            sigma_AB = (sigmaA*sigmaB)^0.5;

            dba = B*T^1.5*(1/M_A+1/M_B)^0.5/(P*sigma_AB*omega);
        end

        function dmix = D_Amix(fractionsArray, specieA, T, P)
            % fractionsArray = y_H2, y_CO, y_CO2, y_H2S, y_SO2, y_H2O, y_S6, y_COS, y_CS2, y_NH3, y_N2
            fractionsArrayNames = ['H2 '; 'CO '; 'CO2'; 'H2S'; 'SO2'; 'H2O'; 'S6 '; 'COS'; 'CS2'; 'NH3'; 'N2 '];
            fractionsArrayNames_cellstr = cellstr(fractionsArrayNames);

            specieA_index = strfind(fractionsArrayNames_cellstr, specieA);
            for i=1:length(fractionsArray)
               if specieA_index{i}==1
                   specieA_index=i;
                   break
               end
            end

            % Blanc's Law, Poling p 11.20
            Dam = 0;

            for j = 1:length(fractionsArray)
               if j ~= specieA_index
                   Dam = fractionsArray(j) / Diffusivity.D_BA(specieA, fractionsArrayNames_cellstr{j}, T, P);
               end
            end

            dmix=1/Dam;
        end
    end

end