TITLE: 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. You will be reading this file NASA thermodynamic data(Tip: use fopen and fgetl). The formulae are shown below.
Here R is the universal gas constant, T is the temeprature
EQUATIONS:
cp/R=a1+a2.T+a3.T^2+a4.T^3+a5.T^4;
h/RT=a1+a2.(T/2)+a3.(t^2/3)+a4.(T^3/4)+a5.(T^4/5)+a6/T;
s/R=a1.logT+a2.T+a3.(T^2/2)+a4.(T^3/3)+a5.(T^4/4)+a7;
where,
cp= specific heat
R=universal gas constant
a1 to a14 are the data values for the coefficients of the polynomials.
T1,T2,T3 are global temperatures
s=entropy
h=enthalpy
# MAIN CODE
clear all
close all
clc
% getting data file
f=fopen('THERMO.dat','r');
fgetl(f);
% reading header information
tline=fgetl(f)
A=str2num(tline);
global_low_temp=A(:,1);
global_mid_temp=A(:,2);
global_high_temp=A(:,3);
% skipping comments lines
for i=1:3
fgetl(f);
end
% reading species data
for i=1:53
line1=fgetl(f);
l=strsplit(line1);
m=length(l)
sp=l(1);
species=char(sp);
atom=['H','C','O','N','Ar','S']
atomic_weight=[1.008 12.001 15.999 14.007 39.948 32.065]
low_temp=str2double(l(m-3));
mid_temp=str2double(l(m-2));
high_temp=str2double(l(m-1));
% storing the values of second line
line2=fgetl(f);
a=strfind(line2,'E');
a1=str2double(line2(1:(a+3)));
a2=str2double(line2((a(1)+4):(a(2)+3)));
a3=str2double(line2((a(2)+4):(a(3)+3)));
a4=str2double(line2((a(3)+4):(a(4)+3)));
a5=str2double(line2((a(4)+4):(a(5)+3)));
% storing value of 3rd line
line3=fgetl(f)
b=strfind(line3,'E');
a6=str2double(line3(1:(b(1)+3)));
a7=str2double(line3((b(1)+4):(a(2)+3)));
a8=str2double(line3((b(2)+4):(a(3)+3)));
a9=str2double(line3((b(3)+4):(a(4)+3)));
a10=str2double(line3((b(4)+4):(a(5)+3)));
% storing value of 4th line
line4=fgetl(f)
c=strfind(line3,'E');
a11=str2double(line4(1:(c(1)+3)));
a12=str2double(line4((c(1)+4):(a(2)+3)));
a13=str2double(line4((c(2)+4):(a(3)+3)));
a14=str2double(line4((c(3)+4):(a(4)+3)));
% creating temp values for lower and higher temp
TL=linspace(low_temp,mid_temp,100);
TH=linspace(mid_temp,high_temp,100);
% values of universal gas constant
R=8.314;
%calculating specific heat,enthalpy ,entropy for lower and higher temp
for k=1:100
% for lower temp
cp_low(k)=(a8+a9*TL(k)+a10*TL(k)^2+a11*TL(k)^3+a12*TL(k)^4)*R;
h_low(k)=(a8+a9*TL(k)/2+a10*(TL(k)^2+a11*(TL(k)^3)/4)+a12*(TL(k)^4)/5)+a13/TL(k)*R;
s_low(k)=(a8*log(TL(k))+a9*TL(k)+a10*(TL(k)^2)/2+a11*(TL(k)^3)/3+a12*(TL(k)^4)/4)+a13*(TL(k)^5)/5+a14*(TL(k)^6)/6
% for higher temp
cp_high(k)= (a1+a2*TH(k)+a3*TH(k)^2+a4*TH(k)^3+a5*TH(k)^4)*R;
h_high(k)=(a1+a2*TH(k)/2+a3*(TH(k)^2+a4*(TH(k)^3)/4)+a5*(TH(k)^4)/5)+a6/TH(k)*R;
s_high(k)=(a1*log(TH(k))+a2*TH(k)+a3*(TH(k)^2)/2+a4*(TH(k)^3)/3+a5*(TH(k)^4)/4)+a6*(TH(k)^5)/5+a7*(TH(k)^6)/6
end
%creating folder
nasa_file=species
mkdir(nasa_file);
current=pwd;
cd(nasa_file)
% creating figure of specific heat
figure(1)
plot(TL,cp_low,'linewidth',3,'color','r');
hold on
plot(TH,cp_high,'linewidth',3,'color','g');
hold off
title(['temperature vs specific heat,species'])
xlabel('temperature')
ylabel('specific heat')
legend('Low temperature','High temperature')
%saving .png file of figure
fig=sprintf('specific heat of %s.png',species)
saveas(gca,fig)
% creating figure of enthalpy
figure(2)
plot(TL,h_low,'linewidth',3,'color','r')
hold on
plot(TH,h_high,'linewidth',3,'color','g');
hold off
title(['temperature vs enthalpy,species']);
xlabel('temperature')
ylabel('enthalpy');
legend('low temperature','high temperature');
% saving .png file of fig
fig=sprintf('enthalpy of %.png,species');
saveas(gca,fig);
% creating fig of entropy
figure(3)
plot(TL,s_low,'linewidth',3,'color','r');
hold on
plot(TH,s_high,'linewidth',3,'color','r');
hold off
title(['temperature vs entropy,species']);
xlabel('temperature');
ylabel('entropy');
legend('low temp','high temp')
% saving .png file of figure
fig=sprintf('entropy of %s.png',species)
saveas(gca,fig);
cd(current);
% calculation of molecular weight
atom=['H','C','O','N','Ar','S']
atomic_weight=[1.008 12.011 15.999 14.007 39.948 32.065]
mw=0
for p=1:length(species);
for q=1:length(atomic_weight);
if strcmp(species(p),atom(q));
mw=mw+atomic_weight(q);
position=q;
end
end
% used for multiplication into the number of atoms
n=str2double(species(p));
if n>1
mw=mw+atomic_weight(position)*(n-1);
end
end
%display molecular weight
fprintf('molecular weight of %s:%f\n',species,mw);
end
#EXPLANATION:
* used fopen command to open the given file format in the starting of the code.
* using'r' gave permission to read the file.
* after opening the file used fgetl command to red the first line.
* again fgetl command is used to read the second line of the file as temperature range which contains global low,mid and high temperature.
* next using strsplit command to store the temperature in cell array.
* str2num command is used to convert the temperature from string to number.
* temperature range is defined using linspace command hich in this codes takes place 100 values between global low and high temperature.
* the next 3 comments are skipped using fgetl command which runs under for loop.
* there are 53 species in the given file and for each species or element there are four lines in which first line given information about elementnae and under of atoms and temperature range .second,third and fourth line contains 14 coefficients in which first 7 are high temperatures cofficients and next 7 are low temperature coefficients.fgetl command is used to read 4 lines.
* next the character E is located using findstr command which shows the position of the coefficients in that line they are then stored in the variables a,b and c.
* the low temperature coefficients are stored in the varibles a1,a2,a3,a4,a5,a6,a7 and high temperature coefficients are stored in the variables b1,b2,b3,b4,b5,b6,b7. the coefficients run for 53 times under a for loop.
* the coefficients for each elements get stored in the variable for each run.
* calculating specific heat, enthalpy and entropy,molecular weight.
* the formula for calculation are as follows.
cp/R=a1+a2.T+a3.T^2+a4.T^3+a5.T^4;
h/RT=a1+a2.(T/2)+a3.(t^2/3)+a4.(T^3/4)+a5.(T^4/5)+a6/T;
s/R=a1.logT+a2.T+a3.(T^2/2)+a4.(T^3/3)+a5.(T^4/4)+a7;
* REFERANCEs used for the projects:
1. Google
2. skill-lync project website.
3. help forum of matlab.