Project 1 - Parsing NASA thermodynamic data

MAIN CODE:-

clear all
close all
clc

% Opening of NASA's Thermodynamic data

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

% Obtaining of global max, min and medium tempratures 

global_temp = fgetl(f1);
temp_range = strsplit(global_temp,' ');
temp_min = str2num(temp_range{2})
temp_med = str2num(temp_range{3})
temp_max = str2num(temp_range{4})

% Range of temprature
T = linspace(temp_min,temp_max,1000);
R = 8.314;

% Skipping the next three comments in the file 
for i = 1:3
    comments = fgetl(f1);
end

% Getting the coeffecient of species

for i = 1:53
% Extracting line 1 and the name of the species 
line1 = fgetl(f1);
name = strtok(line1);

% Extracting higher and lower temprature coeffecient 

line2 = fgetl(f1);
line3 = fgetl(f1);
line4 = fgetl(f1);

A = strfind(line2,'E');
B = strfind(line3,'E');
C = strfind(line4,'E');

% Extracting higher temprature coefficient
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:B(2)+3));

% Extracting lower temprature coeffecient
b1 = str2num(line3(B(2)+4:B(3)+3));
b2 = str2num(line3(B(3)+4:B(4)+3));
b3 = str2num(line3(B(4)+4:B(5)+3));
b4 = str2num(line4(1:C(1)+3));
b5 = str2num(line4(C(1)+4:C(2)+3));
b6 = str2num(line4(C(2)+4:C(3)+3));
b7 = str2num(line4(C(3)+4:C(4)+3));


% Specific heat calculation

cp = specific_heat(a1,a2,a3,a4,a5,b1,b2,b3,b4,b5,T,R,temp_med);

% Enthalpy calculation

h = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,temp_med);

% Entropy calculation

s = entropy_calc(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,temp_med);

% Calculatinng molecular weight of species

Molecular_weight(i) = molecular_weight(name)

% Creating idividual folder for each species and plotting
plotting(cp,h,s,T,name)
end

MOLECULAR WEIGHT:-

% Fuction to calculate molecular weight
function MW = molecular_weight(name)
atoms = ['H','C','O','N','R'];
fprintf('Molecuar 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 = str2num(name(1));
    if n > 1
        MW = MW + atomic_weight(position)*(n-1);
    end
end
fprintf('%d',MW);
disp(' ');
end

PLOTTING:-

function plotting (cp,h,s,T,name)
cd(pwd)
dir_name = name;
mkdir(dir_name);

figure(1)
plot(T,cp,'LineWidth',2,'Color','b')
title(sprintf('Specific heat Vs Temprature for %s',name))
xlabel('TEMPRATURE')
ylabel('SPECIFIC HEAT')

figure(2)
plot(T,h,'LineWidth',2,'Color','r')
title(sprintf('Enthalpy Vs Temprature for %s',name))
xlabel('TEMPRATURE')
ylabel('ENTHALPY')

figure(3)
plot(T,s,'LineWidth',2,'Color','g')
title(sprintf('Entropy Vs Temprature for %s',name))
xlabel('TEMPRATURE')
ylabel('ENTROPY')

cd(dir_name)

saveas(1,'cp_vs_t.png')
saveas(2,'h_vs_t.png')
saveas(3,'s_vs_t.png')

cd ..
end

ENTROPY:-

%Calculation of Etropy
function s = entropy_calc(a1,a2,a3,a4,a5,a7,b1,b2,b3,b4,b5,b7,T,R,temp_med)
if T <= temp_med
    s = R*(b1.*log(T)+a2.*T+a3.*(T.^2)./2+a4.*(T.^3)./3+a5.*(T.^4)./4+a7);
else
    s = R*(b1.*log(T)+b2.*T+b3.*(T.^2)./2+b4.*(T.^3)./3+b5.*(T.^4)./4+b7);
end

ENTHALPY:-

% Calculating Enthalpy
function h = enthalpy(a1,a2,a3,a4,a5,a6,b1,b2,b3,b4,b5,b6,T,R,temp_med)
if T >= temp_med
    h = R.*(a1.*T+(a2.*T.^2)./2+(a3.*T.^3)./3+(a4.*T.^4)./4+(a5.*T.^5)./5+a6);
else
    h = R.*(b1.*T+(b2.*T.^2)./2+(b3.*T.^3)./3+(b4.*T.^4)./4+(b5.*T.^5)./5+b6);

end

SPECIFIC HEAT:-