Project 1 - Parsing NASA thermodynamic data
Aim:
Write a function that extracts the 14 co-efficients and calculates the enthalpy, entropy and specific heats for all the species in the data file.
The FIRST 7 coefficients are HIGH-temperature coefficients and the SECOND 7 coefficients are LOW-temperature coefficients.
Also use the LOCAL TEMPERATURE RANGES which is unique for each individual species, not the GLOBAL TEMPERATURE RANGE.
2.Calculate the molecular weight of each species and display it in the command window.
3. Plot the Cp, Enthalpy and Entropy for the local temperature range (low temperature : high temperature) specific for each species.
4. Save the plots as images with appropriate names and dump them in separate folders for each species with suitable name. All these should be generated automatically.
5. The main program that you will submit will just take the "THERMO.dat" file as an input and generate Cp, H and S plots.
Main code:
clear all
close all
clc
%Definig the constants
R=8.314; %R= Molar gas constant
f1=fopen('THERMO.dat','r');
fgetl(f1);
global_temp=fgetl(f1);%Reading comment line
temp_range=strsplit(global_temp,' '); % Splitting line with space as seperator
% Defining the Temperature ranges
low_temp=str2num(temp_range{2});
med_temp=str2num(temp_range{3});
high_temp=str2num(temp_range{4});
comments1=fgetl(f1); %Reading comment line
comments2=fgetl(f1); %Reading comment line
comments3=fgetl(f1); %Reading comment line
for k=1:53
% Reading next line
A=fgetl(f1);
c=strsplit(A,' ');
species_name=c{1};
% Defining Local Temperature
local_low_temp=c{end-3};
local_low_temp=str2num(local_low_temp);
local_high_temp=c{end-2};
local_high_temp=str2num(local_high_temp);
local_mid_temp=c{end-1};
local_mid_temp=str2num(local_mid_temp);
T=linspace(local_low_temp,local_high_temp,1000);
% Reading the coefficients
line1=fgetl(f1);% Reading next line to get the first set of coefficients.
a=strfind(line1,'E');
a1=line1(1:a(1)+3);
% a1 to a7-High temperature coefficients.
a1=str2num(a1);
a2=line1(a(1)+4:a(2)+3);
a2=str2num(a2);
a3=line1(a(2)+4:a(3)+3);
a3=str2num(a3);
a4=line1(a(3)+4:a(4)+3);
a4=str2num(a4);
a5=line1(a(4)+4:a(5)+3);
a5=str2num(a5);
line2=fgetl(f1);%Reading next line to get the second set of coefficients.
a=strfind(line2,'E');
a6=line2(1:a(1)+3);
a6=str2num(a6);
a7=line2(a(1)+4:a(2)+3);
a7=str2num(a7);
% a8 to a14-Low temperature coefficients.
a8=line2(a(2)+4:a(3)+3);
a8=str2num(a8);
a9=line2(a(3)+4:a(4)+3);
a9=str2num(a9);
a10=line2(a(4)+4:a(5)+3);
a10=str2num(a10);
line3=fgetl(f1);%Reading next line to get the third set of coefficients.
a=strfind(line3,'E');
a11=line3(1:a(1)+3);
a11=str2num(a11);
a12=line3(a(1)+4:a(2)+3);
a12=str2num(a12);
a13=line3(a(2)+4:a(3)+3);
a13=str2num(a13);
a14=line3(a(3)+4:a(4)+3);
a14=str2num(a14);
% Calculating Specific Heat(Cp),Enthalpy(H) and Entropy(S)
Cp=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,local_mid_temp);
H=enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,T,R,local_mid_temp);
S=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,local_mid_temp);
mkdir(species_name);
cd (species_name);
figure(1)
plot(T,Cp,LineWidth=2,Color='r')
xlabel('Temperature(T)')
ylabel('SpecificHeat(Cp)')
title(sprintf('Temperature vs Specific Heat of %s',species_name))
figure(2)
plot(T,H,Linewidth=2,Color='b')
xlabel('Temperature(T)')
ylabel('Enthalpy(H)')
title(sprintf('Temperature vs Enthalpy of %s',species_name))
figure(3)
plot(T,S,Linewidth=2,color='g')
xlabel('Temperature(T)')
ylabel('Entropy(S)')
title(sprintf('Temperature vs Entropy of %s',species_name))
saveas(1,'Specific Heat','jpg');
saveas(2,'Enthaply','jpg');
saveas(3,'Entropy','jpg');
cd ..
element=['O','N','C','H','A','S'];
atomic_weight=[16 14 12 1 40 32];
W=0;
for i=1:length(species_name)
for j=1:length(element)
if strcmp(species_name(i),element(j))
W=W+atomic_weight(j);
pos=j;
end
end
n=str2num(species_name(i));
if n>1
W=W+atomic_weight(pos)*(n-1);
end
end
fprintf('Molecular weight of %s = %f',species_name,W);
fprintf('\n');
end
Function calling for specific heat
%Function for calling Specific Heat
function[Cp]=specific_heat(a1,a2,a3,a4,a5,a8,a9,a10,a11,a12,T,R,local_mid_temp)
if T<=local_mid_temp
Cp=(a8+(a9*T)+(a10*T.^2)+(a11*T.^3)+(a12*T.^4))*R;
else
Cp=(a1+(a2*T)+(a3*T.^2)+(a4*T.^3)+(a5*T.^4))*R;
end
end
Function calling for Enthapy
%Function for calling Enthalpy
function[H]=enthalpy(a1,a2,a3,a4,a5,a6,a8,a9,a10,a11,a12,T,R,local_mid_temp)
if T<local_mid_temp
H=R.*((a8.*T)+(((a9.*T).^2)/2)+(((a10.*T).^3)/3)+(((a11.*T).^4)/4)+(((a12.*T).^5)/5)+a13);
else
H=R.*((a1.*T)+(((a2.*T).^2)/2)+(((a3.*T).^3)/3)+(((a4.*T).^4)/4)+(((a5.*T).^5)/5)+a6);
end
end
Function calling for entropy
% Function for calling Entropy
function[S]=entropy(a1,a2,a3,a4,a5,a7,a8,a9,a10,a11,a12,a14,T,R,local_mid_temp)
if T<local_mid_temp
S=R.*((a8.*log(T))+(a9.*T)+((a10.*T.^2)/2)+((a11.*T.^3)/3)+((a12.*T.^4)/4)+a14);
else
S=R.*((a1.*log(T))+(a2.*T)+((a3.*T.^2)/2)+((a4.*T.^3)/3)+((a5.*T.^4)/4)+a7);
end
end
Step-wise explanation of code.