Parsing NASA Thermodynamic Data
- Aim : Matlab program to parse NASA thermodynamic data, evaluating thermodynamic properties such as C_p, H, S for different elements and molecules, plotting and saving the evaluated properties.
- Equations & Expressions : 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 different species. The polynomials are the following:
- `C_p/R = a_1 + a_2*T + a_3*T^2 + a_4*T^3 + a_5*T^4`
- `H/(RT) = a_1 + a_2*T/2 + a_3*T^2/3 + a_4*T^3/4 + a_5*T^4/5 + a_6/T`
- `S/R = a_1*lnT + a_2*T + a_3*T^2/2 + a_4*T^3/3 + a_5*T^4/4 + a_7`
The data contains 14 cofficients, the first seven for high temperature range and the second seven for low temperature range. We calculate the thermodynamic properties for the whole temperature range and then plot these properties with the respect to the temperatue range they are calculated for.
- Objectives & Workflow : In this project we will be fetching the cofficients by parsing NASA thermodynamic data to evaluate the C_p, H, S for different species. We will be plotting these evaluated properties and saving them for further analysis.
- The workflow of the project is as follows :
- We open the file using fopen command. We then evaluate the number of species for which we gonna evaluate and plot the thermodynamic properties.
- We read the header of the file containing the global temperature range so over that we can move to fetch the cofficients. To calculate the molecular weight for species we need to create variable with elements name and atomic mass for comparison.
- We fetch the local temperature range and all fourteen cofficients for each species and then call the function 'cal' for giving the required output. We take the properties for over the whole temperature range.
- We plot the C_p, H, S for each species with respective to their local temperature range and save them in separate folder with respective names of the species.
- At last we calculate the molecular mass of each species and display it in the command window.
- Matlab Program :
The programming language used is Matlab scripting language.
%Program to parse NASA thermodynamic data, plotting and saving the plots
clear all;
clc;
file = fopen('THERMO.dat','r'); %opening and indexing the given data file
%Calculating the total number of species in the data
tot_lines = 0; %count variable for counting total number of lines
check = ''; %compare variable for comparing lines to the last line
while strcmp('END',check)==0 %creating infinite loop with termination condition
check = fgetl(file); %fetching each line from the data line by line
tot_lines = tot_lines+1; %incrementing line count
end
tot_lines = (tot_lines - 6)/4; %the number of species calculated
frewind(file); %setting the file position indicator to the begining of the file
%Retrieving the header data of the data file
for m=1:5 %using the loop for indexing the header data lines
if m==2 %condition for specific line of global temperature range
line=fgetl(file); %fetching the line
cell=strsplit(line,' '); %splitting the string
global_low_temp = str2double(cell{2}); %assigning the temperature to their respective variable
global_mid_temp = str2double(cell{3});
global_high_temp = str2double(cell{4});
else
fgetl(file); %fetching the remaining lines
end
end
cell_name={'H','C','N','O','S','A','2','3','4','5','6','7','8'};
mass_value=[1,12,14,16,32,40];
%fetching the species cofficients, calculating H,S,C_p, plotting and saving
for l = 1:tot_lines %indexing each species
count=1; %variable for incrementing the cofficient variable
for n=1:4 %indexing each species's cofficients
if n==1 %condition for fetching species name and temperature range
s_line = fgetl(file); %fetching the line
s_cell = strsplit(s_line,' '); %splitting the string
species = string(s_cell(1)); %assigning the species name to variable
loc_l_temp = str2double(s_cell{end-3}); %assigning the temp. to respective variable
loc_m_temp = str2double(s_cell{end-1});
loc_h_temp = str2double(s_cell{end-2});
else
s_line = fgetl(file); %fetching the further lines containing cofficients
mk = strfind(s_line,'E'); %using strfind to create a marker for a cofficient
mk = [-3 mk]; %completing the marker variable for further calculation
for k = 1:length(mk)-1 %indexing each cofficient for one species in a line
s_coff = s_line(mk(k)+4:mk(k+1)+3); %picking up the cofficient from the string
coff(count) = str2double(s_coff); %assigning the cofficients to a variable
count=count+1; %incrementing the count variable
end
end
end
%calculating the C_p, H, S for each species, calling the function cal
l_temp = (loc_l_temp:loc_m_temp-1); %creating low and high temperature range
h_temp = (loc_m_temp:loc_h_temp);
[C_p_h, H_h, S_h, C_p_l, H_l, S_l] = cal(l_temp,h_temp,coff); %calling function
C_p=[C_p_l, C_p_h]; %C_p calculated for both temperature range
H=[H_l, H_h]; %H ans S calculated for both temperature range
S=[S_l, S_h];
temp=[l_temp, h_temp];
%plottting the calculated C_p,S,H in same figure window for each species
figure('Name',species,'NumberTitle','off'); %creating figure window specific for each species
set(gcf,'Position',[350 80 700 600]); %changing window size
subplot(3,1,1); %creating subplot in the window
plot(temp,C_p,'linewidth',1); %plotting C_p with temperature range
title('C-p(Specific Heat)'); %providing title to the plot
subplot(3,1,2);
plot(temp,H,'linewidth',3); %plotting H and S in the futher subplots and giving title
title('H(Enthalpy)');
subplot(3,1,3);
plot(temp,S,'linewidth',3);
title('S(Entropy)');
%saving each plot in respective species folder by first creating it
mkdir('figures',species); %making new directory for each species
saveas(gcf,fullfile('E:Matlabfigures',species,strcat(species,'.png')),'png'); %saving plots of species in their respective folders
%calculating each species molecular mass/weight
p=1; %count variable
cmpr=num2cell(char(species)); %breaking the species name and creating cell
for i = 1:length(cmpr) %indexing the newly created cell
logic=strcmp(cmpr(i),cell_name); %comparing the cell with elements
for j = 1:length(logic) %indexing the logic variable
if logic(j)==1 %condition for finding the match
if j>=7 %condition for repeated element
values(p)=values(p-1)*(str2double(cell_name{j})-1); %assigning the value of repeated element
p=p+1; %incrementing the count variable
else
values(p)=mass_value(j); %assigning the mass value of element
p=p+1; %incrementing the count variable
end
break
end
end
end
mass=sum(values); %adding the mass values of each element of species
clear values; %clearing the variable for next iteration/species
fprintf('%s have %d molecular mass.n',species,mass)
end
- Explanation/Description :
- First of all we need to open the file for parsing the data using fopen command. We are reading this file so in fopen we pass the file name and 'r'. We index the data file by variable file which return a integer greater than 3 after opening the file. If the variable return -1 then it failed to open the file.
- We calculate the number of species for which we evaluate and plot. We create two variables, one which will keep the count of the lines 'tot_lines' and second to store each line as string.
while strcmp('END',check)==0 %creating infinite loop with termination condition check = fgetl(file); %fetching each line from the data line by line tot_lines = tot_lines+1; %incrementing line count end - Our data ends with a keyword 'END' so we use a while loop with the condition of string comparison of every fetched line and keyword 'END'. The loop will perform iterations while the conditions is true which is the comparison of check variable and string 'END' is false. The check variable fetch the next line on each iteration using fgetl command and the count variable will get incremented. In this way we get the total number of lines.
- Now from the total number of lines we know that the first five lines are the header data further which each species takes four lines so doing the math gives us the number of species. Also at the last iteration when the check variable fetch the 'END' the count variable get incremented that's why we subtract six from total number of lines. The frewind command will set the file position indicator to the begining of the file.
for m=1:5 %using the loop for indexing the header data lines if m==2 %condition for specific line of global temperature range line=fgetl(file); %fetching the line cell=strsplit(line,' '); %splitting the string global_low_temp = str2double(cell{2}); %assigning the temperature to their respective variable global_mid_temp = str2double(cell{3}); global_high_temp = str2double(cell{4}); else fgetl(file); %fetching the remaining lines end end - We use the for loop for indexing the header data lines. Then we put if condition specifically for line number two in which there is global temperature range. We use the fgetl command to store the whole line in variable line as a string. So to fetch the temperature we need to split the string with the delimiter being ' ' and store it into variable cell.
- After storing we just need to assign the respective temperature to the variables by converting them into intergers. We need to take care of using the curly bracket for indexing the cell. Otherwise the condition we just need to move down to further lines using the fgetl command.
- We will be using for loop to index each species, on each iteration from fetching the cofficients, evaluating thermodynamic properties to plotting, saving and calculating molecular mass will be done for one specie.
- For calculating the mass we need to define some variables before the loop containing the elements and their mass which we will use for comparison with different species.
-
for n=1:4 %indexing each species's cofficients if n==1 %condition for fetching species name and temperature range s_line = fgetl(file); %fetching the line s_cell = strsplit(s_line,' '); %splitting the string species = string(s_cell(1)); %assigning the species name to variable loc_l_temp = str2double(s_cell{end-3}); %assigning the temp. to respective variable loc_m_temp = str2double(s_cell{end-1}); loc_h_temp = str2double(s_cell{end-2}); else s_line = fgetl(file); %fetching the further lines containing cofficients mk = strfind(s_line,'E'); %using strfind to create a marker for a cofficient mk = [-3 mk]; %completing the marker variable for further calculation for k = 1:length(mk)-1 %indexing each cofficient for one species in a line s_coff = s_line(mk(k)+4:mk(k+1)+3); %picking up the cofficient from the string coff(count) = str2double(s_coff); %assigning the cofficients to a variable count=count+1; %incrementing the count variable end end end - There are four lines for each species containing the name of species and local temperature range in first line and cofficients in further three lines. So we use the for loop to index each line, further in the loop we use if condition to retrieve the name and temperature range.
- We use the fgetl command to read the line and store it in s_line as a character vector. Then we use the string split command to break it with delimiter being ' ' and save the cell output in variable s_cell. Now we just need to assign the name and temperature to respective variables after converting to required datatype. We use string command to convert the species name from cell to string before saving it in variable 'species'. We use the str2double to change the string into interger.
- After first iteration, for n=2 the if condition will fail and the else statements will be executed. Then using the fgetl command we read the next line and overwrite it on variable s_line. There are four to five cofficients in each line. In each cofficient we have common is 'E', before and after of it consist one cofficient.
- So we use strfind command to find at what places in the character vector E is so that we can use it as a mark to fetch the cofficient. We save these values in variable mk. So to retrieve one cofficient we just need to find at what places the cofficient start and end in the character vector s_line. Each cofficient ends after 3 places of E but also at the fourth place next cofficient start so using a for loop of the length of mk variable we can fetch all cofficients in a line.
- We need to add -3 in the variable mk to make the above logic work. The first value of mk will give us the starting point of the second cofficient, so to make that value 1 for the very first cofficient of the line we add -3 in mk which after addition with +4 gives 1.
- We then convert it into integer and save it in variable coff which will save all the 14 cofficients for one specie as the count variable will get incremented after each iteration. The above loop will iterate for 14 times as it is nested in another loop.
-
l_temp = (loc_l_temp:loc_m_temp-1); %creating low and high temperature range h_temp = (loc_m_temp:loc_h_temp); [C_p_h, H_h, S_h, C_p_l, H_l, S_l] = cal(l_temp,h_temp,coff); %calling function C_p=[C_p_l, C_p_h]; %C_p calculated for both temperature range H=[H_l, H_h]; %H ans S calculated for both temperature range S=[S_l, S_h]; temp=[l_temp, h_temp]; - Now we have all 14 cofficients for the specie in the variable coff. We will be evaluating porperties for both low and high temperature range. So we create low and high temperature range using range operator.
- We call the function 'cal' which evaluates the polynomials taking the inputs of low and high temperature range and cofficients. The function gives the output of C_p, H, S for both low and high temperature which then we combine to get results over the whole temperature range.
-
figure('Name',species,'NumberTitle','off'); %creating figure window specific for each species set(gcf,'Position',[350 80 700 600]); %changing window size subplot(3,1,1); %creating subplot in the window plot(temp,C_p,'linewidth',1); %plotting C_p with temperature range title('C-p(Specific Heat)'); %providing title to the plot - Now comes the plotting of evaluated properties of the specie. We create a figure window with specific characterstics using the command figure. The species name will be appear as the title of the window and the number will not be displayed.
- We change the dimensions of the window using the set position command. We use the subplot command to create partitions in the same window for plotting all three properties. Then we plot the properties with providing specific titles using the title command.
-
mkdir('figures',species); %making new directory for each species saveas(gcf,fullfile('E:Matlabfigures',species,strcat(species,'.png')),'png'); %saving plots of species in their respective folders - Now we will save this figure window specific for each species. We create a new specific directory for each species in which we save the plot using the command mkdir. In the command the figure is the parent folder in which further folders will be created specific for each species.
- After creating the folder now comes the saving part for which we use the saveas command. In the command itself we use fullfile command to specify the address and name of the file. The saveas command will save the current figure window in current directory in which the program is saved.
- So to save the figure in their respective folder we use fullfile command also to combine the species name and file extension use the strcat command.
-
for i = 1:length(cmpr) %indexing the newly created cell logic=strcmp(cmpr(i),cell_name); %comparing the cell with elements for j = 1:length(logic) %indexing the logic variable if logic(j)==1 %condition for finding the match if j>=7 %condition for repeated element values(p)=values(p-1)*(str2double(cell_name{j})-1); %assigning the value of repeated element p=p+1; %incrementing the count variable else values(p)=mass_value(j); %assigning the mass value of element p=p+1; %incrementing the count variable end break end end end - Now we are calculating the molecular mass of the species. We break the specie's name string into single characters using the command num2cell and save those in variable cmpr. We will compare each character in the variable cmpr with the already created variable cell_name. The cell_name contains the elements name and numeric values that can be found in any species name. There is another variable called mass_value which is created in such way that mass value of each element is placed at the same position where its element is placed in the variable cell_name.
- We then use the string compare command to compare each character of variable cmpr with the cell array cell_name giving us the output with a logical vector containing 1 and 0. The output will be one when the comparison is true so if we can just find out at what position in variable logic one is we will be able to know which element it is. Then we can also take that same position for mass_value to know its mass.
- The problem arises when it is a numeric value then we need to add the mass value of the previous element for which we use an if condition because in the cell_name the numerics are after a certain position. So if condition is true then we add the previous element mass, one time less the value of the numeric itself. Thus in this way we save all the mass values for each element in the variable 'values'.
- After adding the mass in the variable values, we display it in the command window. Further which we need to clear the variable values so that in the next iteration for next species there is no mass value is remaining in it.
- Cal Function :
function [C_p_h, H_h, S_h, C_p_l, H_l, S_l] = cal(l_temp,h_temp,coff)
R = 8.3145;
C_p_h = (coff(1)+coff(2).*h_temp+coff(3).*(h_temp).^2+coff(4).*(h_temp).^3+coff(5).*(h_temp).^4).*R;
H_h = (coff(1)+ coff(2).*h_temp/2+ coff(3)*(h_temp.^2)/3+ coff(4)*(h_temp.^3)/4+ coff(5)*(h_temp.^4)/5+ coff(6)./h_temp).*(R*h_temp);
S_h = (coff(1)*log(h_temp)+ coff(2)*h_temp+ coff(3)*(h_temp.^2)/2+ coff(4)*(h_temp.^3)/3+ coff(5)*(h_temp.^4)/4+ coff(7)).*R;
C_p_l = (coff(8)+coff(9).*l_temp+coff(10).*(l_temp).^2+coff(11).*(l_temp).^3+coff(12).*(l_temp).^4).*R;
H_l = (coff(8)+ coff(9).*l_temp/2+ coff(10)*(l_temp.^2)/3+ coff(11)*(l_temp.^3)/4+ coff(12)*(l_temp.^4)/5+ coff(13)./l_temp).*(R*l_temp);
S_l = (coff(8)*log(l_temp)+ coff(9)*l_temp+ coff(10)*(l_temp.^2)/2+ coff(11)*(l_temp.^3)/3+ coff(12)*(l_temp.^4)/4+ coff(14)).*R;
end
- Explanation :
In this function we take the input of both low - high temperature range and cofficients of the species. We write the equations for calculating the C_p, H, S for both high and low temperture range and take the output in respective variables. To take the output over the whole temperature range we then club the ouputs of high and low temperature in the main program.
- Outputs/Attachments :
- References :
- https://in.mathworks.com/matlabcentral/answers/392620-how-to-save-multiple-matlab-figures-in-separate-folders
- https://in.mathworks.com/help/matlab/ref/mkdir.html?searchHighlight=mkdir&s_tid=srchtitle
- https://in.mathworks.com/matlabcentral/answers/541832-operator-not-supported-for-operands-of-type-cell?s_tid=answers_rc1-1_p1_Topic
- https://in.mathworks.com/help/matlab/ref/saveas.html?searchHighlight=saveas&s_tid=srchtitle#bup2b2t-1
- https://in.mathworks.com/help/matlab/ref/fullfile.html