Parsing NASA thermodynamic data using MATLAB

Introduction:

This project deals with extracting relevant information from a data file and using that information to solve the following challenges.

1. Extracting 14 coefficients from the file and using them in order to calculate the gas properties of various elements using the local temperature ranges.
   
   
   
2. Using the gas properties namely specific heat, enthalpy and entropy and plotting them against temperature and also automatically saving the plots to a directory.
   
   
   
3. Calculating the molecular weight of all the compounds in the file and displaying it in the command window.
   
   
   

All the challenges above were solved using different function definitions and ultimately calling these functions in the main program. The first step was to extract the coefficients from the data file. This was achieved after carefully observing the '.dat' file and deducing that the coefficients were organized in a systematic order ie 5 coefficients each in the first 2 lines and the last 4 coefficients in the 3rd line. After this observation was made, simple commands like fgetl were used along with the concept from the tutorial.

Function for extracting coefficients:

function [F]=coefficients(f)

%Reading the dat file line by line
f1=fgetl(f);

%looking for E in the coefficient line by line and using it to evaluate the
%coefficient
A=strfind(f1,'E');

%line1
a1=f1(1:A(1)+3);
A1=str2num(a1);
a2=f1(A(1)+5:A(2)+3);
A2=str2num(a2);
a3=f1(A(2)+5:A(3)+3);
A3=str2num(a3);
a4=f1(A(3)+5:A(4)+3);;
A4=str2num(a4);
a5=f1(A(4)+5:A(5)+3);
A5=str2num(a5);


%line2
f2=fgetl(f);
B=strfind(f2,'E');

b1=f2(1:B(1)+3);
B1=str2num(b1);
b2=f2(B(1)+5:B(2)+3);
B2=str2num(b2);
b3=f2(B(2)+5:B(3)+3);
B3=str2num(b3);
b4=f2(B(3)+5:B(4)+3);
B4=str2num(b4);
b5=f2(B(4)+5:B(5)+3);
B5=str2num(b5);


%line3
f3=fgetl(f);
C=strfind(f3,'E');

c1=f3(1:C(1)+3);
C1=str2num(c1);
c2=f3(C(1)+5:C(2)+3);
C2=str2num(c2);
c3=f3(C(2)+5:C(3)+3);
C3=str2num(c3);
c4=f3(C(3)+5:C(4)+3);
C4=str2num(c4);

F=[A1 A2 A3 A4 A5 B1 B2 B3 B4 B5 C1 C2 C3 C4];

end

After extracting the coefficients successfully, the next step was to use them along with the local temperature ranges in order to calculate the following gas properties

1) Specific heat(Cp)


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

2) Enthalpy(H)


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

3) Entropy(S)


`S=(a1*log10(T)+a2*T+(a3*T^2)/2+(a4*T^3)/3+(a5*T^4)/4+a7)*R`

The main challenge here was to make sure that we are using the right coefficients for the right temperature ranges. Therefore first seven coefficients were used for the low range temperature ie lower bound to mid and the next seven were used for high temperature ie mid to high using the function below.


Function to calculate the gas properties:

function [Cp,S,H]=gas_prop(F,T_mid,T,R)

%Function to return the Specific heat(Cp),Entropy(S),and Enthalpy(H) using
%the formulas mentioned for different temperature ranges

for i=1:length(T)
    
    if T(i)>=T_mid
        
        Cp=(F(1)+F(2).*T+F(3).*T.^2+F(4).*T.^3+F(5).*T.^4).*R;
        S=(F(1).*log10(T)+F(2).*T+(F(3).*T.^2)/2+(F(4).*T.^3)/3+(F(5).*T.^4)/4+F(7)).*R;
        H= (F(1).*T + F(2).*(T.^2)/2 +F(3).*(T.^3)/3 + F(4).*(T.^4)/4 + F(5).*(T.^5)/5 + F(6).*(T.^2)).*R.*T;
        
    else
        Cp=(F(8)+F(9).*T+F(10).*T.^2+F(11).*T.^3+F(12).*T.^4).*R;
        S=(F(8).*log10(T)+F(9).*T+(F(10).*T.^2)/2+(F(11).*T.^3)/3+(F(12).*T.^4)/4+F(14)).*R;
         H= (F(8).*T + F(9).*(T.^2)/2 +F(10).*(T.^3)/3 + F(11).*(T.^4)/4 + F(12).*(T.^5)/5 + F(13).*(T.^2)).*R.*T;
        
    end
    
end

end

After calculating the gas properties the next step was to plot them against the local temperature ranges provided in the data file and also saving them automatically in a directory. This was achieved easily as all the data was already extracted using the two functions above


Function to plot the gas properties against temperature:

function [plots]=plotting(T,folder,Cp,H,S,elementname)

     %Making a directory where all the plots will be saved
     c = pwd;
     cd(folder);
     data = sprintf('%s',elementname);
     mkdir(data);
     
     %Plotting CpvsT
     figure(1)
     plot(T,Cp,'linewidth',3,'color','c');
     title(elementname)
     xlabel('Temperature [K]');
     ylabel('Specific heat Cp [J/K]');
     legend('specific heat');

     %plotting HvsT
     figure(2)
     plot(T,H,'linewidth',3,'color','k');
     title(elementname)
     xlabel('Temperature [K]')
     ylabel('Enthalpy H [J/Kg]')
     legend('Enthalpy');
     
     %plotting SvsT
     figure(3)
     plot(T,S,'linewidth',3,'color','m');
     title(elementname)
     xlabel('Temperature [K]')
     ylabel('Entropy S [J/K]')
     legend('Entropy');
     
     %Saving the plots with relevant names
     cd(data)
  saveas(1,'Specific_heat.png')
  saveas(2,'Enthalpy.png')
  saveas(3,'Entropy.png')
  cd(c)
   
end

After saving all the plots to a directory the next challenge was to calculate the molecular weight of all the compounds and display them in the command window. This was done after carefully observing all the compound and noticing that all of them were made up of 6 basic elements namely H, C, N, O, Ar and S. The molecular weight of these 6 fundamental elements became the basis to calculate the weight of all the other compounds. This was achieved using the code below


Function to calculate molecular weight:

function [M_W] = molecularweight(elementname)

R_elem= {'H','C','N','O','AR','S'};
A_W = [1 12 14 16 40 32];

M_W = 0;
for i = 1:length(elementname)
    for j = 1:length(R_elem)
        if strcmp(elementname(i),R_elem(j))
            M_W = M_W+A_W(j);
            c = j;
        end
    end
    n = str2num(elementname(i));
    if n>1
        M_W = M_W+A_W(c)*(n-1);
    end
end
end

Once all the function definitions were ready, it was time to call them in the main function. Smaller tasks like calculating the total number of compunds and getting their individual names were also done in the main program.


Code for the main program:

%clearing the command window and workspace
clear all
close all 
clc


%Reading the NASA file provided
f = fopen('THERMO.dat','r')

%universal gas constant R, 
R=8.314;

%fgetl is used to read the document line by line
fgetl(f);

%Calculating the temperature ranges given in the dat file 
temp=fgetl(f);

%using strsplit and ' spacebar ' as reference to differentiate between
%temperatures
T_A=strsplit(temp,' ');


T_lb=str2num(T_A{2});
T_mid=str2num(T_A{3});
T_ub=str2num(T_A{4});

%100 different values of temperature to get the plots
T=linspace(T_lb,T_ub,200);

%creating the folder where all the elements along with data will be saved
folder='NASA_elements';
%Making a new folder in the directory
mkdir(folder)

%Skipping the first 3 lines as they are comments
    
    skip1=fgetl(f);
    skip2=fgetl(f);
    skip3=fgetl(f);


%calculating total no. of elements in the dat file
%Starting with line=0
line=0;

for i=1:Inf    
    L1=fgetl(f);
    
    if L1=="END"
        break
        
    end
    
    line=line+1;
    
end

Totalelements=(line)/4

%command to come back on top of the dat file
frewind(f)

%skipping the next five lines
for i=1:5
   skip=fgetl(f);
      
end

%Calling all the required functions for every element

for j=1:Totalelements
 
%find element name    
elementname=name_func(f)
     
%extracting coefficients
F=coefficients(f);

%Function to calculate gas properties
[Cp, S, H]=gas_prop(F,T_mid,T,R);


%Molecular weight function
MW=molecularweight(elementname)

%Plotting function
plotting(T,folder,Cp,H,S,elementname)

end

fclose(f)

Results:

1. Screenshot of the directory where all the plots were saved automatically
   
   
   
   
   
   
2. Plots for O2 gas properties:
   
   
   a) Specific heat (Cp)
   
   
   
   
   
   
   
   b) Enthalpy (H)
   
   
   
   
   
   
   c) Entropy (S)
   
   
   
   
   
   
   
   
3. Plots for N2 gas properties:
   
   
   a) Specific heat (Cp)
   
   
   
   
   
   b) Enthalpy (H)
   
   
   
   
   
   c) Entropy (S)
   
   
   
   
   
   
4. Plots for CO2 gas properties:
   
   
   a) Specific heat (Cp)
   
   
   
   
   
   b) Enthalpy (H)
   
   
   
   
   
   
   c) Entropy (S)
   
   
   
   
   
   
   
5. Screenshot of command window with the molecular weights:
   
   
   
   
   
   
   

Errors:

• The initial value for f was -4 (negative) which was fixed by changing the file name
  
  
• the molecular function didn't include sulfur (S) initially but was later accommodated after checking the weight results
  
  
• Naming error in one of the functions was fixed after the prompt.
  
  
  
  
  

References:

1. coming back on top of the file- https://in.mathworks.com/help/matlab/ref/frewind.html
   
   
   
2. saving plots to a directory- https://in.mathworks.com/matlabcentral/answers/46200-how-to-save-a-figure-file-into-a-subfolder
   
   
   
3. molecular weight function help- https://in.mathworks.com/matlabcentral/answers/106320-make-a-loop-to-calculate-molecular-weight