Project 1 - Parsing NASA thermodynamic data
close all
clear all
clc
thermo = fopen('THERMO.dat','r')
fgetl(thermo);
tempdat = fgetl(thermo);
A = strsplit(tempdat,' ');
%Extracting the global temperature range
LowGt = str2num(A{2});
MidGt = str2num(A{3});
HighGt = str2num(A{4});
%Defining variable to store the input from end-user
get = input(' Enter Species Name','s');
while ~feof(thermo)
data = fgetl(thermo);
a = strsplit(data,' ');
b = a(1);
%Comparing input and data of molecule
match = strcmpi(b , get)
if match == 1
disp('Species Data Available')
%Extracting coefficients of desired moecule
Line1 = fgetl(thermo);
SL1 = findstr(Line1,'E');
Coef1 = Line1(1:SL1(1)+3);
Coef2 = Line1(SL1(1)+4:SL1(2)+3);
Coef3 = Line1(SL1(2)+4:SL1(3)+3);
Coef4 = Line1(SL1(3)+4:SL1(4)+3);
Coef5 = Line1(SL1(4)+4:SL1(5)+3);
a1 = str2num(Coef1);
a2 = str2num(Coef2);
a3 = str2num(Coef3);
a4 = str2num(Coef4);
a5 = str2num(Coef5);
Line2 = fgetl(thermo);
SL2 = findstr(Line2,'E');
Coef6 = Line2(1:SL2(1)+3);
Coef7 = Line2(SL2(1)+4:SL2(2)+3);
Coef8 = Line2(SL2(2)+4:SL2(3)+3);
Coef9 = Line2(SL2(3)+4:SL2(4)+3);
Coef10 = Line2(SL2(4)+4:SL2(5)+3);
a6 = str2num(Coef6);
a7 = str2num(Coef7);
a8 = str2num(Coef8);
a9 = str2num(Coef9);
a10 = str2num(Coef10);
Line3 =fgetl(thermo);
SL3 = findstr(Line3,'E');
Coef11 = Line3(1:SL3(1)+3);
Coef12= Line3(SL3(1)+4:SL3(2)+3);
Coef13 = Line3(SL3(2)+4:SL3(3)+3);
Coef14 = Line3(SL3(3)+4:SL3(4)+3);
a11 = str2num(Coef11);
a12 = str2num(Coef12);
a13 = str2num(Coef13);
a14 = str2num(Coef14);
%Extracting the local temperature range of the input molecule
mlt = a(end-1);
midlt = str2double(mlt);
hlt = a(end-2);
highlt = str2double(hlt);
llt = a(end-3);
lowlt = str2double(llt);
break
elseif match == 0
disp('Species Data Not Available')
else
continue
end
end
f = upper(get);
%Displaying Molecular Weight of the input molecule.
moleWeight = MW(get);
disp('Molecule Entered')
disp(f)
disp('Molecular Weight')
disp(moleWeight)
TRange = linspace(100,10000,9000);
R = 8.314;
%Obtaining value of SPECIFIC HEAT
cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,midlt,TRange,R);
%Obtaining value of ENTROPY
s = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,midlt,TRange,R);
%Obtaining value of ENTHALPY
h = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,midlt,TRange,R);
%Plotting Cp v T
figure(1)
plot(TRange,cp,'linewidth',3,'color','r')
grid on
xlabel('Temperature')
ylabel('Specific Heat (J/K)')
title(f)
%Plotting S v T
figure(2)
plot(TRange,s,'linewidth',3,'color','b')
grid on
xlabel('Temperature')
ylabel('Entropy (J/K)')
title(f)
%Plotting h v T
figure(3)
plot(TRange,h,'linewidth',3,'color','g')
grid on
xlabel('Temperature')
ylabel('Enthalpy (J/kg)')
title(f)
fclose('all');
%SPECIFIC HEAT FUNCTION
function cp_val = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,midlt,TRange,R)
if (TRange>midlt)
cp_val = ( a1+ (a2*TRange) + (a3 *(TRange.^2)) + (a4*(TRange.^3)) + (a5 * (TRange.^4)))*R;
else
cp_val = ( a8+ (a9*TRange) + (a10 *(TRange.^2)) + (a11*(TRange.^3)) + (a12 * (TRange.^4)))*R;
end
end
%ENTROPY FUNCTION
function s_val = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,midlt,TRange,R)
if (TRange>midlt)
s_val= ((a1*log(TRange)) + (a2 *TRange) + ((a3*(TRange.^2))/2) + ((a4*(TRange.^3)) /3) + ((a5*(TRange.^4))/4) + (a7))*R;
else
s_val= ((a8*log(TRange)) + (a9 *TRange) + ((a10*(TRange.^2))/2) + ((a11*(TRange.^3)) /3) + ((a12*(TRange.^4))/4) + (a14))*R;
end
end
%Molecular Weight Function
function Mol_Weight = MW(get)
Weight = 0;
% input as molecule name
% atoms dealt with H C N O Ar
% atomic mass 1 12 14 16 40 respectively
%to find molecular weight, we have to compare the atomsbfrom input
atoms = ['H','C','N','O','A'];
%Only A taken for Ar since A&R will be contained two positions in get
mass = [1,12,14,16,40];
%Since we have to compare each characters for mass, using 2 for loop
for i = 1:length(get)
for j = 1:length(atoms)
if strcmpi(get(i),atoms(j))
Weight = Weight + mass(j);
y = j;
end
end
%for more number of same atoms in the molecule
number = str2num(get(i));
if (number>1)
% for c2 we have to add mass of c 2 times
%initially since for C, mass is added to mol_weight
% for extra molecule we use (number -1) so adding one more weight
Weight = Weight + (mass(y)*(number - 1));
end
end
Mol_Weight = Weight;
end
%Enthalpy Function
function h_val = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,midlt,TRange,R)
if (TRange>midlt)
h_val = ( a1+ ((a2 .* TRange)./2) + ((a3.*(TRange.^2))./3) +((a4 .*(TRange.^3))./4) + ((a5 .* (TRange.^4))./5) + (a6./TRange)).* (R.*TRange);
else
h_val = ( a8+ ((a9 .* TRange)./2) + ((a10.*(TRange.^2))./3) +((a11 .*(TRange.^3))./4) + ((a12 .* (TRange.^4))./5) + (a13./TRange)).* (R.*TRange);
end
end