Project 1 - Parsing NASA thermodynamic data
the main block of code is as follows:
clear all
close all
clc
% first we read the data
f1=fopen('THERMO.txt','r');
f=fgetl(f1);
%get global temp range in char form
global_temprature_range=fgetl(f1);
% we split the tempratures into high/low
Temp=global_temprature_range;
A=strsplit(Temp,' ');
% we convert the char data into numbers
global_low_temprature=str2num(A{2});
global_high_temprature=str2num(A{4});
global_medium_temprature=str2num(A{3});
%we skip the header info
for i=1:3
fgetl(f1)
end
%we need to differentiate the following data. we use fgetl
%for line 1
for i=1:53 % as there are 53 species
tline_header=fgetl(f1);
%for getting species name
B=strsplit(tline_header,' ');
Species_Name=(B{1});
%for getting local temprature
C=strsplit(tline_header,'G ');
D=str2num(C{2});
T_local_low=D(1);
T_local_high=D(2);
T_local_medium=D(3);
%we create a vector space taking inputs
T=linspace(T_local_low,T_local_high,100);
tline_1=fgetl(f1);
a=findstr(tline_1,'E'); % here position of E is stored in a cell array
a1=tline_1(1:a(1)+3);
a1=str2num(a1);
a2=tline_1(a(1)+4:a(2)+3);
a2=str2num(a2);
a3=tline_1(a(2)+4:a(3)+3);
a3=str2num(a3);
a4=tline_1(a(3)+4:a(4)+3);
a4=str2num(a4);
a5=tline_1(a(4)+4:a(5)+3);
a5=str2num(a5);
%for line 2
tline_2=fgetl(f1);
b=findstr(tline_2,'E');
b1=tline_2(1:b(1)+3);
b1=str2num(b1);
b2=tline_2(b(1)+4:b(2)+3);
b2=str2num(b2);
b3=tline_2(b(2)+4:b(3)+3);
b3=str2num(b3);
b4=tline_2(b(3)+4:b(4)+3);
b4=str2num(b4);
b5=tline_2(b(4)+4:b(5)+3);
b5=str2num(b5);
%for line 3
tline_3=fgetl(f1);
c=findstr(tline_3,'E');
c1=tline_3(1:c(1)+3);
c1=str2num(c1);
c2=tline_3(b(1)+4:c(2)+3);
c2=str2num(c2);
c3=tline_3(c(2)+4:c(3)+3);
c3=str2num(c3);
c4=tline_3(c(3)+4:c(4)+3);
c4=str2num(c4);
%Universal gas constant in KJ/moleK
R = 8.3143;
%calculating Specific heat,Enthalpy,Entropy
Cp=Specific_heat(a1,a2,a3,a4,a5,b5,c1,c2,c3,c4,T,R,T_local_medium);
H =Specific_enthalpy(a1,a2,a3,a4,b4,a5,b1,b5,c1,c2,c3,c4,T,R,T_local_medium);
S=Entropy(a1,a2,a3,a4,a5,b2,b3,b4,b5,c1,c2,c4,T,R,T_local_medium);
%calculating species molecular weight
Molecularweight = Molecular_weights(Species_Name);
%saving the above into folders
cur_dir=pwd; % pwd returns the current folder
%here we create a directory using input arg as our species name which
% will iterate and change with each loop until 53 names are saved in
% the directory
mkdir(Species_Name);
cd(Species_Name);
% plotting the above and labelling title and x and y axis
figure (1)
plot(T,Cp,"LineWidth",3,"Color",'r')
ylabel('Specific Heat J⋅kg−1⋅K−1')
xlabel('Temprature')
grid on
title('Specific Heat vs Temprature',['Species is ',num2str(Species_Name),3])
saveas(gcf,'Cp vs T.png') %gcf gives the current plot and it's saved as png with the 2nd argument
%as name for the plot
figure (2)
plot(T,H,"LineWidth",3,"Color",'b')
ylabel('Specific Enthalpy kJ/kg')
xlabel('Temprature')
grid on
title('Enthalpy vs Temprature',['Species is ',num2str(Species_Name),3])
saveas(gcf,'H vs T.png')
figure (3)
plot(T,S,"LineWidth",3)
ylabel('Entropy J⋅K−1')
xlabel('Temprature')
grid on
title('Entropy vs Temprature',['Species is ',num2str(Species_Name),3])
saveas(gcf,'S vs T.png')
end
the function definition is as follows
1. Function Specific_heat
% from main code
% create a function Cp using the create function option and take inputs
% from main code
function Cp = Specific_heat(a1,a2,a3,a4,a5,b5,c1,c2,c3,c4,T,R,T_local_medium)
if T >= T_local_medium
Cp=R*((a1)+(a2*T)+(a3*T.^2)+(a4*T.^3)+(a5*T.^4));
else
Cp=R*((b5)+(c1*T)+(c2*T.^2)+(c3*T.^3)+(c4*T.^4));
end
end
2. Function Specific Enthalpy
% create a function Cp using the create function option and take inputs
% from main code
function H = Specific_enthalpy(a1,a2,a3,a4,b4,a5,b1,b5,c1,c2,c3,c4,T,R,T_local_medium)
% we divide the temprature range vector into 2 parts high or low with
% inputs from main code
if T >= T_local_medium
H=R.*T.*((a1)+(a2.*T/2)+(a3*T.^2/3)+(a4.*T^3/4)+(a5*T.^4/5)+(b1./T));
else
H=R.*T.*((b4)+(b5.*T/2)+(c1*T.^2/3)+(c2*T.^3/4)+(c3*T.^4/5)+(c4./T));
end
end
3.Function Entropy
% we do the same for entropy function
function S=Entropy(a1,a2,a3,a4,a5,b2,b3,b4,b5,c1,c2,c4,T,R,T_local_medium)
if T>=T_local_medium
S=R.*((a1*log(T))+(a2*T)+(a3*T.^2)+(a4*T.^3/3)+(a5*T.^4/4)+(b2));
else
S=R.*((b3*log(T))+(b4*T)+(b5*T.^2)+(c1*T.^3/3)+(c2*T.^4/4)+(c4));
end
end
4. Function Molecular Weights
% Molecular Weight Calculation function
function Molecularweight = Molecular_weights(Species_Name)
% Standard - Atomic Name
atomic_name = ['H','C','O','N','A'];
% Hydrogen, Carbon, Oxygen, Nitrogen, Argon
fprintf('Molecular Weight of ');
fprintf(Species_Name);
fprintf(' : ');
% Standard Atomic Number
Atomic_Weight = [1.008 12.011 15.999 14.007 39.948];
% We start the counter with 0 for molecular weight
Molecularweight = 0;
% Adding Atomic number to obtain particular value
for i = 1:length(Species_Name)
for j = 1:length(atomic_name)
if strcmpi(Species_Name(i),atomic_name(j))
%if above comparision is true value of molecular mass is
%calculated
Molecularweight = Molecularweight + Atomic_Weight(j);
position = j;
end
end
n = str2num(Species_Name(i));
if n>1
%comparing if the value of k is greater than executed the further statements
Molecularweight = Molecularweight + Atomic_Weight(position)*(n-1);
end
end
fprintf('molecular weight %s :',Species_Name);
fprintf('%d',Molecularweight);
disp(' ');
end
Molecular weights screenshots
Mistakes made
Here the argument contains scalar arguments hence curly braces are not required
couple of mistakes here the brackets were not closed and T needs a elementwise operator(.)
after a few iterations i found out that i have to use LOCAL TEMPRATURE VARIABLES and not the global temprature from the second line of the thermo.dat file. Hence i redid the code in a seperate instance shown in following screenshot
then i applied the same logic and inserted in into my main code shown below with a little tweaks.also had to change the global temp code in functions to local temp.
Referrences
https://de.mathworks.com/help/matlab/ref/pwd.html
https://de.mathworks.com/help/matlab/ref/strsplit.html
https://de.mathworks.com/help/matlab/matlab_prog/matlab-operators-and-special-characters.html
https://de.mathworks.com/help/matlab/ref/strcmpi.html
https://drive.google.com/file/d/1J_imYOhAVYRKRttYDUu4XY0q5sOu26rL/view. (The matlab pdf provided by skill-Lync)
https://de.mathworks.com/matlabcentral/answers/106320-make-a-loop-to-calculate-molecular-weight