Project 1 - Parsing NASA thermodynamic data
%PROJECT 1 : NASA FILE PARSING
clear all
close all
clc
f1 = fopen('THERMO.DAT','r')
%Reading The first line from data
first_line = fgetl(f1)
%Reading The second line from data
global_temp = fgetl(f1)
%Splitting the second line using a blank space as a delimiter
A = strsplit(global_temp,' ')
%Obtain temperature value from the splitting data
global_low_temp = str2num(A{1})
global_mid_temp = str2num(A{2})
global_high_temp = str2num(A{3})
%Creating a vector of 400 evenly spaced temperature points
R = 8.314 %Univrsal gas constent
%Skipping the next 3 lines
for i = 1:3
fgetl(f1)
end
for j = 1:53
l1 = fgetl(f1);
line1 = strsplit(l1,' ');
T_mid = str2double(line1(length(line1)-1));
T_high = str2double(line1(length(line1)-2));
T_low = str2double(line1(length(line1)-3));
T = linspace(T_low,T_high,400);
element_name = strtok(l1);
line2 = fgetl(f1);
X = strfind(line2,'E');
line3 = fgetl(f1);
Y = strfind(line3,'E');
line4 = fgetl(f1);
Z = strfind(line4,'E');
a1 = str2num(line2(1:X(1)+3));
a2 = str2num(line2(X(1)+4:X(2)+3));
a3 = str2num(line2(X(2)+4:X(3)+3));
a4 = str2num(line2(X(3)+4:X(4)+3));
a5 = str2num(line2(X(4)+4:X(5)+3));
a6 = str2num(line3(1:Y(1)+3));
a7 = str2num(line2(Y(1)+4:Y(2)+3));
a8 = str2num(line2(Y(2)+4:Y(3)+3));
a9 = str2num(line2(Y(3)+4:Y(4)+3));
a10 = str2num(line2(Y(4)+4:Y(5)+3));
a11 = str2num(line2(1:Z(1)+3));
a12 = str2num(line2(Z(1)+4:Z(2)+3));
a13 = str2num(line2(Z(2)+4:Z(3)+3));
a14 = str2num(line2(Z(3)+4:Z(4)+3));
cp = zeros(1,400);
% for v = 1:length(T)
%Specific heat
cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,T,R,global_mid_temp);
%enthalpy
h = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,T,R,global_mid_temp);
%entropy
s = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,R,global_mid_temp);
%caluculatting molecular weight of species
molecular_wt(j) = molecular_weight(element_name);
%Saving the data
cd('C:\Users\manoj\TEST')
dir_name = element_name;
mkdir(dir_name);
%Plotting the data
figure(1)
plot(T,cp,'linewidth',2,'color','b')
title(sprintf('Specific heat vs Temperature for %s',element_name))
xlabel('TEMPERATURE')
ylabel('SPECIFIC HEAT')
grid on
figure(2)
plot(T,h,'linewidth',2,'color','r')
title(sprintf('Enthalpy vs Temperature for %s',element_name))
xlabel('TEMPERATURE')
ylabel('ENTHALPY')
grid on
figure(3)
plot(T,s,'linewidth',2,'color','g')
title(sprintf('Entropy vs Temperature for %s',element_name))
xlabel('TEMPERATURE')
ylabel('ENTROPY')
grid on
%Pause the loop so we can see each species plot in real time
pause(0.0001)
cd(dir_name);
saveas(1,'Specific heat.png');
saveas(2,'Enthalpy.png');
saveas(3,'Entropy.png');
end
fclose(f1)
>>The fucntions I have used in the main program to calculate Specific heat(cp),Enthalpy(h),Entropy(s) and molecular weight(MW).
%Function to caluculate specific heat(cp)
function cp = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,T,R,global_mid_temp)
%comparing with global temp
if T > global_mid_temp
cp = R*(a1+a2.*T+a3.*(T).^2+a4.*(T).^3+a5.*(T).^4);
else
cp = R*(a6+a7.*T+a8.*(T).^2+a9.*(T).^3+a10.*(T).^4);
end
end
%Function to caluculate enthalpy(h)
function h = enthalpy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,T,R,global_mid_temp)
if T > global_mid_temp
h = R.*(a1.*T+(a2.*T.^2)./2+(a3.*T.^3)./3+(a4.*T.^4)./4+(a5.*T.^5)./5+a6);
else
h = R.*(a7.*T+(a8.*T.^2)./2+(a9.*T.^3)./3+(a10.*T.^4)./4+(a11.*T.^5)./5+a12);
end
%Function to caluculate entropy(s)
function s = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,T,R,global_mid_temp)
if T < global_mid_temp
s = R*(a1.*log(T)+a2.*T+a3.*(T.^2)+a4.*(T.^3)+a5.*(T.^4)+a6.*(T.^5)+a7);
else
s = R*(a8.*log(T)+a9.*T+a10.*(T.^2)+a11.*(T.^3)+a12.*(T.^4)+a13.*(T.^5)+a14);
end
- Function for molecular weight(MW)
> I have created a function 'molecular_weight'. To calculate molecular wieght.
> This function takes the weight of hydrogen,carbon,nitrogen and argon as inputs. All the elemnts consists of all these atoms.
> I have used 'strcmp' to compare the elements with species. if they both matches then the molecular weight will be caluculated
> Take the string of species name and position to caluculate the weight. Here 'n' is number of molecules of elements present in the species.
The function commands are given below :
%Function to caluculate molecular weight
function MW = molecular_weight(element_name)
atoms = ['H','C','N','O','A'];
fprintf('molecular_weight of');
fprintf(element_name);
fprintf(' : ');
atomic_weight = [1 12 14 15 39];
MW = 0;
for j = 1:length(element_name)
for i = 1:length(atoms)
if strcmp(element_name(j),atoms(i))
MW = MW + atomic_weight(i)
position = i;
end
end
n = str2num(element_name(j));
if n > 1
MW = MW + atomic_weight(position)*(n-1);
end
end
fprintf('%d',MW);
disp(' ');
end
>>Molecular weights (MW)
- I have edited the molecular weights data to be clearly visible.(I have not changed the molecular weights of the species)
molecular_weight ofO : 15
molecular_weight ofO2 : 30
molecular_weight ofH : 1
molecular_weight ofH2 : 2
molecular_weight ofOH : 16
molecular_weight ofH2O : 17
molecular_weight ofHO2 : 31
molecular_weight ofH2O2 : 32
molecular_weight ofC : 12
molecular_weight ofCH : 13
molecular_weight ofCH2 : 14
molecular_weight ofCH2(S) : 14
molecular_weight ofCH3 : 15
molecular_weight ofCH4 : 16
molecular_weight ofCO : 27
molecular_weight ofCO2 : 42
molecular_weight ofHCO : 28
molecular_weight ofCH2O : 29
molecular_weight ofCH2OH : 30
molecular_weight ofCH3O : 30
molecular_weight ofCH3OH : 31
molecular_weight ofC2H : 25
molecular_weight ofC2H2 : 26
molecular_weight ofC2H3 : 27
molecular_weight ofC2H4 : 28
molecular_weight ofC2H5 : 29
molecular_weight ofC2H6 : 30
molecular_weight ofCH2CO : 41
molecular_weight ofHCCO : 40
molecular_weight ofHCCOH : 41
molecular_weight ofH2CN : 28
molecular_weight ofHCN : 27
molecular_weight ofHNO : 30
molecular_weight ofN : 14
molecular_weight ofNNH : 29
molecular_weight ofN2O : 43
molecular_weight ofNH : 15
molecular_weight ofNH2 : 16
molecular_weight ofNH3 : 17
molecular_weight ofNO : 29
molecular_weight ofNO2 : 44
molecular_weight ofHCNO : 42
molecular_weight ofHOCN : 42
molecular_weight ofHNCO : 42
molecular_weight ofNCO : 41
molecular_weight ofCN : 26
molecular_weight ofHCNN : 41
molecular_weight ofN2 : 28
molecular_weight ofAR : 39
molecular_weight ofC3H8 : 44
molecular_weight ofC3H7 : 43
molecular_weight ofCH3CHO : 43
molecular_weight ofCH2CHO : 42
>> Refernce links.
https://www.mathworks.com/help/matlab/matlab_prog/parse-function-inputs.html
2. https://www.mathworks.com/help/matlab/ref/fscanf.html
3. https://www.mathworks.com/help/instrument/scanstr.html
4. https://www.mathworks.com/help/matlab/ref/mkdir.html
>> Step wise explanation of the code:
The main objective of the project is to write a function to extract 14 coefficients (The first seven coefficients are high-temperature coefficients and the second seven coefficients are low-temperature coefficients) with the input information of given species and calculate the Molecular weight (MW), Specific heat (Cp), Entropy (s) and Enthalpy (h) for each species with help of NASA thermodynamic polynomial data ‘THERMO’.
I have used ‘fopen’ command to open NASA’s experimental data (THERMO).
The first line consists text (THERMO), ‘fgetl’ command is used to read the first line in the data.
The second line consist of three temperature values, here also ‘fgetl’ command is used to read the second line and stored the values in the variable called global_temp.
The ‘strsplit’ command is used to split the data using given delimiter. The ‘blank space’ is a delimiter, the temperature range consist of three values in a vector.
Obtained the global Low, mid and high temperatures from the splitted data, here ‘str2num’ command is used to convert character array of temperature into numeric array of the temperature.
Now created a vector of 400 evenly spaced temperature points by using ‘linspace’ command and defined universal gas constant value used in the calculations.
I have used ‘fgetl’ command in a for loop to skip the third, fourth and fifth lines which are text comments in the data file THERMO, to be read in text form.
From line 6 to 118, there is a set of 4 lines, first line consists of the name of the species and second, third and fourth lines will have experimental data.
Line 1, I have used ‘fgetl’ command to read line 6 and stored the line 6 data in the “l1” variable.
‘strtok’ command is used to get the selected part of the string, used to select the specie which have ‘O’ and store the respective string (Specie) in line1 variable.
I have used ‘fgetl’ command to read line 7 and stored the data in variable line2
To find the places of ‘E’ in line 7, ‘strfind’ command is used and stored respective species in an array named ‘X’.
Similarly, I have repeated the same procedure for next two lines (line three and line four) and saved them in the arrays called ‘Y’ and ‘Z’.
With the help of ‘str2num’ command, extracted all the coefficients (a1, a2, a3, a4, a5, a6, a7,a8, a9, a10, a11, a12, a13 and a14).
This entire process needs to repeat from line 6 to line 118 for all 53 species, a ‘for loop’ is used for above process.
By calling the functions [specific heat (cp), enthalpy (h), entropy (s) and molecular weight (MW)] into the main code to find out the values of cp, h, s and MW for each species.
‘plot’ command is used to plot the
Temperature and specific heat
Temperature and entropy
Temperature and enthalpy for all the species
Used ‘cd’ command to create a new folder for every iteration of loop and it creates a folder with the species name and pasted its directory in cd
‘mkdir’ command is used to provide name for the folder which is created in the earlier step
The below commands are used to save all the plots in a separate folder
cd(C:\Users\manoj\TEST)
dir_name = element_name;
mkdir(dir_name)
‘saveas’ command is used to save the plots in ‘png’ format, the commands are given below
cd(dir_name);
saveas(1,'Specific heat.png');
saveas(2,'Enthalpy.png');
saveas(3,'Entropy.png');
At the end, closed the NASA’s thermodynamic experimental data file (THERMO) with the help of ‘fclose’ command.