Project 1 - Parsing NASA thermodynamic data

NASA came up with polynomials that can be used to evaluate thermodynamic properties such as Cp, H and S. They have also documented the co-efficients that are required to evaluate these polynomials for 1000+ species.

To guide you we have broken down the problem into meaningful bits that you can work on individually.

1. Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file. You will be reading this file NASA thermodynamic data(Tip: use fopen and fgetl). The formulae are shown below.

Here R is the universal gas constant, T is the temeprature

function Enthalp

function [H] = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,tmid,R,T)

% finding enthalpy with low and temperature variables
if T <=tmid
H = (a8+(a9*T/2)+(a10*(T.^2)/3)+(a11*(T.^3)/4)+(a12*(T.^4)/5)+(a13./T))*R.*T;
else
H = (a1+(a2*T/2)+(a3*(T.^2)/3)+(a4*(T.^3)/4)+(a5*(T.^4)/5)+(a6./T))*R.*T;
end

end

function of Entropy

function [S] = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,tmid,T,R)

% finding entropy with low and temperature variables
if T <=tmid
S = ((a8*log(T))+(a9*T)+(a10*(T.^2)/2)+(a11*(T.^3)/3)+(a12*(T.^4)/4)+a14)*R;
else
S = ((a1*log(T))+(a2*T)+(a3*(T.^2)/2)+(a4*(T.^3)/3)+(a5*(T.^4)/4)+a7)*R;
end

function [W] = molecular_weight(species)

atomic_names = ['H', 'C', 'O', 'N', 'Ar'];
atomic_weights = [1.00784, 12.0107, 15.999, 14.0067, 39.948];
W=0;

for i=1:length(species)
for j = 1:length(atomic_names)
if species(i) == atomic_names(j)
W=W+atomic_weights(j);
position = j;
end
end
n=str2num(species(i));
if n >=1
W=W+atomic_weights(position)*(n-1);
end
end
fprintf('Molecular weight of %s = %f',species,W)
disp(' ')
end

end

function [cp] = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,tmid,T,R)

% finding specific heat with low and temperature variables
if T <=tmid
cp = (a8+(a9*T)+(a10*T.^2)+(a11*T.^3)+(a12*T.^4))*R;
else
cp = (a1+(a2*T)+(a3*T.^2)+(a4*T.^3)+(a5*T.^4))*R;
end

end

Mainfunction

% File Parsing - NASA thermodynamic data

clear all
close all
clc

% opening file
f1 = fopen('THERMO.dat','r');

% skipping first line
fgetl(f1);

% global temperatures
temps = fgetl(f1);
temps = str2num(temps);

global_low_temp = temps(1);
global_mid_temp = temps(2);
global_high_temp = temps(3);

% creating array for temperature
T = linspace(global_low_temp,global_high_temp,500);

% universal gas constant
R = 8.314;

% skipping next 3 lines
for i = 1:3
fgetl(f1);
end

%storing starts
for i = 1:53

% creating array for species
tline = fgetl(f1)
A = strsplit(tline,' ')
% species
species = A{1}

%local temperatures
tlow = str2num(A{end-3});
thigh = str2num(A{end-2});
tmid = str2num(A{end-1});

T = linspace(tlow,thigh,500);

% finding coefficients
%line 1
tline1 = fgetl(f1);
pos1 = findstr(tline1,'E');

% higher coefficients
a1 = str2num(tline1(1:pos1(1)+3));
a2 = str2num(tline1(pos1(1)+4:pos1(2)+3));
a3 = str2num(tline1(pos1(2)+4:pos1(3)+3));
a4 = str2num(tline1(pos1(3)+4:pos1(4)+3));
a5 = str2num(tline1(pos1(4)+4:pos1(5)+3));

% line 2
tline2 = fgetl(f1);
pos2 = findstr(tline2,'E');
a6 = str2num(tline2(1:pos2(1)+3));
a7 = str2num(tline2(pos2(1)+4:pos2(2)+3));

% lower coefficients
a8 = str2num(tline2(pos2(2)+4:pos2(3)+3));
a9 = str2num(tline2(pos2(3)+4:pos2(4)+3));
a10 = str2num(tline2(pos2(4)+4:pos2(5)+3));

% line 3
tline3 = fgetl(f1);
pos3 = findstr(tline3,'E');
a11 = str2num(tline3(1:pos3(1)+3));
a12 = str2num(tline3(pos3(1)+4:pos3(2)+3));
a13 = str2num(tline3(pos3(2)+4:pos3(3)+3));
a14 = str2num(tline3(pos3(3)+4:pos3(4)+3));

% specific heat
cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,tmid,T,R);

% enthalpy
H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,tmid,R,T );

% entropy
S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,tmid,T,R );

% changing directory

% creating newfolder for species
mkdir(species)

% changing directory
cd(species)

% plotting

figure(1)
plot(T,cp)
xlabel('Temperature')
ylabel('Specific heat')
title('Specific heat vs Temperature of ',sprintf('%s',species))
saveas(figure(1),'cp vs T.jpg')

figure(2)
plot(T,H)
xlabel('Temperature')
ylabel('Enthalpy')
title('Enthalpy vs Temperature of ',sprintf('%s',species))
saveas(figure(2),'H vs T.jpg')

figure(3)
plot(T,S)
xlabel('Temperature')
ylabel('Entropy')
title('Entropy vs Temperature of ',sprintf('%s',species))
saveas(figure(3),'S vs T.jpg')

% coming back to program directory
cd ..
% calculating molecular weight of species
M = molecular_weight(species);

end