BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

Parsing NASA thermodynamic data

SShabarish P R20 Jul 2018

The Function that extracts the 14 Coefficients for a given species is as below

LOGIC :

  • In this Function we use \"fopen\" to open a data file,
  • \"fgetl\" is used to read a particular line, we can either store the line read into some variable or we may just read and skip the lines in the data file
  • \"strsplit\", is used to split the data read, we must also specify the condition or a parameter based on which the string is split.
  • \"strfind\" and \"strcmp\"are used to find the data searched and compare the data with a reference.
    Using these tools we build a function that segregates the data, from a given file.

CODE:

function [a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14] = Data(user_string)
F1 = fopen(\'NASA_thermodynamic_data.dat\', \'r\');
fgetl(F1);
B = fgetl(F1);
a = strsplit(B, \' \');
global_low_temperature = str2double(a{2});
global_mid_temperature = str2double(a{3});
global_high_temperature = str2double(a{4});
F2 = fopen (\'DATA.txt\', \'r\');
for q = 1:5
    q = fgets(F2);
end
Y = fopen(\'Data2.txt\',\'w\');
for i = 1:209    
    i = fgetl(F2);
    C = strsplit(i, \' \');
    string = C{1};
    U = strcmp(string, user_string);
    if U == 1 
    for k = 1:3
        k = fgetl(F2);
        fprintf(Y, \'%s\\n\', k);
    end
    end
end
H = fopen(\'Data2.txt\', \'r\');
X = fgetl(H);
Y = fgetl(H);
Z = fgetl(H);
x = strfind(X, \'E\');
y = strfind(Y, \'E\');
z = strfind(Z, \'E\');
A1 = X(1:x(1)+3);
a1 = str2double(A1);
A2 = X(x(1)+4:x(2)+3);
a2 = str2double(A2);
A3 = X(x(2)+4:x(3)+3);
a3 = str2double(A3);
A4 = X(x(3)+4:x(4)+3);
a4 = str2double(A4);
A5 = X(x(4)+4:x(5)+3);
a5 = str2double(A5);
A6 = Y(1:y(1)+3);
a6 = str2double(A6);
A7 = Y(y(1)+4:y(2)+3);
a7 = str2double(A7);
A8 = Y(y(2)+4:y(3)+3);
a8 = str2double(A8);
A9 = Y(y(3)+4:y(4)+3);
a9 = str2double(A9);
A10 = Y(y(4)+4:y(5)+3);
a10 = str2double(A10);
A11 = Z(1:z(1)+3);
a11 = str2double(A11);
A12 = Z(z(1)+4:z(2)+3);
a12 = str2double(A12);
A13 = Z(z(2)+4:z(3)+3);
a13 = str2double(A13);
A14 = Z(z(3)+4:z(4)+3);
a14 = str2double(A14);
end

The Function to calculate the Specific Heat, Entropy and Enthalpy is as given below

LOGIC:

In this function we define the formula which is given below:

The 14 Coefficients obtained by using the above function is given as inputs to this function, however we see that the 3 formulas shown in the figure use only 7 coefficients.
Therefore,

  • we use the first 7 coefficients for High Temperature Range Cp, H and S calculation.
  • we use the second 7 coefficients for low Temperature Range Cp, H and S calculation.

CODE:

function [Cp1,H1,S1,Cp2,H2,S2] = CpHandScalculation(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14)
%Defining the value of Universal Gas Constant%
R = 8.314; % in J/mol.k %
%Calculation of Cp, H and S%
ct1 = 1;
for temperature_range = 300:1000
    lt = temperature_range;
    Cp1(ct1) = (R*((a8) + (a9*lt) + (a10*(lt)^2) + (a11*(lt)^3) + (a12*(lt)^4)));
    H1(ct1) = (R*lt*((a8) + (a9*(lt/2)) + (a10*(((lt)^2)/3)) + (a11*(((lt)^3)/4)) + (a12*(((lt)^4)/4))+ (a13/lt)));
    S1(ct1) = (R*((a8*(log(lt))) + (a9*lt) + (a10*(((lt)^2)/2)) + (a11*(((lt)^3)/3)) + (a12*(((lt)^4))/4)) + a14);
    ct1 = ct1+1;
end
ct2 = 1;
for temperature_range = 1001:5000
        ht = temperature_range;
        Cp2(ct2) = (R*((a1) + (a2*ht) + (a3*(ht)^2) + (a4*(ht)^3) + (a5*(ht)^4)));
        H2(ct2) = (R*ht*((a1) + (a2*(ht/2)) + (a3*(((ht)^2)/3)) + (a4*(((ht)^3)/4)) + (a5*(((ht)^4)/4))+ (a6/ht)));
        S2(ct2) = (R*((a1*(log(ht))) + (a2*ht) + (a3*(((ht)^2)/2)) + (a4*(((ht)^3)/3)) + (a5*(((ht)^4))/4)) + a7);
        ct2 = ct2+1;
end
end

The Function to Calculate the Molecular weight of a given species is as given below

LOGIC:

In this Function we first take the user_string which is the Species Name as an Input,
One can observe that the Species Name is an array, thus we take each letter of the array (Which is an Element name) and assign to it, its molecular weight.
Meanwhile if there is a number after an element name, we multiply the number with its Element molecular weight and subtract one of the element molecular weight values.
Finally we add the molecular weight of each of the Element to get the final Molecular weight of the species.
[Please look into the code to understand the calculations Explained above more clearly]

CODE:

function[MolecularWeight] = Molecularweight(user_string)
if strcmp(user_string(1), \'C\') 
    Q=12.0107;
elseif strcmp(user_string(1), \'H\') 
    Q=1.00794;
elseif strcmp(user_string(1), \'O\')
    Q=15.999;
elseif strcmp(user_string(1), \'N\') 
    Q=14.0067;
elseif strcmp(user_string(1), \'A\') 
    Q=39.948;
elseif strcmp(user_string(1), \'R\') 
    Q=0;
else
end
if length(user_string) < 2;
    V=0;
elseif strcmp(user_string(2), \'C\')
    V=12.0107;
elseif strcmp(user_string(2), \'H\') 
    V=1.00794;
elseif strcmp(user_string(2), \'O\')
    V=15.999;
elseif strcmp(user_string(2), \'N\') 
    V=14.0067;
elseif strcmp(user_string(2), \'1\') 
    V=(1*Q)-Q;
elseif strcmp(user_string(2), \'2\') 
    V=(2*Q)-Q;
elseif strcmp(user_string(2), \'3\') 
    V=(3*Q)-Q;
elseif strcmp(user_string(2), \'4\') 
    V=(4*Q)-Q;
elseif strcmp(user_string(2), \'5\') 
    V=(5*Q)-Q;
elseif strcmp(user_string(2), \'6\') 
    V=(6*Q)-Q;
elseif strcmp(user_string(2), \'7\') 
    V=(7*Q)-Q;
elseif strcmp(user_string(2), \'8\') 
    V=(8*Q)-Q;
elseif strcmp(user_string(2), \'A\') 
    V=39.948;
elseif strcmp(user_string(2), \'R\') 
    V=0;
else
end

if length(user_string) < 3;
    W=0;
elseif strcmp(user_string(3), \'C\') 
    W=12.0107;
elseif strcmp(user_string(3), \'H\') 
    W=1.00794;
elseif strcmp(user_string(3), \'O\')
    W=15.999;
elseif strcmp(user_string(3), \'N\') 
    W=14.0067;
elseif strcmp(user_string(3), \'1\') 
    W=(1*V)-V;
elseif strcmp(user_string(3), \'2\') 
    W=(2*V)-V;
elseif strcmp(user_string(3), \'3\') 
    W=(3*V)-V;
elseif strcmp(user_string(3), \'4\') 
    W=(4*V)-V;
elseif strcmp(user_string(3), \'5\') 
    W=(5*V)-V;
elseif strcmp(user_string(3), \'6\') 
    W=(6*V)-V;
elseif strcmp(user_string(3), \'7\') 
    W=(7*V)-V;
elseif strcmp(user_string(3), \'8\') 
    W=(8*V)-V;
elseif strcmp(user_string(3), \'A\') 
    W=39.948;
elseif strcmp(user_string(3), \'R\') 
    W=0;
else
end

if length(user_string) < 4;
    Y=0;
elseif strcmp(user_string(4), \'C\') 
    Y=12.0107;
elseif strcmp(user_string(4), \'H\') 
    Y=1.00794;
elseif strcmp(user_string(4), \'O\')
    Y=15.999;
elseif strcmp(user_string(4), \'N\') 
    Y=14.0067;
elseif strcmp(user_string(4), \'1\') 
    Y=(1*W)-W;
elseif strcmp(user_string(4), \'2\') 
    Y=(2*W)-W;
elseif strcmp(user_string(4), \'3\') 
    Y=(3*W)-W;
elseif strcmp(user_string(4), \'4\') 
    Y=(4*W)-W;
elseif strcmp(user_string(4), \'5\') 
    Y=(5*W)-W;
elseif strcmp(user_string(4), \'6\') 
    Y=(6*W)-W;
elseif strcmp(user_string(4), \'7\') 
    Y=(7*W)-W;
elseif strcmp(user_string(4), \'8\') 
    Y=(8*W)-W;
elseif strcmp(user_string(4), \'A\') 
    Y=39.948;
elseif strcmp(user_string(4), \'R\') 
    Y=0;
else
end

if length(user_string) < 5;
    Z=0;
elseif strcmp(user_string(5), \'C\') 
    Z=12.0107;
elseif strcmp(user_string(5), \'H\') 
    Z=1.00794;
elseif strcmp(user_string(5), \'O\')
    Z=15.999;
elseif strcmp(user_string(5), \'N\') 
    Z=14.0067;
elseif strcmp(user_string(5), \'1\') 
    Z=(1*Y)-Y;
elseif strcmp(user_string(5), \'2\') 
    Z=(2*Y)-Y;
elseif strcmp(user_string(5), \'3\') 
    Z=(3*Y)-Y;
elseif strcmp(user_string(5), \'4\') 
    Z=(4*Y)-Y;
elseif strcmp(user_string(5), \'5\') 
    Z=(5*Y)-Y;
elseif strcmp(user_string(5), \'6\') 
    Z=(6*Y)-Y;
elseif strcmp(user_string(5), \'7\') 
    Z=(7*Y)-Y;
elseif strcmp(user_string(5), \'8\') 
    Z=(8*Y)-Y;
elseif strcmp(user_string(5), \'A\') 
    Z=39.948;
elseif strcmp(user_string(5), \'R\') 
    Z=0;
else
end

if length(user_string) < 6;
    U=0;
elseif strcmp(user_string(6), \'C\') 
    U=12.0107;
elseif strcmp(user_string(6), \'H\') 
    U=1.00794;
elseif strcmp(user_string(6), \'O\')
    U=15.999;
elseif strcmp(user_string(6), \'N\') 
    U=14.0067;
elseif strcmp(user_string(6), \'1\') 
    U=(1*Z)-Z;
elseif strcmp(user_string(6), \'2\') 
    U=(2*Z)-Z;
elseif strcmp(user_string(6), \'3\') 
    U=(3*Z)-Z;
elseif strcmp(user_string(6), \'4\') 
    U=(4*Z)-Z;
elseif strcmp(user_string(6), \'5\') 
    U=(5*Z)-Z;
elseif strcmp(user_string(6), \'6\') 
    U=(6*Z)-Z;
elseif strcmp(user_string(6), \'7\') 
    U=(7*Z)-Z;
elseif strcmp(user_string(6), \'8\') 
    U=(8*Z)-Z;
elseif strcmp(user_string(6), \'A\') 
    U=39.948;
elseif strcmp(user_string(6), \'R\') 
    U=0;
else
end
MolecularWeight = Q+V+W+Y+Z+U;
end

The function to plot Specific Heat Vs Temperature, Enthalpy Vs Temperature and Entropy Vs Temperature and to store them in a folder and with relevant names is as given below

LOGIC:

In this Function we first define the X axis which is the Temperature, we define this using linspace command, we define the temperature range from Low Temperature to Mid Temperature, Later we define from Mid Temperature to High Temperature Range.
Using the Cp, S and H obtained from the function defined above, we plot them against Temperature, we use \"Hold on\" and \"Hold off\" command to generate the complete Plot between the Low and High Temperature Ranges.

Using \"mkdir\", we create the folder, later using \'saveas\' we will store the plots generated in the folder Created

CODE:

function[figure1,figure2,figure3] = Plots(user_string,Cp1,H1,S1,Cp2,H2,S2)
%First Plot%
figure(1);
figure1 = figure(1);
TemperatureRange1A = linspace(300,1000,701);
plot(TemperatureRange1A,Cp1, \'linewidth\', 1, \'color\', \'r\');
title(\'Specific Heat (Cp) V/S Temperature [K]\');
xlabel(\'Temperature [K]\');
ylabel(\'Specific Heat (Cp)[J/molK]\');
hold on 
TemperatureRange1B = linspace(1001,5000,4000);
plot(TemperatureRange1B,Cp2, \'linewidth\', 1, \'color\', \'g\');
hold off;
%Second Plot%
figure(2);
figure2 = figure(2);
TemperatureRange2A = linspace(300,1000,701);
plot(TemperatureRange2A,S1, \'linewidth\', 1, \'color\', \'r\');
hold on;
TemperatureRange2B = linspace(1001,5000,4000);
plot(TemperatureRange2B,S2, \'linewidth\', 1, \'color\', \'g\');
hold off;
title(\'Entropy (S)[J/mol/K] V/S Temperature [K]\');
xlabel(\'Temperature [K]\');
ylabel(\'Entropy (S)[J/mol/K]\');
%Third Plot%
figure(3);
figure3 = figure(3);
TemperatureRange3A = linspace(300,1000,701);
plot(TemperatureRange3A,H1, \'linewidth\', 1, \'color\', \'r\');
hold on;
TemperatureRange3B = linspace(1001,5000,4000);
plot(TemperatureRange3B,H2, \'linewidth\', 1, \'color\', \'g\');
hold off;
title(\'Enthalpy (H)[KJ/mol] V/S Temperature [k]\');
xlabel(\'Temperature [K]\');
ylabel(\'Enthalpy (H) [KJ/mol]\');
mkdir ResultsOfProgram;
mkdir(\'ResultsOfProgram\', user_string)
saveas(figure1, [pwd \'/ResultsOfProgram/SpecificHeatVsTemperature.jpg\']);
saveas(figure2, [pwd \'/ResultsOfProgram/EntropyVsTemperature.jpg\']);
saveas(figure3, [pwd \'/ResultsOfProgram/EnthalpyVsTemperature.jpg\']);
end

The Results are stored within the folder \'ResultsofProgram\'.

THE MAIN PROGRAM TO PERFORM ALL THE INSTRUCTIONS GIVEN IN THE QUESTION IS AS BELOW

Description:

In this Code we call the functions defined above we use them to read the data, calculate the Cp, H and S, Plot them and find the Molecular weight of each. 

%Program for Parsing NASA thermodynamic data%
clear all;
close all;
clc;
%We ask for a Species name%
user_string = input(\'Enter the Species name: \',\'s\');
%We generate the Coefficients for the given Species%
[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14] = Data(user_string);
%We calculate the Cp, H and S values for the given Species%
[Cp1,H1,S1,Cp2,H2,S2] = CpHandScalculation(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14);
%We plot the respective Graphs for the same%
[figure1,figure2,figure3] = Plots(user_string,Cp1,H1,S1,Cp2,H2,S2);
sprintf(\'The Species given as input is: %s\\n\', user_string)
%We calculate the Molecular Weight of each Species%
[MolecularWeight] = Molecularweight(user_string);
fprintf(\'The Molecular weight of the species is: %f\\n\', MolecularWeight)

Some Results:

  1. For Species Name O2:

     Figure: Molecular Weight

Figure: CP Vs Temperature

Figure: S Vs Temperature

Figure: H Vs Temperature

  1. For Species Name N2:

Figure: Molecular Weight

Figure: Cp Vs Temperature

Figure: H Vs Temperature

Figure: S Vs Temperature

  1. For Species Name CO2:

Figure: Molecular Weight

Figure: Cp Vs Temperature

Figure: H Vs Temperature

Figure: S Vs Temperature

Observation:

We see the discontinuities in the plots obtained, the reason for these discontinuities can be due to the difference between the values of 1st set of 7 Coefficients and 2nd set of 7 coefficients extracted from the data file, which was used to calculate Cp, H and S.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.