function f = funkeval(x)
%
%  these two lines force a column vector of length n
%
n = max(size(x));
f = zeros(n,1);
%
%  enter the functions here
%
%
%  adiabatic chemical reaction equilibria calculation
%  nitrogen plus hydrogen to ammonia 
%  (Sandler Illustration 9.1-4 page 509)
%
%
%  Checklist.
%		(1)  Make sure you have entered the correct number of reactions (nr)
%		(2)  Make sure you have entered the correct number of components (nc)
%     (3)  Set adiabatic = 1 for adiabatic system, set adiabatic = 0 for isothermal system
%     (4)  If the system is isothermal, set the isothermal temperature , (Tiso)
%     (5)  Enter initial mole amounts of each component (flowin)
%     (6)  Enter initial temperature of each component (Tin)
%     (7)  Enter reactor pressure in atmospheres (pressure)
%     (8)  Enter the stoichiometric matrix
%     (9)  Enter the heats of formation  (Kcal/mol)
%     (10)  Enter the free energies of formation (Kcal/mol)
%     (11)  Enter the heat capacity constants 
%
% STEP ONE.  ENTER PROBLEM SPECIFICATIONS
%
% 	(1)  number of independent chemical reactions
%
nr = 1;
%
% 	(2)number of components
%
nc = 3;
%
%  (3)  determines whether the system is adiabatic or isothermal
%
adiabatic = 0;
%adiabatic = 1;
%
%  (4)	isothermal temperature
%
Tiso = 520;
%
%  (5)  inlet flowrates
%
flowin=zeros(nc,1);
flowin(1) = 0.5;
flowin(2) = 1.5;
%
% (6)  inlet temperatures (K)
%
Tin=zeros(nc,1);
TempC = 25;
Tin(1) = TempC+273.1;
Tin(2) = TempC+273.1;
%
%  (7)  reactor pressure in atmospheres
%
pressure = 1;  % put this in atmospheres
%
%  (8)  stoichiometric matrix
%  the order is 
%  benzene ethylne ethyl benzene
%
nu=[-1/2	-3/2	1];
%
%  (9) heats of formation at the reference temperature
%  from Sandler (Kcal/mol)
%  benzene ethylne ethyl benzene
%
Hf = [0;	0;	-10.960];
%
%  (10)  free energies of formation at the reference temperature
%  from Sandler (Kcal/mol)
%  benzene ethylne ethyl benzene
%
Gf = [0.0;	0.0;	-3.903];
%
%  (11)  heat capacity constants from Sandler Appendix
%  Each column contains parameters for a given component.
%  The four rows contain the four parameters a, b, c, d for 
%  Cp = 4.184*(a + b*T/10^2 + c*T^2/10^5 +d*T^3/10^9) [Joules/mole/K]
%  benzene ethylne ethyl benzene
%
Cpcon = [6.903		6.952			6.5846	
   	-0.03753		-0.04576		0.61251
      0.1930	 	0.09563		0.23663
      -0.6861		-0.2079		-1.5981];
%
%     you won't need to touch anything below here.   
%  
%
%  STEP TWO:  LET THE PROGRAM DO ITS WORK
%
%
%  STEP TWO A:  Identify unknowns
%
%
% nr extents of reaction
%
X=zeros(nr,1);
for i = 1:nr
   X(i) = x(i);
end
%
% 1 Temperature
%
if (adiabatic == 1) 
   T = x(nr+1);
else
   T=Tiso;
end
%
%  STEP TWO B.  DEFINE KNOWNS
%
Hf = Hf*4.184*1000; % (Joules/mol)
Gf = Gf*4.184*1000; % (Joules/mol)
%
%  reference temperature
%
Tref = 25 + 273.1;
%
%  constants
%
R = 8.314; % [J/mol/k]
RT = R*T;
%
%
%
%  heats of reaction at the reference temperature
%
DHrref = nu*Hf;
%
%  free energies of reaction at at the reference temperature
%
DGrref = nu*Gf;
%
%  scale heat capacity constants 
%  
for k = 1:nc
   Cpcon(2,k) = Cpcon(2,k)*10^-2;   
   Cpcon(3,k) = Cpcon(3,k)*10^-5;   
  	Cpcon(4,k) = Cpcon(4,k)*10^-9;   
end
Cpcon = 4.184*Cpcon;
%Cp = 4.184*(Cpcon(1,2) + Cpcon(2,2)*T/10^2 + Cpcon(3,2)*T^2/10^5 +Cpcon(4,2)*T^3/10^9)
%
%  enthalpies at arbitrary T
%  (need this for flow out of the reactor term in energy balance)
%
for k = 1:nc
   Cpint2(k) = Cpint(T,Cpcon,k) - Cpint(Tref,Cpcon,k);
end
%
%  heats of reaction at arbitrary T
%  (need this for energy balance)
%
for i = 1:nr
   DHr(i) = DHrref(i); 
   for k = 1:nc
      DHr(i) = DHr(i) + nu(i,k)*Cpint2(k);
   end
end
%
%  calculate equilibrium constants at the reference temperature 
%
for i = 1:nr
   Karef(i) = exp(-DGrref(i)/(R*Tref));
end
%
% integrate the heat of reaction over RT^2 from Tref to T
%  (need this to calculate Ka as a function of Temperature)
%
for i = 1:nr
   term1 = DHroRT2intfunk(T,Tref,Cpcon,nu,DHrref,R,i,nc);
   term2 = DHroRT2intfunk(Tref,Tref,Cpcon,nu,DHrref,R,i,nc);
	DHroRT2int2(i) = term1 - term2;
end
%
%  calculate equilibrium constants at arbitrary  temperature 
%
for i = 1:nr
   Ka(i) = Karef(i)*exp(DHroRT2int2(i));
end
%
%  define molar composition based on extent of reactions
%
xE = flowin;
for i = 1:nr
   for k = 1:nc
      xE(k) = xE(k) + nu(i,k)*X(i);
   end
end
E = sum(xE);
xE = xE/E;
%
%  stream enthalpies
%
heatin = zeros(nc,1);
for k = 1:nc
   heatin(k) = flowin(k)*(Cpint(Tin(k),Cpcon,k) - Cpint(Tref,Cpcon,k));
end   
heatintot = sum(heatin);
heatout = zeros(nc,1);
for k = 1:nc
   heatout(k) = E*xE(k)*(Cpint(T,Cpcon,k) - Cpint(Tref,Cpcon,k));
end   
heatouttot = sum(heatout);
%
%  Step Three.  Write down nr+1 equations
%
%
% nr equilibrium contraints
%
for i = 1:nr
   f(i) = 1;
   for k = 1:nc
      f(i) = f(i)*(xE(k)^nu(i,k));
   end
   f(i) = f(i) - (pressure^(-sum(nu(i,:))))*Ka(i);
end
%
% nc material balances
%
%fmole =  flowin + nu'*X - E*xE;
%for k = 1:nc
%   f(k+nr) = fmole(k);
%end
%
% 1 mole fraction constraint
%
%f(nr+nc+1) = sum(xE) - 1;
%
% 1 energy balance
%
convaid = 0.01;
heatrxn = 0.0;
for i = 1:nr
   heatrxn = heatrxn + DHr(i)*X(i);
end
if (adiabatic == 1) 
   f(nr+1) = convaid*(heatintot - heatouttot - heatrxn);
end

%%%fprintf (1,'xE = %15e  %15e  %15e \n', xE(1), xE(2), xE(3))








%
%  integrated heat capacity of component k
%
function y = Cpint(T,Cpcon,k)
%Cp = 4.184*(a + b*T/10^2 + c*T^2/10^5 +d*T^3/10^9) [Joules/mole/K]
y = Cpcon(1,k)*T + Cpcon(2,k)*T^2/2 + Cpcon(3,k)*T^3/3 + Cpcon(4,k)*T^4/4;
%[Joules/mole]

%
% integrate the heat of reaction over RT^2 from Tref to T
%
function y = DHroRT2intfunk(T,Tref,Cpcon,nu,DHrref,R,i,nc)
term1a = -DHrref(i)/(R*T);
term1b = 0.0;
for k = 1:nc
   t1b1 = Cpcon(1,k)*log(T) + Cpcon(2,k)/2*T + Cpcon(3,k)/6*T^2 + Cpcon(4,k)/12*T^3;
   t1b2 = Cpcon(1,k)*Tref + Cpcon(2,k)*Tref^2/2 + Cpcon(3,k)*Tref^3/3 + Cpcon(4,k)*Tref^4/4;
   term1b = term1b + nu(i,k)/R*(t1b1 + t1b2/T);
end
y = term1a + term1b;