Project 1 - Parsing NASA thermodynamic data
Objective
The objective of this challenge is to parse the data from a file 'THERMO.dat' and calculate specific heat, enthalpy, entropy and molecular weight.
For this project I've written two user-defined functions. One of the function 'thermo_properties' calculates the specific heat, enthalpy and entropy of the species. The other function 'molecular_weight' calculates molecular weight of the species. These functions will be called inside the main script. After that plotting is done inside the main script using the returned values of specific heat, enthalpy and entropy. The molecular weight will be printed as output in the command window.
Main Script
clear
close all
clc
thermo_data = fopen('THERMO.dat','r'); % opening file for reading and assigning the file to variable thermo_data
file_name = fgetl(thermo_data); % reading line from the file and assigning it to the variable
temp_line = fgetl(thermo_data); % reading line from the file and assigning it to the variable
temp = strsplit(temp_line,' '); % splitting the string vector 'temp_line' into cells. 'Space' is the delimiter
global_low_temp = str2double(temp{2}); % converting the string to double values and assigning to the variables
global_medium_temp = str2double(temp{3}); % converting the string to double values and assigning to the variables
global_high_temp = str2double(temp{4}); % converting the string to double values and assigning to the variables
for i = 1:3 % initializing the for loop variable i to read next three lines
eval(sprintf('Comment%d = fgetl(thermo_data);',i)) % using eval and sprintf to read the comments from the file and assigning the comments to variables
end % ending the for loop
for i = 1:53 % initializing the for loop variable i to read next 53 lines from the file.
species_header = (strsplit(fgetl(thermo_data),' ')); % reading next line and splitting the string and then assigning cells to the variable
species = species_header{1}; % assgining the string from the first cell to the variable
% converting the strings from the header line to double and assigning to the variable
local_low_temp = str2double(species_header{end-3});
local_mid_temp = str2double(species_header{end-1});
local_high_temp = str2double(species_header{end-2});
for j = 1:3 % initializing for loop variable j to read next three line
eval(sprintf('line%d = fgetl(thermo_data);',j)) % using eval to execute string created using sprintf as code
eval(sprintf('p%d = strfind(line%d,''E'');',j,j)) % using strfind function to get the location of character 'E' and assigning the indices to a variable
end % end the loop
% converting the coefficients from first three lines of the species to double and assigning to the variables.
% a1 to a7 are high temperature coefficients and a8 to a14 are low temperature coefficients
a1 = str2double(line1(1:p1(1)+3));
a2 = str2double(line1(p1(1)+4:p1(2)+3));
a3 = str2double(line1(p1(2)+4:p1(3)+3));
a4 = str2double(line1(p1(3)+4:p1(4)+3));
a5 = str2double(line1(p1(4)+4:p1(5)+3));
a6 = str2double(line2(1:p2(1)+3));
a7 = str2double(line2(p2(1)+4:p2(2)+3));
a8 = str2double(line2(p2(2)+4:p2(3)+3));
a9 = str2double(line2(p2(3)+4:p2(4)+3));
a10 = str2double(line2(p2(4)+4:p2(5)+3));
a11 = str2double(line3(1:p3(1)+3));
a12 = str2double(line3(p3(1)+4:p3(2)+3));
a13 = str2double(line3(p3(2)+4:p3(3)+3));
a14 = str2double(line3(p3(3)+4:p3(4)+3));
T = linspace(local_low_temp,local_high_temp,500); % using linspace to create temperature array of 500 elements between local low temp and local high temp
[cp,h,s] = thermo_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,local_mid_temp); % calling thermo_properties to calculate values specific heat, enthalpy and entropy
mol_wt = molecular_weight(species); % calling molecular_weight function to calculate molecular weight
current_directory = pwd; % assigning path to current directory to variable using pwd
cd parsing_challenge % changing current working directory to folder parsing_challenge
mkdir(species) % creating new folder for species
cd(species) % changing current working directory
% plotting specific heat vs temperature graph
figure(1) % first figure for the specific heat vs temperature graph
plot(T,cp,'linewidth',3) % plotting with default color and linewidth 3
xlabel('Temperature')
ylabel('Specific Heat')
title(sprintf('Temperature vs Specific Heat: %s',species)) % using sprintf to create title for the first figure
% plotting temperature vs enthalpy graph
figure(2) % second figure for the temperature vs enthalpy graph
plot(T,h,'color','k','linewidth',3) % plotting with black color and linewidth 3
xlabel('Temperature')
ylabel('Enthalpy')
title(sprintf('Temperature vs Enthalpy: %s',species)) % using sprintf to create title for the 2nd figure
% plotting temperature vs entropy graph
figure(3) % third figure for temperature vs entropy graph
plot(T,s,'color','c','linewidth',3) % plotting with cyan color and linewidth 3
xlabel('Temperature')
ylabel('Entropy')
title(sprintf('Temperature vs Entropy: %s',species)) % using sprintf to create title for the 3rd figure
saveas(1,'Specific_Heat.png') % saving the figure 1 as png file
saveas(2,'Enthalpy.png') % saving the figure 2 as png file
saveas(3,'Entropy.png') % saving the figure 3 as png file
cd(current_directory) % changing the current working directory initial directory
end % ending the loop
fclose(thermo_data); % closing the file
User-defined Function 'thermo_properties'
function [cp, h, s] = thermo_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,local_mid_temp)
R = 8.314; % specifying the universal gas constant
if T > local_mid_temp % using if statement to specify conditions for high temp coefficient
cp = R.*(a1 + a2.*T + a3.*T.^2 + a4.*T.^3 + a5.*T.^4); % calculating specific heat
h = R.*T.*(a1 + a2.*T./2 + a3.*T.^2./3 + a4.*(T.^3)./4 + a5.*(T.^4)./5 + a6./T); % calculating enthalpy
s = R.*(a1.*log(T) + a2.*T + a3.*T.^2./2 + a4.*T.^3./3 + a5.*T.^4./4 + a7); % calculating entropy
else % low temperature coefficients
cp = R.*(a8 + a9.*T + a10.*T.^2 + a11.*T.^3 + a12.*T.^4); % calculating specific heat
h = R.*T.*(a8 + a9.*T./2 + a10.*T.^2./3 + a11.*T.^3./4 + a12.*T.^4./5 + a13./T); % calculating enthalpy
s = R.*(a8.*log(T) + a9.*T + a10.*T.^2./2 + a11.*T.^3./3 + a12.*T.^4./4 + a14); % calculating entropy
end
endUser-defined Funciton 'molecular_weight'
function [mol_wt] = molecular_weight(species)
elements = {'O','H','N','C','AR'}; % assigning the symbol of the elements to a cell array
atomic_wt = [15.999,1.00794,14.0067,12.0107,39.948]; % assigning atomic weights of the corresponding elements
mol_wt = 0; % initializing variable mol_wt
for i = 1:length(species) % initializing for loop variable i
atoms = str2double(species(i)); % converting the character from the string to double
for j = 1:length(elements) % initializing loop variable j
if strcmpi(species(i),elements{j}(1)) % using conditional if statement and using strcompi to compare the string characters
mol_wt = (mol_wt + atomic_wt(j)); % calculating molecular weight
element_location = j; % assigning the location of matched element to variable
elseif ~isnan(atoms) % using the elseif statement and ~isnan logical funciton to set the condition
mol_wt = mol_wt + atomic_wt((element_location))*(atoms-1); % calculating molecular weight
break % using break statement to exit the loop
end % ending the if-elseif statement
end % ending the for loop
end % ending the for loop
fprintf('Molecular weight of %s is: %.f amu \n',species,round(mol_wt)) % printing the molecular weight.
end Explanation for the main script
First I've opened the file for reading using fopen command and assigned it to variable thermo_data. Then I've used fgetl function to read the first and second line from the file and these lines are assigned to variables. Then I've used strspilt function on the string array 'temp_line' to split the string at delimeter 'space'. This will give a cell array. Now to convert the temperature values from character type to number I've used function str2double. The location of the cells in the cell array is provided using braces{}. Then the converted values are assigned to the variables. After this I've used for loop to read the next three lines which are comments. I've used fgetl, eval and sprintf functions inside the loop to read the lines and assign them to the variables. sprintf function writes the formatted data to string and eval function will execute the string or character as regular code.
thermo_data = fopen('THERMO.dat','r');
file_name = fgetl(thermo_data);
temp_line = fgetl(thermo_data);
temp = strsplit(temp_line,' ');
global_low_temp = str2double(temp{2});
global_medium_temp = str2double(temp{3});
global_high_temp = str2double(temp{4});
for i = 1:3
eval(sprintf('Comment%d = fgetl(thermo_data);',i))
end
After this to read the next lines I've again used the for loop. There are 53 species and each species has four lines of data. First line contains the species name and local temperature ranges and the next lines have the temperature coefficients ( total 14 in number). The first seven coefficients are high temp coefficients and the next seven coefficients are low temp coefficients. I've initialized the for loop variable i from 1 to 53. Then I've used fgetl and strsplit together to read the line and to split it and assigned it to a variable. Then the location of the name of the species is specified and assigned to the variable. After this the local temperature ranges are located and converted from string to double using str2double and are assigned to the variables. Now I've used another for loop inside the for loop to read the next three lines containing temperature coefficients. For this I've initialized the loop variable i from 1 to 3 and then again I've used fgetl, eval and sprintf functions inside the loop to read the lines and assign them to the variables 'line1', 'line2' and 'line3'. After this I've used strfind function with eval and sprintf to locate the character 'E' and assign the location to a variables 'p1', 'p2', 'p3'. Then the loop is ended. Now the coefficients are extracted by specifying the location and are converted to double using str2double. These coefficients are assigned to the variables ( from a1 to a14). After this I've specified the temperature array using linspace command and range is set between the local low temperature and local high temperature with 500 data points. Then I've called the used-defined function thermo properties which will return the values of specific heat, enthalpy and entropy and store them in the variables 'cp', 'h' and 's' respectively. Then I've called the second user-defined function molecular_weight which will return the value of molecular weight of the species and assign it to the variable 'mol_wt'. Now I've used pwd command to assign the current working directory to a variable 'current_directory'. Then I've changed the current working directory to the folder 'parsing_challenge' using command cd. Inside this folder I've created a folder 'species' using command mkdir. Then using cd command the current directory is changed to this new folder named 'species'. Inside this folder the plotting will be done and the graphs will be stored. Figure 1 has temperature vs specific heat graph, figure 2 has temperature vs enthalpy graph and figure 3 has temperature vs entropy graph for the species. Plotting is done using plot command. Linewidth is set to 3 for each plot. Temperature vs specific heat plot has default line color. For temperature vs enthalpy plot the color is set to black and for the temperature vs entropy plot the color is set to cyan. labelling for each plot is done using xlabel and ylabel. For the title of the plot I've used sprintf command along with title command to provide name of the species in the title. Now the graph from each figure is saved to the folder using saveas command. At last I've changed the current working directory to the initial directory 'current_directory' by using cd command. After this I've ended the for loop and then I've closed the file using fclose
for i = 1:53
species_header = (strsplit(fgetl(thermo_data),' '));
species = species_header{1};
local_low_temp = str2double(species_header{end-3});
local_mid_temp = str2double(species_header{end-1});
local_high_temp = str2double(species_header{end-2});
for j = 1:3
eval(sprintf('line%d = fgetl(thermo_data);',j))
eval(sprintf('p%d = strfind(line%d,''E'');',j,j))
end
a1 = str2double(line1(1:p1(1)+3));
a2 = str2double(line1(p1(1)+4:p1(2)+3));
a3 = str2double(line1(p1(2)+4:p1(3)+3));
a4 = str2double(line1(p1(3)+4:p1(4)+3));
a5 = str2double(line1(p1(4)+4:p1(5)+3));
a6 = str2double(line2(1:p2(1)+3));
a7 = str2double(line2(p2(1)+4:p2(2)+3));
a8 = str2double(line2(p2(2)+4:p2(3)+3));
a9 = str2double(line2(p2(3)+4:p2(4)+3));
a10 = str2double(line2(p2(4)+4:p2(5)+3));
a11 = str2double(line3(1:p3(1)+3));
a12 = str2double(line3(p3(1)+4:p3(2)+3));
a13 = str2double(line3(p3(2)+4:p3(3)+3));
a14 = str2double(line3(p3(3)+4:p3(4)+3));
T = linspace(local_low_temp,local_high_temp,500);
[cp,h,s] = thermo_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,local_mid_temp);
mol_wt = molecular_weight(species);
current_directory = pwd;
cd parsing_challenge
mkdir(species)
cd(species)
figure(1)
plot(T,cp,'linewidth',3)
xlabel('Temperature')
ylabel('Specific Heat')
title(sprintf('Temperature vs Specific Heat: %s',species))
figure(2)
plot(T,h,'color','k','linewidth',3)
xlabel('Temperature')
ylabel('Enthalpy')
title(sprintf('Temperature vs Enthalpy: %s',species))
figure(3)
plot(T,s,'color','c','linewidth',3)
xlabel('Temperature')
ylabel('Entropy')
title(sprintf('Temperature vs Entropy: %s',species))
saveas(1,'Specific_Heat.png')
saveas(2,'Enthalpy.png')
saveas(3,'Entropy.png')
cd(current_directory)
end
fclose(thermo_data);
Explanation for the user-defined function 'thermo_properties'
I've written this function to calculate the values of specific heat, enthalp and entropy for a species. Temperature coefficients (a1 to a14) and temperature values from array T and local_mid_temperature value for the species are the inputs specified for the function. Three values will be returned by the function. The value of universal gas constant R is also specified inside the function. If statement is used to specify the condition. If the input temperature is greater then local_mid_temp then high temperature coefficients (a1 to a7) are used for calculations and if input temperature is lower then the local_mid_temp then low temperature coefficients (a8 to a14) are used. Then formula are written in vectorized form.
function [cp, h, s] = thermo_properties(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,local_mid_temp)
R = 8.314;
if T > local_mid_temp
cp = R.*(a1 + a2.*T + a3.*T.^2 + a4.*T.^3 + a5.*T.^4);
h = R.*T.*(a1 + a2.*T./2 + a3.*T.^2./3 + a4.*(T.^3)./4 + a5.*(T.^4)./5 + a6./T);
s = R.*(a1.*log(T) + a2.*T + a3.*T.^2./2 + a4.*T.^3./3 + a5.*T.^4./4 + a7);
else
cp = R.*(a8 + a9.*T + a10.*T.^2 + a11.*T.^3 + a12.*T.^4);
h = R.*T.*(a8 + a9.*T./2 + a10.*T.^2./3 + a11.*T.^3./4 + a12.*T.^4./5 + a13./T);
s = R.*(a8.*log(T) + a9.*T + a10.*T.^2./2 + a11.*T.^3./3 + a12.*T.^4./4 + a14);
end
end
Explanation for the user-defined function 'molecular_weight'
This function is defined in order to calculate the molecular weight. The function will require one input (species) and return one output (molecular weight). The elements symbols are specified using a cell array 'elements'. Atomic weights of the elements are stored in a variable named atomic_wt. Then a variable 'mol_wt' is initialized with value equal to zero. Then for loop variable i is initialized from 1 to lenght of species vector. Afterwards the character from the string is converted to double using str2double and is assigned to a variable named 'atoms'. Then another for is used to with loop variable j. j is specied from 1 to length of vector 'elements'. Inside this for loop if statement is used. The condition specifies that if the characters from vectors species and elements match then use the defined formula. The function strcmpi is used to compare the characters of the string (case of the letter is ignored by strcmpi function). This function compares the strings and returns logical true(1) if there is match or false(0) if there is no match. The formula will just add the atomic weight of the elements that match. The atomic weight will be specified by the loop variable j. The location of the character from the vector 'elements' which has a match is then assigned to a variable named 'element_location' by using j. Then elseif is used to specify a condition that if the value of atoms is not 'not a number' (~isnan) which means if the converted value of vector 'atom' is a number then use the specified formula. The formula will add the atomic weight to the existing cumulative molecular weight by multiplying the atomic weight with number of atoms minus 1 (atoms - 1). This is done because one of the atom is already counted in the if statement where the elements matched. Atomic weight of the elements with more than one atoms will be calculate with this formula. I've used break statement to exit the loop after the calculation of the molecular weight and move on to the next element.
Then finally I've used fprintf to print the molecular weight of the species.
function [mol_wt] = molecular_weight(species)
elements = {'O','H','N','C','AR'};
atomic_wt = [15.999,1.00794,14.0067,12.0107,39.948];
mol_wt = 0;
for i = 1:length(species)
atoms = str2double(species(i));
for j = 1:length(elements)
if strcmpi(species(i),elements{j}(1))
mol_wt = (mol_wt + atomic_wt(j));
element_location = j;
elseif ~isnan(atoms)
mol_wt = mol_wt + atomic_wt((element_location))*(atoms-1);
break
end
end
end
fprintf('Molecular weight of %s is: %.f amu \n',species,round(mol_wt))
end
Following are the outputs returned after executing the program
Screenshots of the plots
Screenshots of the command window
Screenshot of speceis folder
