function [Tout, Pout, rhoOut, gammaOut, VsoundOut, Mout, noEvapOut, Urs, X, phi, aeroLoadFrac, dropVolFrac, absCrossR5] = aerofrosh(varargin)
% 
% When using this code for research, please cite:
% M. F. Campbell, D. R. Haylett, D. F. Davidson, and R. K. Hanson. 
% AEROFROSH: A shock condition calculator for multi-component-fuel aerosol-laden flows. 
% Shock Waves 26:4 (2016) 429-447.
% 
% "AEROsol FROsen SHock solver"
% PURPOSE: 
%   * Calculate conditions behind incident and reflected shocks in the
%       aerosol shock tube
% INPUTS:
%   * Uis = incident shock velocity at endwall (m/s)
%   * T1 = temperature in region 1 (K)
%   * P1 = pressure in region 1 (torr)
%   * absorbR2 = absorbance (-ln(I/Io) in region 2 @ wavelength of absCross
%   * absorbR5 = absorbance (-ln(I/Io) in region 5 @ wavelength of absCross - NOT REQUIRED
%   * bathGasSpec = cell array of strings of gases in bath gas
%         e.g.: {'AR','N2','O2'}
%   * XbathGas = vector of mole fractions of gases in bath gas
%         e.g.: [0.01,0.79,0.20]
%   * fuelSpec = cell array of strings of liquids in fuel mixture
%         e.g.: {'C11H22O2','C12H26'}
%   * Xfuel = vector of mole fractions of liquids in fuel mixture
%         e.g.: [0.995, 0.005]
%   * vpfVals = vector containing vapor pressures of fuel components @ T1 (Torr)
%   * rfVals = vector containing liquid-phase densities of fuel components @ T1,P1 (g/cm3)
%   * hfVals = vector containing liquid-phase enthalpies of fuel components @ T1,P1 (J/kg)
%   * absCrossPolyMatrix = Matrix containing 4th order polynomial of
%         absorption cross sections for fuels.  Each row in the matrix
%         represents one component of the fuel mixture.  
%         Units produced by evaluation of polynomial should be m^2/mol.  
%         e.g. [x1*T2^3 + x2*T2^2 + x3*T2 + x4   --> C11H22O2
%               x1*T2^3 + x2*T2^2 + x3*T2 + x4]; --> C12H26
%   * nasaFile = string name of *.dat file containing NASA polynomial 
%         coefficients for molecules in shock tube gas mixture.  
%         Should be in the same folder as "aerofrosh.m" function 
%         e.g.: 'sandia.dat'
%   * Ltube = optical path length through shock tube (shock tube width) (meters)
%         e.g.: 0.1000
% OUTPUTS:
%   * Tout = [T2 T5]
%        * T2 = temperature in region 2 (K)
%        * T5 = temperature in region 5 (K)  
%   * Pout = [P2 P5]
%        * P2 = pressure in region 2 (Pa)   
%        * P5 = pressure in region 5 (Pa)
%   * rhoOut = [rho1 rho2 rho5]
%        * rho1 = density in region 1 (kg/m3)
%        * rho2 = density in region 2 (kg/m3)
%        * rho5 = density in region 5 (kg/m3)
%   * gammaOut = [gamma1 gamma2 gamma5]
%        * gamma1 = cp1/cv1 in region 1
%        * gamma2 = cp2/cv2 in region 2
%        * gamma5 = cp5/cv5 in region 5
%   * VsoundOut = [Vsound1 Vsound2 Vsound5]
%        * Vsound1 = speed of sound in region 1 (m/s)
%        * Vsound2 = speed of sound in region 2 (m/s)
%        * Vsound5 = speed of sound in region 5 (m/s)
%   * Mout = [Mis Mrs]
%        * Mis = Mach number of incident shock
%        * Mrs = Mach number of reflected shock
%   * noEvapOut = [P2noEvap T2noEvap rho2noEvap]
%        * P2noEvap = pressure without droplet evap in region 2 (Pa)
%        * T2noEvap = temperaure without droplet evap in region 2 (K)
%        * rho2noEvap = density without droplet evap in region 2 (kg/m3)
%   * Urs = velocity of reflected shock at endwall (m/s)
%   * X = fully evaporated gas mole fractions
%   * phi = equivalence ratio of mixture
%   * aeroLoadFrac = molar percent of fuel loading from aerosol droplets
%   * dropVolFrac = fraction of volume in region 1 occupied by droplets
%   * absCrossR5 = absorption ((Io-I)/Io) cross section of fuel mixture 
%      (m^2/mol) at T5, assuming complete evaporation in Region 2.  
%      To verify complete evaporation in Region 2, a visible laser should 
%      be directed through the tube.  Its Mie Scattering will go to zero 
%      as the droplets evaporate.
% DEVELOPMENT HISTORY: 
%   * Originally written by Dr. David Davidson using Mie Scattering in 
%       Region 1 for fuel loading measurement.
%   * Updated by Daniel Haylett to use laser absorption in Region 2 
%       for fuel loading measurement.
%   * Improved and transcribed into MATLAB code by Matthew Campbell.
% VERSION NUMBER:
%   * 0.5: 11/01/2010 - Demo release
%   * 1.0: 05/01/2011 - Initial release; Mis accounts for droplets
%   * 1.1: 05/28/2012 - Check for convergence
%   * 1.2: 05/29/2012 - Changed criteria for viable shock from Mrs to Urs
%   * 1.3: 06/02/2012 - Better SANDIA.DAT file reading routine
%   * 1.4: 05/11/2016 - Fix phi calculation. Change only affects oxygenated fuels.

%% PARSE INPUT 
    
% where are we?
    fprintf('Parsing Input\n')
% Grab information from "varargin"
    for i=1:2:length(varargin)
        switch varargin{i}
            case 'Uis'
                Uis = varargin{i+1}; % m/s
            case 'P1'
                P1 = varargin{i+1}; % torr
            case 'T1'
                T1 = varargin{i+1}; % K
            case 'absorbR2'
                absorbR2 = varargin{i+1}; % -ln(I/Io) in Region 2 (just before reflected shock)
            case 'absorbR5'
                absorbR5 = varargin{i+1}; % -ln(I/Io) in Region 5 (just after reflected shock)
            case 'bathGasSpec'
                bathGasSpec = varargin{i+1}; % cell array of species in bath gas
            case 'XbathGas'
                XbathGas = varargin{i+1}; % vector of bath gas mole fractions
            case 'fuelSpec'
                fuelSpec = varargin{i+1}; % cell array of species in fuel mixture
            case 'Xfuel'
                Xfuel = varargin{i+1}; % vector of fuel component mole fractions
            case 'vpfVals'
                vpfVals = varargin{i+1}; % vector of fuel component vapor pressures @ T1 (Torr)
            case 'rfVals'
                rfVals = varargin{i+1}; % vector of fuel component liquid-phase densities @ T1,P1 (g/cm3)
            case 'hfVals'
                hfVals = varargin{i+1}; % vector of fuel component liquid-phase enthalpies @ T1,P1 (J/kg)
            case 'absCrossPolyMatrix'
                absCrossPolyMatrix = varargin{i+1}; % matrix of absorption cross section polynomial; gives units of m^2/mol
            case 'nasaFile'
                nasaFile = varargin{i+1}; % string name of *.dat file containing NASA polynomial coefficients
            case 'Ltube'
                Ltube = varargin{i+1}; % optical path length through shock tube in meters
            otherwise
                % do nothing
        end
    end
    
%% CHECK INPUTS

% where are we?
    fprintf('Checking Input\n')
% if vector is column, invert it so that it is a row vector.
    if size(bathGasSpec,1)>size(bathGasSpec,2)
        bathGasSpec=bathGasSpec';
    end
    if size(XbathGas,1)>size(XbathGas,2)
        XbathGas=XbathGas';
    end
    if size(fuelSpec,1)>size(fuelSpec,2)
        fuelSpec=fuelSpec';
    end
    if size(Xfuel,1)>size(Xfuel,2)
        Xfuel=Xfuel';
    end
    if size(vpfVals,1)>size(vpfVals,2)
        vpfVals=vpfVals';
    end
    if size(rfVals,1)>size(rfVals,2)
        rfVals=rfVals';
    end
    if size(hfVals,1)>size(hfVals,2)
        hfVals=hfVals';
    end
% check absorption cross section matrix
    if size(absCrossPolyMatrix,2) ~= 4
        error('Absorption cross section polynomial matrix should have 4 columns and each row should represent one fuel component');
    end
% check gas vectors
    if (size(bathGasSpec,1) ~= size(XbathGas,1)) || (size(bathGasSpec,2) ~= size(XbathGas,2))
        error('Bath gas species and mole fraction vectors must have same size')
    end
% check liquid vectors
    if (size(fuelSpec,1) ~= size(Xfuel,1)) || (size(fuelSpec,2) ~= size(Xfuel,2))
        error('Fuel species and mole fraction vectors must have same size')
    end
% check temperature
    if T1<273.15
        error('Temperature should be in Kelvin, not Celsius!');
    end
% ensure that mole fractions sum to 1.0
    if sum(XbathGas)-1.0 > 0.000001
        error('Bath gas mole fractions do not sum to 1.0');
    end
    if sum(Xfuel)-1.0 > 0.000001
        error('Fuel component mole fractions do not sum to 1.0');
    end

%% INITIAL COMPUTATIONS

% where are we?
    fprintf('Performing Initial Computations\n')
% create a cell array of all the species names (bath gases, then liquids)
    speciesVec = [bathGasSpec fuelSpec];
% make it UPPERCASE
    speciesVec = upper(speciesVec);
% remember how many components there are in the fuel mixture
    numFuelComps = size(fuelSpec,2);
% Transform pressure units
    P1torr = P1; % Torr
    P1 = P1torr*(1/760)*(101325/1); % Torr --> Pa
% The velocity of the test gas relative to the shock in region 1 (u1) is
% equal to the incident shock velocity (Uis)
    u1 = Uis; % m/s
    
%% CONSTANTS

% where are we?
    fprintf('Defining Constants\n')
% elemental molecular weights
    eleMolWts = { % g/mol (cell array)
        'AC'	227
        'AG'	107.8682
        'AL'	26.9815386
        'AM'	243
        'AR'	39.948
        'AS'	74.92160
        'AT'	210
        'AU'	196.966569
        'B'     10.811
        'BA'	137.327
        'BE'	9.012182
        'BH'	262
        'BI'	208.98040
        'BK'	247
        'BR'	79.904
        'C'     12.0107
        'CA'	40.078
        'CD'	112.411
        'CE'	140.116
        'CF'	251
        'CL'	35.453
        'CM'	247
        'CO'	58.933195
        'CR'	51.9961
        'CS'	132.9054519
        'CU'	63.546
        'DB'	262
        'DY'	162.500
        'ER'	167.259
        'ES'	252
        'EU'	151.964
        'F'     18.9984032
        'FE'	55.845
        'FM'	257
        'FR'	223
        'GA'	69.723
        'GD'	157.25
        'GE'	72.64
        'H'     1.00794
        'HE'	4.002602
        'HF'	178.49
        'HG'	200.59
        'HO'	164.93032
        'HS'	265
        'I'     126.90447
        'IN'	114.818
        'IR'	192.217
        'K'     39.0983
        'KR'	83.798
        'LA'	138.90547
        'LI'	6.941
        'LR'	262
        'LU'	174.9668
        'MD'	258
        'MG'	24.305
        'MN'	54.938045
        'MO'	95.96
        'MT'	266
        'N'     14.0067
        'NA'	22.98976928
        'NB'	92.90638
        'ND'	144.242
        'NE'	20.1797
        'NI'	58.6934
        'NO'	259
        'NP'	237
        'O'     15.9994
        'OS'	190.23
        'P'     30.973762
        'PA'	231.03588
        'PB'	207.2
        'PD'	106.42
        'PM'	145
        'PO'	209
        'PR'	140.90765
        'PT'	195.084
        'PU'	244
        'RA'	226
        'RB'	85.4678
        'RE'	186.207
        'RF'	261
        'RH'	102.90550
        'RN'	222
        'RU'	101.07
        'S'     32.065
        'SB'	121.760
        'SC'	44.955912
        'SE'	78.96
        'SG'	263
        'SI'	28.0855
        'SM'	150.36
        'SN'	118.710
        'SR'	87.62
        'TA'	180.94788
        'TB'	158.92535
        'TC'	98
        'TE'	127.60
        'TH'	232.03806
        'TI'	47.867
        'TL'	204.3833
        'TM'	168.93421
        'U'     238.02891
        'V'     50.9415
        'W'     183.84
        'XE'	131.293
        'Y'     88.90585
        'YB'	173.054
        'ZN'	65.38
        'ZR'	91.224};
% universal gas constant
    Ru = 8.314471; % J/mol*K
    
%% READ SANDIA THERMODYNAMIC DATA FILE, COMPUTE SPECIES MOLECULAR WEIGHTS

% where are we?
    fprintf('Reading Sandia Thermo Data\n')
% Read file (see subfunction below)
    [coefVec, nameEle, numEle] = readSandia(speciesVec, nasaFile);
% where are we?
    fprintf('Computing Species Molecular Weights\n')
% Find molecular weight of each species, gas and liquid
    specMolWts = getSpecMolecWts(nameEle, numEle, eleMolWts); % g/mol
    
%% DETERMINE MOLAR AND MASS AVERAGE PROPERTIES FOR FUEL MIXTURE

% where are we?
    fprintf('Computing Average Properties Of Fuel Mixture\n')
% Molecular Weight
    fuelMolWtAvg = specMolWts((end-numFuelComps+1):end)*Xfuel'; % (g/mol,i)*(mol,i/mol,total) = g/mol
% Vapor Pressure
% use Raoult's Law of ideal mixtures: P,a = x,a*Psat,a
% compute the vapor pressure for each component in solution
% then sum these partial pressures up to find the total vapor pressure
% according to Dalton's law of partial pressures.
    vpfPartialPressure = vpfVals.*Xfuel; % Torr
    vpfAvg = sum(vpfPartialPressure); % Torr
% Density
% assume 1 mol of liquid fuel.  the total mass of this 1 mol of liquid fuel
% is the fuel's average molecular weight, and the mass of each component is
% the mole fraction of that component multiplied by its molar mass.  the
% total volume is the mass of each component divided by its density. the
% density is found by dividing total mass by total volume.
    totalMass = fuelMolWtAvg*1; % g/mol * 1 mol = g
    massEachComponent = specMolWts((end-numFuelComps+1):end).*Xfuel; % g/mol, note the ".*" operator for element-by-element multiplication.  
    massFracs = massEachComponent/totalMass; % g/g
    rfAvg = 1/(sum(massFracs./rfVals)); % g/mol
% Specific Enthalpy
% assume 1 mol of liquid fuel.  find the mass fraction of each component.
% then do a mass fraction weighted summation to find the specific enthalpy.
    hfAvg = hfVals*massFracs'; % (J/kg)*(g/g) = J/kg
    
%% DETERMINE GAS PHASE MOLE FRACTIONS IN REGION ONE

% region 1: bath gas + liquid vapor pressure;
% region 2: bath gas + liquid vapor pressure + evaporated droplets;
% region 5: assume to be same as region 2 (assume no fuel decomposition in region 2); 

% where are we?
    fprintf('Determining Region 1 Gas Mole Fracs\n')
% assemble mole fractions, accounting for vapor pressure of fuel mole fraction of gas species
% the measured P1 includes the vapor pressure; therefore the partial pressures of bath gas 
% species are less than simply their bath gas mole fraction times total pressure. 
    X1onlyGas = XbathGas*(P1torr-vpfAvg)/P1torr;
% add in fuel component mole fractions (derived from vapor pressure)
% note that X1onlyGas is the mole fractions IN THE GAS PHASE, NOT the
% overall mole fractions in region 1 (which account for both liquid and gas)
    X1fuelMoleFracs = vpfPartialPressure/P1torr;
    X1onlyGas = [X1onlyGas X1fuelMoleFracs];
% the ratio of vapor-phase fuel components due to vapor pressure is NOT the same as the 
% molar fuel ratio in the liquid.  This composition is computed below, though the program 
% uses another method to solve the mass balance equations.  
    fuelFracsVapor = X1fuelMoleFracs/sum(X1fuelMoleFracs); % for reference only

%% DETERMINE GAS-ONLY PROPERTIES IN REGION ONE

% Weighted sum by the mole fractions in the gas; this is the region 1 value
% multiply then add mole fractions by species molec weights 
    gasMolWtNoDrops = X1onlyGas*specMolWts'; % g/mol
% compute the enthalpy of gas without drops
    [CpNoDrops, hNoDrops] = getCph_T(T1, X1onlyGas, gasMolWtNoDrops, Ru, coefVec, speciesVec); % J/kg*K, J/kg
% compute ideal gas constant for gas-part of mixture in region 1
    Rgas1 = (Ru/gasMolWtNoDrops)*1000; % (J/mol*K)*(mol/g)*(1000 g/1 kg) = J/kg*K
% compute gas specific volume in region 1
    v1noDrops = Rgas1*T1/P1; % (J/kg*K)*K/Pa = m3/kg

%% SOLVE FOR INCIDENT SHOCK CONDITIONS

% At this point we have a vector of bath gas with the vapor pressure of the
% fuel mixed in.  However, we also need to account for the change in mole
% fraction which occurs when the aerosolus fuel evaporates.  This loop uses
% the absorption in Region 2 to compute the new mole fractions.  It also
% contains an inner loop which uses Newton-Raphson to solve for T2 and P2.

% Cm is the mole fraction of moles of molecules in the liquid phase, 
% ie, in the droplets.  Cm = mol_droplets / (mol_droplets + mol_gas). 

% where are we?
    fprintf('Solving for Incident Shock Conditions\n')
% First guess is that the aerosol added no fuel, ie just vapor pressure
% ***In general, Xguess is the mole fractions after complete evaporation***
    Xguess = X1onlyGas; 
% Set Cm initially to be zero, ie there were no moles in liquid phase
    Cm = 0; 
% Iteration parameters for Mole Fraction loop
    maxTimeMF = 10; % second
    maxIterMF = 1000; % iterations
    tolMF = 1e-6; % convergence tolerance
    stopperMF = false;
% Iteration loop
    iterMF = 0; 
    startTimeMF = tic;
    while stopperMF == false;
        % inc the counter
            iterMF = iterMF+1;
        % CONDITIONS FOR MIXTURE
            % compute the molecular weight of the mixture once it is fully evaporated
                gasMolWt = Xguess*specMolWts'; % g/mol
            % gas constant for fully evaporated mixture
            % this is defined in regions 2 and 5 where there is only gas.
            % it is undefined for the total mixture in region 1 where there are both gas and droplets.
            % however, it IS fully defined for the gas-portion of the region 1 mixture (computed above).  
                Rgas2 = (Ru/gasMolWt)*1000; % (J/mol*K)*(mol/g)*(1000 g/1 kg) = J/kg*K
            % Assume 1 mole of liquid and gas combined in the pre-evaporation mixture in region 1.
            % Then the mass and volume of each can be determined.
                massLiq = Cm*1*fuelMolWtAvg*(1/1000); % (mol frac)*(1 mol)*(g/mol)*(1 kg/1000 g) = kg
                massGas = (1-Cm)*1*gasMolWtNoDrops*(1/1000); % kg
                volLiq = Cm*1*fuelMolWtAvg*(1/rfAvg)*(1/1e6); % (mol frac)*(1 mol)*(g/mol)*(cm3/g)*(1 m3/1e6 cm3) = m3
                volGas = (1-Cm)*1*Ru*T1/P1; % mol*(J/mol*K)*K/Pa = m3
            % compute density and specific volume of total mixture in region 1
                rho1 = (massLiq+massGas)/(volLiq+volGas); % kg/m3
                v1 = 1/rho1; % m3/kg
            % Adjust the enthalpy for the gas in region 1 by subtracting out the enthalpy 
            % of vaporization and adding in the enthalpy that the evaporated fuel would have at T1.
            % This is equivalent to h_liq: [h_liq = h_vap - h_fg]
            % Add in the enthalpy via mass fractions
            % Recall that "gasMolWt" is the molecular weight of the gas including evaporated liquid
                massFracLiq = (Cm*fuelMolWtAvg)/gasMolWt; % (molfrac*g/mol)/(molfrac*g/mol) = no units
                massFracGas = ((1-Cm)*gasMolWtNoDrops)/gasMolWt;
                if abs(massFracLiq+massFracGas-1) > 0.000001
                    error('Mass fractions not summing')
                end
                h1 = hNoDrops*massFracGas + hfAvg*massFracLiq; % J/kg
            % The sound speed required is the FROZEN sound speed, discussed in Temkin 2000.
            % It is derived from the gaseous sound speed.     
                % specific heat at constant volume
                CvNoDrops = CpNoDrops-Rgas1; % J/kg*K
                % ratio of specific heats
                gamma1 = CpNoDrops/CvNoDrops;
                % density of gas only, without drops
                rho1noDrops =1/v1noDrops; % kg/m3
                % speed of sound in gas only, without drops
                Vsound1gas = sqrt(gamma1*P1/rho1noDrops); % m/s
                % volume fraction of droplets (droplet volume concentration)
                dropVolFrac = volLiq/(volLiq+volGas);
                % droplet mass concentration
                dropMassConc = dropVolFrac*(rfAvg/rho1noDrops)/(1+(1-rho1noDrops/rfAvg)*dropVolFrac*(rfAvg/rho1noDrops));
                % frozen sound speed (sound speed corrected for droplets)
                Vsound1 = Vsound1gas*(1-dropMassConc)^(1/2)/(1-dropVolFrac);
            % compute mach number of incident shock
                Mis = u1/Vsound1; % no units
        % SOLVE IDEAL SHOCK EQNS 
                % Before Newton-Rhapson guess using ideal relations:
                %     P2/P1 = ((2*gamma1*Mis^2-(gamma1-1))/(gamma1+1))
                %     T2/T1 = (gamma1*Mis^2-(gamma1-1)/2)*((gamma1-1)/2*Mis^2+1)/(((gamma1+1)/2)^2*Mis^2)
                %     gamma = cp1/(cp1-1) --> assumes driven gas is mostly Ar with Cp ~ 2.5
                %     Mis = U1/Vsound1
            % compute ideal pressure in region 2
                P2ideal = P1*((2*gamma1*Mis^2-(gamma1-1))/(gamma1+1));
            % compute ideal temperature in region 2
                T2ideal = T1*(gamma1*Mis^2-(gamma1-1)/2)*((gamma1-1)/2*Mis^2+1)/(((gamma1+1)/2)^2*Mis^2); % K
            % specific volume
                v2ideal = Rgas2*T2ideal/P2ideal; % m3/kg
        % SOLVE NORMAL SHOCK EQNS 
            % give it the ideal values to help pinpoint the solution
                [P2, T2, v2, Cp2, h2] = incidentShockNR(...
                    P2ideal, T2ideal, v2ideal, Rgas2, ...
                    P1, v1, h1, u1, ...
                    Xguess, gasMolWt, coefVec, speciesVec, Ru);
        % ACCOUNT FOR OBSERVED ABSORPTION
            % compute absorption coefficient at T2 for all fuel components
                for i=1:numFuelComps
                    absCrossPoly = absCrossPolyMatrix(i,:); % grab the coefficients in the row
                    crossR2vals(i) = polyval(absCrossPoly,T2); % (m^2/mol) evaluate the polynomial at T2
                end
            % compute the molar average absorption coefficient
            % for the first iteration, base this on the mole fractions in
            % the liquid fuel so if the vapor pressure is ~0 the algorithm
            % works.  After that, base it on the vaporized mole fractions
                if iterMF == 1;
                    fuelEvapMoleFracs = Xfuel;
                else
                    tempTot = sum(Xguess((end-numFuelComps+1):end)); % sum of mole fractions of fuel
                    fuelEvapMoleFracs = Xguess((end-numFuelComps+1):end)/tempTot; % divide to find fuel component mole fractions
                end
                crossR2avg = crossR2vals*fuelEvapMoleFracs'; % m^2/mol
            % compute molar concentration of fuel (this is for all evaporated fuel combined)
                nFuel = absorbR2/(crossR2avg*Ltube); % mol/m^3 =(no units)/(m^2/mol*m)
            % compute molar concentration of gas in shock tube at T2, P2
                nTot2 = P2/(Ru*T2); % mol/m^3
            % compute mole fraction of fuel in region 2
                XevapFuelNew = nFuel/nTot2;
            % compute the old mole fraction of fuel total 
            % this is the sum of the mole fractions of individual fuel components
                XevapFuelOld = sum(Xguess((end-numFuelComps+1):end));
        % CHECK RESULTS AND DECIDE WHAT TO DO
            % check variable
                checkMF = abs((XevapFuelOld-XevapFuelNew)/XevapFuelNew);
            % decide what to do
                if checkMF < tolMF;
                    stopperMF = true; % exit the loop
                elseif toc(startTimeMF) > maxTimeMF
                    error('Timed out on MF loop')
                elseif iterMF > maxIterMF
                    error('Itered out on NR loop')
                else
                    % update Cm (mole fraction of molecules in liquid)
                        % Assume 1 mol total of gaseous mixture prior to evaporation, 
                        % then realize that the total moles of gas mixture increase only 
                        % due to the addition of evaporated fuel droplets.  
                        XnonFuelPreEvap = sum(X1onlyGas(1:end-numFuelComps)); % = 1-sum(X1onlyGas(end-numFuelComps+1))
                        XnonFuelPostEvap = 1-XevapFuelNew;
                        molAdded = (XnonFuelPreEvap-XnonFuelPostEvap)/(XnonFuelPostEvap);
                        % Cm is the fraction of molecules (liquid and gas
                        % combined) that are in the liquid phase in region 1
                        Cm = molAdded/(1+molAdded); % mol evaporated divided by total mols (1 mol gas mixture + mol evaporated)
                    % solve for new post-evaporation mole fractions
                        % Again assume 1 mol total of gaseous mixture prior to evaporation.
                        % Find the moles of non-fuel bath gas components present.  
                        molesBathGasComponents = X1onlyGas(1:end-numFuelComps)*(1.0); % (mol/mol)*mol = mol  
                        % Find the moles of each fuel component prior to evaporation 
                        % in the gas phase due to their vapor pressures.
                        molesFuelComponentsVapor = X1onlyGas((end-numFuelComps+1):end)*(1.0); % (mol/mol)*mol = mol
                        % Find the moles of each fuel component added via aerosol evaporation
                        % assume that the fuel nebulizes in proportion to its mole fractions
                        % ie mole fractions in the aerosol drops are identical to the mole fractions in the liquid
                        molesFuelComponentsAerosol = Xfuel*molAdded;
                        % we now know how many moles of each species there are, and from which source.
                        % add the vapor and aerosol fuel mole fractions
                        molesFuelComponents = molesFuelComponentsVapor+molesFuelComponentsAerosol;
                        % assemble new mole factions vector
                        Xguess = [molesBathGasComponents molesFuelComponents]/(1+molAdded);
                        % check to make sure mole fractions are consistent
                        if abs(XevapFuelNew-sum(Xguess(end-numFuelComps+1:end))) > 0.000001
                            error('Mole fractions not summing')
                        end
                end
    end
% ensure that mole fraction of fuel did not decrease
% note that minimum mole fraction of fuel is that of its vapor pressure
    if sum(Xguess((end-numFuelComps+1):end)) < sum(X1fuelMoleFracs)
        error('Xfuelabs(%4.3s) < Xfuelvp(%4.3s)',sum(Xguess(end-numFuelComps+1:end)),sum(X1fuelMoleFracs))
    end
% compute gas velocity in region 2
    u2 = u1*(v2/v1); % m/s
% compute gamma in region 2
    Cv2 = Cp2 - Rgas2; % J/kg*K
    gamma2 = Cp2/Cv2;
% compute density in region 2
    rho2 = 1/v2; % kg/m3
% assign mole fractions
    X = Xguess;
% compute speed of sound in region 2
    Vsound2 = sqrt(gamma2*P2/rho2); % m/s
% check to make sure shock conservations equations are solved
    tolCheck = 1e-6;
    eq1 = ((rho1*u1)-(rho2*u2))/(rho1*u1); % (kg/m3)*(m/s)
    eq2 = ((P1+rho1*u1^2)-(P2+rho2*u2^2))/(P1+rho1*u1^2); % Pa + (kg/m3)*(m2/s2)*(1 Pa / 1 kg/ms2)
    eq3 = ((h1+(1/2)*u1^2)-(h2+(1/2)*u2^2))/(h1+(1/2)*u1^2); % J/kg*(1 kgm2/s2 / 1 J) + m2/s2
    if abs(eq1)>tolCheck
        error('Continuity not holding over incident shock')
    elseif abs(eq2)>tolCheck
        error('Momentum not holding over incident shock')
    elseif abs(eq3)>tolCheck
        error('Energy not holding over incident shock')
    end
    
%% SOLVE INCIDENT SHOCK CONDITIONS WITHOUT EVAPORATION

% where are we?
    fprintf('Solving for Incident Shock with No Evaporation\n')
% It is often helpful to know the shock conditions in Region 2 prior to
% evaporation in order to estimate the effect that the aerosol had on the
% shock wave.  None of the values computed in this section are necessary
% for the computation of the reflected shock conditions.  
% Assume that the aerosol is not present and essentially compute the
% properties using only the test gas mixture (including fuel vapor)
% Initial guess:
%   * Use P2ideal, T2ideal, v2ideal, because these will be closer to
%   frozen shock values (see Guha 1992)
% Gas properties:
%   * Use Rgas1, P1, v1noDrops, h1noDrops, X1onlyGas, gasMolWtNoDrops
    [P2noEvap, T2noEvap, v2noEvap, dummy1, dummy2] = incidentShockNR(...
        P2ideal, T2ideal, v2ideal, Rgas1, ...
        P1, v1noDrops, hNoDrops, u1, ...
        X1onlyGas, gasMolWtNoDrops, coefVec, speciesVec, Ru);
    rho2noEvap = 1/v2noEvap; 

%% SOLVE REFLECTED SHOCK CONDITIONS

% where are we?
    fprintf('Solving for Reflected Shock Conditions\n')
% Before Newton-Rhapson guess using ideal relations:
%     alpha=(gamma+1)/(gamma-1)
%     P5/P2=(alpha+2-P1/P2)/(1+alpha*P1/P2)
%     T5/T1=(P5/P2)*(alpha+P5/P2)/(1+alpha*P5/P2)
%     gamma=cp1/(cp1-1) --> assumes driven gas is mostly Ar with Cp ~ 2.5
% ideal gas constant in region 5 is the same as for region 2 because
% molecular species and mole fractions are the same.
    Rgas5 = Rgas2;    
% compute alpha
    alpha = (gamma2+1)/(gamma2-1);
% compute ideal pressure in region 2
    P5ideal = P2*(alpha+2-P1/P2)/(1+alpha*P1/P2); % Pa
% compute ideal temperature in region 2
    T5ideal = T2*(P5ideal/P2)*(alpha+P5ideal/P2)/(1+alpha*P5ideal/P2); % K
% specific volume
    v5ideal = Rgas5*T5ideal/P5ideal; % m3/kg
% Solve the Reflected Shock Equations  
    [P5, T5, v5, Cp5, h5] = reflectedShockNR(...
        P5ideal, T5ideal, v5ideal, Rgas5, ...
        P2, v2, h2, u2, u1, ...
        X, gasMolWt, coefVec, speciesVec, Ru);
% compute density in region 5
    rho5 = 1/v5; % kg/m3
% compute gas velocity in region 5
% given the assumption that the gas is stagnated after the reflected shock
% in region 5, u5 (the velocity of the test gas in region 5 relative to the 
% reflected shock) is equal and opposite to the velocity of the reflected 
% shock (Urs).
    u5 = v5*(u1-u2)/(v2-v5); % m/s
    Urs = u5; % m/s
% compute gamma in region 5
    Cv5 = Cp5 - Rgas5; % J/kg*K
    gamma5 = Cp5/Cv5;
% compute density in region 5
    rho5 = 1/v5; % kg/m3
% compute speed of sound in region 5
    Vsound5 = sqrt(gamma5*P5/rho5); % m/s
% Mach number of reflected shock
    u2dash = u5+u1-u2; % m/s
    Mrs = u2dash/Vsound2; % no units
% check answers
    % T5 will always be greater than T2, no mater how much evaporation
    % P5 will always be greater than P2
    % Urs will always be less than Uis
    % however, Mrs may not be less than Mis; this is because as droplet
    % volume fraction increases, T2, R2 (often), and gamma2 decrease, 
    % leading to a decrease in Vsound2 and an increase in Mrs.  
    if (T5<=T2) || (P5<=P2) || (Urs>=Uis)
        error('Impossible reflected shock conditions')
    end
% check to make sure shock conservations equations are solved
    tolCheck = 1e-6;
    eq1 = ((rho2*u2dash)-(rho5*u5))/(rho2*u2dash); % (kg/m3)*(m/s)
    eq2 = ((P2+rho2*u2dash^2)-(P5+rho5*u5^2))/(P2+rho2*u2dash^2); % Pa + (kg/m3)*(m2/s2)*(1 Pa / 1 kg/ms2)
    eq3 = ((h2+(1/2)*u2dash^2)-(h5+(1/2)*u5^2))/(h2+(1/2)*u2dash^2); % J/kg*(1 kgm2/s2 / 1 J) + m2/s2
    if abs(eq1)>tolCheck
        error('Continuity not holding over reflected shock')
    elseif abs(eq2)>tolCheck
        error('Momentum not holding over reflected shock')
    elseif abs(eq3)>tolCheck
        error('Energy not holding over reflected shock')
    end
    
%% COMPUTE ALF, DVF, PHI
    
% where are we?
    fprintf('Solving for ALF, DVF, and Phi\n')
% compute molar fraction of fuel loading due to aerosol
% Cm = mole fraction just liquid fuel
% sum(X(...)) = total mole fraction of fuel (gas + liquid)
    aeroLoadFrac = Cm/sum(X((end-numFuelComps+1):end));
% function to compute stoichiometry (see below)
    phi = getStoich(X, speciesVec, nameEle, numEle);
    
%% COMPUTE ABSORPTION CROSS SECTION IN REGION FIVE

% where are we?
    fprintf('Solving for Absorption Coefficient in Region Five\n')
% decide whether to do this part or not
    if exist('absorbR5','var') % if the user gave an absorption value in R5
        % Assume that we have complete evaporation in R2 (this 
        % can be verified by seeing Mie scattering of a visible laser go to
        % zero in Region 2).  Then, we can compute the absorption cross 
        % section for T5 using this mole fraction.  Note that the accuracy
        % of the R5 absorption cross section depends on the accuracy of the
        % R2 absorption cross section! 
            % compute molar concentration in shock tube at T5, P5
                nTotR5 = P5/(Ru*T5); % mol/m3
            % Grab the fuel mole fraction in R5 computed via iteration
                XfuelR5iter = sum(X((end-numFuelComps+1):end));
            % compute fuel particle density in R5 (mole frac * total)
                nFuelR5iter = XfuelR5iter*nTotR5; % mol/m^3
            % compute absorption cross section
            % this is the molar average cross section for the fuel mixture
                absCrossR5 = absorbR5/(nFuelR5iter*Ltube); % m^2/mol; Ltube defined above
    else
        % output an obviously wrong answer
            absCrossR5 = -999;
    end
    
%% ASSEMBLE OUTPUT VECTORS

% where are we?
    fprintf('Assembling Output Vectors\n')
% output vectors
    Tout = [T2 T5]; % K
    Pout = [P2 P5]; % Pa
    rhoOut = [rho1 rho2 rho5]; % kg/m3
    gammaOut = [gamma1 gamma2 gamma5]; 
    VsoundOut = [Vsound1 Vsound2 Vsound5]; % m/s
    Mout = [Mis Mrs]; % m/s
    noEvapOut = [P2noEvap T2noEvap rho2noEvap]; % Pa, K, kg/m3
   
    
    
function specMolWts = getSpecMolecWts(nameEle, numEle, eleMolWts)
% This function computes the species molecular weights

% Find molecular weight of each species (gas and liquid) (g/mol)
    % rows of nameEle represent species
    % columns of nameEle are the elements composing that species
    specMolWts = zeros(1,size(nameEle,1)); % row vector; each entry is mol wt of the species
    gotIt = false; % used to determine if an element was found
    for i=1:size(nameEle,1) % outer loop over all rows of nameEle
        for j=1:size(nameEle,2) % inner loop over all columns of nameEle 
            tempEl = nameEle{i,j};
            if isempty(tempEl)
                continue; % if nothing there, skip it. 
            else
                for k=1:size(eleMolWts,1) %loop over rows of eleMolWts
                    if strcmp(eleMolWts{k,1},tempEl)
                        % add to the total weight; 
                        % the amount to add is the molecular weight of the 
                        % element times the number of that element present 
                        % in the species
                        specMolWts(i) = specMolWts(i)+eleMolWts{k,2}*numEle(i,j);
                        gotIt = true; % we found the element
                        break; % stop looking, we got it;
                    end
                end
                if gotIt == false
                    error('Atomic weight of %1s not found',tempEl);
                end
            end
            gotIt = false; 
        end
    end
    
    
function [Cp, h] = getCph_T(T, X, gasMolWt, Ru, coefVec, speciesVec)
% This finds the specific heat and enthalpy at the specified temperature

% loop over all species to compute for each species first
    HoverRT = zeros(1,size(coefVec,3));
    CpoverR = zeros(1,size(coefVec,3));
    for i=1:size(coefVec,3) % loop over 3rd dimension of coefVec 
        % grab temperature limits
            Tlow = coefVec(3,1,i);
            Tcom = coefVec(3,2,i);
            Thigh = coefVec(3,3,i);
        % check temperature limits
            if T < Tlow
                %fprintf('%3s: T (%1.2f K) < Tlow (%1.0f K)\n',speciesVec{i},T,Tlow);
            elseif T > Thigh
                %fprintf('%3s: T (%1.2f K) > Thigh (%1.0f K)\n',speciesVec{i},T,Thigh);
            end
        % decide which set of data to use
            if T < Tcom % low temp set
                j = 2;
            else % high temp set for T >= Tcom
                j = 1;
            end
        % grab temperature coefficients
            a1 = coefVec(j,1,i);
            a2 = coefVec(j,2,i);
            a3 = coefVec(j,3,i);
            a4 = coefVec(j,4,i);
            a5 = coefVec(j,5,i);
            a6 = coefVec(j,6,i);
            a7 = coefVec(j,7,i);
        % compute Cp/R
            CpoverR(i) = a1 + a2*T + a3*(T^2) + a4*(T^3) + a5*(T^4);
        % compute h/RT
            HoverRT(i) = a1 + a2*T/2 + a3*(T^2)/3 + a4*(T^3)/4 + a5*(T^4)/5 + a6/T;
    end
% now compute mixture properties
    Cp = (X*CpoverR')*Ru*(1000/gasMolWt); % (no units)*(J/mol*K)*(1000 g/kg)*(mol/g) = J/kg*K
    h = (X*HoverRT')*T*Ru*(1000/gasMolWt); % J/kg
    
    
function [P2, T2, v2, Cp2, h2] = incidentShockNR(P2ideal, T2ideal, v2ideal, Rgas2, P1, v1, h1, u1, Xguess, gasMolWt, coefVec, speciesVec, Ru)
% This function solves for the conditions behind the incident shock

% Solve the Incident Shock Equations
%     P*v = Rgas*T
%     rho1*u1 = rho2*u2
%     P1+rho1*u1^2 = P2+rho2*u2^2
%     h1+(1/2)*u1^2 = h2+(1/2)*u2^2
% by minimizing the following pair of equations (f1, f2)
%     f1 = (P2/P1-1)+(u1^2/(P1*v1))*(v2/v1-1)
%     f2 = (h2-h1)/((1/2)*u1^2)+(v2^2/v1^2-1)
% use Newton-Raphson Method

% Assign ideal to first guess
% Note that we do not use the P2 and T2 from the last iteration
% of the mole fractiion loop, because this prevents the new
% Newton-Raphson loop from being ill-conditioned as we changed
% the mole fractions.  A conservative method of solving the
% equations is to change the mole fractions, but start from the
% ideal P2 and T2 each time.  
    P2guess = P2ideal;
    T2guess = T2ideal; 
    v2 = v2ideal;
% NR iteration parameters
    maxTimeNR = 10; % seconds
    maxIterNR = 100; % iterations
    tolNR = 1e-6; % convergence tolerance
    stopperNR = false;
% NR loop
    iterNR = 0; 
    startTimeNR = tic;
    while stopperNR == false;
        % inc the counter
            iterNR = iterNR+1;
        % compute cp and h for gas in region 2 (see subfunction below)
            [Cp2, h2] = getCph_T(T2guess, Xguess, gasMolWt, Ru, coefVec, speciesVec); % J/kg*K, J/kg
        % compute f1
            f1 = (P2guess/P1-1)+(u1^2/(P1*v1))*(v2/v1-1);
        % compute f2
            f2 = (h2-h1)/((1/2)*u1^2)+((v2/v1)^2-1);
        % compute derivatives for jacobian
        % note that dh/dP = 0 (ideal gas)
        % note that dT/dP = 0 (independent variables)
            df1dP2 = (1/P1)+(u1^2/(P1*v1^2))*(-v2/P2guess); % df1/dP2
            df1dT2 = (1/P1)*(0)+(u1^2/(P1*v1^2))*(v2/T2guess);
            df2dP2 = (2/u1^2)*(0)+(2*v2/v1^2)*(-v2/P2guess);
            df2dT2 = (2/u1^2)*(Cp2)+(2*v2/v1^2)*(v2/T2guess);
        % form Jacobian matrix [df1/dP2 df1/dT2; df2/dP2 df2/dT2]
            J = [df1dP2 df1dT2; ...
                 df2dP2 df2dT2];
        % matrix multiplication to solve for new T2 and P2 guesses
        % using matlab backslash operator rather than inverse saves time!
        % general method: [P2new T2new]' = [P2old T2old]' - inv(J)*[f1old f2old]'
            PTold = [P2guess T2guess]';
            fvec = [f1 f2]';
            PTnew = PTold - J\fvec;
            P2new = PTnew(1);
            T2new = PTnew(2);
        % check results - the T2 and P2 we guess should give the same T2,P2
            checkP = abs((P2guess-P2new)/P2new); % percent change from new
            checkT = abs((T2guess-T2new)/T2new);
        % decide whether to keep going or quit
            if (checkP < tolNR) && (checkT < tolNR)
                stopperNR = true; % exit the loop
            elseif toc(startTimeNR) > maxTimeNR;
                error('Timed out on NR loop')
            elseif iterNR > maxIterNR
                error('Itered out on NR loop')
            else
                % reassign the variables for the next iteration
                P2guess = P2new;
                T2guess = T2new;
                % recompute specific volume
                v2 = Rgas2*T2guess/P2guess; % m3/kg
                % prevent the algorithm from going the wrong way
                while P2guess < 0
                    P2guess = mean([P2ideal P2guess]); % average back to ok range 
                end
                while T2guess < 0
                    T2guess = mean([T2ideal T2guess]);
                end
            end
    end
% save results
    P2 = P2guess;
    T2 = T2guess;
    
    
function [P5, T5, v5, Cp5, h5] = reflectedShockNR(P5ideal, T5ideal, v5ideal, Rgas5, P2, v2, h2, u2, u1, X, gasMolWt, coefVec, speciesVec, Ru)
% This function solves for the conditions behind the reflected shock

% Solve the Reflected Shock Equations
%     P*v = Rgas*T
%     rho2*u2dash=rho5*u5
%     P2+rho2*u2dash^2=P5+rho5*u5^2
%     h2+(1/2)*u2dash^2=h5+(1/2)*u5^2
%     u2dash=u5+u1-u2
%     u2dash=(u1-u2)*(rho5/(rho5-rho2))
%     U5=0 (gas relative to lab frame behind reflected shock is stagnant)
% by solving the following pair of equations
%     f3=(P5/P2-1)+(u1-u2)^2/(P2*(v5-v2))
%     f4=(h5-h2)/((1/2)*(u1-u2)^2)+(v5+v2)/(v5-v2)
% use Newton-Raphson Method

% Assign ideal to first guess
    P5guess = P5ideal;
    T5guess = T5ideal; 
    v5 = v5ideal;
% NR iteration parameters
    maxTimeNR = 10; % second
    maxIterNR = 100; % iterations
    tolNR = 1e-6; % convergence tolerance
    stopperNR = false;
% NR loop
    iterNR = 0; 
    startTimeNR = tic;
    while stopperNR == false;
        % inc the counter
            iterNR = iterNR+1;
        % compute cp and h for gas in region 5 (see subfunction below)
            [Cp5, h5] = getCph_T(T5guess, X, gasMolWt, Ru, coefVec, speciesVec); % J/kg*K, J/kg
        % compute f3
            f3 = (P5guess/P2-1)+(u1-u2)^2/(P2*(v5-v2));
        % compute f4
            f4 = (h5-h2)/((1/2)*(u1-u2)^2)+(v5+v2)/(v5-v2);
        % compute derivatives for jacobian
        % note that dh/dP = 0 (ideal gas)
        % note that dT/dP = 0 (independent variables)
            df3dP5 = (1/P2)+(-(u1-u2)^2/(P2*(v5-v2)^2))*(-v5/P5guess); % df3/dP5
            df3dT5 = (1/P2)*(0)+(-(u1-u2)^2/(P2*(v5-v2)^2))*(v5/T5guess);
            df4dP5 = (1/(1/2*(u1-u2)^2))*(0)+(-2*v2/(v5-v2)^2)*(-v5/P5guess);
            df4dT5 = (1/(1/2*(u1-u2)^2))*(Cp5)+(-2*v2/(v5-v2)^2)*(v5/T5guess);
        % form Jacobian matrix [df1/dP2 df1/dT2; df2/dP2 df2/dT2]
            J = [df3dP5 df3dT5; ...
                 df4dP5 df4dT5];
        % matrix multiplication to solve for new T5 and P5 guesses
        % using matlab backslash operator rather than inverse saves time!
        % general method: [P5new T5new]' = [P5old T5old]' - inv(J)*[f3old f4old]'
            PTold = [P5guess T5guess]';
            fvec = [f3 f4]';
            PTnew = PTold - J\fvec;
            P5new = PTnew(1);
            T5new = PTnew(2);
        % check results - the T2 and P2 we guess should give the same T2,P2
            checkP = abs((P5guess-P5new)/P5new); % percent change from new
            checkT = abs((T5guess-T5new)/T5new);
        % decide whether to keep going or quit
            if (checkP < tolNR) && (checkT < tolNR)
                stopperNR = true; % exit the loop
            elseif toc(startTimeNR) > maxTimeNR;
                error('Timed out on NR loop')
            elseif iterNR > maxIterNR
                error('Itered out on NR loop')
            else
                % reassign the variables for the next iteration
                P5guess = P5new;
                T5guess = T5new;
                % recompute specific volume
                v5 = Rgas5*T5guess/P5guess; % m3/kg
                % prevent the algorithm from going the wrong way
                while P5guess < 0
                    P5guess = mean([P5ideal P5guess]); % average back to ok range 
                end
                while T5guess < 0
                    T5guess = mean([T5ideal T5guess]);
                end
            end
    end
% save results
    P5 = P5guess; % Pa
    T5 = T5guess; % K


function phi = getStoich(X, speciesVec, nameEle, numEle)
% This function computes the stoichiometry of the shock
% phi = (nFuel/nO2)_actual / {nFuel/nO2)_stoich
%     = sum{xFueli*(Nci+(1/4)Nhi+Nsi-(1/2)Noi)}/xO2
% Stoich = C-->CO2, H-->H2O, S-->SO2

% Inert species list: these are not included in the calculation of phi
% N2, AR, etc. need not be listed here because they have no C, H, or S.
    inertReactantSpecies = {'CO2','H2O','SO2'}; 

% Find the index of oxygen in the species vector.  
% If oxygen is not present, or if it is listed more than once, phi = -999
    tempArray = strncmp(speciesVec,'O2',2); % vector of 0's with a 1 where O2 is.
    if sum(tempArray)==0 % oxygen not listed
        phi = -999;
        return
    elseif sum(tempArray)>1 % oxygen listed more than once
        fprintf('Error: O2 listed more than once in test gas species vector\n')
        phi = -999;
        return
    end
    oxygenIndex = find(tempArray); % gives index where 1 is present
    
% Compute numerator of phi fraction
    phiNumerator = 0;
    for i=1:length(X)
        % check if this is an inert species
        speciesName = speciesVec{i};
        tempArray = strncmp(inertReactantSpecies,speciesName,length(speciesName));
        if sum(tempArray)>0 % this species is inert and does not contribute to phi
            continue
        end
        % check if this is O2
        if i==oxygenIndex % oxygen is not considered to be a fuel, it is an oxidizer
            continue
        end
        % determine number of C, H, S, and O atoms in this species
        nC = 0;
        nH = 0;
        nS = 0;
        nO = 0;
        for j=1:size(nameEle,2) % looping over columns of nameEle (rows are species, columns are elements)
            % grab the name of the element (ie, 'H', 'C', 'AR', 'HE')
            tempEl = nameEle{i,j};
            % decide what to do
            if isempty(tempEl)
                continue; % if nothing there, skip it. 
            elseif strcmp(tempEl,'C')
                nC = numEle(i,j);
            elseif strcmp(tempEl,'H')
                nH = numEle(i,j);
            elseif strcmp(tempEl,'O')
                nO = numEle(i,j);
            elseif strcmp(tempEl,'S')
                nS = numEle(i,j);
            else
                continue % not one of the atoms we care about ('AR', 'HE', etc.)
            end
        end
        % Compute this part of the numerator and add it to the sum
        phiNumerator = phiNumerator + X(i)*(nC + (1/4)*nH + nS - (1/2)*nO);
    end
% Compute denominator
    phiDenominator = X(oxygenIndex);
% Compute phi
    phi = phiNumerator/phiDenominator; % if rich, phi>1, if lean, phi<1      
    
    
function [coefVec, nameEle, numEle] = readSandia(speciesVec, nasaFile)
% This function reads the thermodynamic data file

% See the accompanying document explaining the SANDIA.DAT file
% there are two fits for data: a high and a low temperature fit
% Each fit has 7 coefficients; we need to reassemble these. 
% Note: we are not reading in the date, and also there is a value in
% the 5th position in line 4 which is the enthalpy of formation divided
% by the gas constant R; this can be computed using the equations from
% the coefficients already read into the file.
% for more information, see the following resources:
%    http://garfield.chem.elte.hu/burcat/THERM.DAT
%    http://www.galcit.caltech.edu/EDL/public/formats/chemkin.html
%    http://www.galcit.caltech.edu/EDL/public/formats/nasaold.html
%    http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19780009781_1978009781.pdf

% make a checklist vector to ensure we got everything
    checkList = zeros(1,length(speciesVec));
% open the file
    fileID = fopen(nasaFile);
% read the first line
    firstLine = fgets(fileID);
% is this the right starter file?
    tempIndex = strfind(firstLine,'THERMO'); % index where 'THERMO' begins
    if isempty(tempIndex) % if 'THERMO' does not occur in first line
        error('%1s is empty',nasaFile)
    end
% read the next line
    secondLine = fgets(fileID);
% find the general temperature limits
    genTlow = str2num(secondLine(1:10)); % general low temp limit (K)
    genTcom = str2num(secondLine(11:20)); % general common temp (K)
    genThigh = str2num(secondLine(21:30)); % general high temp limit (K)
% read the file
    closeNow = false;
    while (closeNow ~= true)
        % reset variables
            elements = {''};
            numAtoms = [];
            Tlow = 0;
            Tcom = 0;
            Thigh = 0;
            TlowWarn = false;
            TcomWarn = false;
            ThighWarn = false;
        % get the first line of the species
            L1 = fgets(fileID);
        % check for ending
            endText = upper(strtrim(L1(1:3)));
            if strcmp(endText,'END')
                closeNow = true;
        % have we found everything?
            elseif sum(checkList) == length(checkList); % all 1's so this sums to the length if full! 
                closeNow = true; 
        % if not at the end, keep reading the file
            else
                % grab the species name
                    speciesName = upper(strtrim(L1(1:18))); % remove white space, make UPPERCASE
                % check the name; sometimes there is other text beside it
                    refIndex = strfind(speciesName,'REF'); % find where 'REF' occurs
                    if ~isempty(refIndex) % if 'REF' occured, then... 
                        speciesName = strtrim(speciesName(1:refIndex-1)); % take a substring
                    end
                % grab the elements and number of atoms of each
                    elpos = [25, 30, 35, 40, 74];
                    for i=1:5
                        % generate indicies
                            tmpIdx1 = elpos(i)+0; %element char 1
                            tmpIdx2 = elpos(i)+1; %element char 2
                            tmpIdx3 = elpos(i)+2; %number of elem digit 1
                            tmpIdx4 = elpos(i)+4; %number of elem digit 3
                        % grab the element name and the number of atoms
                            tmpEl = upper(strtrim(L1(tmpIdx1:tmpIdx2))); % UPPERCASE
                            tmpNum = str2num(L1(tmpIdx3:tmpIdx4));
                        % make sure the element is there 
                        % not all compounds have 5 elements
                        % if it's there, store it.  if not, throw it away. 
                        check = strcmp(tmpEl,''); % true if nothing, false if something
                        % some files have a zero in column 74
                        % this throws off the algorithm so check for it 
                        check2 = strcmp(tmpEl,'0'); 
                        if (check == false) && (check2 == false)
                            elements{1,i} = tmpEl;
                            numAtoms(1,i) = tmpNum;
                        end
                    end
                % grab the phase
                    phase = upper(L1(45)); % UPPERCASE letter
                % grab the low, middle, and high temperatures (K)
                    Tlow = str2num(L1(46:55));
                    Thigh = str2num(L1(56:65));
                    Tcom = str2num(L1(66:73));
                % check the inputs
                    if sum(size(Tlow))>2 || sum(size(Tcom))>2 || sum(size(Thigh))>2
                        error('SANDIA.DAT file format issue for %7s; Check for tabs instead of spaces.',speciesName)
                    end
                % if we did not see Tlow etc, use the general value
                    if isempty(Tlow); 
                        Tlow = genTlow;
                        TlowWarn = true; % flag that temperature was replaced
                    end
                    if isempty(Tcom);
                        Tcom = genTcom;
                        TcomWarn = true; % flag that temperature was replaced
                    end
                    if isempty(Thigh);
                        Thigh = genThigh;
                        ThighWarn = true; % flag that temperature was replaced
                    end
                % get line two
                    L2 = fgets(fileID);
                % grab the temperature coefficients here
                    Ha1 = str2num(L2(1:15)); % high temperature coefficient a1
                    Ha2 = str2num(L2(16:30));
                    Ha3 = str2num(L2(31:45));
                    Ha4 = str2num(L2(46:60));
                    Ha5 = str2num(L2(61:75));
                % get line three
                    L3 = fgets(fileID);
                % grab the temperature coefficients here
                    Ha6 = str2num(L3(1:15));
                    Ha7 = str2num(L3(16:30));
                    La1 = str2num(L3(31:45)); % low temperature coefficient a1
                    La2 = str2num(L3(46:60));
                    La3 = str2num(L3(61:75));            
                % get line four
                    L4 = fgets(fileID);
                % grab the temperature coefficients here
                    La4 = str2num(L4(1:15));
                    La5 = str2num(L4(16:30));
                    La6 = str2num(L4(31:45));
                    La7 = str2num(L4(46:60));              
                % if we want this species, store it, otherwise disregard
                    for i=1:length(speciesVec);
                        % if the same name as one of our wanted species, 
                        % and if it is a gas
                        % and if we have not grabbed it yet
                        if strcmp(speciesName,speciesVec{i}) && strcmp(phase,'G') && checkList(i)==0; 
                            % coefficeint and temperature matrix
                            coefVec(1,:,i) = [Ha1 Ha2 Ha3 Ha4 Ha5 Ha6 Ha7];
                            coefVec(2,:,i) = [La1 La2 La3 La4 La5 La6 La7];
                            coefVec(3,:,i) = [Tlow Tcom Thigh 0 0 0 0];
                            % element name and number vectors
                            for j=1:length(elements)
                                nameEle(i,j) = elements(1,j);  % use parentheses not curly brackets
                                numEle(i,j) = numAtoms(1,j);
                            end
                            % display warnings about changing temperatures
                            if TlowWarn
                                fprintf('WARNING: Defaulting Tlow to %4.2f K for %7s\n',genTlow,speciesName);
                            end
                            if TcomWarn
                                fprintf('WARNING: Defaulting Tcom to %4.2f K for %7s\n',genTcom,speciesName);
                            end
                            if ThighWarn
                                fprintf('WARNING: Defaulting Thigh to %4.2f K for %7s\n',genThigh,speciesName);
                            end
                            % check off that we got it
                            checkList(i) = 1; 
                        end
                    end
            end
    end
% close the file
    fclose(fileID);
% did we get everything?
    if sum(checkList) ~= length(checkList)
        for i=1:length(checkList)
            if checkList(i) == 0
                error('Species not in %5s: %5s',nasaFile,speciesVec{i})
            end
        end
    end
    
