Indian Institute of Technology Madras

Department of Aerospace Engineering

Assignment Report
(Physical Gas Dynamics)

Course No. AS5300

Submitted by:

Akash B

AE17M031
BACKGROUND
The aim of this assignment is to calculate the equilibrium properties of air at different
temperatures ranging from 300 K to 8000 K. We need two intensive properties to fix the
equilibrium state of a gas mixture. In this case, since the temperature is being varied we will
keep the pressure as a constant. Thus, what is obtained is the variation in properties of air kept
at a constant pressure as the temperature is increased from 300 K to 8000 K.

PROCEDURE
We consider air at temperatures of upto 8000 K as being composed of essentially five chemical
species : O2, N2,O,N and NO. We also consider the ionization of NO. NO is the first species to
undergo ionization as the temperature is increased since it has the lowest ionization energy.
Thus the equilibrium is governed by the following four reactions:

𝑂2 ⇌ 𝑂 + 𝑂

𝑁2 ⇌ 𝑁 + 𝑁

𝑁𝑂 ⇌ 𝑁 + 𝑂

𝑁𝑂 ⇌ 𝑁𝑂+ + 𝑒
Thus, there are 7 species in total. The concentration of each of these species at a particular
temperature is what we need to solve for to find the equilibrium properties at that
temperature. Instead of molar concentration we use a more convenient variable which is the
number of moles per unit mass of the mixture, called the mole-mass ratio denoted by the
symbol ǹ. This is related to the concentration by [χ] =  ǹ. Now, the law of mass action in terms
of mole-mass ratio can be written for each of the above four equations as follows:

ǹ2𝑂 1
= 𝐾𝑐,1 (𝑇)
ǹ𝑂2 𝜌

ǹ2𝑁 1
= 𝐾𝑐,2 (𝑇)
ǹ𝑁2 𝜌

ǹ𝑁 ǹ𝑂 1
= 𝐾𝑐,3 (𝑇)
ǹ𝑁𝑂 𝜌
ǹ𝑒 ǹ𝑁𝑂+ 1
= 𝐾𝑐,4 (𝑇)
ǹ𝑁𝑂 𝜌
The above equations can be written in terms of mole fractions and partition functions from
which we can derive the expressions for the equilibrium constants.

[𝑂]2 (𝑄0 )2
= exp⁡(−𝐷02 /𝐾𝑇)
[𝑂2 ] 𝑁̂ 𝑉𝑄00

[𝑁]2 (𝑄𝑁 )2
= exp⁡(−𝐷𝑁2 /𝐾𝑇)
[𝑁2 ] 𝑁̂ 𝑉𝑄𝑁𝑁

[𝑂2 ][𝑁2 ] 𝑄0 𝑄𝑁
= exp⁡(−𝐷𝑁𝑂 /𝐾𝑇)
[𝑁𝑂] ̂ 𝑉𝑄𝑁𝑂
𝑁
+
[𝑒][𝑁𝑂+ ] 𝑄𝑒 𝑄𝑁𝑂
= exp⁡(−𝐼/𝐾𝑇)
[𝑁𝑂] ̂ 𝑉𝑄𝑁𝑂
𝑁
The dissociation energy for the first three reactions and the ionization energy for the fourth
reaction are as follows:

𝐷𝑂2 =8.2x10−19 J

𝐷𝑁2 =1.56x10−18 J

𝐷𝑁𝑂 =1.04x10−18 J

𝐼 =1.482x10−18 J

The partition functions are obtained as follows:

Using the factorization property of partition functions:

𝑎 𝑎
𝑄𝑎 = 𝑄𝑡𝑟 𝑄𝑒𝑙
𝑎𝑎 𝑎𝑎 𝑎𝑎 𝑎𝑎
𝑄𝑎𝑎 = 𝑄𝑡𝑟 𝑄𝑟𝑜𝑡 𝑄𝑣𝑖𝑏 𝑄𝑒𝑙
The translational partition function for A atoms is given by:
3
𝑎
2𝜋𝑚𝑘𝑇 2
𝑄𝑡𝑟 = 𝑉 ( )
ℎ2
In 𝑄 𝑎𝑎 m will be replaced by 2m.

Rotational partition function is given by :

𝑎𝑎
1 𝑇
𝑄𝑟𝑜𝑡 = ( )
𝜎 𝜃𝑟
Where σ=1 for heteronuclear and σ=2 for homonuclear molecules and 𝜃𝑟 is the characteristic
rotational temperature

⁡𝜃𝑟 = 2.1⁡𝐾, 2.9⁡𝐾⁡𝑎𝑛𝑑⁡2.5 K for O2,N2and NO respectively. NO+ has the same value as NO.

Vibrational partition function is given by:

𝑎𝑎
1
𝑄𝑣𝑖𝑏 = −𝜃𝑣
1−𝑒 𝑇

Where 𝜃𝑣 is the characteristic vibrational temperature

𝜃𝑣 = 2270⁡𝐾, 3390⁡𝐾⁡𝑎𝑛𝑑⁡2740 K for O2,N2 and NO respectively. NO+ has the same value as

NO.

Electronic partition functions for the different species are as follows:

228 326
(− ) (− )
𝑄𝑒𝑙 𝑂 = 5 + 3𝑒 𝑇 +𝑒 𝑇

174
𝑁𝑂 (− )
𝑄𝑒𝑙 =2+ 2𝑒 𝑇
11390
(− )
𝑄𝑒𝑙 𝑂2 = 3 + 2𝑒 𝑇

𝑄𝑒𝑙 𝑁2 = 1

𝑄𝑒𝑙 𝑁 = 4

𝑄𝑒𝑙 𝑒− = 2
In addition to the above four equations we need three more equations to find the seven
unknowns. Two of these are given by the conservation of oxygen and nitrogen nuclei and the
third one is given by the condition of charge neutrality.
At 300 K, the mole fraction of N2 and O2 in air can be approximately taken to be 0.788 and 0.21
𝑃𝑡𝑜𝑡 𝜒
respectively. The mole-mass ratio at 300 K can then be found using ǹ =
𝜌𝑅𝑢 𝑇

Where 𝑃𝑡𝑜𝑡 = 1.01325 × 105 Pa

𝜌 = 1.173⁡𝑘𝑔/𝑚3

And 𝑅𝑢 = 8.314⁡𝐽/𝑚𝑜𝑙𝐾
𝑚𝑜𝑙 𝑚𝑜𝑙
This gives (ǹ𝑂2 ) = 7.33 𝑎𝑛𝑑⁡(ǹ𝑁2 ) = 27.3 ⁡
0 𝑘𝑔 0 𝑘𝑔

2ǹ𝑂2 +ǹ𝑂 +ǹ𝑁𝑂 + ǹ𝑁𝑂+ = 2(ǹ𝑂2 )

0

2ǹ𝑁2 +ǹ𝑁 +ǹ𝑁𝑂 + ǹ𝑁𝑂+ =2⁡(ǹ𝑁2 )

0

ǹ𝑁𝑂+ = ǹ𝑒

The above set of 7 non linear equations can be solved numerically. fsolve function in MATLAB
has been used to solve the equations at different temperatures. On solving, the mole-mass
ratio of the seven species at different temperatures is obtained. These can then be used to
compute various thermodynamic properties of air at different temperatures. The property
values thus calculated are tabulated and plotted below.
 RESULTS AND DISCUSSION

a)

Fig 1. Variation of specific heat at constant pressure with temperature for three different pressures 1
bar, .01 bar and 10 bar
 b)

Fig 2. Variation of specific heat at constant volume with temperature for three different pressures 1
bar, .01 bar and 10 bar
 c)

Fig 3. Variation of Helmholtz free energy with temperature for three different pressures 1 bar, .01 bar
and 10 bar
d)

Fig 4. Variation of specific energy with temperature for three different pressures 1 bar, .01 bar and 10
bar
e)

Fig 5. Variation of specific enthalpy with temperature for three different pressures 1 bar, .01 bar and
10 bar
 f)

Fig 6. Variation of specific entropy with temperature for three different pressures 1 bar, .01 bar and 10
bar
g)

Fig 7. Variation of density with temperature for three different pressures 1 bar, 5 bar and 10 bar
h)
Fig 8. Variation of concentration in terms of mole fraction of different species with temperature
at 1 bar, .01 bar and 10 bar respectively
Temperature Density Specific Specific Specific Specific Specific Ratio of
(in Kelvin) (in Kg/m^3) Energy Enthalpy Entropy heat at heat at Specific
(in Kj/kg) (in Kj/kg) (in Kj/kg) Constant Constant Heats
Pressure Volume
(in Kj/kgK) (in
Kj/kgK)
300 1.157754 212.4111 298.7852 25.0572 1.009711 0.721797 1.398885
400 0.868315 284.9214 400.0869 25.26867 1.018409 0.730495 1.394135
500 0.694652 358.6325 502.5894 25.43483 1.034315 0.746402 1.385736
600 0.578877 434.1907 606.939 25.57368 1.055385 0.767471 1.375146
700 0.49618 511.9774 713.5171 25.69435 1.078642 0.790728 1.364112
800 0.434158 592.0983 822.4294 25.80191 1.101748 0.813834 1.353775
900 0.385918 674.469 933.5915 25.89936 1.123335 0.835421 1.344633
1000 0.347326 758.9057 1046.82 25.98868 1.142804 0.85489 1.336785
1100 0.315751 845.188 1161.893 26.0712 1.160025 0.872111 1.330134
1200 0.289438 933.095 1278.592 26.14794 1.175118 0.887204 1.324518
1300 0.267174 1022.424 1396.711 26.21965 1.188306 0.900392 1.319765
1400 0.24809 1112.995 1516.074 26.28695 1.199845 0.91193 1.31572
1500 0.231551 1204.654 1636.525 26.35035 1.209976 0.922062 1.31225
1600 0.217079 1297.272 1757.934 26.41026 1.21892 0.931006 1.30925
1700 0.20431 1390.737 1880.191 26.46705 1.22686 0.938946 1.306635
1800 0.192959 1484.957 2003.202 26.52102 1.233951 0.946037 1.304337
1900 0.182803 1579.853 2126.889 26.57243 1.240319 0.952406 1.302302
2000 0.173663 1675.357 2251.185 26.62151 1.246069 0.958155 1.300488
2100 0.165393 1771.411 2376.031 26.66846 1.251284 0.96337 1.298861
2200 0.157876 1867.966 2501.376 26.71345 1.256032 0.968118 1.297395
2300 0.150978 1964.973 2627.323 26.76239 1.260639 0.972415 1.2964
2400 0.144576 2062.375 2754.054 26.82287 1.264831 0.976253 1.295597
2500 0.138788 2160.167 2880.691 26.86372 1.268244 0.979895 1.294265
2600 0.133448 2258.309 3007.664 26.90258 1.271509 0.983243 1.29318
2700 0.128113 2356.716 3137.278 27.0295 1.287418 0.984733 1.307378
2800 0.123283 2455.162 3266.304 27.11093 1.295188 0.985783 1.313867
2900 0.1187 2553.764 3396.223 27.21147 1.304314 0.986169 1.322606
3000 0.11432 2652.38 3527.114 27.33119 1.314782 0.985911 1.333571
3100 0.11011 2750.917 3659.103 27.47474 1.326245 0.984509 1.347113
3200 0.106035 2849.243 3792.326 27.64409 1.339154 0.981749 1.364049
3300 0.102093 2947.263 3926.763 27.84248 1.348931 0.978328 1.378813
3400 0.09825 3044.82 4062.63 28.06951 1.362636 0.973113 1.400286
3500 0.094493 3141.799 4200.075 28.32751 1.406609 0.963612 1.459726
3600 0.090838 3238.12 4338.985 28.61472 1.41502 0.957106 1.478435
3700 0.087293 3333.717 4479.286 28.92767 1.409573 0.951941 1.480736
3800 0.08391 3428.637 4620.395 29.25425 1.419174 0.944923 1.501894
3900 0.080688 3522.908 4762.251 29.58902 1.419652 0.939838 1.510528
4000 0.077661 3616.678 4904.32 29.91939 1.420429 0.935571 1.518248
4100 0.074843 3710.111 5046.249 30.23611 1.416903 0.933237 1.518267
4200 0.072241 3803.401 5187.656 30.53032 1.410414 0.932732 1.512131
4300 0.069853 3896.718 5328.305 30.79684 1.401427 0.933964 1.500515
4400 0.067668 3990.211 5468.02 31.03269 1.391949 0.936271 1.486696
4500 0.065667 4083.967 5606.797 31.23878 1.382829 0.939197 1.472352
4600 0.063829 4178.025 5744.711 31.418 1.375023 0.942218 1.459347
4700 0.062129 4272.377 5881.921 31.57487 1.368869 0.945029 1.448494
4800 0.060547 4366.989 6018.603 31.7143 1.3647 0.947343 1.440557
4900 0.059061 4461.805 6154.959 31.84131 1.362508 0.949035 1.435677
5000 0.057656 4556.76 6291.183 31.96053 1.362241 0.950047 1.433867
5100 0.056317 4651.784 6427.462 32.07622 1.363796 0.950347 1.435051
5200 0.055031 4746.806 6563.976 32.19216 1.367084 0.949927 1.439146
5300 0.053788 4841.753 6700.894 32.31171 1.372039 0.948782 1.446105
5400 0.052581 4936.553 6838.381 32.43793 1.378623 0.946903 1.455928
5500 0.051402 5031.131 6976.599 32.57361 1.386828 0.944272 1.468674
5600 0.050244 5125.411 7115.71 32.72136 1.39666 0.940866 1.484441
5700 0.049102 5219.315 7255.879 32.88366 1.408145 0.936654 1.503377
5800 0.047973 5312.76 7397.273 33.06287 1.421315 0.931601 1.525669
5900 0.046857 5404.858 7539.02 33.25764 1.435421 0.925036 1.551747
6000 0.045741 5497.077 7683.318 33.47732 1.452529 0.918546 1.581335
6100 0.044627 5588.604 7829.404 33.72061 1.471359 0.911081 1.614961
6200 0.043514 5679.338 7977.449 33.98946 1.491868 0.902587 1.652881
6300 0.042401 5769.176 8127.619 34.28565 1.513999 0.89304 1.695332
6400 0.041287 5858.012 8280.072 34.61075 1.537665 0.882441 1.742513
6500 0.040173 5945.741 8434.956 34.96607 1.562735 0.870821 1.794555
6600 0.03906 6032.265 8592.403 35.35259 1.589017 0.85824 1.851483
6700 0.037957 6115.555 8750.106 35.76451 1.617212 0.845454 1.912832
6800 0.036851 6199.478 8913.088 36.21518 1.645308 0.8315 1.978722
6900 0.035748 6283.742 9081.065 36.70255 1.672074 0.81592 2.049312
7000 0.034664 6364.66 9249.522 37.21414 1.699716 0.800905 2.122244
7100 0.033595 6444.072 9420.705 37.75383 1.726396 0.785864 2.196811
7200 0.032547 6521.992 9594.487 38.3188 1.751437 0.77112 2.271289
7300 0.031524 6598.462 9770.674 38.90536 1.774123 0.757028 2.343536
7400 0.03053 6673.567 9948.997 39.50896 1.793729 0.743958 2.411063
7500 0.029571 6747.424 10129.12 40.12431 1.809576 0.732276 2.471166
7600 0.02865 6820.211 10310.66 40.74554 1.820986 0.722241 2.521299
7700 0.02777 6892.06 10493.1 41.36618 1.82774 0.714311 2.558743
7800 0.026934 6963.214 10675.99 41.97979 1.829466 0.708577 2.58189
7900 0.026144 7033.899 10858.83 42.57999 1.826174 0.705154 2.589751
 8000 0.025401 7104.344 11041.12 43.1608 1.818021 0.703981 2.582487

Table 1: Variation of different properties of air with temperature at 1 bar

 The variations in specific heat at constant pressure and specific heat at constant
volume are sharper at lower temperatures.
 Both internal energy and enthalpy increases almost linearly with temperature at all
pressures.
 Entropy increases in steps with the entropy remaining nearly constant for certain
temperature ranges.
 Density decreases with increase in temperature. Higher density is obtained at higher
pressures for a given temperature
 O2 and N2 are stable at lower temperatures of upto 2500 K.
 The concentration of O2 decreases to trace values at around 4000 K indicating the
dissociation of O2 into O is near completion.
 Also at 4000 K, the concentration of N starts increasing appreciably. The dissociation
of N2 into N approaches completion near 8000 K.
 The concentrations of NO, e- and NO+ do not increase appreciably in the
temperature range considered. This indicates that ionization is not predominant in
this temperature range.

CONCLUSION
The properties of equilibrium air at different temperatures was calculated by solving the
reaction rate equations and the conservation of nuclei equations simultaneously. The
equilibrium constants were obtained in terms of partition functions. The properties of air were
then calculated using the respective equations in terms of partition functions and the results
were plotted an tabulated fe different pressures. The temperature ranges at which dissociation
of O2 and N2 became significant were found out. Also, it was found that ionization does not
have a significant impact on species concentration for the temperature range considered.
APPENDIX

a) Main program

clc;
clear all;
pconst=[1e5 1e3 1e6]; %three different pressures used

T=300:10:8020; %range of temperatures

Y=transpose(T); %transposed temperature vector for plotting
Y1=Y(1:end-1);
c=zeros(773,7,3); %concentration of different species for each pressure
et=zeros(772,3); %specific energy for each pressure
cp=zeros(771,3); %specific heat at constant pressure for each pressure
st=zeros(772,3); %specific entropy for each pressure
ht=zeros(772,3); %specific enthalpy for each pressure
ft=zeros(772,3); %specific Helmholtz free energy for each pressure
cv=zeros(771,3); %specific heat at constant volume for each pressure
g=zeros(771,3); %ratio of specific heats for each pressure
rho=zeros(772,3);

for j=1:3

P=pconst(j);
[n, c,et,st,cv,cp,ht,ft,g,rho ] = concentrationsolver( P,T,j,et,cp,cv,st,ht,c,ft,g,rho );

figure(1)
plot(Y1,et(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Specific Energy(in kJ/kg)')
hold on;
linkdata on

figure(2)
plot(Y1,ht(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Specific Enthalpy(in kJ/kg)')
hold on;
linkdata on

figure(3)
plot(Y1,st(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Specific Entropy(in kJ/kgK)')
hold on;
linkdata on
figure(4)
plot(Y1(1:end-1),cp(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Specific heat at constant pressure(in kJ/kgK)')
hold on;
linkdata on

figure(5)
plot(Y1(1:end-1),cv(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Specific heat at constant volume(in kJ/kgK)')
hold on;
linkdata on

figure(6)
plot(Y1(1:end-1),g(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Ratio of specific heats')
hold on;
linkdata on

figure(7)
plot(Y1,ht(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Helmholtz free energy(in kJ/kg)')
hold on;
linkdata on

figure(8)
plot(Y1,rho(:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Density(in kg/m^3)')
hold on;

figure(j+8)
plot(Y,c(:,:,j))
xlim([300 8000])
xlabel('Temperature(in degree K)')
ylabel('Concentration in terms of mole fraction(%)')

end
%to tabulate the values corresponding 1 bar
vect1=transpose(Y(1):100:Y(end-2));to
vect2=real(rho(1:10:end-1,1));
vect3=real(et(1:10:end-1,1));
vect4=real(ht(1:10:end-1,1));
vect5=real(st(1:10:end-1,1));
vect6=real(cp(1:10:end,1));
vect7=real(cv(1:10:end,1));
vect8=real(g(1:10:end,1));

table=[vect1(:),vect2(:),vect3(:),vect4(:),vect5(:),vect6(:),vect7(:),vect8(:)];

b) Function to solve for mole mass ratios using fsolve

function [n, c,et,st,cv,cp,ht,ft,g,rho ] = concentrationsolver( P,T,j,et,cp,cv,st,ht,c,ft,g,rho)

Navo=6.022*10^23; %Avogadro's number

k=1.38*10^-23; %Boltzmann's constant
hp=6.626*10^-34; %Planck's constant
r=8.314; %Universal gas constant

%molecular masses of species in the order O,N,O2,N2,NO,e,NO+

m=[2.65*10^-26 2.33*10^-26 2*2.65*10^-26 2*2.33*10^-26 2.65*10^-26+2.33*10^-26 9.1*10^-31
2.65*10^-26+2.33*10^-26];

thetav=[inf inf 2270 3390 2740 inf 2740]; %vibrational characteristic temperature of species
x0=[1e-100 1e-100 .00733 .0273 1e-100 1e-100 ]; %initial condition
%dissociation energies
Doo=8.2*10^-19;
Dnn=1.56*10^-18;
Dno=1.04*10^-18;

%ionization energy
I=1.482*10^-18;

qe=zeros(7);
q=zeros(771,7);
n=zeros(771,7);

for i=1:1:length(T)

thetar=[T(i) T(i) 2.1/2 2.9/2 2.5 T(i) 2.5]; %rotational characteristic temp.
qt=(((2*pi*k*T(i))/(hp^2))^1.5)*(m.^1.5); %translational partition function
qr=T(i)./thetar; %rotational partition function
qv=1./(1-exp(-thetav/T(i))); %vibrational partition function
%electronic partition functions
qe(1)=5+3*exp(-228/T(i))+exp(-326/T(i));
qe(2)=4;
qe(3)=3+2*exp(-11390/T(i));
qe(4)=1;
qe(5)=2+2*exp(-174/T(i));
qe(6)=2;
qe(7)=2+2*exp(-174/T(i));
%total partition functions
q(i,1)=qt(1)*qr(1)*qv(1)*qe(1);
q(i,2)=qt(2)*qr(2)*qv(2)*qe(2);
q(i,3)=qt(3)*qr(3)*qv(3)*qe(3);
q(i,4)=qt(4)*qr(4)*qv(4)*qe(4);
q(i,5)=qt(5)*qr(5)*qv(5)*qe(5);
q(i,6)=qt(6)*qr(6)*qv(6)*qe(6);
q(i,7)=qt(7)*qr(7)*qv(7)*qe(7);

ro=P/(T(i)*r); %not the actual density

A=(1/(Navo*ro))*((q(i,1)^2)/(q(i,3)))*exp(-Doo/(k*T(i))); %equilibrium constants
B=(1/(Navo*ro))*((q(i,2)^2)/(q(i,4)))*exp(-Dnn/(k*T(i)));
C=(1/(Navo*ro))*((q(i,1)*q(i,2))/(q(i,5)))*exp(-Dno/(k*T(i)) );
D=(1/(Navo*ro))*((q(i,6)*q(i,7))/(q(i,5)))*exp(-I/(k*T(i)));

options=optimset('TolFun',1e-50,'TolX',1e-50,'MaxFunEvals',5000000,'MaxIter',500000);
x = fsolve(@(x)rootfinder( x, A,B,C,D),x0,options); %finding mole mass ratios using fsolve
x(7)=x(6); %charge neutrality
n(i,:)=x;

end
[n, c,et,st,cv,cp,ht,ft,g,rho ] = propertycalculator(
Navo,k,q,n,T,j,et,cp,cv,st,ht,c,ft,g,rho,P); %calculating properties

end

c) Function to calculate properties

d) function [n, c,et,st,cv,cp,ht,ft,g,rho ] = propertycalculator(

Navo,k,q,n,T,j,et,cp,cv,st,ht,c,ft,g,rho,P )

Y=transpose(T);
Y1=Y(1:end-1);
e=zeros(772,7);
s=zeros(772,7);
f=zeros(772,7);
ctemp=zeros(773,7);

for i=1:7
ctemp(:,i)=(n(:,i)./sum(n,2))*100;
end;
c(:,:,j)=ctemp; %concentration in percentage mole fraction
for i=1:7

dq=diff(log(q(:,i)))./diff(Y); %differentiating q with respect to temp.

e(:,i)=Navo*k*(n(1:end-1,i).*(Y1.^2.*dq));
lq=log(q(1:end-1,i)./n(1:end-1,i)); %log q/N
s(:,i)=Navo*k*(n(1:end-1,i).*(lq+1+Y1.*dq));
f(:,i)=-Navo*k*(n(1:end-1,i).*(Y1.*(lq+1)));
 end
et(:,j)=sum(e,2);
st(:,j)=sum(s,2);
cv(:,j)=diff(et(:,j))./diff(Y1);
ft(:,j)=sum(f,2);
nt=sum(n,2);
rho(:,j)=(P*(8.314*(nt(1:end-1).*Y1)).^(-1))/1000;
ht(:,j)=et(:,j)+8.314*(nt(1:end-1).*Y1);
cp(:,j)=diff(ht(:,j))./diff(Y1);
g(:,j)=cp(:,j)./cv(:,j);
end

e) function containing set of equations for fsolve

function F =rootfinder( x, A,B,C,D)

F(1) =x(1)^2-sum(x)*A*x(3);
F(2) =x(2)^2-sum(x)*B*x(4);
F(3) =x(1)*x(2)-sum(x)*C*x(5);
F(4)=x(6)^2-sum(x)*D*x(5);
F(5)=2*x(3)+x(1)+x(5)+x(6)-2*.00733;
F(6)=2*x(4)+x(2)+x(5)+x(6)-2*.0273;

end