Parsing NASA Thermodynamic data-Matlab

Aim:

1. Write a function that extracts the 14 coefficients and calculates the enthalpy, entropy and specific heats for all the species in the data file.
2. Calculate the molecular weight of each species and display it in the command window.
3. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
4. Save the plots as images with appropriate names and dump them in separate folders for each species with a suitable name. All these should be generated automatically.
5. The main program should take the "THERMO.dat" file as an input and generate Cp, H and S plots.

Objective:

File parsing of 'THERMO.dat'

Red-global temperatures

Yellow-Local temperatures

Blue-Elements or Species

Formulae:

This is the formulae to calculate the specific heat, enthalpy and entropy.

Here R is the universal gas constant, T is the temeprature

This is the formula to calculate the molecular mass. It is sum of Atomic Number times the number of atoms.

Procedure:

1. Initially function program is made with parameters such as Gas constant, coefficients and Temperature values.

2. In main program, Gas constant is defined and the file THERM.dat is opened using fopen command.

3. fgetl is used further to read each and every line as per the requirement for parsing.

4. strsplit is often used to split the values after reading through fgetl.

5. str2double is used to convert the str to numbers with double precision.

6. A for loop is used to skip the comments.

7. An another for loop is used to read the values from as per the convinience.

8. mkdir is used to create new directory.

9. Temperature is discretized using linspace and the function call is made

10. Plotting is done for cp,enthalpy and entropy.

11. Now, Molecular weight is calculated for each species

Program:

Function Program:

function [cp H S]=coeff(R,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,mid_T)

for i=1:length(T)
    if T(i)>mid_T
      cp=(a1+a2.*T+a3.*T.^2+a4.*T.^3+a5*T.^4).*R;
      H=(a1+a2.*T./2+a3.*T.^2./3+a4.*T.^3./4+a5.*T.^4./5+a6./T).*R.*T;
      S=(a1.*log(T)+a2.*T+a3.*T.^2./2 + a4.*T.^3./3+a5.*T.^4./4+a7).*R;
    else
      cp=(a8+a9.*T+a10.*T.^2+a11.*T.^3+a12*T.^4).*R;
      H=(a8+a9.*T./2+a10.*T.^2./3+a11.*T.^3./4+a12.*T.^4./5+a13./T).*R.*T;
      S=(a8.*log(T)+a9.*T+a10.*T.^2./2 + a11.*T.^3./3+a12.*T.^4./4+a14).*R;  
    end
    
end

Main Program:

clear all
close all
clc

R=8.314; %Universal Gas constant

file=fopen('THERMO.dat','r'); %opens the file and read mode is on
title=fgetl(file);  %reads the first line and it repeats the process again and again for next consecutive lines
tline1=fgetl(file);
gtemp=strsplit(tline1,' '); %global temperature

%global temperatures 
g_low_temp=str2double(gtemp{2});
g_mid_temp=str2double(gtemp{3});
g_high_temp=str2double(gtemp{4});

%skipping of comments
for i=1:3 
    comment=fgetl(file);
end

%computation loop
for j=1:53
      tline2=fgetl(file);
      header=strsplit(tline2,' ');
      species=header{1};      %header of the file
      g_find=strfind(tline2,'G ');
      loctemp=tline2(g_find+4:length(tline2));
      ltemp=strsplit(loctemp,' ');
      
      l_low_temp=str2double(ltemp{1});
      l_mid_temp=str2double(ltemp{2});
      l_high_temp=str2double(ltemp{3});
      
      line1=fgetl(file);
      l1=strfind(line1,'E');
      a1=str2double(line1(1:(l1(1)+3)));
      a2=str2double(line1((l1(1)+4:l1(2)+3)));
      a3=str2double(line1((l1(2)+4:l1(3)+3)));
      a4=str2double(line1((l1(3)+4:l1(4)+3)));
      a5=str2double(line1((l1(4)+4:l1(5)+3)));
      
      line2=fgetl(file);
      l2=strfind(line2,'E');
      a6=str2double(line2(1:(l2(1)+3)));
      a7=str2double(line2((l2(1)+4:l2(2)+3)));
      a8=str2double(line2((l2(2)+4:l2(3)+3)));
      a9=str2double(line2((l2(3)+4:l2(4)+3)));
      a10=str2double(line2((l2(4)+4:l2(5)+3)));
      
      line3=fgetl(file);
      l3=strfind(line3,'E');
      a11=str2double(line3(1:(l3(1)+3)));
      a12=str2double(line3((l3(1)+4:l3(2)+3)));
      a13=str2double(line3((l3(2)+4:l3(3)+3)));
      a14=str2double(line3((l3(3)+4:l3(4)+3)));
      
      %creating a new directory to save the plots
      file_parsing=species;
      mkdir(file_parsing);
      dir=pwd;
      cd(file_parsing);
      
      %discretisation of temps to separate low and high temps      
      T=linspace(l_low_temp,l_high_temp,50);
      mid_T=l_mid_temp;
      
      %function call
      [cp H S]=coeff(R,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,mid_T);
      
      
      %Plotting
      figure(1)
      plot(T,cp,'linewidth',2,'color','r');
      xlabel('Temperature[K]');
      ylabel('Specific heat [J/mol-K]');
          
      
      figure(2)
      plot(T,H,'linewidth',2,'color','g');
      xlabel('Temperature[K]');
      ylabel('Enthalpy [J/kg]');
      
      
      figure(3)
      plot(T,S,'linewidth',2,'color','b');
      xlabel('Temperature[K]');
      ylabel('Entropy [J/K]');
      
      directory=('C:\Users\91790\Desktop\');
      
      % saving the plots inside folders
           
      folder1=fullfile(directory,species,'cp.jpg');
      folder2=fullfile(directory,species,'H.jpg');
      folder3=fullfile(directory,species,'S.jpg');
      
      saveas(figure(1),folder1);
      saveas(figure(2),folder2);
      saveas(figure(3),folder3);
      
      
      %Calculation of Molecular mass
      
      element=['O','H','C','N','A']; %elements in the species
      atomic_wt=[16 1 12 14 40]; %atomic weight corresponding the elements
      molecular_mass=0; %initial value of molecular mass
      for i=1:length(species)
          for j=1:length(atomic_wt)
              if strcmp(species(i),element(j))
                  molecular_mass=molecular_mass+atomic_wt(j);
                  location=j;
              end
          end
          
      n=str2double(species(i));
      if n>1
          molecular_mass=molecular_mass+atomic_wt(location)*(n-1);
      end
      
      end
      
      fprintf('molecular_weight of %s:',species);
      fprintf('%d',molecular_mass);
      
 end
      
      fclose(file);
      
    
      
      
      
      

 

Result:

Plots for O:

For O2:

Note: All other plots are attached in the document .zip file

Result of Molecular weights: