Parsing NASA thermodynamic data

1. Aim: File parsing the NASA thermodynamic data file to evaluate the thermodynamic properties such as specific Specific Heat (Cp), Enthalpy (H) and Entropy (S) for various species.

2. Objectives:

• To read the coefficients from the thermodynamic data file.
• To evaluate the polynomial equations using these coefficients.
• To plot the graphs between Cp, H and S with respect to the temperature values and dump these plots indenpendently in the corresponding species folder (directory).
• To calculate the molecular weight automatically of each species.

3. Theory:

1. File parsing is a technique to read the data from a particular file. In MATLAB, fopen(fileID), fclose(fileID) and fgetl(obj) reads one line of text from the instrument connected to obj. Here, the file ID is THERMO_new.txt file.
2. The formulae to calculate Cp, H and S are: `Cp = (a1 + (a2 * T) + (a3 *(T^2)) + (a4*(T^3)) + (a5*(T^4)))*R`, `H = (a1 + ((a2*T) /2) + ((a3*(T^2))/3) + ((a4*(T^3))/4) + ((a5*(T^4))/5) + (a6)/T)*R*T` and `S = ((a1*log(T)) + (a2*T) + ((a3*(T^2))/2) + ((a4*(T^3))/3) + ((a5*(T^4))/4) + a7)*R`, where a1, a2, a3, a4, a5, a6 and a7 are the high temperature coefficients.
3. Whereas, `Cp = (a8 + (a9*T) + (a10*(T^2)) + (a11*(T^3)) + (a12*(T^4)))*R`, `H = (a8 + ((a9*T)/2) + ((a10*(T^2))/3) + ((a11*(T^3))/4) + ((a12*(T^4))/5) + (a13)/T)*R*T` and `S = ((a8*log(T)) + (a9*T) + ((a10*(T^2))/2) + ((a11*(T^3))/3) + ((a12*(T^4))/4) + a14)*R`, where a8, a9, a10, a11, a12, a13 and a14 are the low temperature coefficients.

4. MATLAB program to extract the 14 coefficients and calculate the Cp, H and S polynomials:

• The following is the coded data_extract function along with the comments, which calculates the Cp, H and S polynomials of each species.
• It takes 3 output arguments as Cp, H and S and 17 input arguments as all 14 coefficients, universal gas constant (R in J/mol), local temperature values from the linspace command in the main program and local middle temperature value (1000 K).
• if and else condition is used to divide the 14 temperature coefficients.

% A program to 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

% A function named data_extract is used 
function [Cp, H, S] = data_extract (a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,local_mid_temp)
for i = length(T)
    if (T(i) > local_mid_temp)
      Cp = (a1 + (a2.*T) + (a3.*(T.^2)) + (a4.*(T.^3)) + (a5.*(T.^4))).*R;
      H = (a1 + ((a2.*T)./2) + ((a3.*(T.^2))./3) + ((a4.*(T.^3))./4) + ((a5.*(T.^4))./5) + (a6./T)).*R.*T;
      S = ((a1.*log(T)) + (a2.*T) + ((a3.*(T.^2))./2) + ((a4.*(T.^3))./3) + ((a5.*(T.^4))./4) + a7).*R;
    else
      Cp = (a8 + (a9.*T) + (a10.*(T.^2)) + (a11.*(T.^3)) + (a12.*(T.^4))).*R;
      H = (a8 + ((a9.*T)./2) + ((a10.*(T.^2))./3) + ((a11.*(T.^3))./4) + ((a12.*(T.^4))./5) + (a13./T)).*R.*T;
      S = ((a8.*log(T)) + (a9.*T) + ((a10.*(T.^2))./2) + ((a11.*(T.^3))./3) + ((a12.*(T.^4))./4) + a14).*R;
    end
end
end

5. MATLAB program to calculate the molecular weight of each species:

• The following is the coded mol_weight along with the comments, which calculates the molecular weight of each species.
• It takes 1 output argument as molecular_weight and 1 input argument as species_name.
• Elements like H, C, N, O, N, AR and S are stored in an array named elements.
• The atomic weights of these elements are stored in an array named atomic_weight.
• The molecular weight of these elements will be equal to their atomic weights as they consist of a single atom.
• The molecular weights of other species like O2, N2, CO2, CH3CHO etc. are calculated using the formulae given in the if and elseif conditions with their respective names.
• Lastly, the molecular weight of each species is displayed using disp(variable_name) function.

% A function to extract the molecular weight of the species using
% mol_weight

function [molecular_weight] = mol_weight(species_name)

% assigning the elements in string array
elements = ['H','C','O','N','AR','S'];

% assigning the atomic weights of the elements in atomic_weight array
atomic_weight = [1.00794 12.0107 15.999 14.0067 39.948 32.065];

for i = 1:length(species_name)
    for j = 1:length(atomic_weight)
        
        % string comparison to match the species name and the element's
        % name
        if strcmp(species_name(i),elements(j))
            
            % calculating the molecular weight
            molecular_weight = atomic_weight(j);
            
            % breaking the for loop when the if conditions becomes true
            break;
        end
    end
   
   % convering the cell array to matrix array using the cell2mat function 
   n = cell2mat(species_name(i));
   
   % calculating the size of the n array 
   sz = size(n);

   % using if and else if conditions to calculate the molecular weight of
   % the species other than that of the elements and different molecular
   % weight formulae has been formulated as per the species number of atoms
   % and molecules respectively
    if (sz(2) > 1)
        if (n == "H2")
           molecular_weight = (atomic_weight(1))*2;
        elseif (n == "C2")
            molecular_weight = (atomic_weight(2))*2;
        elseif (n == "O2")
            molecular_weight = (atomic_weight(3))*2;
        elseif (n == "N2")
            molecular_weight = (atomic_weight(4))*2;
        elseif (n == "AR")
            molecular_weight = (atomic_weight(5));
        elseif (n == "OH")
            molecular_weight = atomic_weight(3) + atomic_weight(1);
        elseif (n == "H2O")
            molecular_weight = ((atomic_weight(1))*2) + atomic_weight(3);
        elseif (n == "HO2")
            molecular_weight = ((atomic_weight(3))*2) + atomic_weight(1);
        elseif (n == "H2O2")
            molecular_weight = ((atomic_weight(1))*2) + ((atomic_weight(3))*2);
        elseif (n == "CH")
            molecular_weight = atomic_weight(2) + atomic_weight(1);
        elseif (n == "CH2")
            molecular_weight = ((atomic_weight(1))*2) + atomic_weight(2);
        elseif (n == "CH2(S)")
            molecular_weight = ((atomic_weight(1))*2) + atomic_weight(2) + atomic_weight(6);
        elseif (n == "CH3")
            molecular_weight = ((atomic_weight(1))*3) + atomic_weight(2);
        elseif (n == "CH4")
            molecular_weight = ((atomic_weight(1))*4) + atomic_weight(2);
        elseif (n == "CO")
            molecular_weight = atomic_weight(2) + atomic_weight(3);
        elseif (n == "CO2")
            molecular_weight = ((atomic_weight(3))*2) + atomic_weight(2);
        elseif (n == "HCO")
            molecular_weight = atomic_weight(1) + atomic_weight(2) + atomic_weight(3);
        elseif (n == "CH2O")
            molecular_weight = ((atomic_weight(1))*2) + atomic_weight(2) + atomic_weight(3);
        elseif (n == "CH2OH" || n == "CH3O")
             molecular_weight = ((atomic_weight(1))*3) + atomic_weight(2) + atomic_weight(3);
        elseif (n == "CH3OH")
            molecular_weight = ((atomic_weight(1))*4) + atomic_weight(2) + atomic_weight(3);
        elseif (n == "C2H")
             molecular_weight = ((atomic_weight(2))*2) + atomic_weight(1);
        elseif (n == "C2H2")
            molecular_weight = ((atomic_weight(1))*2) + (atomic_weight(2))*2;
        elseif (n == "C2H3")
            molecular_weight = ((atomic_weight(1))*3) + (atomic_weight(2))*2;
        elseif (n == "C2H4")
            molecular_weight = ((atomic_weight(1))*4) + (atomic_weight(2))*2;
        elseif (n == "C2H5")
            molecular_weight = ((atomic_weight(1))*5) + (atomic_weight(2))*2;
        elseif (n == "C2H6")
            molecular_weight = ((atomic_weight(1))*6) + (atomic_weight(2))*2;
        elseif (n == "CH2CO" || n == "HCCOH")
            molecular_weight = ((atomic_weight(1))*2) + ((atomic_weight(2))*2) + atomic_weight(3);
        elseif (n == "HCCO")
            molecular_weight = atomic_weight(1) + ((atomic_weight(2))*2) + atomic_weight(3);
        elseif (n == "H2CN")
            molecular_weight = ((atomic_weight(1))*2) + atomic_weight(2) + atomic_weight(4);
        elseif (n == "HCN")
             molecular_weight = atomic_weight(1) + atomic_weight(2) + atomic_weight(4);
        elseif (n == "HNO")
            molecular_weight = atomic_weight(1) + atomic_weight(3) + atomic_weight(4);
        elseif (n == "NNH")
            molecular_weight = ((atomic_weight(4))*2) + atomic_weight(1);
        elseif (n == "N2O")
            molecular_weight = ((atomic_weight(4))*2) + atomic_weight(3);
        elseif (n == "NH")
             molecular_weight = atomic_weight(1) + atomic_weight(4);
        elseif (n == "NH2")
             molecular_weight = ((atomic_weight(1))*2) + atomic_weight(4); 
        elseif (n == "NH3")
             molecular_weight = ((atomic_weight(1))*3) + atomic_weight(4);  
        elseif (n == "NO")
             molecular_weight = atomic_weight(3) + atomic_weight(4);
        elseif (n == "NO2")
             molecular_weight = ((atomic_weight(3))*2) + atomic_weight(4);
        elseif (n == "HCNO" || n == "HOCN" || n == "HNCO")
            molecular_weight = atomic_weight(1) +  atomic_weight(2) + atomic_weight(3) + atomic_weight(4);
        elseif (n == "NCO")
            molecular_weight = atomic_weight(2) + atomic_weight(3) + atomic_weight(4);
         elseif (n == "CN")
             molecular_weight = atomic_weight(2) + atomic_weight(4);
        elseif (n == "HCNN")
            molecular_weight = ((atomic_weight(4))*2) + atomic_weight(1) + atomic_weight(2);
        elseif (n == "C3H7")
            molecular_weight = ((atomic_weight(1))*7) + (atomic_weight(2))*3;
        elseif (n == "C3H8")
            molecular_weight = ((atomic_weight(1))*8) + (atomic_weight(2))*3;
        elseif (n == "CH2CHO")
            molecular_weight = ((atomic_weight(1))*3) + ((atomic_weight(2))*2) + atomic_weight(3);
        elseif (n == "CH3CHO")
            molecular_weight = ((atomic_weight(1))*4) + ((atomic_weight(2))*2) + atomic_weight(3);   
        end 
    end
    
    % to display the molecular weight of the species with their names, an
    % cell array is created having strings 
    X = ['Molecular weight of the',' ', species_name(i),' ', 'is', ' ', num2str(molecular_weight)];
    
    % converting the cell array to the matrix array
    Y = cell2mat(X);
    
    % displaying the Y array using disp function
    disp(Y);
    
end
end

6. Main MATLAB program:

• The following is the main MATLAB program along with the comments, which plots the graphs of Cp, H and S with respect to the local temperature values.
• It calls the data_extract function and mol_weight function.
• The plots of each species are dumped into the respective species folder along with the axis labels and titles.
• pwd command is used to get the current working directory and cd command is used to the change the directory.

% A file parsing program to plot the graphs of the calculated enthalpy, entropy and 
% specific heat from the thermodynamic data file.

% clearing the workspace
clear all

% closing the plots
close all

% clearing the command window
clc

% to read the data file
file = fopen('THERMO_new.txt','r');

% to get the first line from the line
f_line = fgetl(file);

% for global temperature values
g_temp = fgetl(file);

% converting string array to number array
g_temp_num = str2num(g_temp);

% using for loop and if, elseif and else conditions to store the low,
% middle and high global temperatures
for i = 1:length(g_temp_num)
   if (g_temp_num(i) < 1000)
       g_low_temp = g_temp_num(i);
   elseif (g_temp_num(i) > 1000)
       g_high_temp = g_temp_num(i);
   else
       g_mid_temp = g_temp_num(i); 
   end
end

% skipping the lines
for j = 1:3
    tline = fgetl(file);
end

% number of species
num_species = 53;

% getting the lines from the file one by one and storing the coefficient
% values using the for loop 
for k = 1:num_species
    
  % getting the line
  r1_tline = fgetl(file);
  
  % splitting the tline string array
  tline_split = strsplit(r1_tline, ' ');
  
  % calculating the size of the tline_split
  len = size(tline_split);
  
  % extracting the species name
  species_name = tline_split(1);
  
  % converting the tline_split string to double precision values using 
  % str2double for the local low, middle and high temperatures respectively
  % local low, middle and high temperatures respectively
  local_low_temp = str2double(tline_split(len(2) - 3));
  local_high_temp = str2double(tline_split(len(2) - 2));
  local_mid_temp = str2double(tline_split(len(2) - 1));  
  
  % getting the first line of the coefficients 
  r2_tline = fgetl(file);
  
  % converting character array or string array to numeric array using
  % str2num
  r2_tline_array = str2num(r2_tline);
  
  % storing the 5 coeffcients as a1,a2,a3,a4,a5 respectively
  a1 = r2_tline_array(1);
  a2 = r2_tline_array(2);
  a3 = r2_tline_array(3);
  a4 = r2_tline_array(4);
  a5 = r2_tline_array(5);
  
  % getting the second line of the coefficients 
  r3_tline = fgetl(file);
  
  % converting character array or string array to numeric array using
  % str2num
  r3_tline_array = str2num(r3_tline);
  
  % storing the next 5 coeffcients as a6,a7,a8,a9,a10 respectively
  a6 = r3_tline_array(1);
  a7 = r3_tline_array(2);
  a8 = r3_tline_array(3);
  a9 = r3_tline_array(4);
  a10 = r3_tline_array(5); 
  
  % getting the second line of the coefficients 
  r4_tline = fgetl(file);
  
  % converting character array or string array to numeric array using
  % str2num
  r4_tline_array = str2num(r4_tline);
  
  % storing the next 4 coeffcients as a11,a12,a13,a14 respectively
  a11 = r4_tline_array(1);
  a12 = r4_tline_array(2);
  a13 = r4_tline_array(3);
  a14 = r4_tline_array(4);  
  
  % to get the current directory using pwd command
  current_dir = pwd;
  
  % storing the species name in thermo_data to get the folder with the
  % species name using sprintf function
  thermo_data = sprintf('%s',species_name{1});
  
  % to create a folder with the species name 
  mkdir(thermo_data);
     
  % dividing the local temperature values into 100 divisions from local low
  % temperature to local high temperature using linspace function
  T = linspace(local_low_temp,local_high_temp,100);

  % univeral gas constant in J/mol
  R = 8.314;

  % calculating specific heat, enthalpy, entropy using data_extract function
  [Cp,H,S] = data_extract(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,T,local_mid_temp);

  % plotting the graphs
  % opening the Figure 1 window
  figure(1)
  
  % plotting the graph between specific heat and local temperature using
  % red colour and linewidth as 2
  plot(Cp,T,'color','r','linewidth',2);
  
  % adding title as Specific heat Vs Temperature
  title(['Specific heat Vs Temperature plot of',' ',species_name{1}]);
  
  % labelling the X - axis
  xlabel('Temperature');
  
  % labelling the Y - axis
  ylabel('Specific heat');
  
  % turning the grid on
  grid on

  % opening the Figure 2 window
  figure(2)
  
  % plotting the graph between enthalpy and local temperature using
  % green colour and linewidth as 2
  plot(H,T,'color','g','linewidth',2);
  
  % adding title as Enthalpy Vs Temperature
  title(['Enthalpy Vs Temperature plot of',' ',species_name{1}]);
  
  % labelling the X - axis 
  xlabel('Temperature');
  
  % labelling the Y - axis
  ylabel('Enthalpy');
  
  % turning the grid on 
  grid on

  % opening the Figure 3 window
  figure(3)
  
  % plotting the graph between entropy and local temperature using
  % blue colour and linewidth as 2
  plot(S,T,'color','b','linewidth',2);
  
  % adding title as Entropy Vs Temperature
  title(['Entropy Vs Temperature plot of',' ',species_name{1}]);
  
  % labelling the X - axis
  xlabel('Temperature');
  
  % labelling the Y - axis
  ylabel('Entropy');
  
  % turning the grid on
  grid on
 
  % changing the directory to save the plots in the corresponding species
  % folder
  cd(thermo_data);
  
  % saving the plots
  saveas(1,'Specific heat.jpg');
  saveas(2,'Enthalpy.jpg');
  saveas(3,'Entropy.jpg');
  
  % change the current folder back to the original folder and using only cd
  % command original folder can be displayed
  cd(current_dir);
   
  % calculating the molecular weight of species using mol_weight function 
  molecular_weight = mol_weight(species_name);
end

% closing the file
fclose(file);

7. Results:

The following are the Cp, H and S plots with respect to the local temperature values of O2:

The following are the Cp, H and S plots with respect to the local temperature values of N2:

The following are the Cp, H and S plots with respect to the local temperature values of CO2:

The following images shows the species folder in which the plots are dumped:

The following images shows the molecular weight of each species displayed in the command window of MATLAB:

8. Conclusion:

• The plots of Cp, H and S with respect to the local temperature values are plotted and dumped into the respective species folder.
• The nature of the plots will vary from species to species according to the coefficients given to the polynomials.
• The molecular weights of all species are calculated and displayed in the command window.

9. References:

1. fgetl - https://www.mathworks.com/help/instrument/fgetl.html
2. cell2mat - https://www.mathworks.com/help/matlab/ref/cell2mat.html
3. str2num - https://www.mathworks.com/help/matlab/ref/str2num.html
4. strsplit - https://www.mathworks.com/help/matlab/ref/strsplit.html
5. strcmp - https://www.mathworks.com/help/matlab/ref/strcmp.html
6. fopen - https://www.mathworks.com/help/matlab/ref/fopen.html
7. fclose - https://www.mathworks.com/help/matlab/ref/fclose.html
8. pwd - https://www.mathworks.com/help/matlab/ref/pwd.html
9. cd - https://www.mathworks.com/help/matlab/ref/cd.html
10. MATLAB Answers - https://www.mathworks.com/matlabcentral/answers/74938-mkdir-with-variable-name
11. NASA Thermodynamic data - https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19940013151.pdf