(* chemBalance, a package for balancing chemical reactions,
by Frank Kampas *)
(* Version 1.0, January, 2000 *)
(* This package requires Mathematica version 3.0 or later *)
(* The chemical reaction is entered as an argument of the function balance.
The reaction must be in the form of a string. Inside the string,
the reactants and products are separated by the "equilibrium sign" entered
as < escape > equi < escape > . The molecules on each side must have plus signs between them.
Each atom in the molecule must have a subscript indicating the number
of atoms in the molecule, even if there is only one of that atom.
There should be no parentheses in the molecules.
If there are atoms on one side of the equation that do not appear on
the other, an error message "atoms don't match" will be produced.
See chemBalanceExamples.nb for examples *)
(* Please send comments to fkampas@msn.com *)
BeginPackage["chemBalance`"]
balance::usage = "balance[eq_String] balances chemical reactions"
Begin["`Private`"]
stringDivide[s1_String, s2_String] :=
Module[
{p = Flatten[
StringPosition[s1,
s2] /. {x_Integer, y_Integer} -> {x - 1, y + 1}],
l = StringLength[s1]},
StringTake[s1, #] & /@ Partition[Join[{1}, p, {l}], 2, 2]
]
(* divide string s1 into lists of strings at symbol s2 *)
parse\[UnderBracket]equation[s_String] :=
stringDivide[#, "+"] & /@ stringDivide[s, "\[Equilibrium]"] (*
divide equation into left - hand and right -
hand sides and at plus signs *)
convert[s_List] := (ToExpression /@ StringReplace[s, "I" -> "i"]) /.
Times -> Plus /. Subscript[x_, n_] -> n*x
(* convert I to i to avoid confusion with square root of minus 1,
convert molecule to sum of atoms *)
atoms[a_] := List @@ (Plus @@ a) /. n_Integer*x_ -> x
(* find the atoms in a molecule *)
stringPlus[x__List] :=
Insert[x, " + ", Range[2, Length[x]] /. y_?NumberQ -> {y}] (*
insert plus signs between elements in a list *)
stringDot[a_List, b_List] :=
StringJoin @ Inner[StringJoin[##] &, a, b, stringPlus[{##}] &];(*
inner product for two lists of strings *)
balance[eq_String] :=
Module[
{eq1, reactants, products, numr, nump, rcoeff, pcoeff, rc, pc, r1, p1, \
r2,
p2, soln, coeff, mult, lhs, rhs},
eq1 = parse\[UnderBracket]equation[eq]; (* parse equation *)
{reactants, products} = eq1; (* find reactants and products *)
{numr, nump} = Length /@ eq1; (*
find number of reactants and number of products *)
rcoeff = Array[rc, numr];
pcoeff = Array[pc, nump]; (* create arrays for coefficients *)
{r1, p1} = convert /@ {reactants, products};(*
convert molecules to sums of atoms *)
r2 = rcoeff.r1; p2 = pcoeff.p1;(* multiply coefficients by atoms *)
{ratms, patms} = atoms /@ {r1, p1}; (*
find atoms for reactants and products *)
If[ratms =!= patms, Return["atoms don't match"]]; (*
check if atoms match *)
soln = SolveAlways[{r2 == p2, rc[1] == 1}, ratms][[1]];(*
solve for coefficients *)
coeff = {rcoeff, pcoeff} /. soln; (* set coefficients from solution *)
mult = LCM @@ Denominator /@ Flatten[coeff];(*
find least common multiple of denominators of coefficients *)
coeff = mult*coeff; (* clear denominators *)
coeff = Map[ToString, coeff, {-1}]; (* convert coefficients to strings *)
coeff = Map[StringReplace[#, "1" -> ""] &, coeff, {-1}]; (*
remove unit coefficients *)
lhs = stringDot[coeff[[1]], reactants]; (* left - hand - side of result \
*)
rhs = stringDot[coeff[[2]], products];(* right - hand - side of result *)
StringJoin[lhs, " \[Equilibrium] ", rhs](* final answer *)
]
End[]
EndPackage[]