Project 1 - Parsing NASA thermodynamic data

File Parsing :

         The given NASA polynomials have the form: 

where a1, a2, a3, a4, a5, a6, and a7 are the numerical coefficients supplied in the NASA thermodynamic file. The first 7 numbers starting on the second line of each species entry (five of the second line and the first two of the third line) are the seven coefficients (a1 through a7, respectively) for the high-temperature range (above 1000 K, the upper boundary is specified on the first line of the species entry). The following seven numbers are the coefficients (a1 through a7, respectively) for the low-temperature range (below 1000 K, the lower boundary is specified on the first line of the species entry)

 Main program:    

1. The fgetl() command is used to read lines from a file. The syntax is given as follows:

                       Tline = fgetl(FileID)

2. The str2num() command is used to convert character array or string to numeric array.

           Syntax:  x = str2num(chr)    

3. The strsplit() is used to split a character or string vector with a specified delimiter.

           Syntax:  C = strsplit(str,delimiter)

4. The findstr() is used to find a string within another string.

           Syntax:  K = findstr(str1,str2)

5. mkdir is used make a new folder. 

           Syntax:  mkdir foldername  

We use the fgetl() to read each line of the Thermo.dat file and use a for loop to read the coeffecients of all the 53 species.

%File parsing NASA thermodynamic data
clear all
close all
clc

%Reading NASA Thermodynamic data  
f1 = fopen('THERMO.dat','r');

%Gas constant value
R = 8.314;

%Skipping the first line 
fgetl(f1);

%Fetching Global temperature from second line 
A = fgetl(f1);
Global_temp_ranges = strsplit(A,' ');
Global_low_temp = str2num(Global_temp_ranges{1});
Global_mid_temp = str2num(Global_temp_ranges{2});
Global_high_temp = str2num(Global_temp_ranges{3});
%Skipping the comments
for i = 1:3
    fgetl(f1);
end 
%Looping the program to fetch data from all species
 for i = 1:53  
     %Getting first line from the species
       B = fgetl(f1);
       C = strsplit(B,' ');
       %Getting species name
       Species = C{1};
       %Getting local temperature ranges of the species
       Loc_temp = strsplit(B,' ');
       Local_low_temp = str2num(Loc_temp{end-3});
       Local_mid_temp = str2num(Loc_temp{end-1});
       Local_high_temp = str2num(Loc_temp{end-2});
       
       %Temperature value
       temp = linspace(Local_low_temp,Local_high_temp,1000);
       
       %Fetching the coefficient values
       line_1 = fgetl(f1);
       b = findstr(line_1,'E'); 
       %First 7 coeffecients(for temp > 1000 K)
       a1 = line_1(1:b(1)+3);
       a_1 = str2num(a1);
       a2 = line_1(b(1)+4:b(2)+3);
       a_2 = str2num(a2);
       a3 = line_1(b(2)+4:b(3)+3);
       a_3 = str2num(a3);
       a4 = line_1(b(3)+4:b(4)+3);
       a_4 = str2num(a4);
       a5 = line_1(b(4)+4:b(5)+3);
       a_5 = str2num(a5);
       line_2 = fgetl(f1);
       c = findstr(line_2,'E');
       a6 = line_2(1:c(1)+3);
       a_6 = str2num(a6);
       a7 = line_2(c(1)+4:c(2)+3);
       a_7 = str2num(a7);
       
       %Last 7 coeffecients(for temp < 1000 K)
       a8 = line_2(c(2)+4:c(3)+3);
       a_8 = str2num(a8);
       a9 = line_2(c(3)+4:c(4)+3);
       a_9 = str2num(a9);
       a10 = line_2(c(4)+4:c(5)+3);
       a_10 = str2num(a10);
       line_3 = fgetl(f1);
       d = findstr(line_3,'E');
       a11 = line_3(1:d(1)+3);
       a_11 = str2num(a11);
       a12 = line_3(d(1)+4:d(2)+3);
       a_12 = str2num(a12);
       a13 = line_3(d(2)+4:d(3)+3);
       a_13 = str2num(a13);
       a14 = line_3(d(3)+4:d(4)+3);
       a_14 = str2num(a14);
     for j = 1:length(temp)   
      %Calculating Specific heat,Enthalpy & Entropy values 
      Cp(j) = Specific_heat(a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9,a_10,a_11,a_12,a_13,a_14,R,Local_mid_temp,temp(j));
      H(j) = Enthalpy(a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9,a_10,a_11,a_12,a_13,a_14,R,Local_mid_temp,temp(j));
      S(j) = Entropy(a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9,a_10,a_11,a_12,a_13,a_14,R,Local_mid_temp,temp(j));
     end
     
     %Calculating molecular weight of sepcies
      mol_weight(i) = molecular_weight(Species);
      
      %Creating new folder
      mkdir(Species)
      
      %Changing directory
      cd(Species)
      
      %Plotting Specific heat w.r.t Temperature range 
      figure(1)  
      plot(temp,Cp,'linewidth',2,'color','r')
      xlabel('Temperatue(in K)')
      ylabel('Specific heat(in KJ)')
      title(sprintf('Specific heat vs Temperature range of %s',Species))
      grid on
      saveas(1,'Specific heat.jpg')
      
      %Plotting Enthalpy w.r.t Temperature range
      figure(2)
      plot(temp,H,'linewidth',2,'color','g')
      xlabel('Temperatue(in K)')
      ylabel('Enthalpy(in KJ/mol)')
      title(sprintf('Enthalpy vs Temperature range of %s',Species))
      grid on
      saveas(2,'Enthalpy.jpg')
      
      %Plotting Entropy w.r.t Temperature range
      figure(3)
      plot(temp,S,'linewidth',2,'color','b')
      xlabel('Temperatue(in K)')
      ylabel('Entropy(in KJ/mol-K)')
      title(sprintf('Entropy vs Temperature range of %s',Species))
      grid on
      saveas(3,'Entropy.jpg')
      
      cd('C:\Users\Admin\Documents\MATLAB')
 end   
   

 Specific Heat:

             Specific heat, the quantity of heat required to raise the temperauture of one gram of a substance by one degree Celcius.

function Cp = Specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Local_mid_temp,t)

      if t <= Local_mid_temp                  
          Cp = R*(a8+a9.*t+a10.*t.^2+a11.*t.^3+a12.*t.^4);
      else 
         Cp = R*(a1+a2.*t+a3.*t.^2+a4.*t.^3+a5.*t.^4);
      end
end

Enthalpy:

            The enthalpy H  of a thermodynamic system is defined as the sum of its internal energy U and the work required to achieve its pressure and volume.

function H = Enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Local_mid_temp,t)

      if t <= Local_mid_temp                  
        H = R.*t.*(a8+a9.*t./2+a10.*(t.^2)./3+a11.*(t.^3)./4+a12.*(t.^4)/5+a13./t);
      else 
          H = R.*t.*(a1+a2.*t./2+a3.*(t.^2)./3+a4.*(t.^3)./4+a5.*(t.^4)/5+a6./t); 
      end
end

 Entropy:

          Its is a thermodynamic quantity representing the unavailability of a system's thermal energy for conversion into mechanical work, often interpreted as the degree of disorder or randomness in the system.

function S = Entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,R,Local_mid_temp,t)

      if t <= Local_mid_temp                  
        S = R*(a8.*log(t)+a9.*t+a10.*(t.^2)./2+a11.*(t.^3)./3+a12.*(t.^4)./4+a14);
      else 
        S = R*(a1.*log(t)+a2.*t+a3.*(t.^2)./2+a4.*(t.^3)./3+a5.*(t.^4)./4+a7)   
      end
end

Molecular weight: 

                The molecular mass is the mass of the given molecule. It is measured in daltons.Different molecule of the same compound may have different molecular mass because they contain different isotope of the element. 

               The atomic weight of few atoms are gathered from the internet. Then we use two for loops to compare the strings of the species with the atomic names and therefore calculating the molecular weight of the species. 

function mol_weight = molecular_weight(Species)
     Atomic_name = ['H','C','O','N','A','S']
     Atomic_weight = [1.008,12.011,15.99,14.007,39.948,32.065]
     mol_weight = 0;
     
     for i = 1:length(Species)
        for j = 1:length(Atomic_name)
             if strcmpi(Species(i),Atomic_name(j))
                 mol_weight = mol_weight + Atomic_weight(j);
                pos = j;
             end
        end
        n = str2num(Species(i));
        if n>1
             mol_weight = mol_weight+Atomic_weight(pos)*(n-1);
        end
     end
      sprintf('Molecular weight of %s is %f u \n',Species,mol_weight)
end

Errors:

       1. Simple errors such as missing paranthesis where faced.      

References:

 1. http://combustion.berkeley.edu/gri-mech/data/nasa_plnm.html

 2. https://in.mathworks.com/help/matlab/ref/mkdir.html

 3. https://in.mathworks.com/help/matlab/ref/fgetl.html

 4. https://in.mathworks.com/help/matlab/ref/str2num.html

 5. https://in.mathworks.com/help/matlab/ref/strsplit.html

 6. https://in.mathworks.com/help/matlab/ref/findstr.html

 7. https://link.springer.com/chapter/10.1007/978-1-4684-0186-8_8