%Function to calculate Entropy
function S=Entropy_cal(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,temp,R,global_med_temp)
%Comparre with global temperature
%If temperature is less than global temperature then minimum coefficients are taken else maximum coefficients are taken
if temp<=global_med_temp
S=R*(b1.*log(temp)+(b2.*temp)+(b3*((temp).^2)/2)+(b4*((temp).^3)/3)+(b5*((temp).^4)/4)+b7);
else
S=R*(a1.*log(temp)+(a2.*temp)+(a3*((temp).^2)/2)+(a4*((temp).^3)/3)+(a5*((temp).^4)/4)+a7);
end
%Function to calculate specific heat
function Cp=specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,temp,R,global_med_temp)
%Compare woth global temperature
%If temperature is less than global temeperature then minimum coefficients are taken else maximum coefficients are taken
if temp<=global_med_temp
Cp=R*(b1+(b2*temp)+(b3*(temp).^2)+(b4*(temp).^3)+(b5*(temp).^4));
else
Cp=R*(a1+(a2*temp)+(a3*(temp).^2)+(a4*(temp).^3)+(a5*(temp).^4));
end
%Function to calculate Enthalpy
function H=Enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,temp,R,global_med_temp)
%Compare with global temperature
%If temperature is less than global temperature then minimum coefficients are taken else maximum coefficients are taken
if temp<=global_med_temp
H=R*temp.*(b1+((b2*temp)/2)+((b3*(temp).^2)/3)+((b4*(temp).^3)/4)+((b5*(temp).^4)/5)+((b6./temp)));
else
H=R*temp.*(a1+((a2*temp)/2)+((a3*(temp).^2)/3)+((a4*(temp).^3)/4)+((a5*(temp).^4)/5)+((a6./temp)));
end
%Function to calculate Molecular weight
function Molecularweight=Molecular_weight(Species_SP)
atomic_name=['H' , 'c' , 'O' , 'N' , 'A' ];
%Hydrogen,Carbon,Oxygen,Nitrogen,Argon
fprintf('Molecular weight of');
fprintf(Species_SP);
fprintf( ' : ' );
%Standard Atomic Number
Atomic_weight=[1.008 12.011 15.999 14.007 39.948];
Molecularweight=0;
%Add atomic weight to obtain particular value
for i=1:length(Species_SP)
for j=1:length(atomic_name)
if strcmpi(Species_SP(i),atomic_name(j))
Molecularweight=Molecularweight+Atomic_weight(j);
position=j;
end
end
n=str2num(Species_SP(i));
if n>1
Molecularweight=Molecularweight+Atomic_weight(position)*(n-1);
end
end
fprintf('Molecular weight %s :',Species_SP);
fprintf('%d',Molecularweight);
disp(' ');
end
clear all
close all
clc
%Make directory
%Reading data file-NASA Thermodynamic
f1=fopen('THERMO.dat','r')
gt=fgetl(f1);
%Obtaining temperatures by reader header information
global_temp=fgetl(f1);
%Splitting
t=strsplit(global_temp);
%Converting string cells to numbers using str2num
global_low_temp=str2num(t{2});
global_med_temp=str2num(t{3});
global_high_temp=str2num(t{4});
%Temperature range & Gas Constant
%Temperature Range assumed
temp=linspace(global_low_temp,global_high_temp,1000);
%Universal gas constant J/(mol-K)
R=8.314;
%Skipping three comment lines in the file
for i=1:3
gt=fgetl(f1);
end
for i=1:53
gt=fgetl(f1);
A=strsplit(gt,' ');
Species_SP=(A{1})
line_1=fgetl(f1);
a=strfind(line_1,'E');
%Extracting higher and lower temperature coefficients
%Extracting first 7 coefficients-High temperature coefficients
%First coefficient
a1=line_1(1:a(1)+3);
a1=str2num(a1);
%second coefficient
a2=line_1(a(1)+4:a(2)+3);
a2=str2num(a2);
%Third coefficient
a3=line_1(a(2)+4:a(3)+3);
a3=str2num(a3);
%Fourth coefficient
a4=line_1(a(3)+4:a(4)+3);
a4=str2num(a4);
%Fifth coefficient
a5=line_1(a(4)+4:a(5)+3);
a5=str2num(a5);
line_2=fgetl(f1);
b=strfind(line_1,'E');
%Sixth coefficient
a6=line_2(1:b(1)+3);
a6=str2num(a6);
%Seventh coefficient
a7=line_2(a(1)+4:b(2)+3);
a7=str2num(a7);
%Extracting second 7 coefficients-Low temperature coefficients
%Eighth coefficient
b1=line_2(b(2)+4:b(3)+3);
b1=str2num(b1);
%Ninth coefficient
b2=line_2(b(3)+4:b(4)+3);
b2=str2num(b2);
%Tenth coefficient
b3=line_2(b(4)+4:b(5)+3);
b3=str2num(b3);
line_3=fgetl(f1);
c=strfind( line_1,'E');
%Eleventh coefficient
b4=line_3(1:c(1)+3);
b4=str2num(b4);
%Twelth coefficient
b5=line_3(c(1)+4:c(2)+3);
b5=str2num(b5);
%Thirteenth coefficient
b6=line_3(c(2)+4:c(3)+3);
b6=str2num(b6);
%Fourteenth coefficient
b7=line_3(c(3)+4:c(4)+3);
b7=str2num(b7);
%Specific heat calculation
Cp=specificheat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,temp,R,global_med_temp);
%Enthalpy calculation
H=Enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,temp,R,global_med_temp);
%Entropy calculation
S=Entropy_cal(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,temp,R,global_med_temp);
%Calculating molecular weight of species
Molecularweight(i)=Molecular_weight(Species_SP);
%Read & Write Molecular list
Molecular_list=fopen('Molecular weight.txt','a+');
fprintf(Molecular_list,'Molecular weight of %s is %d n',Species_SP,Molecularweight);
fclose(Molecular_list);
%plotting(CP,H,S,Temp,Species_SP)
%Saving as a folder
cur_dir=pwd;
mkdir(Species_SP);
cd(Species_SP);
%Plot 1,Specific heat(CP) vs Temperature Range (Temp)
figure(1)
plot(temp,Cp,'linewidth',2,'color','r');
xlabel('Temperature [k]');
ylabel('Specific Heat [KJ]');
title(sprintf('Varying specific heat with temperature range of %s',Species_SP));
grid on
saveas(figure(1),'Specific Heat.jpg');
%Plot 2,Enthalpy (H) vs Tmperature Range(Temp)
figure(2)
plot(temp,H,'linewidth',2,'color','g');
xlabel('Temperature[K]');
ylabel('Enthalpy in [KJ/(mol)]');
title(sprintf('Varying Enthalpy with Temperature Range of %S',Species_SP));
grid on
saveas(figure(2),'Enthalpy.jpg');
%Plot 3,Entropy (S) vs Temperature Range (Temp)
figure(3)
plot(temp,S,'linewidth',2,'color','b');
xlabel('Temperature[k]');
ylabel('Entropy in [KJ/(mol-K)]');
title(sprintf('varying Entropy with Temperature Range of %S',Species_SP));
grid on
saveas(figure(3),'Entropy.jpg');
cd(cur_dir);
%Molecular weight'
M=Molecular_weight(Species_SP)
f2=fopen('molecular mass','w');
fprintf('%s n %s n %f n''species','Molecular mass',M)
fclose(f2)
end