Project 1 - Parsing NASA thermodynamic data

clear all; close all; clc;

f1 = fopen('THERMO.dat','r');

% Getting the Global Temperatures
l1 = fgetl(f1);
l2 = fgetl(f1);

Global_temp = strsplit(l2,' ');
LT_global = str2num(Global_temp{2});
AT_global = str2num(Global_temp{3});
HT_global = str2num(Global_temp{4});

% Calculating the Temperature range
T = linspace(LT_global,HT_global,1000);

% Molar Gas Constant from given data
R = 8.314;

% Skipping the next 3 lines
l3 = fgetl(f1);
l4 = fgetl(f1);
l5 = fgetl(f1);

% Getting the local temperature
for i = 1:53
l6 = fgetl(f1);
a = strsplit(l6,' ');
LT_local = str2num(a{end-3});
AT_local = str2num(a{end-1});
HT_local = str2num(a{end-2});

species_name = a{1};

% Getting the Higher & Lower Coefficients
l7 = fgetl(f1);
A = strfind(l7,'E');
a1 = str2num(l7(1:A(1)+3));
a2 = str2num(l7(A(1)+4:A(2)+3));
a3 = str2num(l7(A(2)+4:A(3)+3));
a4 = str2num(l7(A(3)+4:A(4)+3));
a5 = str2num(l7(A(4)+4:A(5)+3));

l8 = fgetl(f1);
A1 = strfind(l8,'E');
a6 = str2num(l8(1:A1(1)+3));
a7 = str2num(l8(A1(1)+4:A1(2)+3));
a8 = str2num(l8(A1(2)+4:A1(3)+3));
a9 = str2num(l8(A1(3)+4:A1(4)+3));
a10 = str2num(l8(A1(4)+4:A1(5)+3));

l9 = fgetl(f1);
A2 = strfind(l9,'E');
a11 = str2num(l9(1:A2(1)+3));
a12 = str2num(l9(A2(1)+4:A2(2)+3));
a13 = str2num(l9(A2(2)+4:A2(3)+3));
a14 = str2num(l9(A2(3)+4:A2(4)+3));

% Calculating Specific Heat (CP), Enthalfy (H) & Entropy (S)
[CP] = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,T,R,AT_global);
[H] = enthalfy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,T,R,AT_global);
[S] = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a14,T,R,AT_global);

% Calculating Molecular Weight
Molecular_weight = wt(species_name);
fprintf('molecular_weight of %s = %f', species_name, Molecular_weight)
disp(' ');
% Plotting the CP, H & S WRT Temperature (T)
% Saving the plots in a seperate folder
pwd = 'C:NASA_project_1';
mkdir(species_name);
cd(species_name);

% Plotting CP vs T
figure(1)
plot(T, CP);
title(['CP vs T graph for ',species_name]);
xlabel('Temperature T');
ylabel('CP');

% Plotting H vs T
figure(2)
plot(T, H);
title(['H vs T graph for ',species_name]);
xlabel('Temperature T');
ylabel('Enthalfy H');

% Plotting S vs T
figure(3)
plot(T, S);
title(['S vs T graph for ',species_name]);
xlabel('Temperature T');
ylabel('Entropy S');
cd ..
end

FUNCTIONS

function CP = specific_heat(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,T,R,AT_global)

% Calaculating the CP for higher & lower Temperatures
if (T>= AT_global)
CP = R*(a1+a2*T+a3*T.^2+a4*T.^3+a5*T.^4);
else
CP = R*(a8+a9*T+a10*T.^2+a11*T.^3+a12*T.^4);
end
end

function H = enthalfy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,T,R,AT_global)

% Calaculating the H for higher & lower Temperatures
if (T>= AT_global)
H = R*(a1+(a2*T)/2+(a3*T.^2)/2+(a4*T.^3)/2+(a5*T.^4)/2);
else
H = R*(a8+(a9*T)/2+(a10*T.^2)/2+(a11*T.^3)/2+(a12*T.^4)/2);
end
end

function S = entropy(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a14,T,R,AT_global)

% Calaculating the S for higher & lower Temperatures
if (T>= AT_global)
S = R*(a1.*log(T)+a2*T+(a3*T.^2)/2+(a4*T.^3)/3+(a5*T.^4)/4);
else
S = R*(a8.*log(T)+a9*T+(a10*T.^2)/2+(a11*T.^3)/3+(a12*T.^4)/4);
end
end

function Molecular_weight = wt(species_name)

Elements = ['H','C','O','N','A'];
Atomic_weight = [1 12 16 14 40];

Molecular_weight = 0;

for i = 1:length(species_name)
for j = 1:length(Atomic_weight)
if strcmp(species_name(i),Elements(j))
Molecular_weight = Molecular_weight + Atomic_weight(j);
k=j;
end
end
n = str2double(species_name(i));
if (n>1)
Molecular_weight = Molecular_weight + Atomic_weight(k) * (n-1);
end
end
end

https://docs.google.com/document/d/1Xd26ob6r_yFxYpaEJJsbrzFpNOLDo1LeIXMyjcS0ZI8/edit?usp=sharing