Parsing NASA thermodynamic data using matlab

The formulaes for calculating enthalpy, entropy, and specific heat are as follows:-

`(Cp)/R = a1 + a2 T + a3 T^2 + a4 T^3 + a5 T^4`

`(h)/(RT) = a1 + a2 T/2 + a3 T^2/3 + a4 T^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`

where, R =universal gas constant

T=temperature

The main program is given below:-

lear all
close all
clc

%Inputs
R= 8.314;

%Reading the data file
file = fopen('THERMO.dat','r');

%Reading the first line of the file 
l1 =fgetl(file);

%Getting the global minimum, maximum and mid temperatures
l2=fgetl(file);
temperatures = strsplit(l2);
global_min = str2num(temperatures{2});
global_mid = str2num(temperatures{3});
global_max = str2num(temperatures{4});

%Defining the temperature range 
temp = linspace (global_min,global_max,500);

%Ignoring the comments in the file 
for i=1:3
fgetl(file);
end

%Creating a loop to go through all the 53 species of data in the file 
for i=1:53
    %Reading the first line which includes the information for species
    l=fgetl(file);
    
    %Reading the name of the species
    element_name=strtok(l);
    
    %Reading the three lines which contains 7 high temp coefficients & 7 low temp 
    %coefficients
    l1=fgetl(file);
    l2=fgetl(file);
    l3=fgetl(file);
    
    %Finding the position of E so that we can extract the coefficients
    a=strfind(l1,'E');
    b=strfind(l2,'E');
    c=strfind(l3,'E');

    %extracting high temp coeff
    a1=str2num(l1(1:a(1)+3));
    a2=str2num(l1(a(1)+3:a(2)+3));
    a3=str2num(l1(a(2)+4:a(3)+3));
    a4=str2num(l1(a(3)+4:a(4)+3));
    a5=str2num(l1(a(4)+4:a(5)+3));
    a6=str2num(l2(1:b(1)+3));
    a7=str2num(l2(b(1)+4:b(2)+3));

    %Extracting low temp coeff
    b1=str2num(l2(b(2)+4:b(3)+3));
    b2=str2num(l2(b(3)+4:b(4)+3));
    b3=str2num(l2(b(4)+4:b(5)+3));
    b4=str2num(l3(1:c(1)+3));
    b5=str2num(l3(c(1)+4:c(2)+3));
    b6=str2num(l3(c(2)+4:c(3)+3));
    b7=str2num(l3(c(3)+4:c(4)+3));

    %Calculating specific heat
    c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid);
    
    %Calculalting enthalpy
    H = enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid);
    
    %calculating entropy
    S = entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid);
    
    %Plotting graphs
    plotting(c_p,H,S,temp,element_name)
    
    %Calculating molecular weights
    mol_wt(i) = molecular_weight(element_name)
    
    %Function to write molecular weight to a file 
    weight = mol_wt(i);
    writing_to_file(element_name, weight)
    
end

fclose (file)

Explaination of the code:-

• The given data is then formatted as a thermo.dat file.
• the global maximum, global mid and global minimum values which is read bby using fgetl command and split into the three values by using the strsplit command and converted into numbers.
• Then a for loop is created to loop through the 53 species data given in the file and clculate the required values.
• The low temperatures values are stored in variables a1,a2,a2,a4,a5,a6 and a7 and the high temperature values are stored invariables b1,b2,b3,b4,b5,b6 and b7.
• The next line is the reference to a function called specific heat which calculates and returns the value of specific heat using the high and low temperature values as inputs.
• The next line is the reference to a function called enthalpy which calculates and returns the value of enthalpy using the high and low temperature values as inputs.
• The next line is the reference to a function called entropy which calculates and returns the value of entropy using the high and low temperature values as inputs.
• The next line is the reference to a function that plots graphs using the given inputs.
• The next line is the reference to a function that calculates the molecular weights of the different elements.
• The next line is the reference to a function that writes the molecular weights of diferent elements to a file.

The NASA thermdynamic data file is closed using the fclose command.

1. Function to calculate the Specific heat

%Function to caculate specific heat
function c_p = specific_heat(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid)
    if temp<=global_mid
        c_p = R*(a1 + a2.*temp + a3.*(temp.^2) + a4.*(temp.^3) + a5.*(temp.^4))
    else
        c_p = R*(b1 + b2.*temp+b3.*(temp.^2) + b4.*(temp.^3) + b5.*(temp.^4))
end

2. Function to calculate the Enthalpy

%Function to calculate enthalpy
function H = enthalpy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid)
    if temp <=global_mid
        H = R.*temp.*(a1 + a2.*(temp./2) + a3.*((temp.^2)./3) + a4.*((temp.^3)./4) + a5.*((temp.^4)./5) + a6./temp)
    else
        H = R.*temp.*(b1 + b2.*(temp./2) + b3.*((temp.^2)./3) + b4.*((temp.^3)./4) + b5.*((temp.^4)./5) + b6./temp)
end

3. Function to clculate the Entropy

%Function to calculate entropy
function S =entropy(a1,a2,a3,a4,a5,a6,a7,b1,b2,b3,b4,b5,b6,b7,R,temp,global_min,global_max,global_mid)
    if temp<=global_mid
        S = R*(a1.*log(temp) + a2.*temp + a3.*(temp.^2)./2 + a4.*(temp.^3)./3 + a5.*(temp.^4)./4 + a7)
    else
        S = R*(b1.*log(temp) + b2.*temp + b3.*(temp.^2)./2 + b4.*(temp.^3)./3 + b5.*(temp.^4)./4 + b7)

end

4. Function to plot the graphs

%Function to plot graphs for all elements
function plotting(c_p,H,S,temp,element_name)
cd(pwd)
dir_name=element_name;
mkdir(dir_name);

%Plotting temperature vs specific heat
figure(1)
plot(temp,c_p,'lineWidth',3)
title(sprintf('Specific Heat vs Temperature rn for %s:',element_name))
xlabel('Temperature')
ylabel('Specific Heat')

%Plotting temperature vs specific heat
figure(2)
plot(temp,H,'linewidth',3)
title(sprintf('Enthalpy vs Temperature rn for %s:',element_name))
xlabel('Temperature')
ylabel('Enthalpy')

figure(3)
plot(temp,S,'linewidth',3)
title(sprintf('Entropy vs Temperature rn for %s:',element_name))
xlabel('Temperature')
ylabel('Entropy')

cd(dir_name)
saveas(1,'c_p_vs_temp.png')
saveas(2,'H_vs_temp.png')
saveas(3,'S_vs_temp.png')

cd ..

end

Explaination while plotting

• Using the mkdir command a new folder is created in the current folder pwd with the same name as the element name.
• The first figure(1) represents the plots the temperature vs specific heat values.
• The first figure(2) represents the plots the temperature vs enthalpy values.
• The first figure(3) represents the plots the temperature vs entropy values.
• Using the saveas command the three figures are saved in the new folder which is created by the element name.
• Using the cd.. command the folder is changed back to the original folder.
• The function runs for 53 times through the for loop and it creates seperate folders for each species and saves all the plots in the folders.

5. Function to calculate the molecular weight of all elements

%Function for clculating molecular weight
function [w] = molecular_weight(element_name)

s = {'H' 'C' 'N' 'O' 'AR'}; 
%atomic weight of each elements
atomic_wt = [1 12 14 16 40];                                                      
w = 0;                                                                       

for i = 1:length(element_name)  
 for j = 1:length(s)
    if strcmp(element_name(i), s(j))                           
    w= w + atomic_wt(j); 
    %storing the position of matched element
    position = j;                                                         
  end 
end 
   
n = str2num(element_name(i));
if n > 1
                                                      
w = w + atomic_wt(position)*(n-1);                                        
end 
end
end

• A string array is created which contains the name of the atoms.
• A variable atomic_wt is created which contains the atomic weight of respective atoms in the string array.
• Then used strcmp command to compare the name of the element with the string array.
• If there is repeatition of any element, str2num is used to find that and (n-1) is multiplied to atomic_wt of element to calculate molecular weight.

6. Function to write the molecular weight data of all elements to a file

%Writing molecular weight of elements to a file
    function writing_to_file(element_name,weight)
    dir_name = element_name
    mol_file = fopen('mol_wt.docx' , 'a+')
    fprintf (mol_file, 'Molecular weight of %s is %f n', dir_name,weight)
    fclose (mol_file)
    end

• The above function takes the element name and the molecular weights as the data inputs and writes the molecular weight along with the name of the element in a file.

• Firstly a new file is created with name 'mol_wt.docx' and it is opened in append mode denoted by 'a+' so that we can write the molecular weight of all elements in a single file.

• Using the fprintf command the name of the element along with their molecular weight is written to the file.

• The file is then closed.

Screenshot of the command window:-

Screenshot of the window where all the folders are saved.

Plots for O2:-

plots for N2:-

plots for CO2:-

Referances: Google , Skill-lync, skill-lync helpdesk.