Parsing NASA thermodynamic data

Objective :-

• Reading the 'NASA thermodynamic data file' and write the function that extracts the 14 coefficient and calculate the enthalphy , entrophy and specific heats for all the species
• Calculate the molecular weight of each species and display it in the command window.
• plot the Cp, Enthalphy, and , Entrophy for the local temperature range.
• Save the plots as image with appropriate names and dump them in separate folders for each species.

Theory :-

Where ,

1. R is the universal gas constant (R = 8.314 J/mol . K).
2. T is the temperature (Local Temperature).
3. Cp is the specific heat (J/Kg-K)
4. H is the enthalphy (J/Kg)
5. S is the entrophy (J/Kg)

The atomic weight for the elements which is used to calculated Molecular Weight (MW) is given below

• H = 1.00784
• C =12.0107
• Ar= 39.948
• N = 14.0067
• O = 15.9999

MATLAB CODE :

Function to calculate the Specific Heat (Cp) :-

% Function to find Specific Heat
function Cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Lower_temp,Higher_temp,mt)
if Higher_temp>=mt
    Cp =(a1+a2.*(Higher_temp)+a3.*(Higher_temp.^2)+a4.*(Higher_temp.^3)+a5.*(Higher_temp.^4)+a6.*(Higher_temp.^5)+a7.*(Higher_temp.^6))*R;
else 
    Cp =(a8+a9.*(Lower_temp)+a10.*(Lower_temp.^2)+a11.*(Lower_temp.^3)+a12.*(Lower_temp.^4)+a13.*(Lower_temp.^5)+a14.*(Lower_temp.^6))*R;
end

Explanation :-

• Function command is determined for Cp with inputs of low & high temperature coefficient Universal gas constant (R) , Temperature (T) and local mid temperature.
• Then if else command is used to check the condition wheater Higher temp is less than or equal to (<=) local mid temp
• Cp is then calculated by writing the formula as code for both low and high temperature coefficients.

Fuction to calculate the Enthalphy (H) :-

function H=enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Lower_temp,Higher_temp,mt)
if Higher_temp>=mt
    H= (a1+a2.*(Higher_temp./2)+a3.*((Higher_temp).^2)./3+a4.*((Higher_temp).^3)./4+a5.*((Higher_temp).^4)./5+a6.*((Higher_temp).^5)./6+a7.*((Higher_temp).^6)./7)*R;
else 
    H=(a8+a9.*(Lower_temp./2)+a10.*((Lower_temp).^2)./3+a11.*((Lower_temp).^3)./4+a12.*((Lower_temp).^4)./5+a13.*((Lower_temp).^5)./6+a14.*((Lower_temp).^6)./7)*R;
end 

Explanation :-

• Function command is determined for inputs of low & high temperature coefficient Universal gas constant (R) , Temperature (T) and local mid temperature.
• Then if else command is used to check the condition wheater Higher temp is less than or equal to (<=) local mid temp
• Low temp coefficient a6 & High temperature b6 is calculated additionally.
• H is then calculated by writing the formula as code for both low and high temperature coefficients.

Fuction to calculate the Entrophy (s) :-

% Funbction to find Entropy
function S=entrophy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Lower_temp,Higher_temp,mt)
if Higher_temp>=mt
    S=(a1.*log(Higher_temp)+a2.*Higher_temp+a3.*(Higher_temp.^2)./2+a4.*(Higher_temp.^3)./3+a5.*(Higher_temp.^4)./4+a7)*R;
else 
    S=(a8.*log(Lower_temp)+a9.*Lower_temp+a10'*(Lower_temp.^2)./2+a11.*(Lower_temp.^3)./3+a12.*(Lower_temp.^4)./4+a14)*R;
end

Explanation :-

• Function command is determined forSinputs of low & high temperature coefficient Universal gas constant (R) , Temperature (T) and local mid temperature.
• Then if else command is used to check the condition wheater Higher temp is less than or equal to (<=) local mid temp
• Low temp coefficient a6 & High temperature b6 is calculated additionally.
• S is then calculated by writing the formula as code for both low and high temperature coefficients.

Main programme :-

% program for Parsing NASA thermodynamic data
clear all
close all
clc

% universal gas constant
R = 8.314;

% Read the data file 
file = fopen('THERMO.dat','r');
header = fgetl(file);                            

%reading the global temperature 
global_temp = fgetl(file);
a = strsplit(global_temp,' ');
% global low temperature
global_low_temp =str2double(a{2});
% global middle temperature
global_mid_temp =str2double(a{3});
% global high temperature
global_high_temp =str2double(a{4});

% skipping the next three line 
fgetl(file);
fgetl(file);
fgetl(file);

%using for loop command to consist 53

for i = 1:53
    line= fgetl(file);% reading the line
    A=strsplit(line,'');% seperating the variable
    compounds=A(1); % naming the first variables as 'compounds'
    Q = strfind(line,'G');
    local_low_temp = line(Q+3:Q+10);
    local_high_temp = line(Q+12:Q+20);
    local_mid_temp = line(Q+22:Q+30);
    Lower_temp = str2double(local_low_temp);
    Higher_temp = str2double(local_high_temp);
    mt = str2double(local_mid_temp);
    
    % line 1-Reading the coefficient data
    
    cfline1= fgetl(file);
    X = strfind(cfline1,'E');
    a1 = cfline1(1:X(1)+3);
    a1 = str2double(a1); % storing the first coefficient
    a2 = cfline1(X(1)+4:X(2)+3);
    a2 = str2double(a2); % storing the second coefficient
    a3 = cfline1(X(2)+4:X(3)+3);
    a3 = str2double(a3); % storing the third coefficient
    a4 = cfline1(X(3)+4:X(4)+3);
    a4 = str2double(a4); % storing the fourth coefficient
    a5 = cfline1(X(4)+4:X(5)+3);
    a5 = str2double(a5); % storing the fifth coefficient
    
    %line 2-Reading the coefficient data
    
    cfline2= fgetl(file);
    G = strfind(cfline2,'E');
    a6 = cfline2(1:G(1)+3);
    a6 = str2double(a6); % storing the sixth coefficient
    a7 = cfline2(G(1)+4:G(2)+3);
    a7 = str2double(a7);% storing the seventh coefficient
    a8 = cfline2(G(2)+4:G(3)+3);
    a8 = str2double(a8);% Storing the eight coefficient
    a9 = cfline2(G(3)+4:G(4)+3);
    a9 = str2double(a9);% storing the ninth coefficient
    a10 = cfline2(G(4)+4:G(5)+3);
    a10 = str2double(a10);% storing the tenth coefficient
    
    %line 3-Reading the coefficient data
    
    cfline3= fgetl(file);
    I = strfind(cfline3,'E');
    a11 = cfline3(1:I(1)+3);
    a11 = str2double(a11); % storing the eleventh coefficient
    a12 = cfline3(I(1)+4:I(2)+3);
    a12 = str2double(a12); % storing the twelfth coefficient
    a13 = cfline3(I(2)+4:I(3)+3);
    a13 = str2double(a13); % storing the thirteen coefficient
    a14 = cfline3(I(3)+4:I(4)+3);
    a14 = str2double(a14);% storing the fourteenth coefficient
    
    % First 7 high temperature & second 7 low temperature
    
    Higher_temp = linspace(mt,Higher_temp,100);
    Lower_temp = linspace(Lower_temp,mt,100);
    
    % specific heat of high and low temperature
    
    Cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Lower_temp,Higher_temp ,mt);
    
    % enthalphy of the given data
    
    H=enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Lower_temp,Higher_temp,mt);
    
    %Entropy
    
    S=entrophy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Lower_temp,Higher_temp,mt);
    
    % saving as a folder
    compound=char(compounds);
    fd = sprintf('%s',compound);
    cur_dir=pwd;
    mkdir(compound);
    cd(compound);
    
    %plot 1, specific heat with high temperature vs low temperature
    figure(1)
    plot(Lower_temp,Cp,'linewidth',3)
    hold on
    plot(Higher_temp,Cp,'linewidth',3)
    hold off
    figure(2)
    plot(Lower_temp,H,'linewidth',3)
    hold on
    plot(Higher_temp,H,'linewidth',3)
    hold off
    figure(3)
    plot(Lower_temp,S,'linewidth',3)
    hold on
    plot(Higher_temp,S,'linewidth',3)
    hold off
    % saving the plot in folder
    saveas(figure(1),'cp.png');
     saveas(figure(2),'H.png');
      saveas(figure(3),'S.png');
      cd(cur_dir);
   % calculating molecular weight
      h=0;c=0;n=0;o=0;ar=0;
   for k=1:length(compounds)
       if strcmp(compounds(k),'H')==1
           if k==length(compounds)
               h=h+1;
           else
               if isstrprop(compounds(k+1),'digit')==1
                   h=h+str2double(compounds(k+1));
               elseif isstrprop(compounds(k+1),'digit')==0
                   h=h+1;
               end
           end
       elseif strcmp(compounds(k),'C')==1
           if k==length(compounds)
           c=c+1;
       else
           if isstrprop(compounds(k+1),'digit')==1
               c=c+str2double(compounds(k+1));
           elseif isstrprop(compounds(k+1),'digit')==0
               c=c+1;
           end
           end
       elseif strcmp(compounds(k),'N')==1
           if k==length(compounds)
               n=n+1;
           
           else
               if isstrprop(compounds(k+1),'digit')==1
                   n=n+str2double(compounds(k+1));
               elseif isstrprop(compounds(k+1),'digit')==0
                   n=n+1;
               end
           end
   
           elseif strcmp(compounds(k),'O')==1
               if k==length(compounds)
                   o=o+1;
               else
                   if isstrprop(compounds(k+1),'digit')==1
                       o=o+str2double(compounds(k+1));
                   elseif isstrprop(compounds(k+1),'digit')==0
                       o=o+1;
                   end
               end
               elseif strcmp(compounds(k),'AR')==1
                   if k==length(compounds)
                       ar=ar+1;
                   else
                       if isstrprop(compounds(k+1),'digit')==1
                           ar=ar+str2double(compounds(k+1));
                       elseif isstrprop(compounds(k+1),'digit')==0
                           ar=ar+1;
                       end
                   end
                   elseif isstrprop(compounds(k),'digit')==1
                       continue;
  
       end
   end

H=1;C=12;N=14;O=16;AR=40;
%displaying the molecular weight for  each compounds
disp(compounds)
molecular_weight=H*h+C*c+N*n+O*o+AR*ar;


end

Explanation for Main programme :-

• The program is started with commands clear all , close all, and clc to establish a fresh start without the manipulation of previous values.
• The given data file 'THERMO.dat' file is read using command fopen which is assigned to variable file .
• The command fgetl is used to navigate to the next line which is the heading of the file.
• Same command is used again to read the next line wgich give us the overal global temperature ranging from min to max with middle value.
• Minimum , mid and max global temperature is split and assigned to variables using command strsplit & str2double.
• then the loop has been continue for 53 times because we given 53 compounds.
• Again use the command fgetl, str2split, str2double we find the local minimum , mid and max temperature as wel as find the compounds data for all 53 loop.
• The hgh temperature coefficient is extracted by finding the position of E in each line
• So this will repeated for 7 high temp coefficient (a1,a2,a3,a4,a5,a6,a7) and 7 low temperature coefficient (b1,b2,b3,b4,b5,b6,b7).
• the function of cp,H,Sand MW the values are inputed automatically for the thermodynamic properties will be calculated.
• then we have to plot for monoatomic gas species for specific heat vs temperature, enthalphy vs temperature and entropy vs temperature
• then molecular function we have to give in for loop command
• After that it will print the molecular weight in command window.
• then loop will ended.

The errors faced while programming.

No error .

Output :-

N2

o2

co2

Molecular Weight :-