Project 1 - Parsing NASA thermodynamic data
#Aim: 1.To write a program in MATLAB that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the NASA Thermodynamic 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.
The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.
The formulae are shown below.
-----------------------------------------------******************************--------------------------------
clear all
close all
clc
%Nasa thermodynamic file challenge program
%No of species in the thermo dat file
n=53;
%Temperature array
T=linspace(200,3500,100);
R=8.314 ; %J/mol-K
f1=fopen('THERMO.dat','r'); %file 'THERMO.dat' is opened for reading
%fgetl -command used to read lines
%First 5 lines contain only header files and global temperature
%coefficients
fgetl(f1);
fgetl(f1);
fgetl(f1);
fgetl(f1);
fgetl(f1);
%For assigning values to characters map variable is used.
map('C')=12;
map('O')=16;
map('H')=1;
map('N')=14;
map('A')=40;
map('R')=0;
%Array containing all the character present in species in the file
Element_array=['C','O','H','N','A','R'];
%Directory data_for_species created
mkdir data_for_species
%for loop for all 53 species
for i=1:n
X=fgetl(f1); %First line of species element is assigned to X
A=strsplit(X,' '); %X is splitted where ' ' is present
All_Species_name(i) = A(1); %Species elements stored in array- All_Species_name
Species_name=All_Species_name(1,i); %Particular species extracted from array All_Species_name for every iteration
Molecular_weight_of_species_initial=zeros(n,1); %Molecular weight of species initialised with array of (n,1) zeros.
if length(char(Species_name))<=1 %Condition for species having single element
Molecular_weight_of_species(i)=map((char(Species_name))); %Character of species element value mapped
fprintf('<strong>Molecular Weight of the Species- %s is %d \n </strong>',char(Species_name),Molecular_weight_of_species(i))
else %Condition for species having element more than 1
for j=1:(length(char(Species_name))-1) %Loop for calculating molecular weight of species
%Condition for species ,having element to be checked as
%character and next elemnet to it is real number
if any(Element_array(:)==(char(Species_name{:}(j)))) && (str2double(Species_name{:}(j+1)))>= 1;
Molecular_weight_of_species_element(i)=map(char(Species_name{:}(j)))*(str2double(Species_name{:}(j+1)));
%Condition for species ,having element to be checked as
%character and next elemnet to it is also a character and next
%element is not last element in the evaluating species
elseif any(Element_array(:)==(char(Species_name{:}(j)))) && any(Element_array(:)==(char(Species_name{:}(j+1)))) && (j+1 ~= length(Species_name{:}));
Molecular_weight_of_species_element(i)=map(char(Species_name{:}(j)));
%Condition for species ,having element to be checked as
%character and next elemnet to it is also a character and next
%element is last element in the evaluating species
elseif any(Element_array(:)==(char(Species_name{:}(j)))) && any(Element_array(:)==(char(Species_name{:}(j+1)))) && (j+1 == length(Species_name{:}));
Molecular_weight_of_species_element(i)=map(char(Species_name{:}(j)))+ map(char(Species_name{:}(j+1)));
%Condition for species ,having element to be checked as real
%number and next element is character and next element is the
%last element in the evaluating species
elseif (str2double(Species_name{:}(j)))>= 1 && any(Element_array(:)==(char(Species_name{:}(j+1)))) && (j+1 == length(Species_name{:}));
Molecular_weight_of_species_element(i)=map(char(Species_name{:}(j+1)));
%Condition for species ,having element to be checked as real
%number and next element is character and next element is not
%last element in the evaluating species
else (str2double(Species_name{:}(j)))>= 1 && (j+1 ~= length(Species_name{:}));
Molecular_weight_of_species_element(i)=0;
end
Molecular_weight_of_species(i)=Molecular_weight_of_species_initial(i)+Molecular_weight_of_species_element(i);
Molecular_weight_of_species_initial=Molecular_weight_of_species;
end
fprintf('<strong>Molecular Weight of the Species- %s is %d\n </strong>',char(Species_name),Molecular_weight_of_species(i))
end
%Middle,Higher and Lower local temp coeff for species extracted
Middle_local_temp_coeff = str2double(A(end-1));
Higher_local_temp_coeff = str2double(A(end-2));
Lower_local_temp_coeff = str2double(A(end-3));
Line_1=fgetl(f1); %First line of species coefficients assigned to Line_1
Line_1_Delimiter_E_index = strfind(Line_1,'E'); %Index of delimiter 'E' in the line calculated
%There are total 14 coefficients for every species
%First 7 high temp coeffs
%Last 7 low temp coeffs
%All the coefficients extracted from the index of 'E'
High_temp_coeff_1 = str2double(Line_1(1:Line_1_Delimiter_E_index(1)+3));
High_temp_coeff_2 = str2double(Line_1(Line_1_Delimiter_E_index(1)+4:Line_1_Delimiter_E_index(2)+3));
High_temp_coeff_3= str2double(Line_1(Line_1_Delimiter_E_index(2)+4:Line_1_Delimiter_E_index(3)+3));
High_temp_coeff_4 = str2double(Line_1(Line_1_Delimiter_E_index(3)+4:Line_1_Delimiter_E_index(4)+3));
High_temp_coeff_5 = str2double(Line_1(Line_1_Delimiter_E_index(4)+4:Line_1_Delimiter_E_index(5)+3));
Line_2=fgetl(f1); %Second line of species coefficients assigned to Line_2
Line_2_Delimiter_E_index = strfind(Line_2,'E');
High_temp_coeff_6 = str2double(Line_2(1:Line_2_Delimiter_E_index(1)+3));
High_temp_coeff_7 = str2double(Line_2(Line_2_Delimiter_E_index(1)+4:Line_2_Delimiter_E_index(2)+3));
Low_temp_coeff_1 = str2double(Line_2(Line_2_Delimiter_E_index(2)+4:Line_2_Delimiter_E_index(3)+3));
Low_temp_coeff_2 = str2double(Line_2(Line_2_Delimiter_E_index(3)+4:Line_2_Delimiter_E_index(4)+3));
Low_temp_coeff_3 = str2double(Line_2(Line_2_Delimiter_E_index(4)+4:Line_2_Delimiter_E_index(5)+3));
Line_3=fgetl(f1); %Third line of species coefficients assigned to Line_3
Line_3_Delimiter_E_index = strfind(Line_3,'E');
Low_temp_coeff_4 = str2double(Line_3(1:Line_3_Delimiter_E_index(1)+3));
Low_temp_coeff_5 = str2double(Line_3(Line_3_Delimiter_E_index(1)+4:Line_3_Delimiter_E_index(2)+3));
Low_temp_coeff_6 = str2double(Line_3(Line_3_Delimiter_E_index(2)+4:Line_3_Delimiter_E_index(3)+3));
Low_temp_coeff_7 = str2double(Line_3(Line_3_Delimiter_E_index(3)+4:Line_3_Delimiter_E_index(4)+3));
%For Temperature between 200 to 999 Kelvin
%Low temeprature coefficients will be used to calculate
%Specific heat capacity(Cp),Enthalpy(H)and Entropy(S)
%For Temperature more than 1000 kelvin
%High temperature coefficients will be used.
for k = 1:length(T)
if (200 <= T(k)) && (T(k) < 1000)
a1 = Low_temp_coeff_1;
a2 = Low_temp_coeff_2;
a3 = Low_temp_coeff_3;
a4 = Low_temp_coeff_4;
a5 = Low_temp_coeff_5;
a6 = Low_temp_coeff_6;
a7 = Low_temp_coeff_7;
else
a1 = High_temp_coeff_1;
a2 = High_temp_coeff_2;
a3 = High_temp_coeff_3;
a4 = High_temp_coeff_4;
a5 = High_temp_coeff_5;
a6 = High_temp_coeff_6;
a7 = High_temp_coeff_7;
end
Cp(k)= (a1 + a2*T(k) + a3*T(k)^2 + a4*T(k)^3 + a5*T(k)^4)*R;
H(k) = (a1 + a2*T(k)/2 + a3*T(k)^2/3 + a4*T(k)^3/4 + a5*T(k)^4/5 + a6/T(k))*R*T(k);
S(k) = (a1*log(T(k)) + a2*T(k) + a3*T(k)^2/2 + a4*T(k)^3/3 + a5*T(k)^4/4 + a7)*R;
end
%Change directory to data_for_species from present working directory
cd data_for_species
%Create directory with species name
mkdir(sprintf('%s',char(Species_name)));
%Assign folder name to species_folder
species_folder=sprintf('%s',char(Species_name));
%Change directory to species_folder
cd(species_folder);
%Array storing plot file name
plot_array=["Cp vs Temp.png","H vs Temp.png","S vs Temp.png"];
%Array storing ylabels
plot_ylabel_array=["Cp(J/kg-K)","H(J)","S(J/K)"];
%Array storing titles
plot_titles_array=["Plot between Cp & T","Plot between H & T","Plot between S & T"];
%Array storing Cp,H,S value
Parameters_array=[Cp;H;S];
for m=1:3
f=figure('visible','off');
plot(T,Parameters_array(m,:));
xlabel('Temperature (kelvin)');
ylabel(plot_ylabel_array(m));
title(plot_titles_array(m));
filename=sprintf('%s',plot_array(m));
fullfile_name=fullfile(pwd,filename);
saveas(f,fullfile_name);
end
cd ..
cd ..
end
fclose(f1);
-------------------------------------------------***********************----------------------------------------
#Detailed report for the program:
1.Number of species in file defined that is 53.
2.Temperature array 'T' defined by linspace.
3.Using 'fopen' command,file 'THERMO.dat' is opened for reading.
4.First 5 lines are not so important ,3 are header files and one global temperature range file.
5.For passing molecular weight to element in element array ,map variable is used to assign values wherever needed.
6.Element array defined which will be used in calculating molecular weight.
7.Make a folder 'data_for_species' using mkdir.
8.A 'for' loop is used for calculating molecular weight.
9.Sixth line is extracted from where species data is starting and is assigned to X variable for every loop.
10.Using 'strsplit' split line 'X' where 'space' is present and is assigned to A.
11.First element in A is species,and is assigned to 'All_Species_name',which is an array containing species name for all species.
12.Variable 'Species_name' is defined which contain species name for particular species for every loop.
13.'Molecular_weight_of_species_initial' defined which is an array of molecular weight initially which is zero.
14.Different conditions are defined using 'if' and 'else' statement.
a).First 'if statement' is for species having single element in the thermo.dat file.
SO whichever species having single element,its molecular weight will be calculated by using 'map' variable.
And that molecular weight is assigned to array 'Molecular_weight_of_species'.
b).In case species having more than one element ,'else' statement will work.
Under 'else' condition a for loop is used for calculating molecular weight.
Different condition under this 'else' statement are as follows:
1).First 'if'--This condition will follow when species element which is being analysed is 'character' and next element to this in that same species is 'real number'.Example:'H2'.-when 'H' will be checked it is a character and next element whcih is '2' is a real number.
2).Next 'elseif'--This condition will follow when species element which is being analysed is 'character' and next element to this in that same species is also character and that next element is not last element in that species..Example:'HO2'-when 'H' will be checked it is a character and next element 'O' which is also a character and 'O' is not last element in the species.
3).Second 'elseif'--This condition will follow when species element which is being analysed is 'character' and next element to this in that same species is also character and that next element is last element in that species.Example:'OH'-when 'O' will be checked,it is a character and next element 'H' which is also a character and it is also a last element.
4).Third 'elseif'--This condition will follow when species element which is being analysed is 'real number' and next element to this in that same species is character and next element is also last element.Example:'H2O'-when '2' will be checked it is a real number and next element 'O' which is a character and is the last element .
5).'else'--This condition will follow when species element which is being analysed is 'real number'
and next element is not last element.Example:'CH2CHO'.
With the help of an example,lets see how Molecular weight is being calculated;
Take 'CH2CHO':as species length greater than 1 so else statement inside for loop for evry species will work.
1).When first element 'C' will be checked it is a character and next element 'H' which is also a character and 'H' is not last element,so first 'elseif' will be satisfied.
So,Molecular_weight_of_species_element=Molecular weight of 'C'.
Since,
Molecular_weight_of_species=Molecular_weight_of_species_initial+Molecular_weight_of_species_element
Molecular_weight_of_species_initial=0(initially)
Therefore,
Molecular_weight_of_species=Molecular weight of 'C'=12
Then,initial molecular weight updated;
Molecular_weight_of_species_initial=Molecular_weight_of_species=12
2).Again
when 'H' will be checked it is a character and next element '2' is a real number .So 'if condition is satisfied.
Molecular_weight_of_species_element=Molecular weight of 'H'*2=2
Therefore,
Molecular_weight_of_species=Molecular_weight_of_species_initial+Molecular_weight_of_species_element
Molecular_weight_of_species=12+2=14
Now updating initial molecular weight,
Molecular_weight_of_species_initial=Molecular_weight_of_species=14.
3).Now when '2' will be checked ,else condition will be satisfied.
therefore,Molecular_weight_of_species_element=0
Therefore,
Molecular_weight_of_species=Molecular_weight_of_species_initial+Molecular_weight_of_species_element
Molecular_weight_of_species=14+0=14
Now updating initial molecular weight,
Molecular_weight_of_species_initial=Molecular_weight_of_species=14.
4).Then,'C' will be checked which is a character and next element 'H' which is also a character and 'H is also not the last element.So first 'elseif' will be satisfied.
Molecular_weight_of_species_element=Molecular weight of 'C'=12
Therefore,
Molecular_weight_of_species=Molecular_weight_of_species_initial+Molecular_weight_of_species_element
Molecular_weight_of_species=14+12=26
Now updating initial molecular weight,
Molecular_weight_of_species_initial=Molecular_weight_of_species=26.
5).When 'H' will be checked which is character and next element 'O' which is also a character and which is also a last character.
Therefore,second 'elseif' will be satisfied .
Molecular_weight_of_species_element=Molecular weight of 'H'+Molecular weight of 'O'
Molecular_weight_of_species_element=1+16=17
Therefore,
Molecular_weight_of_species=Molecular_weight_of_species_initial+Molecular_weight_of_species_element
Molecular_weight_of_species=26+17=43(final molecular weight)
Molecular weight of "CH2CHO"=43
---------------------------------------------***************************------------------------------
#From the first line in every species ,3 local temperature coefficient range will be extracted.Low,Medium and High.
#Below that line ,to 2 next lines coefficients are present,first 2 line having 5 coefficients and last line having 4 coefficient making total 14 coefficients for every species.
#Like this all coefficients are extracted for every species ,first 7 are high and last 7 are low temperature coefficients
#For calculating Cp,H and S ,if temperature range is between 200 to 999 K,low temperature coeff will be used and above 999 high temperature coeff will be used.
#Then enter in directory ' data_for_species' by using 'cd data_for_species'.
#Then create folder by species name.
#Enter in that folder using 'cd(species_folder)'
#Inside that species folder plot 3 graphs by plot title .
1.Cp vs T.
2.H vs T
3. S vs T.
#For every species each species folder will be created and inside that folder 3 plots will be created.
#Wherever the program will be placed,species folder will be created in that directory.
-------------------------------------**************************************-------------------------------
#Screenshot of the window where all the folders are saved.
----------------------------------------*******************************-------------------------------------
#Plots for o2, N2 and Co2 species:
1.O2:
-----------------------------------------------******************************--------------------------------
2.N2:
---------------------------------------------***********************------------------------------------
3.CO2:
-----------------------------------------***************************----------------------------------