function [Cp,H,S]=Shomate(T,molecule_data) % Program to compute thermo stored as Shomate coefficients % W.H. Green, 9/27/02 % Outputs: % 1) constant-pressure heat capacity Cp in J/mol K % 2) enthalpy in kJ/mol % 3) entropy in J/mol K % Note, if T is a vector, these will all be vectors of the same shape. % Inputs: T in K, molecule_data = [A B C D E F G Tmin Tmax]' % A-G, Tmin, Tmax from NIST Webbook http://webbook.nist.gov % in Joule/mole Kelvin units % Cp0= A + B*t + C*t^2 + D*t^3 + E/t^2 % H0= A*t + B*t^2/2 + C*t^3/3 + D*t^4/4 - E/t + F % S0= A*ln(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G % Tmin and Tmax define the range over which % these A-G are accurate % Sometimes there are different parameter sets for different T ranges % if so, make molecule_data a matrix, with a column for each T range. % First: figure out how many parameter sets we have % and how many T's the user has input dim_molecule_data=size(molecule_data); number_of_parameter_sets=dim_molecule_data(2); dim_T=size(T); number_of_T=dim_T(1); Cp=-1e20*ones(dim_T); % if you see any -1e20's in output, something went wrong H=Cp; % most likely your T was outside the T range where the S=Cp; % Shomate fit is valid. for i=1:number_of_parameter_sets Tmin=molecule_data(8,i); Tmax=molecule_data(9,i); for j=1:number_of_T if((T(j)>=Tmin)&(T(j)<=Tmax)) A = molecule_data(1,i); B=molecule_data(2,i); C = molecule_data(3,i); D=molecule_data(4,i); E = molecule_data(5,i); F=molecule_data(6,i); G = molecule_data(7,i); t=0.001*T(j); % NIST t= T/1000 Kelvin Cp(j,1)=A+(B*t)+(C*t*t)+(D*(t^3))+(E/(t*t)); H(j,1)=(A*t)+(0.5*B*t*t)+(C*(t^3)/3)+(0.25*D*(t^4))-(E/t)+F; S(j,1)=(A*log(t))+(B*t)+(0.5*C*t*t)+(D*(t^3)/3)-(0.5*E/(t*t))+G; end end end % if any of the values are -1e20, it means that the Shomate polynomial % was not valid in the requested T range. % % Sample input % For H2 % H2_data=[33.10780 -11.50800 11.60930 -2.844400 -0.159665 -9.991971 172.7880 298 1500]'; % By running [Cp,H,S]=Shomate(300, H2_data); % % Sample output % » Cp % Cp = % 28.8494 % » H % H = % 0.0534 % » S % S = % 130.8586