% Parsing NASA Thermodynamic data.
clear all % to clear the workspace window
close all % to close the all previous figure or plot
clc % to clear the command window
% Reading the file using fopen
f = fopen('THERMO.dat','r');
% Reading the header data
for m=1:5
fgetl(f); % for ignoring the comments in the files
end
% reading all line data
for n = 1:53
% Extracting data from species block
line_1 = fgetl(f);
a = strsplit(line_1,' '); % splitting the string
b = length(a);
% To extract species name
species_name=a{1};
% To extracting local temperature
local_low_temp = a{b-3};
local_high_temp = a{b-2};
local_mid_temp = a{b-1};
% Converting local temperature from string to array
local_low_temp = str2num(local_low_temp);
local_high_temp = str2num(local_high_temp);
local_mid_temp = str2num(local_mid_temp);
% Extraction of coefficient line 1
coeff_line_1 = fgetl(f);
% Reading position of E in coeff_line_1 = fget(l)
posE = strfind(coeff_line_1,'E');
a1 = str2num(coeff_line_1(1:posE(1)+3));
a2 = str2num(coeff_line_1(posE(1)+4:posE(2)+3));
a3 = str2num(coeff_line_1(posE(2)+4:posE(3)+3));
a4 = str2num(coeff_line_1(posE(3)+4:posE(4)+3));
a5 = str2num(coeff_line_1(posE(4)+4:posE(5)+3));
% Extraction of coefficient line 2
coeff_line_2 = fgetl(f);
% Reading position of E in coeff_line_2 = fget(2)
posE = strfind(coeff_line_2,'E');
a6 = str2num(coeff_line_2(1:posE(1)+3));
a7 = str2num(coeff_line_2(posE(1)+4:posE(2)+3));
a8 = str2num(coeff_line_2(posE(2)+4:posE(3)+3));
a9 = str2num(coeff_line_2(posE(3)+4:posE(4)+3));
a10 = str2num(coeff_line_2(posE(4)+4:posE(5)+3));
% Extraction of coefficient line 3
coeff_line_3 = fgetl(f);
% Reading positionn of E in coeff_line_3 = fget(3)
posE = strfind(coeff_line_3,'E');
a11 = str2num(coeff_line_3(1:posE(1)+3));
a12 = str2num(coeff_line_3(posE(1)+4:posE(2)+3));
a13 = str2num(coeff_line_3(posE(2)+4:posE(3)+3));
a14= str2num(coeff_line_3(posE(3)+4:posE(4)+3));
% Take input for temp
T = linspace(local_low_temp,local_high_temp,500);
% Calculation of Cp, H, S
R = 8.314;
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)).*R;
else T(i)<=local_mid_temp;
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).*R;
end
end
% Saving the file and creating the folders
cd('E:\Desktop\skill-lync courses\Matlab')
mkdir(species_name)
cd(species_name)
% plotting and saving of the results
% plotting Temperature vs Specific Heat
figure(1)
plot(T,Cp,'linewidth',2,'color','b')
saveas(figure(1),'Varation of Cp wrt temp.jpg')
xlabel('Temperature in kelvin(k)')
ylabel('specific heat(Cp)')
title(sprintf('Variation of speciffic heat wrt to temp for %s',species_name))
grid on
% plotting Temperature vs Enthalpy
figure(2)
plot(T,H,'linewidth',2,'color','r')
saveas(figure(2),'Varation of H wrt temp.jpg')
xlabel('Temperature in kelvin(k)')
ylabel('Enthalpy(H)')
title('Enthalpy plot')
title(sprintf('Variation of Enthalpy wrt to temp for %s',species_name))
grid on
% plotting Temperature vs Entropy
figure(3)
plot(T,S,'linewidth',2,'color','g')
saveas(figure(3),'Varation of S wrt temp.jpg')
xlabel('Temperature in kelvin(k)')
ylabel('Entropy(S)')
title('Entropy plot')
title(sprintf('Variation of Entropy wrt to temp for %s',species_name))
grid on
% Calculating molecular weight.
atomic_name = ['H','C','O','N','A'];
atomic_weight = [1,12,16,14,39];
molecular_weight = 0;
for j = 1:length(species_name)
for k = 1:length(atomic_name)
if strcmp(species_name(j),atomic_name(k))
molecular_weight = molecular_weight + atomic_weight(k);
position = k;
end
end
n = str2num(species_name(j));
if (n>1)
molecular_weight = molecular_weight + atomic_weight(position)*(n-1);
end
end
fprintf('molecular weight of %s is %f',species_name, molecular_weight)
end
OUTPUT:-
molecular weight of O is 16.000000molecular weight of O2 is 32.000000molecular weight of H is 1.000000molecular weight of H2 is 2.000000molecular weight of OH is 17.000000molecular weight of H2O is 18.000000molecular weight of HO2 is 33.000000molecular weight of H2O2 is 34.000000molecular weight of C is 12.000000molecular weight of CH is 13.000000molecular weight of CH2 is 14.000000molecular weight of CH2(S) is 14.000000molecular weight of CH3 is 15.000000molecular weight of CH4 is 16.000000molecular weight of CO is 28.000000molecular weight of CO2 is 44.000000molecular weight of HCO is 29.000000molecular weight of CH2O is 30.000000molecular weight of CH2OH is 31.000000molecular weight of CH3O is 31.000000molecular weight of CH3OH is 32.000000molecular weight of C2H is 25.000000molecular weight of C2H2 is 26.000000molecular weight of C2H3 is 27.000000molecular weight of C2H4 is 28.000000molecular weight of C2H5 is 29.000000molecular weight of C2H6 is 30.000000molecular weight of CH2CO is 42.000000molecular weight of HCCO is 41.000000molecular weight of HCCOH is 42.000000molecular weight of H2CN is 28.000000molecular weight of HCN is 27.000000molecular weight of HNO is 31.000000molecular weight of N is 14.000000molecular weight of NNH is 29.000000molecular weight of N2O is 44.000000molecular weight of NH is 15.000000molecular weight of NH2 is 16.000000molecular weight of NH3 is 17.000000molecular weight of NO is 30.000000molecular weight of NO2 is 46.000000molecular weight of HCNO is 43.000000molecular weight of HOCN is 43.000000molecular weight of HNCO is 43.000000molecular weight of NCO is 42.000000molecular weight of CN is 26.000000molecular weight of HCNN is 41.000000molecular weight of N2 is 28.000000molecular weight of AR is 39.000000molecular weight of C3H8 is 44.000000molecular weight of C3H7 is 43.000000molecular weight of CH3CHO is 44.000000molecular weight of CH2CHO is 43.000000>>