Project 1 - Parsing NASA thermodynamic data

In this project, I will be parsing a data file prepared by NASA. The contents of the data file can be used to generated thermodynamic properties such as Specific Heat at Constant Pressure 'C_p' (`J//(kg.K)`), Enthalpy `H or Q` (`J`) and Entropy `S` (`J//(kg.mol)`) at various temperatures. The files will be parsed in MATLAB and the molecular weight of each species will be calculated as well.

The data file can be accessed from here: https://drive.google.com/file/d/1fvkbSd3b3ryDMabqAroMgl_9HcOkBhLT/view?usp=sharing

The formulae used to calculate the properties are as below:

`C_p = R(a_1 + a_2T + a_3 T^2 + a_4T^3 + a_5T^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 ln(T) + a_2 T + a_3 T^2 /2 + a_4 T^3 /3 + a_5 T^4 /4 + a_7)`

Where, `R` is the Universal Gas Constant (`J//(mol.K)`), and `a_n` are the coefficients in the file for the particular species.

In the data file, when reading from the top left, there are 14 values. These values are the coefficients required. The first seven coefficients are for higher temperatures, while the rest seven are for lower temperatures. The temperature range is mentioned in the species line (or line 1).

Algorithm with Code and explanation

Algorithm

Step 1: Begin timer, intialise gas constant.

Step 2: Open the file "THERMO.dat" and extract global temperature limits.

Step 3: Close the file.

Step 4: Call function "thermo_data_extractor" .

Step 5: End timer.

Code

clear all
close all
clc

tic
% constants
R = 8.314; % gas constant J/(mol.K)

% openning the file to extract global limits
f1 = fopen('THERMO.dat','r');
fgetl (f1); % skipping file pointer to second line
line2 = fgetl(f1);
A = strsplit (line2,' '); % spliting data separated by a space " "
global_limits(1) = str2num(A{2});
global_limits(2) = str2num(A{3});
global_limits(3) = str2num(A{4});
fclose(f1);

thermo_data_extractor (global_limits, R); 

total_time = toc

Function thermo_data_extractor

Algorithm

Step 1: Assign parameter global temperatures to it's own varaibles.

Step 2: Open the file "THERMO.dat" and skip unwanted lines.

Step 3: Intialise flag as 1 and perform steps 4 to 17 while flag is 1.

Step 4: Read speciesline and split the line.

Step 5: Check if species is "CH2CHO" (which is the last compound), If so set flag to 0.

Step 6: Call function molecular_weight_calculator with species as parameter.

Step 7: Extract lower, middle and higher temperature limits for the current species.

Step 8: Extract five higher coefficients of the current species from line 2 in the file.

Step 9: Extract two higher coefficients and three lower coeffcients of the current species from line 3 in the file.

Step 10: Extract four lower coefficients of the current species from line 4 in the file.

Step 11: Calculate upper and lower temperature bounds from limits calculated in step 7.

Step 12: Calculate lower and upper specific heat, enthalpy and entropy for the current species.

Step 13: Create an export directory for current species.

Step 14: Plot specific heat vs temperature and save the plot to folder.

Step 15: Plot enthalpy vs temperature and save the plot to folder.

Step 16: Plot entropy vs temperature and save the plot to folder.

Step 17: Close open plots and display calculated molecular weight.

Step 18: Close file.

Code

function [] = thermo_data_extractor (global_limits, R)
    
    % assigning global temperatures
    global_low_temp = global_limits(1);
    global_mid_temp = global_limits(2);
    global_high_temp = global_limits(3);
    
    % opening file and skipping unwanted lines 
    f1 = fopen('THERMO.dat','r');
    fgetl(f1);
    fgetl(f1);
    fgetl(f1);
    fgetl(f1);
    fgetl(f1);

    % flag to determine end of file
    flag = 1;

    % this loop runs until flag's value changes
    while (flag == 1)
        % current line is speciesline. splitting line to extract
        % information
        speciesline = fgetl(f1);
        speciesline = strsplit(speciesline);
        % 1st cell contains the species' name
        if (speciesline(1) == "CH2CHO")
            flag = 0;                   % value of flag changes as file pointer is at last species
        end
        % calculating molecular weight of given species
        species = char(speciesline(1));
        [molecular_weight] = molecular_weight_calculator (species);
        
        % extracting lower, midddle and higher temperature limits of
        % current species
        local_low_limit = str2num(speciesline {end-3});
        local_high_limit = str2num(speciesline {end-2});
        local_mid_limit = str2num(speciesline {end-1});
    
        % parsing line 2 of current species. this line contains five high
        % temperature coefficients. these coeffcients are extracted
        % and saved in high_coeff
        t = fgetl(f1);
        a = strfind(t,'E');
        high_coeff(1) = str2num (t(a(1)-11:a(1)+3));
        high_coeff(2) = str2num (t(a(2)-11:a(2)+3));
        high_coeff(3) = str2num (t(a(3)-11:a(3)+3));
        high_coeff(4) = str2num (t(a(4)-11:a(4)+3));
        high_coeff(5) = str2num (t(a(5)-11:a(5)+3));
    
        % parsing line 3 of current species. this line contains two high
        % temperature coefficients and three lower temperaruture coefficients. 
        % these coeffcients are extracted
        % and saved in high_coeff and low_coeff respectively
        t = fgetl(f1);
        a = strfind(t,'E');
        high_coeff(6) = str2num (t(a(1)-11:a(1)+3));
        high_coeff(7) = str2num (t(a(2)-11:a(2)+3));
        low_coeff(1) = str2num (t(a(3)-11:a(3)+3));
        low_coeff(2) = str2num (t(a(4)-11:a(4)+3));
        low_coeff(3) = str2num (t(a(5)-11:a(5)+3));
    
        % parsing line 4 of current species. this line contains four low
        % temperature coefficients. these coeffcients are extracted
        % and saved in low_coeff
        t = fgetl(f1);
        a = strfind(t,'E');
        low_coeff(4) = str2num (t(a(1)-11:a(1)+3));
        low_coeff(5) = str2num (t(a(2)-11:a(2)+3));
        low_coeff(6) = str2num (t(a(3)-11:a(3)+3));
        low_coeff(7) = str2num (t(a(4)-11:a(4)+3));

        % creating temperature arrays from from limits
        t_lower = local_low_limit:local_mid_limit;
        t_upper = local_mid_limit:local_high_limit;

        % Cp/R = a1 + a2 T + a3 T^2 +a4 T^3 + a5 T^4
        % H/RT = a1 + a2 I /2 + a3 T^2 /3 + a4 I^3 /4 + a5 T^4 /5 + a6/T
        % S/R = a1 lnT + a2 T + a3 T^2 /2 + a4 T^3 /3 + a5 T^4 /4 + a7
        
        % calculating specific heat, enthalpy and entropy for lower
        % temperature
        cp_lower = R * (low_coeff(1) + low_coeff(2).*t_lower + low_coeff(3).*t_lower.^2 + low_coeff(4).*t_lower.^3 + low_coeff(5).*t_lower.^4);
        h_lower = R .* t_lower .* (low_coeff(1) + low_coeff(2).*t_lower./2 + low_coeff(3).*t_lower.^2./3 + low_coeff(4).*t_lower.^3./4 + low_coeff(5).*t_lower.^4./5+low_coeff(6)./t_lower);
        s_lower = R * (low_coeff(1).*log(t_lower) + low_coeff(2).*t_lower + low_coeff(3).*t_lower.^2./2 + low_coeff(4).*t_lower.^3./3 + low_coeff(5).*t_lower.^4./4 + low_coeff(7));
        
        % calculating specific heat, enthalpy and entropy for upper
        % temperature
        cp_upper = R * (high_coeff(1) + high_coeff(2).*t_upper + high_coeff(3).*t_upper.^2 + high_coeff(4).*t_upper.^3 + high_coeff(5).*t_upper.^4);
        h_upper = R .* t_upper .* (high_coeff(1) + high_coeff(2).*t_upper./2 + high_coeff(3).*t_upper.^2./3 + high_coeff(4).*t_upper.^3./4 + high_coeff(5).*t_upper.^4./5+high_coeff(6)./t_upper);
        s_upper = R * (high_coeff(1).*log(t_upper) + high_coeff(2).*t_upper + high_coeff(3).*t_upper.^2./2 + high_coeff(4).*t_upper.^3./3 + high_coeff(5).*t_upper.^4./4 + high_coeff(7));
        
        % creating export directory for current species
        mkdir(['O:\MATLAB\NASA_THERMO_OUTPUT\',char(speciesline(1))])
        figuresdir = fullfile('O:\MATLAB\NASA_THERMO_OUTPUT\',char(speciesline(1)));
        
        % plotting specific heat vs temperature
        figure(1)
        plot(t_lower,cp_lower,t_upper,cp_upper,'color','b')
        xlabel('Temperature (K)')
        ylabel('Specific Heat at Constant Pressure (J/(kg.K)')
        title ('Variation of Specific Heat Capacity for ',speciesline(1))
        saveas(gcf, fullfile(figuresdir, 'SpecificHeat'), 'jpeg'); % saving plot in current species' folder with appropiate file name
        
        % plotting enthalpy vs temperature
        figure(2)
        plot(t_lower,h_lower,t_upper,h_upper,'color','b')
        xlabel('Temperature (K)')
        ylabel('Enthalpy (J)')
        title ('Variation of Enthalpy for ',speciesline(1))
        saveas(gcf, fullfile(figuresdir, 'Enthalpy'), 'jpeg'); % saving plot in current species' folder with appropiate file name
        
        % plotting entropy vs temperature
        figure(3)
        plot(t_lower,s_lower,t_upper,s_upper,'color','b')
        xlabel('Temperature (K)')
        ylabel('Entropy (J/(K.mol))')
        title ('Variation of Entropy for ',speciesline(1))
        saveas(gcf, fullfile(figuresdir, 'Entropy'), 'jpeg'); % saving plot in current species' folder with appropiate file name
        close all

        fprintf('Molecular Weight of %s is %f\n',string(speciesline(1)),molecular_weight) % displaying calculated molecular weight
    end
    % closing file
    fclose(f1);
end

Function molecular_weight_calculator

Algorithm

Step 1: Intialise symbols of elements and their respective atomic weights.

Step 2: Traverse the entireity of charecters in the species' name for every element.

Step 3: Check if element is present in current species, if so adjust the molecular weight accordingly.

Step 4: If current position contains a number, adjust weight accordingly.

Step 5: Return molecular weight.

Code

function [molecular_weight] = molecular_weight_calculator (species)
    % this function takes the species as a parameter and sends the
    % molecular weight of the species

    % elements or atoms present in the file and their atomic weights
    atoms = {'H','C','N','O','A'};
    atomic_weight = [1.00784,12.0096,14.00643,15.99903,39.948];
    molecular_weight = 0; % intialising weight as 0 to prevent data corruption

    % for loop to run through each elements in the species
    for i = 1:length(species)
        % for loop to check if atoms are present in the current species
        for j = 1:length(atoms)
            if strcmpi(species(i),atoms(j))
                % if current atom is found in current species, the
                % appropiate atomic weight is called. if found, the
                % position is noted to check for multiple atoms later
                molecular_weight = molecular_weight + atomic_weight(j);
                pos = j;
            end
        end
        % if the current position contains a number, the weight is
        % multiplied with the number of atoms - 1
        n = str2num(species(i));
        if n>1
            molecular_weight = molecular_weight + atomic_weight(pos)*(n-1);
        end
    end
end

Output

Command Window

Molecular Weight of O is 15.999030
Molecular Weight of O2 is 31.998060
Molecular Weight of H is 1.007840
Molecular Weight of H2 is 2.015680
Molecular Weight of OH is 17.006870
Molecular Weight of H2O is 18.014710
Molecular Weight of HO2 is 33.005900
Molecular Weight of H2O2 is 34.013740
Molecular Weight of C is 12.009600
Molecular Weight of CH is 13.017440
Molecular Weight of CH2 is 14.025280
Molecular Weight of CH2(S) is 14.025280
Molecular Weight of CH3 is 15.033120
Molecular Weight of CH4 is 16.040960
Molecular Weight of CO is 28.008630
Molecular Weight of CO2 is 44.007660
Molecular Weight of HCO is 29.016470
Molecular Weight of CH2O is 30.024310
Molecular Weight of CH2OH is 31.032150
Molecular Weight of CH3O is 31.032150
Molecular Weight of CH3OH is 32.039990
Molecular Weight of C2H is 25.027040
Molecular Weight of C2H2 is 26.034880
Molecular Weight of C2H3 is 27.042720
Molecular Weight of C2H4 is 28.050560
Molecular Weight of C2H5 is 29.058400
Molecular Weight of C2H6 is 30.066240
Molecular Weight of CH2CO is 42.033910
Molecular Weight of HCCO is 41.026070
Molecular Weight of HCCOH is 42.033910
Molecular Weight of H2CN is 28.031710
Molecular Weight of HCN is 27.023870
Molecular Weight of HNO is 31.013300
Molecular Weight of N is 14.006430
Molecular Weight of NNH is 29.020700
Molecular Weight of N2O is 44.011890
Molecular Weight of NH is 15.014270
Molecular Weight of NH2 is 16.022110
Molecular Weight of NH3 is 17.029950
Molecular Weight of NO is 30.005460
Molecular Weight of NO2 is 46.004490
Molecular Weight of HCNO is 43.022900
Molecular Weight of HOCN is 43.022900
Molecular Weight of HNCO is 43.022900
Molecular Weight of NCO is 42.015060
Molecular Weight of CN is 26.016030
Molecular Weight of HCNN is 41.030300
Molecular Weight of N2 is 28.012860
Molecular Weight of AR is 39.948000
Molecular Weight of C3H8 is 44.091520
Molecular Weight of C3H7 is 43.083680
Molecular Weight of CH3CHO is 44.049590
Molecular Weight of CH2CHO is 43.041750

total_time =

  118.6513

>> 

Screenshot of Folder Structure

Plots for N2, O2 and CO2

N2

Specific Heat vs Temperature

Enthalpy vs Temperature

Entropy vs Temperature

O2

Specific Heat vs Temperature

Enthalpy vs Temperature

Entropy vs Temperature

CO2

Specific Heat vs Temperature

Enthalpy vs Temperature

Entropy vs Temperature

Folder containing all plots: https://drive.google.com/drive/folders/1MQyxiL5xk2-tPJQEfNa7C9YCTybyHNE7?usp=sharing

References

1. https://www.britannica.com/science/atomic-weight

2. https://www.mathworks.com/help/matlab/ref/disp.html

3. https://www.mathworks.com/matlabcentral/answers/423-cell-array-to-char-array#answer_594

4. https://www.mathworks.com/matlabcentral/answers/141403-cell-array-to-string-array#answer_324468

5. https://www.mathworks.com/help/matlab/ref/mkdir.html

6. https://www.mathworks.com/help/matlab/ref/fullfile.html