Project 1 - Parsing NASA thermodynamic data

OBJECTIVE

1. To 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. To calculate the molecular weight of each species.
3. Tolot the Cp, Enthalpy, and Entropy for the local temperature range (low temperature: high temperature) specific for each species.
4. To save the plots as images with appropriate names and dump them in separate folders for each species with a suitable name.

EQUATIONS USED

`(CP)/R=a1+a2*T+a 3*T^2+a4*T^3+a5*T^4`

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

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

where,

T is the temperature.

CP is the specific heat.

H is the enthalpy.

S is the entropy.

R is the universal gas constant.

PROCEDURE AND THE EXPLAINATION FOR THE CODE

• Will starts with the regular clear commands clear all, close all, clc.
• the first thing here we are going to open the file by using the open command i.e, f1=fopen(f1), and then we are going to read the header information in the first line of the file that will be using fgetl command i.e, gt=fgetl(f1), here we can get the first line i.e, 'THERMO'. And we are going to assign value for k=1 which we will discuss later.
• and then we are going to read the next line which is just containing the temperature values that is also by using the command fgetl i.e, gt1=fgetl(f1). then we are going to split the string into characters by using the command i.e, gt1=strsplit(gt1). And then we are going to convert those characters into numbers which are the global temperatures by using the command i.e, global_low_temp=str2double(gt1{2}), global_mid_temp=str2double(gt1{3}), global_high_temp=str2double(gt1{4}), here we get the temperature values.
• And then we need to skip the next three lines which are actually not containing the user data. To skip those lines we are creating a for loop.
• and again we are going to create another for loop in which we are going to extract the first 14 coefficients and going to calculate for enthalpy, entropy, and specific heat and going to plot for it with respect to temperatures and also we are going to plot it for all the 53 species. And also we are going to calculate the molecular weight of each species.
• in the loop we are going to read the line from which we can get the local temperatures, we will take it as line =fgetl(f1), and next we have to split it i.e, line_split=strsplit(line), so that we can able to convert it into numbers i.e, local_low_temp=str2double(line_split{end-3}), local_high_temp=str2double(line_split{end-2}), local_mid_temp=str2double(line_split{end-1}). These are the local temperature values.
• Now, we need to fetch the coefficients from the next three lines to do that first we are going to read the line i.e, tline=fgetl(f1) and then we are going to find the position of the letter 'E' which is present in the line so that we can get the coefficients to find the position of 'E' we can use the command i.e,
• a=findstr(tline, E) after finding the position of 'E' we can be able to get the values of the coefficients, here, after every 'E' it is having three numbers that should also we need to consider so that can be done as, a1=tline(1:a(1)+3) and the next will be written as a2=tline(a(1)+4:a(2)+3) and so on up to a5=tline(a(4)+4:a(5)+3), along with fetching the coefficients we can convert a1, a2, ....a5 into numbers i.e, a1=str2double(a1) and so on upto a5 can be done. And the next coefficients are present in the next lines and again the same procedure will be repeated for the next two lines.
• After fetching the coefficients we are going to calculate for the specific heat, Enthalpy, and Entropy to calculate these three we need to assign the values for temperatures so, by using linspace command i.e, T=linspace(local_low_temp,local_high_temp,1000) we are assuming for 1000 temperature values in between. And we are taking R=8.314 the universal gas constant in the calculation.
• Then, we can create a 'for loop' in that we are going to calculate for cp, H, and S from 1 to the length of T, by using the suitable equations in that we are using the 'if' condition, i.e T(i)>local_mid_temp in this we are going to get the values for the local high-temperature ranges and the 'else' condition is going to get the values for the local low temperatures.
• Then, we are going to create the folders in the name of individual species that is to save the plots inside the respective folders. That can be done by using the command i.e, mkdir('D:\Bassu', species_name), here we need to get the current path and in which name we are going to create folders that are species names.
• And then, we are going to plot for the temperature vs cp, H, S in figure 1, 2, 3 respectively with respective title and labels required and we are going to assign the figure handling i.e, P1=figure(k), P2=figure(k+1), P3=figure(k+2) as we assigned the value for 'k' as 1 in the beginning that we are using to save the plots as figure number with respective plots of each species and we are incrementing the value of k i.e, k=k+3, so that it can assign the figures of respective species with respective number.
• Next, we are going to create the paths to print the figures into images in the specific format i.e, path1=sprintf('D:/Bassu/%s/temp_vs_cp.png', species_name), here this image will be going to save in the png format that is of temperature vs specific heat-plot and the below two will be for Enthalpy and Entropy respectively.

  path2 =sprintf('D:/Bassu/%s/temp_vs_H.png', species_name)

  path3 =sprintf('D:/Bassu/%s/temp_vs_S.png', species_name)

• then, we need to get the command to save the images to particular paths to do that we are using the 'saveas' command, i.e, savesa(P1,path1)

  saveas(P2,path2)

  saveas(P3,path3), here P1, P2, P3 are the figure handles, and paths are assigned to save the images in the particular folder.

• Then, we are going to calculate for the molecular weight of each species to calculate that we are first going to display all the species name by using the command, disp(species_name) and then going to assign the elements which are present in the data file to get the elements we need to check the data file, and here we have element={'O', 'C', 'H', 'A', 'R', 'N', 'S'} actually here Argon is assigned as two different elements because we are going to compare the elements with the species name that is one by one so it will not be going to consider argon as a single element so we are taking it as different elements and then we are going to get the respective atomic masses i.e, Atomic_mass=[15.999,12.011,1.008,39.948,0,14.007,32.06].
• next, initializing the molecular weight of the array, i.e, m_wt=zeros(1,length(species_name)).
• using 'for loop' to calculate the molecular weight of each species, for p it will be going to iterate from 1 to the length of the species name, and then we are assigning the condition that if the character of the species name is a letter then it will going to take the respective atomic mass of the element. Again we are using 'for loop' for q it is 1 to the length of the element here going to compare with the species name. if the element matches with the species name it will be going to get the respective atomic mass, here we are comparing the species name with the element by using the command, with 'if' condition, strcmp(species_name(p), element(q)) then it will going to calculate the molecular weight as we assigned i.e, m_wt(p)=m_wt(p)+Atomic_mass(q) that is if the element is letter and 'else' condition will be in case it is a number it will be multiplied with the previous letter as that we assigned i.e, m_wt(p-1)=m_wt(p-1)* str2double(species_name(p). and the molecular weight will be the sum of the 'm_wt' that is molecular_weight=sum(m_wt). The molecular weight will be displayed in the command window that is shown in the screenshots.

CODE EXPLAINED WITH THE COMMENTS

clear all
close all
clc




% we are going to open the file with the file handling variable f1.


f1=fopen('THERMO.dat','r'); 

gt=fgetl(f1);           %reading the line here.



% here we are splitting the gt1 to convert it into a number.


k = 1;                   % we are assigning the k value as 1 this we can use in saving the figures.
gt1=fgetl(f1);           % reading the line from where we can obtain the global temperature values.
    
gt1 = strsplit(gt1);     %splitting the string to convert it into character.



% now converting the cell array into a number




global_low_temp=str2double(gt1{2});     % here we are using the str2double command to convert the cell array into numbers.
global_mid_temp=str2double(gt1{3});
global_high_temp=str2double(gt1{4});

% now skipping the next three lines which are not having the useful data
% to skip those we are using the below for loop.



 for l=1:3
     skip_line=fgetl(f1);
 end
    
 
 
 

 for j=1:53
     
     
     
% now we need to extract the lines from line6 in the file which are having the useful data.



line=fgetl(f1);           % reading the sixth line.


% splitting the sixth line in the data file.


line_split=strsplit(line);




% converting the cell array into a number


atomic_num=str2double(line_split{4});              %in the line it is the atomic number.

local_low_temp=str2double(line_split{end-3});      % these are the local temperature values.

local_high_temp=str2double(line_split{end-2});

local_mid_temp=str2double(line_split{end-1});


%line_num=str2double(line_split{9});


species_name = line_split{1};


% now we need to get the 6th line where we are actually going to fetch the
% coefficients.

tline=fgetl(f1);


% we are going to get the locations of the E so that we can find out the
% coefficients of it.


a=findstr(tline,'E');

% we found the locations of E and with the help of these we can find the
% coefficients, we can get 3 numbers after E using that we can find the
% coefficients. The below 7coefficients are for high temperature.

a1=tline(1:a(1)+3);

a1=str2double(a1);

a2=tline(a(1)+4:a(2)+3);

a2=str2double(a2);

a3=tline(a(2)+4:a(3)+3);

a3=str2double(a3);

a4=tline(a(3)+4:a(4)+3);

a4=str2double(a4);

a5=tline(a(4)+4:a(5)+3);

a5=str2double(a5);



% another 2 coefficients are present in the next line so to read the next
% line we are taking it as tline_2.


tline_2=fgetl(f1);



b=findstr(tline_2,'E');          % finding  the locations for E.

a6=tline_2(1:b(1)+3);            % we got the locations of E so here we are fetching the coefficients again.

a6=str2double(a6);

a7=tline_2(b(1)+4:b(2)+3);

a7=str2double(a7);


% the below 7 coefficients are for low temperatures

a8=tline_2(b(2)+4:b(3)+3);

a8=str2double(a8);

a9=tline_2(b(3)+4:b(4)+3);

a9=str2double(a9);

a10=tline_2(b(4)+4:b(5)+3);

a10=str2double(a10);

% the below 4 coefficients are in the next line.

tline_3=fgetl(f1);              % reading the line.

c=findstr(tline_3,'E');         % finding the location of the E.

a11=tline_3(1:c(1)+3);          %we found the locations and again we are going to fetch the remaining coefficients.

a11=str2double(a11);

a12=tline_3(c(1)+4:c(2)+3);

a12=str2double(a12);

a13=tline_3(c(2)+4:c(3)+3);

a13=str2double(a13);

a14=tline_3(c(3)+4:c(4)+3);

a14=str2double(a14);


% assuming the temperature values from low to high global temperature as
% 1000
%  taking the value for a universal gas constant.

R=8.314;


 T=linspace(local_low_temp,local_high_temp,1000);    % we are assigning the temperature values as 1000 intervals between local low temperature and local high temperature.

 
%  calculating the specific heat, enthalpy, entropy 

%we are creating a for loop to perform the calculation for specific heat,
%enthalpy, and entropy by using the equations as follows from 1 to length
%of the temperature. Here the 'if' condition is going to evaluate for high
%temperature and the 'else' condition will be for low temperatures.



 for i=1:length(T)
     
     if T(i)>local_mid_temp
         
       cp(i)=(a1+a2.*T(i)+a3.*T(i)^2+a4.*T(i)^3+a5.*T(i)^4).*R;
       H(i)=(a1+a2.*(T(i)./2)+a3.*(T(i)^2./3)+a4.*(T(i)^3./4)+a5.*(T(i)^4./5)+a6./T(i)).*R.*T(i);
       S(i)=(a1*log(T(i))+a2.*T(i)+a3.*(T(i)^2./2)+a4.*(T(i)^3./3)+a5.*(T(i)^4./4)+a7);
         
      
 
     else 
       cp(i)=(a8+a9.*T(i)+a10.*T(i)^2+a11.*T(i)^3+a12.*T(i)^4).*R;
       H(i)=(a8+a9.*(T(i)./2)+a10.*(T(i)^2./3)+a11.*(T(i)^3./4)+a12.*(T(i)^4./5)+a13./T(i)).*R.*T(i);
       S(i)=(a8*log(T(i))+a9.*T(i)+a10.*(T(i)^2./2)+a11.*(T(i)^3./3)+a12.*(T(i)^4./4)+a14);
         
    
     
     end
 end
 
 
%creating the folders with the name of the species to save the plots
%  within the folder


   mkdir('D:\Bassu',species_name)       % here folders will be generated with respective species names.
   
%  plotting the graph for temperature vs specific heat 

 P1 = figure(k);                        % here we are assigning the figure with the figure handle that we can use in the save as command to save the figures.

 plot(T,cp);                                %plotting temperature vs cp
 
 xlabel('Temperature');                   % labeling x axis as temperature values.
 
 ylabel('Specific heat');                 %labeling y axis as speciifc heat values.
 
 title('Temperature vs Specific heat for ',species_name);       %giving some title to the plot.

%  plotting the graph for temperature vs enthalpy
 

P2 = figure(k+1);                     % here we are assigning the figure with the figure handle that we can use in the save as command to save the figures.

plot(T,H)                             %plotting temperature vs H.

xlabel('Temperature')                     % labeling x axis as temperature values.
 
ylabel('Enthalpy'))                   %labeling y axis as Enthalpy values.

title('Temperature vs Enthalpy for ' ,species_name)      %giving some title to the plot.

% plotting the graph for temperature vs entropy 

P3 = figure(k+2);                   % here we are assigning the figure with the figure handle that we can use in the save as command to save the figures.

plot(T,S);                           %plotting temperature vs Entropy.

xlabel('Temperature');               % labeling x axis as temperature values.

ylabel('Entropy');                   %labeling y axis as Entropy values.

title('Temperature vs Entropy for ' ,species_name);          %giving some title to the plot.
 
 
%  we are going to save the plots as images inside the respective folders
%  of particular species by assigning the path and using the 'saveas' command
%  to save the figure

 k = k+3;                  % here we are giving some increment to k value so that it can get the next figures plotted with the respective number.

path1 =sprintf('D:/Bassu/%s/temp_vs_cp.png', species_name);     %creating the paths to print the images in the particular species folder and in which format it should be.

path2 =sprintf('D:/Bassu/%s/temp_vs_H.png', species_name) ;

path3 =sprintf('D:/Bassu/%s/temp_vs_S.png', species_name);


saveas(P1,path1);           % here we can save the images of plots by assigning the figure handle and path in the 'saveas' command.

saveas(P2,path2);

saveas(P3,path3);

% calculating the molecular weight of each species


 disp(species_name)         % it will display the particular species name to the respective molecular weight.


 element={'O','C','H','A','R','N','S'};     % these are the only elements present in the thermo data so we are using these elements to compare with the species name to calculate the molecular weight.
 
 Atomic_mass=[15.999,12.011,1.008,39.948,0,14.007,32.06];   % these are the respective atomic masses of the element.
 
 
 m_wt=zeros(1,length(species_name));                  %initializing the molecular weight of the array.
  
 %using for loops to calculate the molecular weight of each species.
 
 
 for p=1:length(species_name)              % for p, it will be going to iterate from 1 to the length of the species name.
     
     if isnan(str2double(species_name(p)))  % here assigning the condition that if the character of the species name is letter then it will going to take the respective atomic mass of the element.
         
         for q=1:length(element)            % for 1 to the length of element here going to compare with the species name.
             
             if strcmp(species_name(p),element(q))
            
                 m_wt(p)=m_wt(p)+Atomic_mass(q);         %if the element matches with the species name it will going to get the respective mass.
            
             end
         end
         
     else
         
        m_wt(p-1)=m_wt(p-1)* str2double(species_name(p));   %in case it is a number it will be multiplied with the previous letter.
         
     end
 end
             
   
          molecular_weight=sum(m_wt)         %this will be the molecular weight of the particular species.
 
 end 



 
 

OUTPUT

The following screenshot is of the window where all the folders for particular species are saved or created in the directory.

• The plots for O2 species are as follows.

• 

The following plots for N2.

The following plots for CO2.

The screenshot attached here is of the molecular weights of some of the species displayed in the command window and we have calculated for all the 53 species.

ERRORS MADE

The following error was an incorrect dimension to a raising matrix and was corrected by using the dot before the rising power.

and also I made a lot of mistakes but couldn't be captured.

CONCLUSION

• The program to extract the coefficients from the data file and calculating them to get the Specific heat, Enthalpy, and Entropy values to get the plots of temperature versus CP, H, S has done. And the plots are also made for all the 53 species and saved in the respective created folders.
• Molecular weight for all 53 species also been calculated.