Project 1 - Parsing NASA thermodynamic data
Main Code:
clear all
close all
clc
%opening the thermo.dat file
f1=fopen('THERMO.dat')
%skipping first five lines
for i=1:5
fgetl(f1);
end
%reaching the first element line to access the values of temperature range
a=fgetl(f1);
a=strsplit(a);
%Accessing the local temperature ranges
low_temp=str2num(a{6});
high_temp=str2num(a{7});
mid_temp=str2num(a{8});
%defining temperature
t=linspace(low_temp,high_temp,100);
%setting back the file pointer to the very first line
frewind(f1)
%Skipping the lines
for i=1:5
fgetl(f1);
end
for i=1:53
a= strsplit(fgetl(f1));
element{i}=a{1};
b=fgetl(f1);
c=fgetl(f1);
d=fgetl(f1);
%Finding the position of E in the data prior to getting the
%coefficients
e=strfind(b,'E');
f=strfind(c,'E');
h=strfind(d,'E');
%Obtaining the first 7 coefficients for lower temp range
c1(i)=str2num(b(1:e(1)+3));
c2(i)=str2num(b(e(1)+4:e(2)+3));
c3(i)=str2num(b(e(2)+4:e(3)+3));
c4(i)=str2num(b(e(3)+4:e(4)+3));
c5(i)=str2num(b(e(4)+4:e(5)+3));
c6(i)=str2num(c(1:f(1)+3));
%obtaining coefficients for higher temperature range
c7(i)=str2num(c(f(1)+4:f(2)+3));
C1(i)=str2num(c(f(2)+4:f(3)+3));
C2(i)=str2num(c(f(3)+4:f(4)+3));
C3(i)=str2num(c(f(4)+4:f(5)+3));
C4(i)=str2num(d(1:h(1)+3));
C5(i)=str2num(d(h(1)+4:h(2)+3));
C6(i)=str2num(d(h(2)+4:h(3)+3));
C7(i)=str2num(d(h(3)+4:h(4)+3));
%Obtaining the values of Cp,h,s for local temperature ranges
for j=1:length(t)
if t(j)<mid_temp
cp(j)=(c1(i)+c2(i)*t(j)+c3(i)*t(j)^2+c4(i)*t(j)^3+c5(i)*t(j)^4)*8.314;
h(j)=(c1(i)+(c2(i)*t(j)/2)+(c3(i)/3)*(t(j)^2)+(c4(i)/4)*(t(j)^3)+(c5(i)/5)*(t(j)^4)+(c6(i)/t(j)))*8.314*t(j);
s(j)=(c1(i)*log(t(j))+(c2(i)*t(j))+(c3(i)/2)*(t(j)^2)+(c4(i)/3)*(t(j)^3)+(c5(i)/4)*(t(j)^4)+c7(i))*8.314;
else
cp(j)=(C1(i)+C2(i)*t(j)+C3(i)*t(j)^2+C4(i)*t(j)^3+C5(i)*t(j)^4)*8.314;
h(j)=(C1(i)+(C2(i)*t(j)/2)+(C3(i)/3)*(t(j)^2)+(C4(i)/4)*(t(j)^3)+(C5(i)/5)*(t(j)^4)+(C6(i)/t(j)))*8.314*t(j);
s(j)=(C1(i)*log(t(j))+(C2(i)*t(j))+(C3(i)/2)*(t(j)^2)+(C4(i)/3)*(t(j)^3)+(C5(i)/4)*(t(j)^4)+C7(i))*8.314;
end
end
%creating a directory to save the plot files
mkdir ('C:\Users\arulp\OneDrive\Documents\MATLAB\Examples\R2021a\optim\GetStartedWithOptimizeLiveEditorTaskExample\elements',element{i});
%reaching that directory for saving files
reachthis=(['C:\Users\arulp\OneDrive\Documents\MATLAB\Examples\R2021a\optim\GetStartedWithOptimizeLiveEditorTaskExample\elements\',element{i}]);
cd(reachthis)
%Plotting Cp vs T
figure(1)
plot(t,cp,'-r');
title('CP Vs temperature of',element{i})
xlabel('Temperature');
ylabel('cp values');
%saving Cp vs T in the folder
saveas(figure(1),'cp.jpg')
%Plotting H vs T
figure(2)
plot(t,h,'-g');
title('Enthalpy Vs temperature of',element{i})
xlabel('Temperature');
ylabel('H');
%Saving H vs T
saveas(figure(2),'enthalpy.jpg')
%Plotting S vs T
figure(3)
plot(t,s,'-m');
title('Entropy Vs temperature of',element{i})
xlabel('Temperature');
ylabel('Entropy');
%Saving S vs T
saveas(figure(3),'entropy.jpg')
end
Explanation of main Code:
1)Thermo.dat file was downloaded and saved in the current directory so that it could later be opened by matlab
2)File was opened using 'fopen' command
3)fgetl command to reach the first few lines .Global temps has no relevance to this project. Hence,skipped and frewind command was sometimes used to get back to the first line
4)For loop was used to go over from one element to the other and to repeat the following steps for every element to get its name, coefficients and plot
5)'strsplit' was used to break the first line string of each element into cell array to obtain the species name
6) 'Strfind' was used to find the location of 'E' with reference to which the coefficients would be later segregated
7)Using 'str2num' command the coeffcients were obtained for local temp ranges
8)mkdir command was used to create directory for storing the plots of different species in the respective folder
9)dir command was then used to reach the directory where the files specific to each species was to be stored
10)H vs T, S vs T, Cp vs T graphs were plotted using plot command
11)File was closed using fclose command
Function to find molecular weight:
%%FUNCTION TO CALCULATE MOLECULAR WEIGHT OF MOLECULES
function molecular_weight=MOL_WEIGHT(x)
elements={'N','H','O','C','AR'};
atomicweights=[14,1,16,12,40];
molecular_wt=0;
for i=1:length(x)
if str2num(x(i))>=1
molecular_wt=molecular_wt+(str2num(x(i))-1)*at_wt;
else
for j=1:length(elements)
if strcmpi(x(i),elements{j})
at_wt=atomicweights(j);
molecular_wt=molecular_wt+strcmpi(x(i),elements{j})*at_wt;
end
end
end
end
Z=sprintf('The molecular weight of %s is %f amu',x,molecular_wt);
disp(Z)
end
Explanation of the function code:
1) Distinct elements present in the 53 compounds wee identified. They are just 5 in numbers. So, 2 arrays one for the element-a cell array and the other for corresponding atomic weights were created.
2)strcmp command was used to go over from one element of the species to the other and compared with each element in the previously defined cell array containing elements. If matched, molecular weight was incremented by a value equal to the atomic weight of that element.
3)If it was a number, then the molecular weight was incremented by a value equal to atomic weight of previous element multiplied by one less than this value.
4) Thus, when this process is over, the molecular weight is finally obtained
Output files have been uploaded:
1)Molecular weight -2 files- mol weight1.jpg, mol.weight2.jpg
2)screenshot of folder containing species files- element screenshot1.jpg, elementscreenshot2.jpg
3)Cp, H,S vs T for O2,CO2,N2-(3+3+3) jpg files