clear
clearglobal
funcprot(0)

load(uigetfile('*.sod')); // filename -> therm.sod, contains mlist M 

function r=%hm_e(varargin)
  M=varargin($)
  r=mlist(['hm','Symbol','Temp','Coefficients'],M.Symbol(varargin(1:$-1)),M.Temp(varargin(1:$-1),:),M.Coefficients(varargin(1:$-1),:));
endfunction

function [y]=Therm(species,T)
// species -> species name, ex: 'CH4'
// T -> temperature scalar or vector, K

ind=find(M.Symbol==species);
temp = M.Temp(ind,:);

a = M.Coefficients(ind,8:14); //Low temp region
A = M.Coefficients(ind,1:7); //High temp region

//R = 1.9872041 cal * mol^-1 * K^-1
//R = 8.3144621    J * mol^-1 * K^-1 

//Cp/R -> Specific Heat
y(:,1) = ([ones(T);T;T.^2;T.^3;T.^4]'*a(1:5)').*(T'<=temp(3))+..  //T'>=temp(1) &
            ([ones(T);T;T.^2;T.^3;T.^4]'*A(1:5)').*(T'>temp(3) & T'<=temp(2));

// H/R/T -> Enthalpy
y(:,2) = ([ones(T);T;T.^2;T.^3;T.^4;T.^-1]'*(a(1:6)'.*[1;1/2;1/3;1/4;1/5;1])).*(T'<=temp(3))+.. //T'>=temp(1) & 
            ([ones(T);T;T.^2;T.^3;T.^4;T.^-1]'*(A(1:6)'.*[1;1/2;1/3;1/4;1/5;1])).*(T'>temp(3) & T'<=temp(2));

// S/R -> Entropy
y(:,3) = ([log(T);T;T.^2;T.^3;T.^4;ones(T)]'*([a(1:5),a(7)]'.*[1;1;1/2;1/3;1/4;1])).*(T'<=temp(3))+.. //T'>=temp(1) & 
            ([log(T);T;T.^2;T.^3;T.^4;ones(T)]'*([A(1:5),A(7)]'.*[1;1;1/2;1/3;1/4;1])).*(T'>temp(3) & T'<=temp(2));

endfunction

test=Therm('AL',[300:50:5000]);

scf();
subplot(1,3,1);plot2d([300:50:5000]',test(:,1),style=2);xgrid() // Cp/R
xtitle('','Temperature, K','Specific Heat, Cp/R')
subplot(1,3,2);plot2d([300:50:5000]',test(:,2),style=14);xgrid() // H/R/T
xtitle('','Temperature, K','Enthalpy,H/R/T')
subplot(1,3,3);plot2d([300:50:5000]',test(:,3),style=24);xgrid() // S/R
xtitle('','Temperature, K','Entropy,S/R')
