function  keq=kequil(t2,j)
% kequil	- returns the equilibrium constants for gas phase reactions
% function keq=kequil(t2,j)
% Argument list
% Argument    Gives
%   t2      temperature in units of Tdeg
%    j      indices of the reactions that will be at equilibrium
% Function return the variable: keq,
% for reactions specified and will put a zero for reactions in
% which equil. const. is _not_ solved for. If j is omitted,
% equil. constant for all reactions in stoic will be determined.
% Example : keq1=kequil(550,[1 3])
% LDM

global stoic hcpv Gibb
T2=at(t2);
a=[];  c=[];
szcf = size(stoic);
keq = zeros(szcf(1),1);
if nargin == 1
  j = [1:szcf(1)];
  szj = size(j);
  for w = 1:szcf(1)
     for i = 1:szcf(2)
        a(i,:) = stoic(w,i) .* hcpv(i,:);
     end
     c = [c ; sum(a)];
  end
else
  szj = size(j);
  for w = 1:szj(2)
     for i = 1:szcf(2)
        a(i,:) = stoic(j(1,w),i) .* hcpv(i,:);
     end
     c = [c ; sum(a)];
  end
end

R = 0.00831451;
T = 298.15;

for r = 1:szj(2);
   c1 = c(r,1);c2 = c(r,2);c3 = c(r,3);
   c4 = c(r,4);c5 = c(r,5);c6 = c(r,6);
   G = stoic(j(1,r),:) .* Gibb';
   dG = sum(G);
   lnK = -dG/(R*T);
   cnst298 = (c1/4)*T^4 + (c2/3)*T^3 + (c3/2)*T^2 + c4*T + c5*log(T) - c6/T;
   I = (R*lnK - cnst298 )/R;
   cnstT2 =(c1/4)*T2^4 + (c2/3)*T2^3 + (c3/2)*T2^2 + c4*T2 + c5*log(T2) - c6/T2;
   lnK2 = (cnstT2)/R + I;
   keq(j(1,r),1) = exp(lnK2);
end
