/*
**  NL3.E
**
**
**  Chemical Equilibrium
**
**  This example solves 7 non-linear equations which
**  represent a methane - oxygen reaction.  Equations
**  1-3 represent the atom conservation balances for
**  Oxygen, Hydrogen and Carbon respectively.  Equation 4
**  gives the enthalpy balance, equation 5 the mole fraction
**  constraint, and equations 6 and 7 the equilibrium
**  relationships. (See Carnahan et. al. 1969, for a complete
**  description)
**
**  The task is to solve for the variables:
**
**      x1 - mole fraction of CO in equilibrium mixture
**      x2 - mole fraction of CO2 in equilibrium mixture
**      x3 - mole fraction of H2O in equilibrium mixture
**      x4 - mole fraction of H2 in equilibrium mixture
**      x5 - mole fraction of CH4 in equilibrium mixture
**      x6 - number of moles of O2 per mole of CH4 in feed gases
**      x7 - number of moles of product gases in equilibrium mixture
**                per mole of CH4 in the feed gases.
**
**  The equations to be solved are:
**
**      Eq 1:   0.5*x1 + x2 + 0.5*x3 - x6/x7 = 0
**      Eq 2:   x3 + x4 + 2*x5 -2/x7 = 0
**      Eq 3:   x1 + x2 + x5 - 1/x7 = 0
**      Eq 4:   -28837*x1 - 139009*x2 -78213*x3 + 18927*x4
**              + 8427*x5 + 13492/x7 - 10690*x6/x7  = 0
**      Eq 5:  x1 + x2 + x3 + x4 + x5 - 1 = 0
**      Eq 6:  (P^2)*x1*x4^3 - 1.7837e5*x3*x5 = 0
**      Eq 7:  x1*x3 - 2.6058*x2*x4 = 0;
**
**      where the operating pressure (P) is 20 atmospheres.
**
**  Solution:
**
**      The solution to the problem is
**
**      x1 - mole fraction CO           0.322871
**      x2 - mole fraction CO2          0.009224
**      x3 - mole fraction H2O          0.046017
**      x4 - mole fraction H2           0.618172
**      x5 - mole fraction CH4          0.003717
**      x6 - feed ratio O2/CH4          0.576715
**      x7 - total moles of product     2.977863
**
**
**
**  See Also:
**      Example NL3.E for a solution to this problem utilizing
**      analytic derivatives.
**
**  References:
**      Carnahan, Luther and Wilkes 1969,
**      Applied Numerical Methods (New York: Wiley), Example 5.5
*/


library nlsys;
nlset;

proc  F(x);
    local ff,p,eqns;
    ff = zeros(7,1);
    P = 20;
    ff[1] = 0.5*x[1] + x[2] + 0.5*x[3] - x[6]/x[7];
    ff[2] = x[3] + x[4] + 2*x[5] - 2/x[7];
    ff[3] = x[1] + x[2] + x[5] - 1/x[7];
    ff[4] = -28837*x[1] - 139009*x[2] - 78213*x[3]
            + 18927*x[4] + 8427*x[5] + 13492/x[7]
            - 10690*x[6]/x[7];
    ff[5] = x[1] + x[2] + x[3] + x[4] + x[5] - 1;
    ff[6] = (P^2)*x[1]*x[4]^3 - 1.7837e5*x[3]*x[5];
    ff[7] = x[1]*x[3] - 2.6058*x[2]*x[4];
    retp(ff);
endp;


let x0 = 0.5 0 0 0.5 0 0.5 2.0 ;
_nlalgr = 2;
_nlchpf = 1;
__altnam = { co, co2, h2o, h2, ch4, o2/ch4, total };

output file = nl3.out reset;

__title = "Chemical Equilibrium Problem";
{ x2,fvp,jc,tcode } = nlprt(nlsys(&f,x0));

output off;