Parsing NASA thermodynamic data
AIM:
To write a code that extracts the NASA thermodynamic data from the file and also calculating the molecular weight of the species present in the file.
OBJECTIVES:
- Displaying the molecular weight of the species
- Creating plot for specific heat, enthalpy,entropy againg temperature ranges and saving them in a separate folders having folder name same as species names automatically with the code.
PROGRAM:
% Parsing NASA Thermodynamic data
clear all
close all
clc
% Universal gas constant
R = 8.314;
% Loading the data file
Thermo_data = fopen('THERMO.dat','r'); % fopen opens the THERMO.dat file and r syntax is for reading it
% Reading the header data
header = fgetl(Thermo_data); % Here fgetl reads the first line data and stores in header
% Reading the temperature data
temp = fgetl(Thermo_data); % Here fgel reads the second line data and stores in temp
A = strsplit(temp,' '); % strsplit is used to spilt the string data in temp with ' ' as a delimiter and stores in A
Global_low_temp = str2num(A{2}); % Here str2num is used to convert numerical string data in second element of A and converts to numerical data
Global_mid_temp = str2num(A{3}); % Here str2num is used to convert numerical string data in third element of A and converts to numerical data
Global_high_temp = str2num(A{4}); % Here str2num is used to convert numerical string data in fourth element of A and converts to numerical data
% Reading the comments
for i = 1:3
comments = fgetl(Thermo_data); % Here fgel is running 3 times using for loop and stores in comments
end
% Reading the species data
for i = 1:53
species_header = fgetl(Thermo_data); % Here fgel reads the next line data and stores in species_header
s = strsplit(species_header,' '); % strsplit is used to spilt the string data with ' ' as a delimiter and stores in s
species_name(i) = s(1); % Assigning first cell element of s to species_name and as the loop runs it stores the data
% Reading species name
species = strtok(species_header); % Here strtok extract the string element by ' ' as a delimiter and stores in species
Atoms_name = {'H' 'O' 'C' 'N' 'AR'}; % Assign string element of Atoms_name
Weight = [1,16,12,14,40]; % Assign molecular weight of above Atoms_name
molecular_weight = 0; % Assuming initial molecular weight as zero
for l = 1:length(species) % Running for loop from first element to last element of species
for m = 1:length(Atoms_name) % Running for loop from first element to last element of Atoms_name
if strcmp(species(l) , Atoms_name(m)) % Comparing string elements of lth element of species with mth element of Atoms_name
molecular_weight = molecular_weight + Weight(m); % If above logic is true then mth element of Weight is added with initial element of molecular weight and it is assigned to molecular weight
y = m; % If the above if stateent is true then the m value is stored in y
end
end
n = str2num(species(l)); % Here str2num is used to convert numerical string data in lth element of species if present and converts to numerical data
if n>1 % If the numerical element is present and is greater than 2 the below statement is done
molecular_weight = molecular_weight + Weight(y)*(n-1); % Here y th element of Weight is multiplied with (n-1) and is added to the above new molecular weight and this all is added to new nmolcular weight
end
end
fprintf('Molecular weight of ') % fprintf prints the string data
disp(species) % disp displays the species name that iteration
fprintf('is ')
disp(molecular_weight)
% Reading local temperature ranges
llt(i) = str2double(species_header(49:55)); % str2double converts the string data to double data type it works only if the string data as numerical value
lht(i) = str2double(species_header(58:65));
lmt(i) = str2double(species_header(68:75));
% Reading the coefficients of each species
species_data1 = fgetl(Thermo_data);
s1 =strfind(species_data1,'E'); % strfind finds the position of mentioned string element in the data and stores them in s1
a1(i) = str2num(species_data1(1:s1(1)+3)); % Assigns the data from first to first position number in s1 + 3 after using str2num to it
a2(i) = str2num(species_data1(s1(1)+4:s1(2)+3)); % Assigns the data from first position number of s1 + 4 to second position number in s1 + 3 after using str2num to it
a3(i) = str2num(species_data1(s1(2)+4:s1(3)+3)); % Assigns the data from second position number of s1 + 4 to third position number in s1 + 3 after using str2num to it
a4(i) = str2num(species_data1(s1(3)+4:s1(4)+3)); % Assigns the data from third position of s1 + 4 to fourth position number in s1 + 3 after using str2num to it
a5(i) = str2num(species_data1(s1(4)+4:s1(5)+3)); % Assigns the data from fourth position of s1 + 4 to fifth position number in s1 + 3 after using str2num to it
species_data2 = fgetl(Thermo_data);
s2 =strfind(species_data2,'E');
a6(i) = str2num(species_data2(1:s2(1)+3));
a7(i) = str2num(species_data2(s2(1)+4:s2(2)+3));
a8(i) = str2num(species_data2(s2(2)+4:s2(3)+3));
a9(i) = str2num(species_data2(s2(3)+4:s2(4)+3));
a10(i) = str2num(species_data2(s2(4)+4:s2(5)+3));
species_data3 = fgetl(Thermo_data);
s3 =strfind(species_data3,'E');
a11(i) = str2num(species_data3(1:s3(1)+3));
a12(i) = str2num(species_data3(s3(1)+4:s3(2)+3));
a13(i) = str2num(species_data3(s3(2)+4:s3(3)+3));
a14(i) = str2num(species_data3(s3(3)+4:s3(4)+3));
end
% Calculating Specific heat, Enthalpy, Entropy of the species for local temperature ranges
for j = 1:53 % Using for loop for calculating the cp,h,s
temp_range = linspace(llt(j),lht(j),100); % Creating 100 values having equal intervals between jth element of llt and jth lement of lht using linspace
for k = 1:length(temp_range)
if temp_range(k)>lmt(j) % If the kth elemnt of temp_range is greater than jth element of lmt then below code runs
cp(k) = (a1(j) + a2(j)*temp_range(k) + a3(j)*(temp_range(k)^2) + a4(j)*(temp_range(k)^3) + a5(j)*(temp_range(k)^4))*R;
h(k) = (a1(j) + a2(j)*temp_range(k)/2 + a3(j)*(temp_range(k)^2)/3 + a4(j)*(temp_range(k)^3)/4 + a5(j)*(temp_range(k)^4)/5 + a6(j)/temp_range(k))*R*temp_range(k);
S(k) = (a1(j)*log(temp_range(k)) + a2(j)*temp_range(k) + a3(j)*(temp_range(k)^2)/2 + a4(j)*(temp_range(k)^3)/3 + a5(j)*(temp_range(k)^4)/4 + a7(j))*R;
else
cp(k) = (a8(j) + a9(j)*temp_range(k) + a10(j)*(temp_range(k)^2) + a11(j)*(temp_range(k)^3) + a12(j)*(temp_range(k)^4))*R;
h(k) = (a8(j) + a9(j)*temp_range(k)/2 + a10(j)*(temp_range(k)^2)/3 + a11(j)*(temp_range(k)^3)/4 + a12(j)*(temp_range(k)^4)/5 + a13(j)/temp_range(k))*R*temp_range(k);
S(k) = (a8(j)*log(temp_range(k)) + a9(j)*temp_range(k) + a10(j)*(temp_range(k)^2)/2 + a11(j)*(temp_range(k)^3)/3 + a12(j)*(temp_range(k)^4)/4 + a14(j))*R;
end
end
% Creating separate folders for each species
Species_name = char(species_name(j)); % char is used to create charater array of jth element in species_name
current_directory = pwd; % pwd identifies the current folder MATLAB files are storing
mkdir(Species_name); % Creating new folder with the above charater array as its name uisng mkdir
cd(Species_name); % Making the above folder as curent folder for MATLAB using cd
TvsCp(j) = figure('Visible','off'); % Storing the below plot and using visible off for not appearing the plot on he screen
plot(temp_range,cp,'-','linewidth',5,'color','r');
title(species_name(j));
xlabel('Temperature(K)');
ylabel('Specifc heat');
saveas(TvsCp(j), 'T_vs_Cp.jpeg'); % Saving the above plot to the current folder as the specified format
TvsH(j) = figure('Visible','off');
plot(temp_range,h,'-','linewidth',5,'color','b');
title(species_name(j));
xlabel('Temperature(K)');
ylabel('Enthalpy(H)');
saveas(TvsH(j), 'T_vs_H.jpeg');
TvsS(j) = figure('Visible','off');
plot(temp_range,S,'-','linewidth',5,'color','k');
title(species_name(j));
xlabel('Temperature(K)');
ylabel('Entropy(S)');
saveas(TvsS(j), 'T_vs_S.jpeg');
cd(current_directory); % Redirecting to the first directory folder at the end of the loop
end
% Plotting for O2,N2,CO2
% For O2
figure(TvsCp(2)) % Shows the plot figure in the TvsCp
figure(TvsH(2))
figure(TvsS(2))
% For CO2
figure(TvsCp(16))
figure(TvsH(16))
figure(TvsS(16))
% For N2
figure(TvsCp(48))
figure(TvsH(48))
figure(TvsS(48))
PROCEDURE:
- Defining the input data universal gas constant
- Using the fopen to open the THERMO.dat file
- Using fgetl to read the data line by line and it reads only a single line
- Reading the global temperature ranges line using fgel by first spliting them using strsplit and then converting them to numbers using str2num
- Using fgel 3 times for skipping the comments in data file
- For reading the species data we are using for loop as the species data is in 53 lines we are using for loop for 53 times
- For the first loop run we will do the following:
- Reading the first line a data using fgetl and assigning it to the variable species_header
- Splitting the first line with space as delimiter and storing it in the s variable using strsplit and it is used later
- Assigning the first element of s variable as it is the species name to the variable named species_name
- For calculating molecular weight we use the strtok to split the data in the species_header instead of strsplit as the data containing numerical value and have to be extrated later for calculation to species variable
- Assigning the individual atoms present in the data to Atoms_name
- Assigning the molecular weight of above atoms to Weight
- Assuming initial value of moecular weight as zero
- Using for loop 2 times first for loop with the length of species( i.e., 6) and other for the length of Atoms_name( i.e., 5)
- The logic is that first it take the first species and the comparing each element of the first species with each element of the Atoms_name using strcmp. If the if statement is true then the element number of Atoms_name at which it is true at that element number of Weight is added the initial assumed molecular weight and this new one is added to molecular_weight variable.
- Next for every lth element of species str2num is used. If it is more than 2 then the next code is run
- fprintf and dis is used to display the character and molecular weight data on to the screen
- Likewise for the species it runs automatically.
-
Next in the same loop for extracting the coefficients that are in the next lines after the species name and calculation and plotting of graph all the logic is mentioned in the code
- Formula for calculating the cp,h,s
- `(Cp)/R = a1 + a2 *T + a3 *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`
- Next specific folders are created and the plots are stored in that
- Finally the plot for O2,N2,CO2 are displayed
RESULTS:
1. Folder creation images
2. Species name and the molecular weight displaying
3. Plots for O2,CO2,N2
REFERENCES:
- https://www.mathworks.com/help/matlab/ref/fgetl.html
- https://www.mathworks.com/help/matlab/ref/strsplit.html
- https://www.mathworks.com/help/matlab/ref/str2num.html
- http://matlab.izmiran.ru/help/techdoc/ref/strtok.html
- http://www.ece.northwestern.edu/local-apps/matlabhelp/techdoc/ref/fopen.html
- http://matlab.izmiran.ru/help/techdoc/ref/pwd.html
- https://www.mathworks.com/help/matlab/ref/cd.html
- https://www.mathworks.com/help/matlab/ref/mkdir.html
- https://www.mathworks.com/help/matlab/ref/saveas.html