Parsing NASA thermodynamic data using MATLAB

Parsing NASA thermodynamic data using MATLAB

Main code -

clear all
close all
clc

% Input 
R = 8.314; % Universal gas constant

% skipping the heading
f = fopen(\'Thermo.dat\',\'r\'); % fopen - opens target file for reading and writing

Header = fgetl(f);  % fgetl - reads line from file % reading the header

global_temp_range = fgetl(f);
global_temp_range = str2num(global_temp_range); % str2num - string to numeric array

% skipping the comments
line1 = fgetl(f);
line2 = fgetl(f);
line3 = fgetl(f);

% creating a loop for all 53 chemical elements

for i = 1:53
    
    line0 = fgetl(f);
    chemical_element = strtok(line0); % strtok - to select part of string
    
    % extracting the local temp from each element
    
    a = strfind(line0, \'G   \'); % strfind - to find string / we can use \'findstr\' also
    local_temp = line0(a+1:end-1);
    local_temp = str2num(local_temp);
    
    % local temp range for plotting the graph
    
    t = linspace(local_temp(1),local_temp(2),1000); % linspace - generates arrays
    
    % reading the coefficicents
    
    line1 = fgetl(f);
    line2 = fgetl(f);
    line3 = fgetl(f);
    
    % finding the string \'E\'
    
    a = strfind(line1,\'E\');
    b = strfind(line2,\'E\');
    c = strfind(line3,\'E\');
    
    % finding individual coefficients
    % line1 
    % str2double - converts strings to duble precision values

    a1 = str2double(line2(1:a(1)+3)); 
    a2 = str2double(line2(a(1)+4:a(2)+3));
    a3 = str2double(line2(a(2)+4:a(3)+3));
    a4 = str2double(line2(a(3)+4:a(4)+3));
    a5 = str2double(line2(a(4)+4:a(5)+3));
    
    % line2
    
    a6 = str2double(line2(1:b(1)+3));
    a7 = str2double(line2(b(1)+4:b(2)+3));
    a8 = str2double(line2(b(2)+4:b(3)+3));
    a9 = str2double(line2(b(3)+4:b(4)+3));
    a10 = str2double(line2(b(4)+4:b(5)+3));
    
    % line3
    
    a11 = str2double(line3(1:c(1)+3));
    a12 = str2double(line3(c(1)+4:c(2)+3));
    a13 = str2double(line3(c(2)+4:c(3)+3));
    a14 = str2double(line3(c(3)+4:c(4)+3));
    
    % required calculation
    
    Cp     = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,t,R,local_temp);
    H      = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,t,R,local_temp);
    S      = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,t,R,local_temp);
    m_w(i) = molecular_weight(chemical_element);
    
    % saving the plots in respected folders for elements
    
    e = cd(pwd); % identifies the current folder
    folder = chemical_element;
    mkdir(folder); % makes new directory/folder on target computer
    
    % plotting Cp vs temp
    
    figure(1)
    plot(t,Cp,\'linewidth\',2,\'color\',\'r\')
    % sprintf - formats data into string or character vector
    title(sprintf(\'Specific Heat vs Temperature plot of %s\', chemical_element))
    xlabel(\'Temperature\')
    ylabel(\'Specific Heat\')
    grid on
    
    % plotting H vs temp
    
    figure(2)
    plot(t,H,\'linewidth\',2,\'color\',\'g\')
    title(sprintf(\'Enthalpy vs Temperature plot of %s\', chemical_element))
    xlabel(\'Temperature\')
    ylabel(\'Enthalpy\')
    grid on
    
    % plotting S vs temp
    
    figure(3)
    plot(t,S,\'linewidth\',2,\'color\',\'b\')
    title(sprintf(\'Entropy vs Temperature plot of %s\', chemical_element))
    xlabel(\'Temperature\')
    ylabel(\'Entropy\')
    grid on
    
    % saving the graphs
    
    cd(folder)
    saveas(1,\'Cp_vs_temp.png\')
    saveas(2,\'H_vs_temp.png\')
    saveas(3,\'S_vs_temp.png\')
    cd(e)
    
end

Specific Heat function code -

function Cp = specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,t,R,local_temp)

if t > local_temp(3)
    
    Cp = R*(a1 + a2.*(t) + a3.*(t.^2) + a4.*(t.^3) + a5.*(t.^4));
    
else 
    
    Cp = R*(a8 + a9.*(t) + a10.*(t.^2) + a11.*(t.^3) + a12.*(t.^4));

end

Enthalpy function code -

function H = enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,a13,t,R,local_temp)

if t > local_temp(3)
    
    H = R*(a1.*t + [a2.*(t.^2)./2] + [a3.*(t.^3)./3] + [a4.*(t.^4)./4] + [a5.*(t.^5)./5] + a6);
    
else 
    
    H = R*(a8.*t + [a9.*(t.^2)./2] + [a10.*(t.^3)./3] + [a11.*(t.^4)./4] + [a12.*(t.^5)./5] + a13);  

end

Entropy function code -

function S = entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,t,R,local_temp)

if t > local_temp(3)
    
    S = R*(a1.*log(t) + a2.*(t) + a3.*(t.^2)./2 + a4.*(t.^3)./3 + a5.*(t.^4)./4 + a7);
    
else 
    
    S = R*(a8.*log(t) + a9.*(t) + a10.*(t.^2)./2 + a11.*(t.^3)./3 + a12.*(t.^4)./4 + a14);

end

Molecular Weight function code -

function m_w = molecular_weight(chemical_element)

    element       = [\'H\' \'C\' \'N\' \'O\' \'S\' \'A\'];
    atomic_weight = [1 12 14 16 32 40];
    
    m_w = 0;
    for i = 1:length(chemical_element)
        for j = 1:length(element)
            if strcmpi ( chemical_element(i), element(j) ) % compares strings (case insensitive)
                m_w = m_w + atomic_weight(j);
                position = j;
            end
        end
        
        n = str2num(chemical_element(i));
        
        if n>1
            m_w = m_w + atomic_weight(position)*(n-1);
        end
    end
    
    fprintf(\'Molecular Weigth of %s is : \',chemical_element);
    fprintf(\'%d\\n\', m_w);

end

Explanation -

Have to download the NASA data file, the \'thermo.dat\' file is saved in the same directory where the program code file will be saved. In the main code, necessary inputs such as universal gas constant value is assigned to variable R and the \'fopen\' command is used to open the \'thermo.dat\' file, header information is extracted using \'fgetl\' command. Global temperature range is obtained using \'fgetl\' and stored in the form of numbers using \'str2num\' command and the comments of next 3 lines are skipped using \'fgetl\' again.

For loop is used to read the data for all the 53 chemical elements one at a time, the required data from lineO to line3. LineO describes the name of species and local temperature range for the chemical elements. Line1 to line3 contains 14 coefficients for the chemical elements out of which first 7 are high temperature coefficients and the next 7 are low temperature coefficients which are designated from a1 to a14. The coefficients are extracted and converted from string data to double precision numbers using \'str2double\' command. Then the local temperature range is assigned to variable t for the plots.

The functions for \'specific_heat\', \'enthaply\', \'entropy\' and \'molecular_weight\' are defined using the equations stated. The defined functions are called into main program code to calculate the corresponding properties of the chemical element, then \'enthalpy\' entropy\' and \'specific heat\' are plotted against temperature. These plots are saved in the sub directories for each of the 53 chemical elements labelled with its respective chemical name.

Output -

Screenshot of the window where all the folders are saved

O2 (Specific Heat vs Temperature)

O2 (Enthalpy vs Temperature)

O2 (Entropy vs Temperature)

N2 (Specific Heat vs Temperature)

N2 (Enthalpy vs Temperature)

N2 (Entropy vs Temperature)

CO2 (Specific Heat vs Temperature)

CO2 (Enthalpy vs Temperature)

CO2 (Entropy vs Temperature)

Molecular weight displayed in command window -

Errors -

\'Ar\' was given as input but there where only six values in atomic weight array, So had to change to \'A\'

instead of \';\' had a typo with \':\'