Parsing the NASA thermodynamic data

MAIN CODE:

clear all
close all
clc

%inputs
%Universal gas constant R
R = 8.314;

% Opening the data file 
f = fopen(\'THERMO.dat\',\'r\');

heading = fgetl(f); %reading the header

Global_temp_range = fgetl(f);
Global_temp_range = str2num(Global_temp_range);

% reading comments
line1 = fgetl(f);
line2 = fgetl(f);
line3 = fgetl(f);

% iterating the loop for 53 species

for i = 1:53  
line0 = fgetl(f);
species_name = strtok(line0);

% finding G in line0 

a = findstr(line0,\'G   \');
local_temp = line0(a+1:end-1);
local_temp = str2num(local_temp);

%  local temperature range 

t = linspace(local_temp(1),local_temp(2),500);


line1 = fgetl(f);
line2 = fgetl(f);
line3 = fgetl(f);

 % finding the position of E 
   a = strfind(line1,\'E\');                
   b = strfind(line2,\'E\');
   c = strfind(line3,\'E\');
    
    % high temp coefficients
    h1 = str2double(line1(1:a(1)+3));       
    h2 = str2double(line1(a(1)+3:a(2)+3));
    h3 = str2double(line1(a(2)+4:a(3)+3));
    h4 = str2double(line1(a(3)+4:a(4)+3));
    h5 = str2double(line1(a(4)+4:a(5)+3));
    h6 = str2double(line2(1:b(1)+3));
    h7 = str2double(line2(b(1)+4:b(2)+3));
    
    % low temp coefficients
    l1 = str2double(line2(b(2)+4:b(3)+3));
    l2 = str2double(line2(b(3)+4:b(4)+3));
    l3 = str2double(line2(b(4)+4:b(5)+3));
    l4 = str2double(line3(1:c(1)+3));
    l5 = str2double(line3(c(1)+4:c(2)+3));
    l6 = str2double(line3(c(2)+4:c(3)+3));
    l7 = str2double(line3(c(3)+4:c(4)+3));
    
 
    %Calculating Specific heat, enthalpy , entropy and molecular weight
    
   Cp = specific_heat(l1,l2,l3,l4,l5,h1,h2,h3,h4,h5,t,R,local_temp);                      
    H = enthalpy(l1,l2,l3,l4,l5,l6,h1,h2,h3,h4,h5,h6,t,R,local_temp);                      
    S = entropy(l1,l2,l3,l4,l5,l7,h1,h2,h3,h4,h5,h7,t,R,local_temp);                                                                                        
    mw(i) = molecular_weight(species_name);      
    
    

% saving the plots in multiple folders for each species

e = cd(pwd);
folder = species_name;
mkdir(folder);

% plot specific heat vs temperature

figure(1)
plot(t,Cp,\'linewidth\',1,\'color\',\'g\')
title(sprintf(\'Specific Heat vs Temperature plot for %s\',species_name))
xlabel(\'Temperature\')
ylabel(\'Specific Heat\')

% plot enthalpy vs temperature

figure(2)
plot(t,H,\'linewidth\',1,\'color\',\'r\')
title(sprintf(\'Enthalpy vs Temperature plot for %s\',species_name))
xlabel(\'Temperature\')
ylabel(\'Enthalpy\')

% plot entropy vs temperature

figure(3)
plot(t,S,\'linewidth\',1,\'color\',\'b\')
title(sprintf(\'Entropy vs Temperature plot for %s\',species_name))
xlabel(\'Temperature\')
ylabel(\'Entropy\')


% saving the plots

cd(folder)
saveas(1,\'Cp_vs_temp.png\')
saveas(2,\'H_vs_temp.png\')
saveas(3,\'S_vs_temp.png\')
cd(e)


end

ENTHALPY function code:

function H = enthalpy(l1,l2,l3,l4,l5,l6,h1,h2,h3,h4,h5,h6,t,R,local_temp)

if t > local_temp(3)
    H = R*(l1.*t+l2.*(t.^2)./2+l3.*(t.^3)./3+l4.*(t.^4)./4+l5.*(t.^5)./5+l6);

else
    H=R*(h1.*t+h2.*(t.^2)./2+h3.*(t.^3)./3+h4.*(t.^4)./4+h5.*(t.^5)./5+h6);
end

ENTROPY function code:

function S = entropy(l1,l2,l3,l4,l5,l7,h1,h2,h3,h4,h5,h7,t,R,local_temp)

if t > local_temp(3)
    S = R*(l1.*log(t)+l2.*t+l3.*(t.^2)./2+l4.*(t.^3)./3+l5.*(t.^4)./4+l7);

else
    S = R*(h1.*log(t)+h2.*t+h3.*(t.^2)./2+h4.*(t.^3)./3+h5.*(t.^4)./4+h7);
end

SPECIFIC HEAT function code:

function Cp = specific_heat(l1,l2,l3,l4,l5,h1,h2,h3,h4,h5,t,R,local_temp)

if t > local_temp(3)
    Cp = R*(l1+l2.*(t)+l3.*(t.^2)+l4.*(t.^3)+l5.*(t.^4));

else
    Cp = R*(h1+h2.*(t)+h3.*(t.^2)+h4.*(t.^3)+h5.*(t.^4));
end

MOLECULAR WEIGHT function code:

function mw = molecular_weight(species_name)

species = [\'H\' \'C\' \'N\' \'O\' \'S\' \'R\'];
atomic_weight = [1 12 14 16 32 40];

mw = 0;
for i = 1:length(species_name)
  for j = 1:length(species)
    if strcmpi(species_name(i),species(j))
      mw = mw + atomic_weight(j);
      position = j;
    end
  end
  n = str2num(species_name(i));
  if n>1
    mw = mw + atomic_weight(position)*(n-1);
  end
  end
  fprintf(\'Molecular weight of %s is : \',species_name);
  fprintf(\'%d\\n\',mw); 
end

Explanation to the code:

• The \'thermo.dat\' file is saved in the same directory were the matlab files are saved.
• The functions such as \'enthaply\' , \'entropy\' , \'specific_heat\' and \'molecular_weight\' are defined using the equations stated.
• 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 stored in the form of numbers using\'str2num\' command.And the comments in the next 3 lines are skipped.
• For loop is used to read the data for all the 53 species one at a time. The required data from line0 to line3. Line0 describes the name of species and local temperature range for the species. Line1 to line3 contains 14 coefficients for the species out of which first 7 are high temperature coefficients and the next 7 are low temperature coefficients.
• These coefficients are extracted and converted from string data to numbers using \'str2double\' command.
• Then the local temperature range is assigned to variable t for the plots.
• The defined functions are called to calculate the corresponding properties of the species.
• Finally, \'enthalpy\' \'entropy\' and \'specific heat\' are plotted against temperature. These plots are saved in the sub directories for each of the 53 species labelled with species name.

Output Graphs:

O2:

N2:

CO2:

MOLECULAR WEIGHTS DISPLAYED IN COMMAND WINDOW:

FOLDER CONTAINING THE SUB-DIRECTORIES: