Project 1 - Parsing NASA thermodynamic data

Main Program-

%Thermo-NASA file reading program
clear all; close all;clc

%Reading data file
f1=fopen('THERMO.dat','r')

%Skip header line
fgetl(f1)

%global temperatures
Temp=fgetl(f1);

Temp=str2num(Temp);
global_low_temp=Temp(1);
global_mid_temp=Temp(2);
global_high_temp=Temp(3);

%Universal Gas Constant
R=8.314;

%skip next 3 lines
for i=1:3
    fgetl(f1);
end

%Reading species
for j = 1:4


%In first Species block
tline=fgetl(f1);
a = strsplit(tline,' ');

%First species Monoautomic Oxygen
Species= a{1};

%Defining local temperatures
low_temp=str2num(a{end-3});
mid_temp=str2num(a{end-1});
high_temp=str2num(a{end-2});

%Define temperature range for lower and high temperature coeff.
T_low=linspace(low_temp,mid_temp,100);
T_high=linspace(mid_temp,high_temp,100);

%Defining coefficients from 2nd line
tline1=fgetl(f1);
pos1=findstr(tline1,'E');

%High temperature coefficients
a1=str2num(tline1(1:pos1(1)+3));
a2=str2num(tline1(pos1(1)+4:pos1(2)+3));
a3=str2num(tline1(pos1(2)+4:pos1(3)+3));
a4=str2num(tline1(pos1(3)+4:pos1(4)+3));
a5=str2num(tline1(pos1(4)+4:pos1(5)+3));

%Coefficients from 3rd line
tline2=fgetl(f1);
pos2=findstr(tline2,'E');

a6=str2num(tline2(1:pos2(1)+3));
a7=str2num(tline2(pos2(1)+4:pos2(2)+3));

%Low temperature coefficients
a8=str2num(tline2(pos2(2)+4:pos2(3)+3));
a9=str2num(tline2(pos2(3)+4:pos2(4)+3));
a10=str2num(tline2(pos2(4)+4:pos2(5)+3));


%Coefficients from 4th line
tline3=fgetl(f1);
pos3=findstr(tline3,'E');

a11=str2num(tline3(1:pos3(1)+3));
a12=str2num(tline3(pos3(1)+4:pos3(2)+3));
a13=str2num(tline3(pos3(2)+4:pos3(3)+3));
a14=str2num(tline3(pos3(3)+4:pos3(4)+3));

%Cp,H,S calculation for High temp coeff
Cp_H=R*(a1+(a2*T_high)+(a3*(T_high).^2)+(a4*(T_high).^3)+(a5*(T_high).^4)); %KJ
H_H=R*T_high.*(a1+((a2*T_high)/2)+((a3*(T_high).^2)/3)+((a4*(T_high).^3)/4)+((a5*(T_high).^4)/5)+((a6./T_high))); %KJ/mole
S_H=R*(a1.*log(T_high)+(a2*T_high)+((a3*(T_high).^2)/2)+((a4*(T_high).^3)/3)+((a5*(T_high).^4)/4)+a7); %KJ/mole-K

%Cp,H,S calculation for low temp coeff
Cp_L=R*(a8+(a9*T_low)+(a10*(T_low).^2)+(a11*(T_low).^3)+(a12*(T_low).^4)); %KJ
H_L=R*T_low.*(a8+((a9*T_low)/2)+((a10*(T_low).^2)/3)+((a11*(T_low).^3)/4)+((a12*(T_low).^4)/5)+((a13./T_low))); %KJ/mole
S_L=R*(a8.*log(T_low)+(a9*T_low)+((a10*(T_low).^2)/2)+((a11*(T_low).^3)/3)+((a12*(T_low).^4)/4)+a14); %KJ/mole-K
        
%Creating array of temperature,Cp,S,H
T_array=[T_low T_high];
Cp_array=[Cp_L Cp_H];
S_array=[S_L S_H];
H_array=[H_L H_H];
%plot Cp,H,S with temperature 
figure(1);
plot(T_array,Cp_array);
figure(2);
plot(T_array,H_array);
figure(3);
plot(T_array,S_array);
%making new folder Species to save plots
mkdir(Species);
cd(Species); %Change directory
saveas(figure(1),'Specific Heat','jpg');
saveas(figure(2),'Enthalpy','jpg');
saveas(figure(3),'Entropy','jpg');
 cd ..
%W=molecular_weight(Species); %by function
%fprintf('molecular weight %s = %f',Species,W)
end

Explanation-

First of all, we will open the file 'THERMO.dat' and read this file using fopen command.

Then we will read the header line of the Thermo file using the fgetl command and then

read the temperature file and then convert it to string to number by using

the str2num command.

after that, we will skip the next 3 comment lines by using for loop.

As we have to calculate the molecular weight of 53 species hence we are using the next for loop that is j.

Then we will run the loop to find 53 elements in the data file by assigning fgetl command and then splitting it into a cell array to find the species line, species name and temperature.

Then we read the new line. We will find string 'E' to find the coefficient in line using strfind.

Like that, we will find 14 coefficients and convert them to numbers out of which 7 are higher and 7 are lower coefficient

Now we find Cp, H, S values for Low and high coefficients and by using the array we plot each of them using the plot command.

Then make a new folder using mkdir .

saving all plots in this folder using the saveas command.

At last, we will get molecular weight values in the command window with the help of function program 'W'

cp .. is returned to root folder.

Function Code-

function [W] = molecular_weight(Species)

atomic_name = ['H','C','O','N','AR'];
atomic_weight= [1.00784, 12.0107, 15.999, 14.0067, 39.948];
W=0;
for i = 1:length(Species)
    for j = 1:length(atomic_name)
        if strcmp(Species(i),atomic_name(j))
            W=W+atomic_weight(j)
            pos=j
           
        end
    end
    n=str2num(Species(i))
    if n>1
        W=W+atomic_weight(pos)*(n-1)
end
end

Function Explanation-

We create molecular weight function W and assign species input to it.

Then define the main element of all 53 species in array atomic_name - H C O N Ar.

Then we will create an atomic_weight array with values.

Now we define the molecular weight starting range as zero and then create a for loop condition which runs to the length of the species.

and an other nested for loop is created which runs to the length of the string atomic_name

Then we will use if condition to find whether the species contain more than one element if there is more than 1 element that is n>1

then we will increment the atomic weight by the multiplication of position and n-1 to the atomic _weight array.

SCREENSHOTS OF OUTPUTS -