Project 1 - Parsing NASA thermodynamic data
PROGRAM TO OBTAIN THE SPECIFIC HEAT, ENTHALPY & ENTROPY VS TEMPERATURE PLOTS AND MOLECULAR WEIGHT
AIM :
To write a program to extract the global low,mid and high temperature along with the compound name, local low, local mid, local high temperature and all the low and high temperature co-efficients.
THEORY :
The nasa thermodynamic file contains the global and local temperature along with the polynomials for the calculation of Specific heat, Enthalpy, Entropy and Molecular weight for individual compounds (species).
ACTUAL CODE :
1) MAIN PROGRAM
clear all
close all
clc
% File reading and global temperature extraction
f1 = fopen('THERMO.dat','r');
fgetl(f1);
global_temp = fgetl(f1);
a = findstr(global_temp,' ');
g1 = global_temp(1:10);
g2 = global_temp(12:21);
g3 = global_temp(22:30);
global_min = str2double(g1);
global_mid = str2double(g2);
global_max = str2double(g3);
fgetl(f1);
fgetl(f1);
fgetl(f1);
while ~ feof(f1)
% Extracting the local temperatures
tline1 = strsplit(fgetl(f1),' ');
compound_name = tline1{1};
if strcmpi(compound_name,'END')==1
break
else
r = length(tline1);
local_min_temp = str2double(tline1{r-3});
local_mid_temp = str2double(tline1{r-1});
local_max_temp = str2double(tline1{r-2});
% Extracting the coefficients
tline2 = fgetl(f1);
E_pos = strfind(tline2,'E');
a1 = str2double(tline2(1:E_pos(1)+3));
a2 = str2double(tline2(E_pos(1)+4:E_pos(2)+3));
a3 = str2double(tline2(E_pos(2)+4:E_pos(3)+3));
a4 = str2double(tline2(E_pos(3)+4:E_pos(4)+3));
a5 = str2double(tline2(E_pos(4)+4:E_pos(5)+3));
tline3 = fgetl(f1);
a6 = str2double(tline3(1:E_pos(1)+3));
a7 = str2double(tline3(E_pos(1)+4:E_pos(2)+3));
a8 = str2double(tline3(E_pos(2)+4:E_pos(3)+3));
a9 = str2double(tline3(E_pos(3)+4:E_pos(4)+3));
a10 = str2double(tline3(E_pos(4)+4:E_pos(5)+3));
tline4 = fgetl(f1);
a11 = str2double(tline4(1:E_pos(1)+3));
a12 = str2double(tline4(E_pos(1)+4:E_pos(2)+3));
a13 = str2double(tline4(E_pos(2)+4:E_pos(3)+3));
a14 = str2double(tline4(E_pos(3)+4:E_pos(4)+3));
T = linspace(local_min_temp,local_max_temp,200);
cd('C:\Users\Asus\Desktop\MATLAB\Challenge\Thermal')
[specific_heat,enthalpy,entropy] = cp_h_s(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temp,T);
molecular_wt = mol_wt (compound_name);
mkdir(compound_name)
cd(compound_name)
figure(1)
plot(T,specific_heat,'linewidth',2,'color','b')
xlabel('Temperature (K)');
ylabel('Specific heat (Cp)');
title('Temperature vs Specific heat plot')
grid on
saveas(gcf,'Specific heat (Cp) vs Temperature (T) .png')
figure(2)
plot(T,enthalpy,'linewidth',2,'color','r')
xlabel('Temperature (K)');
ylabel('Enthalpy (H)');
title('Temperature vs Enthalpy')
grid on
saveas(gcf,'Enthalpy (H) vs Temperature (T) .png')
figure(3)
plot(T,entropy,'linewidth',2,'color','g')
xlabel('Temperature (K)');
ylabel('Entropy (S)');
title('Temperature vs Entropy plot')
grid on
saveas(gcf,'Entropy (S) vs Temperature (T) .png')
end
end
cd('C:\Users\Asus\Desktop\MATLAB\Challenge\Thermal')
% Inputs
compound_name = input('Enter the Compound name : ','s')
molecular_wt = mol_wt (compound_name);
sprintf('The molecular weight of %s is %d',compound_name, molecular_wt)2) FUNCTION PROGRAM
function [c_p ,H ,S] = cp_h_s(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temp,T)
R = 8.314;
if (T>=local_mid_temp)
c_p = R.*(a1 + a2*T + a3*T.^2 + a4*T.^3 + a5*T.^4);
H = R.*T.*(a1 + a2*(T./2)+ a3*((T.^2)./3)+ a4*((T.^3)./4)+ a5*((T.^4)./5) + a6./T);
S = R.*(a1*(log(T)) + a2*T + a3*(T.^2)./2 + a4*(T.^3)./3 + a5*(T.^4)./4 + a7);
else
c_p = R.*(a8+a9*T + a10*T.^2+a11*T.^3+a12*T.^4);
H = R.*T.*(a8+a9*(T./2)+a10*((T.^2)./3)+a11*((T.^3)./4)+a12*((T.^4)./5) + a13./T);
S = R.*(a8*(log(T)) + a9*T + a10*(T.^2)./2 + a11*(T.^3)./3 + a12*(T.^4)./4 + a14);
end3) FUNCTION PROGRAM - 2
function [molecular_wt] = mol_wt (compound_name)
elements = ['H','C','N','O','A'];
atomic_wt = [1,12,14,16,40];
molecular_wt = 0;
for i = 1:length(compound_name)
for j = 1:length(elements)
if strcmpi(compound_name(i),elements(j))
molecular_wt = molecular_wt + atomic_wt(j);
position = j;
end
end
f = str2double(compound_name(i));
if f>1
molecular_wt = molecular_wt + atomic_wt(position)*(f-1);
end
end
sprintf('The molecular weight of %s is %d',compound_name, molecular_wt);
end
PROCEDURE :
1) The below command is used to read the 'THERMO.dat' file and to 'fgetl(f1)' to access the individual line of the file. The 'findstr' command is used to find the position of the required character which is then used to extract the global min, mid and max temperature. These characters are converted to numbers
% File reading and global temperature extraction
f1 = fopen('THERMO.dat','r');
fgetl(f1);
global_temp = fgetl(f1);
a = findstr(global_temp,' ');
g1 = global_temp(1:10);
g2 = global_temp(12:21);
g3 = global_temp(22:30);
global_min = str2double(g1);
global_mid = str2double(g2);
global_max = str2double(g3);
fgetl(f1);
fgetl(f1);
fgetl(f1);
2) The command below is used to extract the compound name and is compared with end of the file character 'END'. The command 'while ~feof(f1)' is used to read every line in the file till the end.
while ~ feof(f1)
% Extracting the local temperatures
tline1 = strsplit(fgetl(f1),' ');
compound_name = tline1{1};
if strcmpi(compound_name,'END')==1
break
3) The local min, mid and max temperature are extracted according to the length of the line after which the 'fgetl' command is used to get to the next line.
else
r = length(tline1);
local_min_temp = str2double(tline1{r-3});
local_mid_temp = str2double(tline1{r-1});
local_max_temp = str2double(tline1{r-2});
% Extracting the coefficients
tline2 = fgetl(f1);
4) Since each of the co-efficients contains the letter "E" the command "strfind' is used to obtain the positions of the letter "E" for the first set of co-efficients which is similar to the rest of the co-efficients. Accordingly the command "str2double" is used to convert the characters to numerical format.
E_pos = strfind(tline2,'E');
a1 = str2double(tline2(1:E_pos(1)+3));
a2 = str2double(tline2(E_pos(1)+4:E_pos(2)+3));
a3 = str2double(tline2(E_pos(2)+4:E_pos(3)+3));
a4 = str2double(tline2(E_pos(3)+4:E_pos(4)+3));
a5 = str2double(tline2(E_pos(4)+4:E_pos(5)+3));
tline3 = fgetl(f1);
a6 = str2double(tline3(1:E_pos(1)+3));
a7 = str2double(tline3(E_pos(1)+4:E_pos(2)+3));
a8 = str2double(tline3(E_pos(2)+4:E_pos(3)+3));
a9 = str2double(tline3(E_pos(3)+4:E_pos(4)+3));
a10 = str2double(tline3(E_pos(4)+4:E_pos(5)+3));
tline4 = fgetl(f1);
a11 = str2double(tline4(1:E_pos(1)+3));
a12 = str2double(tline4(E_pos(1)+4:E_pos(2)+3));
a13 = str2double(tline4(E_pos(2)+4:E_pos(3)+3));
a14 = str2double(tline4(E_pos(3)+4:E_pos(4)+3));
5) The Specific heat (Cp), Enthalpy (H) and Entropy (S) is calulated for temperature ranging between Local min and local max temperature which is obtained using the linspace command. The directory is changed to the common directory and below functions are executed to calculate the Specific heat, Enthalpy, Entropy and Molecular weight of the compound.
T = linspace(local_min_temp,local_max_temp,200);
cd('C:\Users\Asus\Desktop\MATLAB\Challenge\Thermal')
[specific_heat,enthalpy,entropy] = cp_h_s(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,local_mid_temp,T);
molecular_wt = mol_wt (compound_name);
6) The functions are calculated according the formula mentioned below
where,
Cp - Specific heat of the compound (species)
H - Enthalpy of the compound
S - Entropy of the compound
R - Universal gas constant
T - Actual temperature (K)
a1, a2, a3, a4, a5, a6 and a7 are the high temperature co-efficients
a8, a9, a10, a11, a12, a13 and a14 are the low temperature co-efficients
The low and high temperature co-efficients are determined based on the local mid temperature, if the actual temperature (T) is above the local mid temperature then the first 7 set of co-efficients are used, if less than the local mid temperature then the remaining 7 set of co-efficients are used.
The molecular weight of the individual compound are also calculated using a different function as shown below.
molecular_wt = molecular_wt + atomic_wt(j);
7) The results obtained are used to plot graphs of different parameters for each compound (species) and new folders (directories) are created for individual compounds. The process is repeated for all the different compounds (species) and the results are stored in individual folders.
mkdir(compound_name)
cd(compound_name)
figure(1)
plot(T,specific_heat,'linewidth',2,'color','b')
xlabel('Temperature (K)');
ylabel('Specific heat (Cp)');
title('Temperature vs Specific heat plot')
grid on
saveas(gcf,'Specific heat (Cp) vs Temperature (T) .png')
figure(2)
plot(T,enthalpy,'linewidth',2,'color','r')
xlabel('Temperature (K)');
ylabel('Enthalpy (H)');
title('Temperature vs Enthalpy')
grid on
saveas(gcf,'Enthalpy (H) vs Temperature (T) .png')
figure(3)
plot(T,entropy,'linewidth',2,'color','g')
xlabel('Temperature (K)');
ylabel('Entropy (S)');
title('Temperature vs Entropy plot')
grid on
saveas(gcf,'Entropy (S) vs Temperature (T) .png')
8) The molecular weight of the required compound is calculated and displayed as per the user input by using the below mentioned commands.
cd('C:\Users\Asus\Desktop\MATLAB\Challenge\Thermal')
% Inputs
compound_name = input('Enter the Compound name : ','s')
molecular_wt = mol_wt (compound_name);
sprintf('The molecular weight of %s is %d',compound_name, molecular_wt)
ERRORS :
1) The above error was obtained since a new directory was created to dump the results into the individual folders, after executing results for the first compound (species) the program would try to calculate the functions from the new directory.
CORRECTIVE ACTION :
1) The above error can be solved changing the file directory to the main directory where the main program and the respective function programs are saved as shown below which eliminated the error.
https://drive.google.com/file/d/1WbzFZXFBsym7-u1B_cLd7n1W41W_a0ub/view?usp=sharing