Parsing NASA Thermodynamics Data

AIM: Write a program in MATLAB to parse the NASA thermodynamic data file and then calculate thermodynamic properties of various gas species.

OBJECTIVES:

1. Write the function that exacts the 14 coefficients. Calculate the Entropy, Enthalpy and Specific Heat with help of formula
2. Calculate the molecular weight of each species and display it on command window.
3. Plot the graph for Entropy, Enthalpy and Specific Heat for given local temperature range.
4. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.

FILE PARSING

File parsing is referred to the read the information from the any type of data file and use this information for the calculate the other things.

From THERMO.dat file we exact the unique 14 coefficients for each species in which the first 7 coefficients are high-temperature coefficients and the second 7 coefficients are low-temperature coefficients. Also use the local temperature ranges which is unique for each individual species.

With the help of this data, we calculate Entropy(S), Enthalpy(H) and Specific Heat (C_p) with help of following formula:

Where

C_p = Specific Heat

H = Enthalpy

S = Entropy

T = Temperature

R = Universal gas constant

a = Coefficients

Here The first 7 coefficients are high-temperature coefficients and the second 7 coefficients are low-temperature coefficients.

STEP-WISE EXPLANATION OF CODE:

1. First, we open the THERMO.dat file with help of the fopen command.
2. We read the main header of the file with help of fgetl command.
3. We split the line2 form the space with help of strsplit command to obtain the temperature.
4. We assign the global low, medium and high temperature with str2num command which convert the string to number.
5. Then we close the file f1.
6. Now we use another fileread command to read the file.
7. We use the regexp command for blocks variable. Syntax of the line is as follows:
   1. If the line begins with whitespace, go on to the next line.
   2. Skip the whitespace that is immediately after that first field.
   3. Scan any number of characters that are not newline characters, and put those characters into a field named header -- in other words header gets the entire rest of the line.
   4. Throw away the newline after the header.
   5. Any number of characters that are not newline characters, along with the following newline character. Repeat this 3 times (that is, read three lines). Put all of those lines in a field named numbers.So, we get like this

8. Then we will try to get only species name in another variable.
9. We take the input of species from command window.
10. We use the upper command for case insensitive.
11. Then we search the input species name in block and take the relevant data for this species.
12. Then we will try to get the coefficient number in relevant data variable.
13. We reshape the matrix which contain the coefficients.
14. Now we exact the local temperature values in xyz matrix.
15. We assign the local low, medium and high temperature value.
16. Then we assign the 14 coefficient values.
17. We calculate the molecular weight with help of mol_wt subfunction.
18. We calculate the Entropy(S), Enthalpy(H) and Specific Heat (C_p) with help of variable subfunction.
19. We change the directory and make new folder with particular species name to save the graph.
20. We plot the graph for Entropy(S), Enthalpy(H) and Specific Heat (C_p) vs local temperature range.
21. We save the graph in species folder with help of saveas command.

SUB FUNCTION FOR MOLECULAR WEIGHT

function [mol_wt]= mwt1(target_species)

elements={'H','C','N','O','AR'};
element_wt = [1,12,14,16,40];
species=upper(target_species);
mol_wt =0;

if (strcmp(species,elements(5)));
    mol_wt=40;
else
  for i=1:length(species)
    for j=1:length(elements)
        if(strcmp(species(i),elements(j)))
            mol_wt=mol_wt+element_wt(j);
            v=j;
        end
    end
    n=str2double(species(i));
    if(n>1)
       mol_wt=mol_wt+element_wt(v)*(n-1) ;
    end
end
end
 fprintf('molecular weight of %s:-',target_species)
 fprintf('%d',mol_wt)
 fprintf(' ')
 fprintf('amu')
end

SUB FUNCTION FOR CALCULATE C_p, H AND S

function [Cp H S] =variable(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_midtemp,T)
R=8.314; %unniversal gas constant

for i = 1:length(T)
    if (i>local_midtemp)
Cp(i) = R*(a1 + (a2.*T(i)) + (a3.*(T(i).^2)) + (a4.*(T(i).^3)) + (a5.*(T(i).^4)));

H(i) = R.*T(i).*(a1 + a2.*(T(i)./2) + a3.*(T(i).^2)./3 + a4.*(T(i).^3)./4+ a5.*(T(i).^4)./5 + a6./T(i));

S(i) = R* (a1.*log(T(i)) + a2.*T(i) + a3.*(T(i).^2)./2 + a4.*(T(i)^3)/3 + a5.*(T(i).^4)/4 + a7 );

    else
        
        Cp(i) = R*(a8 + (a9.*T(i)) + (a10.*(T(i).^2)) + (a11.*(T(i).^3)) + (a12.*(T(i).^4)));

        H(i) = R.*T(i).*(a8 + a9.*(T(i)/2) + a10.*(T(i).^2)./3 + a11.*(T(i).^3)./4+ a12.*(T(i).^4)./5 + a13./T(i));

        S(i) = R* (a8.*log(T(i)) + a9.*T(i) + a10.*(T(i).^2)./2 + a11.*(T(i)^3)./3 + a12.*(T(i).^4)./4 + a14 );
    end
end

MAIN PROGRAM

close all
clear all
clc

f1 = fopen('THERMO.dat','r'); %open the THERMO file
L1 = fgetl(f1);  %reading header 1st line
L2 =fgetl(f1);  %2nd line
temp =strsplit(L2,' '); %converting into cell array

global_lowtemp = str2num(temp{2}); %convert string in number to obtain temperature
global_midtemp = str2num(temp{3});
global_hightemp = str2num(temp{4});

L3 = fgetl(f1); %3rd line
L4 = fgetl(f1); %4th line
L5 = fgetl(f1); %5th line
fclose(f1)   %close the file

filename = 'THERMO.dat';
S = fileread(filename); %read the data file
%extract the species,header and coefficient in different token in blocks
blocks = regexp(S, '^(?<species>\S+)\s*(?<header>[^\n]*)\n(?<numbers>([^\n]+\n){3})', 'names', 'lineanchors');
species_names = {blocks.species}; %extract species name from block

input_species = input('enter species name =')
target_species=upper(input_species);

% match the target species form the species_name
relevant_block = blocks(ismember(species_names, target_species)); 
%store the coefficients in relevant data 
relevant_data = sscanf(relevant_block.numbers, '%f', [1 17]);
relevant_data(18) = relevant_data(17); %17th place is empty so 17th place=18th place
relevant_data(17) = nan; %place 17 is empty
relevant_data = reshape(relevant_data, [6, 3]).'; %reshape the data in matrix form of 3 row and 6 column
relevant_data(:,end) = [];    %just numbering

%extract the header detail of traget species
relevant_header = relevant_block.header;
descr = strtrim(relevant_header(1:10));
shell = strtrim(relevant_header(11:27));

%xyz store the local temperature values
xyz = sscanf(relevant_header(28:end-1), '%g', [1 3]);


local_lowtemp =xyz(1); %obtain local low temperature
local_midtemp =xyz(3); %obtain local medium temperature
local_hightemp =xyz(2); %obtain local high temperature

%High temperature coefficients 
a1= relevant_data(1,1);
a2= relevant_data(1,2);
a3= relevant_data(1,3);
a4= relevant_data(1,4);
a5= relevant_data(1,5);
a6= relevant_data(2,1);
a7= relevant_data(2,2);

%Low temperature coefficients 
a8= relevant_data(2,3);
a9= relevant_data(2,4);
a10= relevant_data(2,5);
a11= relevant_data(3,1);
a12= relevant_data(3,2);
a13= relevant_data(3,3);
a14= relevant_data(3,4);

%molecular weight of the species

[mol_wt]= mwt1(target_species);


T=linspace(local_lowtemp,local_hightemp,1000); %for plotting the graph the temperature 
%find the values of cp,H and S
[Cp H S] =variable(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_midtemp,T);   


cd('C:\Users\Krunal Kabariya\Desktop\NASA1')
%create the new folder
mkdir(target_species)

%plot the graph
figure(1)
plot(T,S,'linewidth',3,'color','r')
xlabel('Tempertaure')
ylabel('Entropy')
title(target_species)

figure(2)
plot(T,H,'linewidth',3,'color','b')
xlabel('Tempertaure')
ylabel('Enthalpy')
title(target_species)

figure(3)
plot(T,Cp,'linewidth',3,'color','k')
xlabel('Tempertaure')
ylabel('specific Heat')
title(target_species)

cd(target_species)

saveas(1,'Entropy.png')
saveas(2,'Enthalpy.png')
saveas(3,'Specific Heat.png')

cd('F:\Skill lync\MAT Lab\NASA')

INPUT IN COMMAND WINDOW

OUTPUT IN COMMAND WINDOW

GRAPH FOR O₂

GRAPH FOR N₂

GRAPH FOR CO₂

MOLECULAR WEIGHT OF EACH SPECIES

SAVING GRAPH OF EACH SPECIES IN SEPARATE FOLDER

ALTERNATE PROGRAM

I also write alternate program for parsing NASA thermodynamic data for plot the graph of each species. In this program, I used for loop which execute 53 times to generate the plot for each species. For this we have to get the line from file than split this line and convert to number to obtain coefficients.We have to do some common changes in subfunction of molecular weight because of some notation changes. Here is the program for this method:

close all
clear all
clc

f1 = fopen('THERMO.dat','r'); %open the THERMO file
L1 = fgetl(f1);  %reading header 1st line
L2 =fgetl(f1);  %2nd line
temp =strsplit(L2,' '); %converting into cell array

global_lowtemp = str2num(temp{2}); %convert string in number to obtain temperature
global_midtemp = str2num(temp{3});
global_hightemp = str2num(temp{4});

L3 = fgetl(f1); %3rd line
L4 = fgetl(f1); %4th line
L5 = fgetl(f1); %5th line

for i=1:53 %for 53 species
    
    species_line=fgetl(f1); %new line which contain species name and temperature range 
    a=strsplit(species_line,' '); %convert in to cell array
    species_name=(a{1}); %obtain species name
   
    
    local_lowtemp =str2num(a{end-3}); %obtain local low temperature
    local_midtemp =str2num(a{end-1}); %obtain local medium temperature
    local_hightemp =str2num(a{end-2}); %obtain local high temperature
    
    coeff_L1=fgetl(f1); %new line which contain coefficients
    b=strfind(coeff_L1,'E'); %find the position of E in the line
    
    %finding first 5 coefficient
    a1=coeff_L1(1:b(1)+3);
    a1=str2num(a1);
    a2=coeff_L1(b(1)+4:b(2)+3);
    a2=str2num(a2);
    a3=coeff_L1(b(2)+4:b(3)+3);
    a3=str2num(a3);
    a4=coeff_L1(b(3)+4:b(4)+3);
    a4=str2num(a4);
    a5=coeff_L1(b(4)+4:b(5)+3);
    a5=str2num(a5);
    
    coeff_L2=fgetl(f1); %new line which contain coefficients
    c=strfind(coeff_L2,'E'); %find the position of E in the line
    
    %finding second 5 coefficient
    a6=coeff_L2(1:c(1)+3);
    a6=str2num(a6);
    a7=coeff_L2(c(1)+4:c(2)+3);
    a7=str2num(a7);
    a8=coeff_L2(c(2)+4:c(3)+3);
    a8=str2num(a8);
    a9=coeff_L2(c(3)+4:c(4)+3);
    a9=str2num(a9);
    a10=coeff_L2(c(4)+4:c(5)+3);
    a10=str2num(a10);
    
    
    
    coeff_L3=fgetl(f1); %new line which contain coefficients
    d=strfind(coeff_L3,'E'); %find the position of E in the line
    
    %finding second 5 coefficient
    a11=coeff_L3(1:d(1)+3);
    a11=str2num(a11);
    a12=coeff_L3(d(1)+4:d(2)+3);
    a12=str2num(a12);
    a13=coeff_L3(d(2)+4:d(3)+3);
    a13=str2num(a13);
    a14=coeff_L3(d(3)+4:d(4)+3);
    a14=str2num(a14);
    
    T=linspace(local_lowtemp,local_hightemp,1000); %for plotting the graph the temperature 
    %find the values of cp,H and S
    [Cp H S] =variable(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_midtemp,T);    
   
    %molecular weight of the species
        [mol_wt]= mwt(species_name);
    
   
     cd('C:\Users\Krunal Kabariya\Desktop\NASA');   
    %create the new folder
    mkdir(species_name);
    
    %plot the graph
    figure(1)
    plot(T,S,'linewidth',3,'color','r')
    xlabel('Tempertaure')
    ylabel('Entropy')
    title(species_name)
    
    figure(2)
    plot(T,H,'linewidth',3,'color','b')
    xlabel('Tempertaure')
    ylabel('Enthalpy')
    title(species_name)
    
    figure(3)
    plot(T,Cp,'linewidth',3,'color','k')
    xlabel('Tempertaure')
    ylabel('specific Heat')
    title(species_name)
    
    cd(species_name)
    saveas(1,'Entropy.png')
    saveas(2,'Enthalpy.png')
    saveas(3,'Specific Heat.png')
    
    cd('F:\Skill lync\MAT Lab\NASA')
    

   
end

REFERENCE

https://in.mathworks.com/help/matlab/ref/regexp.html

https://in.mathworks.com/help/matlab/ref/reshape.html