Project - Parsing NASA thermodynamic data

Q1. Write a function that extracts the 14 co-efficients for a given species.

%PARSING NASA\'S THERMODYNAMIC DATA

clear all
close all
clc

R = 8.314; %boltzmann constant


f1 = fopen(\'thermo30.dat\',\'r\'); %assign a handle to file
title = fgetl(f1); %read the first line \'THERMO\'    

%read the globall thermal range string
gtemprangearray = fgetl(f1); 

gtemprangesplit = strsplit(gtemprangearray,\' \'); %Splits string based on character

gtemplowr = str2num(gtemprangesplit{2}); %use of curly brace for CELL 
gtempmidr = str2num(gtemprangesplit{3});
gtempupr = str2num(gtemprangesplit{4});

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

for i=1:53

%indentifying specie 
sp = fgetl(f1);
spe = strsplit(sp,\' \');
element = spe{1};

%defining local temp range
for k = 1:length(spe) 
    if strcmp(spe(k),\'G\') == 1 
       llowt = str2num(spe{k+1});
       lmidt = str2num(spe{k+3});
       luprt = str2num(spe{k+2});
       llt = llowt:lmidt;
       lut = (lmidt+1):luprt;
    end
end

%reading coefficients
%findstr() can also be used to split on the basis of exponential power
coe = fscanf(f1,\'%E\'); %exponential data format specifier
coe = coe\';

%assigning values to coefficients 
a1u(i) = coe(1);
a2u(i) = coe(2);
a3u(i) = coe(3);
a4u(i) = coe(4);
a5u(i) = coe(5);
a6u(i) = coe(7);
a7u(i) = coe(8);
a1l(i) = coe(9);
a2l(i) = coe(10);
a3l(i) = coe(11);
a4l(i) = coe(13);
a5l(i) = coe(14);
a6l(i) = coe(15);
a7l(i) = coe(16);

cpl = cplcalc(element,a1l(i),a2l(i),a3l(i),a4l(i),a5u(i),llt,R);
cpu = cpucalc(element,a1u(i),a2u(i),a3u(i),a4u(i),a5u(i),lut,R);

h1 = hlcalc(element,a1l(i),a2l(i),a3l(i),a4l(i),a5l(i),a6l(i),llt,R);
hu = hucalc(element,a1u(i),a2u(i),a3u(i),a4u(i),a5u(i),a6u(i),lut,R);

sl = slcalc(element,a1l(i),a2l(i),a3l(i),a4l(i),a5l(i),a7l(i),llt,R);
su = sucalc(element,a1u(i),a2u(i),a3u(i),a4u(i),a5u(i),a7u(i),lut,R);

%disp(\'Molecular Weight of \')
%element
%disp(\' = \');
%molwt = atomicwt(element)


end

This program, firstly, opens the file \'thermo30.dat\' (line 10) using fopen.

Then it reads the title (line 11) , followed by reading the global temperature range (line 14), both using fgetl.

gtemprangesplit = strsplit(gtemprangearray,\' \') splits the global temperature range on the basis of spaces between the temperatures.

Then the global temperature limits are assigned to variables called gtemplowr, gtempmidr and gtempupr, representing the lower, middle and upper global temperature limit respectively.

From line 22 to 24, the 3 comments in the NASA file are skipped using for loop and fgetl.

From line 26, a for loop is entered. Till the end of this for loop, the code determines the species, calculates CP, entropy, enthalpies and molecular weights of the species encountered in the rest of thermo30.dat\' file.

From line 29 to 31, the first line (of 1,2,3,4) is always read, and broken to determine the species. The entire string is stored in a string called sp, which is then split depending upon the spaces and further stored into a cell called spe. The first element of spe i.e spe{1} is the species chemical formula, hence it is stored in a variable called element.

From line 33 to 42, the local temperature range is detected. Since in the entire file, the temperature range always begins after G , a for loop is used to detect where \'G\' lies in element. Once the index (k) is found, k+1 is the index of lower temperature limit, k+3 is the index of middle temperature limit, and k+2 is the index of upper temperature limit.

The lower temperature limit, middle temperature limit and upper temperature limit are assigned to variables called llowt, lmidt, luprt. using which lower range and upper range of temperature is generated using llowt:lmidt and lmidt+1:luprt, called llt and lut respectively.

coe = fscanf(f1,\'%E\');
coe = coe\';

The variable coe reads the coefficients in line 2,3 and 4.

%E is the exponential data format specifier, and works equally well, without requiring the need to split string on the basis of index of E.

Q2. Write functions that calculate the enthalpy, entropy and specific heats using the NASA polynomials.

From line 50 to 63, all coefficients (both upper and lower) are assigned for each species.

cpl = cplcalc(element,a1l,a2l,a3l,a4l,a5u,llt,R);
cpu = cpucalc(element,a1u,a2u,a3u,a4u,a5u,lut,R);

From line 65 to 66 , CP is calculated for lower temp range (cpl) and upper temperature range (cpu). These functions take the element name string, the coefficients, lower temperature range, Boltzmann constant and arguements, calculate the value of CP, and plot graph for each species.

Given is the syntax for these 2 functions :

function [cpl] = cplcalc(ele,a1l,a2l,a3l,a4l,a5l,llt,R)
 %calculation of cp 
%lower range

cpl = (a1l + a2l.*llt + a3l*(llt.^2) + a4l*(llt.^3) + a5l*(llt.^4))*R;


figure(1)
plot(llt,cpl);
xlabel(\'Temperature(lower range)\');
ylabel(\'CP\');
saveas(figure(1),ele,\'jpg\')


end

function [cpu] = cpucalc(element,a1u,a2u,a3u,a4u,a5u,lut,R)
 %calculation of cp 
%upper range

cpu = (a1u + a2u.*lut + a3u*(lut.^2) + a4u*(lut.^3) + a5u*(lut.^4))*R;

figure(1)
plot(lut,cpu);
xlabel(\'Temperature(upper range)\');
ylabel(\'CP\');
saveas(figure(1),element,\'png\')

end

In the line 5 (OF BOTH FUNCTIONS) the value of cp is calculated. From line 8 to 12 the graph is plotted, for which the name is taken as arguement. The graph is saved using saveas command.

Now coming back to the main program,

h1 = hlcalc(element,a1l,a2l,a3l,a4l,a5l,a6l,llt,R);
hu = hucalc(element,a1u,a2u,a3u,a4u,a5u,a6u,lut,R);

From line 68 to 69 , enthalpy is calculated for lower temp range (hl) and upper temperature range (hu). These functions take the element name string, the coefficients, lower temperature range, Boltzmann constant and arguements, calculate the value of enthalpy and plot graph for each species.

Given is the syntax for these 2 functions :

function [hl] = hlcalc(ele,a1l,a2l,a3l,a4l,a5l,a6l,llt,R)
 %calculation of enthalpy
%lower range

hl = (a1l + (a2l.*llt)/2 + a3l*(llt.^2)/3 + a4l*(llt.^3)/4 + a5l*(llt.^4)/5 + a6l./llt)*R.*llt;


figure(1)
plot(llt,hl);
xlabel(\'Temperature(lower range)\');
ylabel(\'Enthalpy\');
saveas(figure(1),ele,\'jpg\')


end



x----------------------------------------------------------------x


function [hu] = hucalc(ele,a1u,a2u,a3u,a4u,a5u,a6u,lut,R)
 %calculation of enthalpy
%upper range

hu = (a1u + (a2u.*lut)/2 + a3u*(lut.^2)/3 + a4u*(lut.^3)/4 + a5u*(lut.^4)/5 + a6u./lut)*R.*lut;


figure(1)
plot(lut,hu);
xlabel(\'Temperature(upper range)\');
ylabel(\'Enthalpy\');
saveas(figure(1),ele,\'png\')


end

In the line 5 (OF BOTH FUNCTIONS) the value of enthalpy is calculated. From line 8 to 12 the graph is plotted, for which the name is taken as arguement. The graph is saved using saveas command.

Now coming back to the main program,

sl = slcalc(element,a1l,a2l,a3l,a4l,a5l,a7l,llt,R);
su = sucalc(element,a1u(i),a2u(i),a3u(i),a4u(i),a5u(i),a7u(i),lut,R);

From line 71 to 72 , entropy is calculated for lower temp range (sl) and upper temperature range (su). These functions take the element name string, the coefficients, lower temperature range, Boltzmann constant and arguements, calculate the value of entropy and plot graph for each species.

Given is the syntax for these 2 functions :

function [sl] = slcalc(ele,a1l,a2l,a3l,a4l,a5l,a7l,llt,R)
 %calculation of entropy
%lower range

sl = (a1l*log(llt) + a2l.*llt + a3l*(llt.^2)/2 + a4l*(llt.^3)/3 + a5l*(llt.^4)/4 + a7l)*R;


figure(1)
plot(llt,sl);
xlabel(\'Temperature(lower range)\');
ylabel(\'Entropy\');
saveas(figure(1),ele,\'jpg\')


end


x-----------------------------------------------------------x


function [su] = sucalc(ele,a1u,a2u,a3u,a4u,a5u,a7u,lut,R)
 %calculation of entropy
%upper range

su = (a1u*log(lut) + a2u.*lut + a3u*(lut.^2)/2 + a4u*(lut.^3)/3 + a5u*(lut.^4)/4 + a7u)*R;


figure(1)
plot(lut,su);
xlabel(\'Temperature(upper range)\');
ylabel(\'Entropy\');
saveas(figure(1),ele,\'png\')


end

In the line 5 (OF BOTH FUNCTIONS) the value of entropy is calculated. From line 8 to 12 the graph is plotted, for which the name is taken as arguement. The graph is saved using saveas command.

Q3. Write a function that calculates the molecular weight of each species.

Returning back to the main program,

disp(\'Molecular Weight of \')
element
disp(\' = \');
molwt = atomicwt(element)

In line 74 molecular weight of species is calculated using the function atomicwt, syntax for which is as follows :

function mass = atomicwt(eqn)
 
 wtmap=zeros(83,1); %make a row vector of 83 zeros
 wtmap( \'CHONS\'/1 )= [12 1 16 14 32]; %assign molwts of CHONS to ASCII eq positions of CHONS
 
 mass=0;
 
 ELE=regexp(eqn,\'[CHONS]\'); % Locate all elements
   for i=1:length(ELE)
       wt=wtmap(eqn(ELE(i))); % ie for ch3cooh ele(1) = 1, eqn(1) = C, wtmap(C) = 12
       
       if i==length(ELE) %last character
          if ELE(end)==length(eqn) % single atom 
             mass=mass+wt;
          else % Case of final numeric eg. CH4
              qty=str2num(eqn(ELE(i)+1:end));
              mass=mass+wt*qty;
          end
       elseif ELE(i)+1==ELE(i+1) % single atom intra eqn
              mass=mass+wt;
       else % Element followed by a qty  allow >9
           qty=str2num(eqn(ELE(i)+1:ELE(i+1)-1));
           mass=mass+wt*qty;
       end
   end
end

eqn is a string passed to this function, and it evaluates and returns its weight back to molwt where it is displayed.

Q4. Write a function that plots the Cp, Enthalpy and Entropy of an input species

The function which calculates the value of Cp, enthalpy and entropy also plots their graphs, given in Q2. The plots are saved with the species name as the image name. The names are generated automatically as the species name after being stored in element is also passed as an arguement to these functions. This string is also used to assign these plots a name.

Q5. Display useful plots (enthalpy, entropy and cp). You should make these plots for o2, N2 and Co2.

Plots of O2 :

CP :

Enthalpy :

Entropy :

For N2 :

CP :

Enthalpy :

Entropy :

For C02 :

CP :

Enthalpy :

Entropy :