Project 1 - Parsing NASA thermodynamic data

AIM:

1. To write a code for extracting the 14 co-efficients and calculating the Specific Heat, Enthalpy and Entropy for all the species in the THERMO.dat file.

2. Calculate the Molecular weight of each species and display in the command window.

3. Plot Specific Heat , Enthalpy and Entropy using Local temperature range specific for each species.

4. Save the plots as images with correct names and put them in separate folders using program code.

5. Attach screenshot of window where all the folders are saved.

6. Attach the plots for species O2 , N2 and CO2.

Main Program Code :

clear all
close all
clc

% NASA data file Reading
f = fopen('THERMO.dat','r');

for i = 1:5
    fgetl(f);
end

for k = 1 : 53

a = fgetl(f);

% Extraction of local temperature
position_dot = strfind(a,'.');
low_temp = a(position_dot(1)-3:position_dot(1)+3);
low_temp = str2double(low_temp);
high_temp = a(position_dot(2)-4:position_dot(2)+3);
high_temp = str2double(high_temp);
mid_temp = a(position_dot(3)-4:position_dot(3)+3);
mid_temp = str2double(mid_temp);

b = strsplit(a,' ');

species_name = b{1}

% Extraction of co-efficients
coeff = fgetl(f);
position = strfind(coeff,'E');

b1 = str2double(coeff(1:(position(1)+3)));
b2 = str2double(coeff((position(1)+4):(position(2)+3)));
b3 = str2double(coeff((position(2)+4):(position(3)+3)));
b4 = str2double(coeff((position(3)+4):(position(4)+3)));
b5 = str2double(coeff((position(4)+4):(position(5)+3)));

coeff1 = fgetl(f);
position1 = strfind(coeff1,'E');

b6 = str2double(coeff1(1:position1(1)+3));
b7 = str2double(coeff1((position1(1)+4):(position1(2)+3)));
b8 = str2double(coeff1((position1(2)+4):(position1(3)+3)));
b9 = str2double(coeff1((position1(3)+4):(position1(4)+3)));
b10 = str2double(coeff1((position1(4)+4):(position1(5)+3)));

coeff2 = fgetl(f);
position2 = strfind(coeff2,'E');

b11 = str2double(coeff2(1:position2(1)+3));
b12 = str2double(coeff2((position2(1)+4):(position2(2)+3)));
b13 = str2double(coeff2((position2(2)+4):(position2(3)+3)));
b14 = str2double(coeff2((position2(3)+4):(position2(4)+3)));

% Temperature array creation

T = linspace(low_temp,high_temp,500);

% for i = 1:length(T)
%     t = T(i);
%     if (t<mid_temp)
%         ltr(i) = t;
%     else
%         htr1(i) = t;
%     end
% end
%   
% d = length(ltr);
% htr = htr1((d+1):length(T));

% ltr = linspace(low_temp,mid_temp,500);
% htr = linspace(mid_temp,high_temp,500);

% Cp = R*(a1 + a2T + a3T^2 + a4T^3 + a5T^4)
% H  = R*T*(a1 + a2T/2 + a3T^2/3 + a4T^3/4 + a5T^4/5 + a6/T)
% S  = R*(a1*lnT + a2T + a3T^2/2 + a4T^3/3 + a5T^4/4 + a7)

% Universal gas constant , R
R = 8.3;

% % Calculation of Cp , H , S (LTR - Low Temp Region)
% 
% Cp1 = R*(b8 + (b9*ltr) + (b10*(ltr).^2) + (b11*(ltr).^3) + (b12*(ltr).^4));
% H1  = R*ltr.*(b8 + ((b9*ltr)/2) + ((b10*(ltr).^2)/3) + ((b11*(ltr).^3)/4) + ((b12*(ltr).^4)/5) + (b13./ltr));
% S1  = R*(b8.*log(ltr) + (b9*ltr) + ((b10*(ltr).^2)/2) + ((b11*(ltr).^3)/3) + ((b12*(ltr).^4)/4) + b14);
% 
% % Calculation of Cp , H , S (LTR - High Temp Region)
% 
% Cp2 = R*(b1 + (b2*htr) + (b3*(htr).^2) + (b4*(htr).^3) + (b5*(htr).^4));
% H2  = R*htr.*(b1 + ((b2*htr)/2) + ((b3*(htr).^2)/3) + ((b4*(htr).^3)/4) + ((b5*(htr).^4)/5) + (b6./htr));
% S2  = R*(b1.*log(htr) + (b2*htr) + ((b3*(htr).^2)/2) + ((b4*(htr).^3)/3) + ((b5*(htr).^4)/4) + b7);

% if(T<=mid_temp)
%     Cp = R*(b8+(b9*T)+(b10*(T).^2)+(b11*(T).^3)+(b12*(T).^4));
% else
%     Cp = R*(b1+(b2*T)+(b3*(T).^2)+(b4*(T).^3)+(b5*(T).^4));
% end
% 
% if(T<=mid_temp)
%     H = R*T.*(b8 + ((b9*T)/2) + ((b10*(T).^2)/3) + ((b11*(T).^3)/4) + ((b12*(T).^4)/5) + (b13./T));
% else
%     H = R*T.*(b1 + ((b2*T)/2) + ((b3*(T).^2)/3) + ((b4*(T).^3)/4) + ((b5*(T).^4)/5) + (b6./T));
% end
% 
% if(T<=mid_temp)
%     S = R*(b8.*log(T)+(b9*T)+((b10*(T).^2)/2)+((b11*(T).^3)/3)+((b12*(T).^4)/4)+b14);
% else
%     S = R*(b1.*log(T)+(b2*T)+((b3*(T).^2)/2)+((b4*(T).^3)/3)+((b5*(T).^4)/4)+b7); 
% end

if(T<=mid_temp)
    Cp = R*(b8+(b9*T)+(b10*(T).^2)+(b11*(T).^3)+(b12*(T).^4));
    H = R*T.*(b8 + ((b9*T)/2) + ((b10*(T).^2)/3) + ((b11*(T).^3)/4) + ((b12*(T).^4)/5) + (b13./T));
    S = R*(b8.*log(T)+(b9*T)+((b10*(T).^2)/2)+((b11*(T).^3)/3)+((b12*(T).^4)/4)+b14);
else
    Cp = R*(b1+(b2*T)+(b3*(T).^2)+(b4*(T).^3)+(b5*(T).^4));
    H = R*T.*(b1 + ((b2*T)/2) + ((b3*(T).^2)/3) + ((b4*(T).^3)/4) + ((b5*(T).^4)/5) + (b6./T));
    S = R*(b1.*log(T)+(b2*T)+((b3*(T).^2)/2)+((b4*(T).^3)/3)+((b5*(T).^4)/4)+b7);
end


% Cp = [Cp1,Cp2];
% H  = [H1,H2];
% S  = [S1,S2];

% Plot for Temperature vs Specific heat
figure(1)
plot(T,Cp,'linewidth',2,'Color','r')
xlabel('Temperature(K)')
ylabel('Specific Heat')
title('Specific heat vs Temperature')

species_name = b{1};
% Creation of Folders

mkdir(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

% changing the folder to Species folder

cd(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

saveas(1,'Cp vs T.jpg')

% Plot for Temperature vs Enthalpy
figure(2)
plot(T,H,'linewidth',2,'Color','g')
xlabel('Temperature(K)')
ylabel('Enthalpy')
title('Enthalpy vs Temperature')

% Creation of Folders

mkdir(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

% changing the folder to Species folder

cd(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

saveas(2,'H vs T.jpg')

% Plot for Temperature vs Entropy
figure(3)
plot(T,S,'linewidth',2,'Color','b')
xlabel('Temperature(K)')
ylabel('Entropy')
title('Entropy vs Temperature')


% Creation of Folders

mkdir(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

% changing the folder to Species folder

cd(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

saveas(3,'S vs T.jpg')

% Calculating molecular weight of different species
Mol_Weight = Molecular_weight(species_name)

end

Explanation: We start the program with clear all , close all and clc. Save the NASA thermodynamic data file in the same folder as the program file. Use fopen command to open and read the input data file and in turn it returns an integer output which is greater than or equal to 3. And after every subsequent execution of the code will increase the integer by 1.

The fgetl command is used 5 times in the for loop to repeat the fgetl command ,which is used to read the first 5 lines from the input data file.

Program for one species :

we use the fgetl command to read the next line(ie the 6th line from the input data file) as the species data starts from this line till 3 more lines after this line.

Next comes the extraction of the local temperatures of the species. Inorder to read the identify the local temperatures we take the 'dot'(ie '.') as the reference point of the tempertures. hence we use strfind command to identify the positons of the '.' . next , we are able to identify and store the value of low_temp in number form using the 'str2double' command. same process is used for high temperature and mid temperature.

The next line shows the use of strsplit(a,' ') command to split the string line stored under the variable named 'a'(ie species data first line). The command will the split and form a cell array made of various strings seperated by a space and stored in the name 'b'.

Next comes the extraction of coefficients of the species. There are 14 coefficients in total for each species which are distributed in the next 3 lines of the species data. After reading the next line using fgetl we use the strfind(coeff,'E') command to identify and locate the coefficients using the 'E' as reference. After locating the coefficients we convert into number and store in a variable using str2double. to find the next coefficient we use the postion of the first coefficient keeping in mind whether the coefficient is a +ve or a -ve coefficient. This way we are able to locate and store the rest of the coefficients of the the species.

Next part is the temperature array creation from the low temperature to the high temperature using the linspace command with 500 values.

Next step is the calculation of the Specific heat , Enthalpy and Entropy values using the inputs as universal gas constant(R) , the coefficients and the mid temprature as the reference. The equations used in the calculation are shown below.

Here we are use the if else command to calculate Cp,H,S. If the temperature is below or equal to the mid temperature then we use the coefficients meant for the lower temperature range. else we use the coefficients meant for the Higher temperature range. We have to be very careful when using the parenthesis and the dot operators. Here we improvised the earlier method which used 3 different if else commands to calculate Cp , H and S separetely. hence we clubbed the 3 equations of Cp , H , S together under if command and the same 3 equations together under else command.

There is an alternate way to calculate the Cp , H, S. It is shown below.

for i = 1:length(T)
    t = T(i);
    if (t<mid_temp)
        ltr(i) = t;
    else
        htr1(i) = t;
    end
end
  
d = length(ltr);
htr = htr1((d+1):length(T));

% Calculation of Cp , H , S (LTR - Low Temp Region)

Cp1 = R*(b8 + (b9*ltr) + (b10*(ltr).^2) + (b11*(ltr).^3) + (b12*(ltr).^4));
H1  = R*ltr.*(b8 + ((b9*ltr)/2) + ((b10*(ltr).^2)/3) + ((b11*(ltr).^3)/4) + ((b12*(ltr).^4)/5) + (b13./ltr));
S1  = R*(b8.*log(ltr) + (b9*ltr) + ((b10*(ltr).^2)/2) + ((b11*(ltr).^3)/3) + ((b12*(ltr).^4)/4) + b14);

% Calculation of Cp , H , S (LTR - High Temp Region)

Cp2 = R*(b1 + (b2*htr) + (b3*(htr).^2) + (b4*(htr).^3) + (b5*(htr).^4));
H2  = R*htr.*(b1 + ((b2*htr)/2) + ((b3*(htr).^2)/3) + ((b4*(htr).^3)/4) + ((b5*(htr).^4)/5) + (b6./htr));
S2  = R*(b1.*log(htr) + (b2*htr) + ((b3*(htr).^2)/2) + ((b4*(htr).^3)/3) + ((b5*(htr).^4)/4) + b7);

Cp = [Cp1,Cp2];
H  = [H1,H2];
S  = [S1,S2];

In this , first we calculate the low temperature range and high temperature range using for loop command and inside the for loop we use the if else command to seperate the lower temperture range and the high temperature range using the mid temperature as reference.

Alternatively we can use the linspace command to calculate the low temperature range and the high temperature range as shown below.

ltr = linspace(low_temp,mid_temp,500);
htr = linspace(mid_temp,high_temp,500);

After the low temperature range and high temperature range arrays are found we calculate seperately the Cp, H and S for the low temperature range first and then we calculate the Cp,H and S for the high temperature range. After that the Cp , H and S for different ranges is combined together in single arrays to find the complete range. This method has a small flaw. For certain species the Plots are not smooth with a very slight dip found in the curve. but the plot is continuous. probably could be because we calculate for the Cp,H and S seperately and then combine later in a single array.

Next step is the plotting and saving of the plots of the species automatically using a program. The code is shown below.

% Plot for Temperature vs Specific heat
figure(1)
plot(T,Cp,'linewidth',2,'Color','r')
xlabel('Temperature(K)')
ylabel('Specific Heat')
title('Specific heat vs Temperature')

species_name = b{1};
% Creation of Folders

mkdir(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

% changing the folder to Species folder

cd(['C:\Users\DELL\Desktop\Matlab\NASA File Parsing challenge\Species folder for NASA file parsing\',species_name])

saveas(1,'Cp vs T.jpg')

We use different figure for Cp vs T , H vs T and S vs T using figure() command. Labels and titles are also added. we mention the species name. Then we form a new folder in which all the plots will be saved. we open the folder and copy the path of the folder and paste it inside brackets of the mkdir command. the mkdir command creates a parent folder and with sub folders of the different species along with the species name. illustration below.

The cd command basically changes the directory of the plots which are saved to the same path as what is mentioned inside the mkdir command brackets. illustration below.

Function Code :

function [MW] = Molecular_weight(species_name)

atomic_name = ['H','C','N','O','Ar']; % Hydrogen,Carbon,Nitrogen,Oxygen,Argon

% Atomic masses of the main species
atomic_weight = [1.008 12.011 14.007 15.990 39.948];

%adding atomic number to obtain particular value
MW = 0;
for i = 1:length(species_name)
    for j = 1:length(atomic_name)
        if strcmp(species_name(i),atomic_name(j))
            MW = MW + atomic_weight(j);
            Position = j;
        end
    end
    n = str2double(species_name(i));
    if n > 1
        MW = MW + atomic_weight(Position)*(n-1);
    end
end

fprintf('molecular weight of %s:',species_name);
fprintf('%f',MW);
disp('');

Here we have first formed arrays for atomic names of the 5 elements used in different species and atomic weight of each element mentioned in the atomic names in the same sequence as that in the atomic names. Then using the initial condition of molecular weight as 0 and a nested for loop we start with the program. in the first for loop we use i for iterations and the length of the species will be the total number of characters and numbers in the species name. Then we use another for loop with iteration as j to find the atomic name. within this for loop we use if command and then compare between the species name with atomic name using strcmp command. in order to find and match the letters of the species with the elements of the atomic names. Next line shows the updating of the molcular weight based on the iteration and atomic name. which keeps updating after every iterations. next we use if command to find the number of molecules of each element. We close all the loops.

The end of the function shows that we use the fprintf command to print the species name and the corresponding moleculer weight.

Inorder to follow the full process for the entire 53 species we use a for loop with length from 1 to 53.

Conclusion:

All the plots are stored in their respective folders. and all the sub folders are stored in the main folder. The screenshot of the main folder showing all the species sub folders is attached. And the zip file containing the plots for O2 , N2 and CO2 are also in the attachments.

Errors Faced:

1. Instead of finding low temp and high temperature , i used ltr(low temp range) and htr(high temp range) as variables to store the low temp and high temp values. low temp and high temp are single element. but instead ltr and htr were added with the single value. hence the use of ltr and htr was not giving the correct output in one of the methods used to calculate cp,H,S.

Corrected this and output was better.

2. for finding the species name I used ( ) brackets instead of the special brackets { } meant for cell arrays . Corrected .

3. In the calculation of Cp , H and S use of proper parenthesis and use of dot operators have to be made.

4. In the alternate method to calculate the Cp , H and S we store and combine the 2 different values using semicolon. corrected and used a comma to make it a row vector.

5. molecular weight calculating function should be placed within the whole for loop and it should not be loop separetely in the for loop as it will give the output of the last species and repeat it 53 times.