Project 1 - Parsing NASA thermodynamic data
clear all
close all
clc
%Data file opening (NASA thermodynamic data)
f = fopen('THERMO.dat','r');
f_line = fgetl(f);
% Global high , low and medium temperatures
global_temperature = fgetl(f);
temperature_range = strsplit(global_temperature, ' ');
low_temp = str2num(temperature_range {1});
mid_temp = str2num(temperature_range{2});
high_temp = str2num(temperature_range {3});
T =linspace(1 , 5, 10);
R = 8.314;
% Skipping the further 3lines
for i = 3:5
Comments = fgetl(f);
end
% Coefficient of species
for i = 1:53
%Extracting line1 and name of species
line1 = fgetl(f);
name= strtok(line1);
%Extracting higher and lower temperature coefficients
line2 = fgetl(f);
line3 = fgetl(f);
line4 = fgetl(f);
a = strfind(line2 ,'E');
b = strfind(line3 ,'E');
c = strfind(line4 ,'E');
%Extracting higher temperature coeffecients
a1 = str2num(line2(1:(a(1)+3)));
a2 = str2num(line2((a(1)+4) : (a(2)+3)));
a3 =str2num(line2((a(2)+4) : (a(3)+3)));
a4 = str2num(line2((a(3)+4) : (a(4)+3)));
a5 = str2num(line2((a(4)+4) : (a(5)+3)));
a6 = str2num(line3(1:(b(1)+3)));
a7 = str2num(line3((b(1)+4) : (a(2)+3)));
%Extracting lower temperature coefficients
a8 = str2num(line3((b(2)+4) : (a(3)+3)));
a9 = str2num(line3((b(3)+4) : (a(4)+3)));
a10 = str2num(line3((b(4)+4) : (a(5)+3)));
a11 = str2num(line3((1:c(1)+3)));
a12 = str2num(line3((c(1)+4) : (a(2)+3)));
a13= str2num(line3((c(2)+4) : (a(3)+3)));
a14= str2num(line3((c(3)+4) : (a(4)+3)));
%cp (specific heat calculation)
cp = specific_heat(a1, a2 , a3 , a4 , a5, a8 , a9 , a10 , a11 , a12 , T,R ,mid_temp);
%h (Enthalpy calculation)
h = enthalpy_final(a1, a2 , a3 , a4 , a5, a6 , a8 , a9 , a10 , a11 , a12 , a13, T, R , mid_temp);
%s (Entropy calculation)
s = Entropy(a1, a2 , a3 , a4 , a5,a7, a8 , a9 , a10 , a11 , a12 , a14, T, R , mid_temp);
%Molecular weight of species calculation
Molecular_wt(i) = molecular_weight(name);
%Plotting
plotting(cp,h,s,T,name)
end
% Function calculating molecular weight
function MW = molecular_weight(name)
atoms = ['H','C','O','N','R'];
fprintf ('Molecular weight of ');
fprintf(name);
fprintf(' : ');
atomic_weight = [1 12 16 14 40];
MW = 0;
for i = 1 : length (name)
for j = 1 : length (atoms)
if strcmp(name(i) , atoms(j))
MW = MW + atomic_weight(j);
position = j;
end
end
n = str2double(name(i));
if n >1
MW = MW + atomic_weight(position)*(n-1);
end
end
fprintf('%d' , MW);
disp(' ');
end
%% Calculation of entropy
function s = Entropy (a1, a2 , a3 , a4 , a5, a7, a8 , a9 , a10 , a11 , a12 , a14, T, R , mid_temp);
if T<= mid_temp
s = (a1.*log(T)+ a2.*T+ a3.*(T.^2)./2 +a4.*(T.^3)./3 + a5.*(T.^4)./4 +a7 ).*R;
else
s = (a8.*log(T) + a9*T+ a10.*(T.^2)./2 + a11.*(T.^3)./3+ a12.*(T.^4)./4+a14).*R;
end
%% Calculation of enthalpy
function h = enthalpy_final(a1, a2 , a3 , a4 , a5, a6 , a8 , a9 , a10 , a11 , a12 , a13, T, R , mid_temp)
if T > mid_temp
h = (a1.*T + (a2.*T.^2)./2 + (a3.*T.^3)./3 +(a4.*T.^4)./4 + (a5.*T.^5)./5 +a6./T).*R.*T;
else
h = (a8.*T + (a9.*T.^2)./2 + (a10.*T.^3)./3 + (a11.*T.^4)./4+ (a12.*T.^5)./5+a13./T).*R.*T;
end
% Calculation of Specific heat
function cp = specific_heat(a1, a2 , a3 , a4 , a5, a8 , a9 , a10 , a11 , a12 , T ,R ,mid_temp)
if T > mid_temp
cp = (a1 + a2.*T + a3.*T.^2 +a4.*T.^3 + a5.*T.^4)*R;
else
cp = (a8 + a9.*T + a10.*T.^2 + a11.*T.^3 + a12.*T.^4)*R;
end
function plotting(cp,h,s,T,name)
nasa_file = name;
mkdir (nasa_file);
figure(1)
plot(T, cp, 'linewidth',2, 'color','r');
title ( sprintf('Temperature vs Specific heat, name'))
xlabel ('Temperature');
ylabel ('Specific heat');
figure(2)
plot(T, h ,'linewidth', 2 , 'color','b');
title( sprintf('Temperature vs Enthalpy for %s',name))
xlabel ('Temperature');
ylabel ('Enthalpy');
figure(3)
plot(T, s ,'linewidth',2, 'color','y');
title( sprintf('Temperature vs Entropy for %s',name))
xlabel ('Temperature');
ylabel ('Entropy');
saveas (1 , 'cp_vs_t.png');
saveas (2 , 'hvs_t.png');
saveas (3 , 'hvs_t.png');
end