LECTURE5_chemical_equilibrium.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "16ba4d47-6066-48c0-af76-1f509a12c559",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "error: package statistics is not installed\n",
      "error: called from\n",
      "    load_packages at line 47 column 7\n",
      "    pkg at line 588 column 7\n"
     ]
    }
   ],
   "source": [
    "% some housekeeping stuff\n",
    "register_graphics_toolkit (\"gnuplot\");\n",
    "available_graphics_toolkits ();\n",
    "graphics_toolkit (\"gnuplot\")\n",
    "clear\n",
    "% load packages\n",
    "pkg load statistics\n",
    "% end of housekeeping"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "725f486a",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Chemical speciation modelling\n",
    "\n",
    " Solving for the equilibrium position of a set of simultaneous reactions subject to the constraints of mass balance and mass action is a common problem in environmental modeling.\n",
    "There exist many available computer programs to solve these types of chemical equilibrium\n",
    "problems, examples include MINEQL (Allison et al., 1991) and PHREEQ (Saini-Eidukat\n",
    "and Yahin, 1999). For engineering and scientific practices though these stand-alone pro-\n",
    "grams are not necessarily completely useful. Often, engineers and scientists will use so-called high-level programming languages such as Matlab or Scilab to implement their larger scale problems. Problems such as reactive transport, or simulation of industrial scale processes can be “coded” into Matlab in a very flexible task-specific manner.\n",
    "\n",
    "Many processes require the solution of simultaneous equilibria and thus the speciation\n",
    "code presented here was developed for Matlab. In the researcher’s lab this code is used\n",
    "to fit experimental data to chemical equilibrium binding constants for metal/organic and\n",
    "metal/phosphate systems but with slight modifications the code could be embedded in Mat-\n",
    "lab code for other applications (such as transport modeling).\n",
    "\n",
    "The Matlab code was written based on the paper by (Carrayrou et al., 2002). An example calculation for the Fe(III) system in water is presented below.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cae293-60d0-42cd-90c1-e7bac1552907",
   "metadata": {},
   "source": [
    "### Iron (III) system\n",
    "\n",
    "Mass balance and mass action for solving chemical equilibrium problems can all be represented in Tableau notation (Smith and Ferris, 2001; Morel and Hering, 1993). For a given defined system it is desirable to determine the equilibrium concentration of all species. Using the iron (III) system as an example a list of species of interest could include H+, OH – , Fe3+, FeOH2+, Fe(OH) +2 , and Fe(OH) –4 . There are other possible species but to keep the discussion simple these are the only ones included here. Also, for dilute solutions water can be assumed as a fixed component. The list includes six species so it is necessary to define six relationships in order to solve this system. The situation can be simplified if it is realized that each of these species are not independent. We can select components from the list of species and use those components to solve the equilibrium problem. For example, if we know H+ (pH) and Fe3+ we can determine the concentration of all other species from their logK values (mass action). In matrix notation we think of pH and Fe3+ as spanning the basis set.\n",
    "\n",
    "Now we need two equations and two unknowns. First relation is mass balance of iron and\n",
    "the second is proton balance (related to electroneutrality). Now we can write the tableau:\n",
    "\n",
    "\n",
    "| H+ | Fe3+ | logK | species |\n",
    "| --- | --- | --- | --- |\n",
    "| 1 | 0 | 0 | H$^+$ |\n",
    "| 0 | 1 | 0 | Fe$^{3+}$ |\n",
    "| -1 | 0 | -14 | OH$^-$ |\n",
    "| -1 | 1 | -2.19 | FeOH$^{2+}$ |\n",
    "| -2 | 1 | -5.67 | Fe(OH)$_2^+$ |\n",
    "| -4 | 1 | -21.6 | Fe(OH)$_4^-$ |\n",
    "||||\n",
    "| TOTH | FeT | |  |\n",
    "\n",
    "Table 1 completely defines the equilibrium problem. The entries in the columns are the\n",
    "stoichiometric coefficients required for formation of each species. For example, there is 1 H+ and 0 Fe3+ in H+. Also, Fe(OH) –4 is formed with one iron and removing 4 protons from\n",
    "water. Given the tableau, all that remains is to determine the values for [H+] and [Fe3+] for specified TOTH and FeT. If you multiply across rows it is possible to determine species concentration and if you sum down columns the total values (mass balance) are recovered. This is best understood by writing out some entries:\n",
    "\n",
    "\\begin{align}\n",
    "{[\\mathrm{H}+]}&={[\\mathrm{H}+]}^1\\times{[\\mathrm{Fe}^{3+}]}^0\\times10^0\\\\\n",
    "{[\\mathrm{Fe}^{3+}]}&={[\\mathrm{H}+]}^0\\times{[\\mathrm{Fe}^{3+}]}^1\\times10^0\\\\\n",
    "{[\\mathrm{OH}^-]}&={[\\mathrm{H}+]}^{-1}\\times{[\\mathrm{Fe}^{3+}]}^0\\times10^{-14} \\label{eqn:Kw}\\\\\n",
    "{[\\mathrm{FeOH}^{2+}]}&={[\\mathrm{H}+]}^{-1}\\times{[\\mathrm{Fe}^{3+}]}^1\\times10^{-2.19} \\label{eqn:KH1}\\\\\n",
    "&\\text{....}\n",
    "\\end{align}\n",
    "\n",
    "If equation Kw is rearranged to solve for $\\mathrm{[OH-]}$ then the equation above for [OH-] s obtained.  As another example, to understand equation formation of FeOH given above, consider that $K_{H1}$ for the first hydrolysis of iron (III) corresponds to the following reaction:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{Fe^{3+}} + \\mathrm{H_2O} \\rightleftharpoons \\mathrm{FeOH^{2+}} + \\mathrm{H+} \\quad K_{H1} = \\dfrac{\\mathrm{[FeOH^{2+}]}\\mathrm{[H+]}}{\\mathrm{[Fe^{3+}]}} \n",
    "\\end{equation}\n",
    "\n",
    "if K is rearranged to solve for $\\mathrm{[FeOH^{2+}]}$ then the equation $\\mathrm{[FeOH^{2+}]}={[\\mathrm{H}+]}^{-1}\\times{[\\mathrm{Fe}^{3+}]}^1\\times10^{-2.19}$  is obtained.  To explain how the summation down the columns is the mass balance.  Consider total iron ....\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathrm{Fe_T}=0\\mathrm{[H+]}+1\\mathrm{[Fe^{3+}]}+0\\mathrm{[OH^-]}+1\\mathrm{[FeOH^{2+}]}+1\\mathrm{[Fe(OH)2^+]}+1\\mathrm{[Fe(OH)4^-]}\n",
    "\\end{equation}\n",
    "\n",
    "notice that the coefficients are the entries down the iron column in the tableau.\n",
    "\n",
    "The problem can easily be expressed in matrix notation.  We'll define the $1 \\times 2$ vector of unknown component concentrations as $\\mathbf{X}$ so we can write $\\mathbf{X}=\\bigl( \\begin{smallmatrix}\n",
    "\\log\\mathrm{[H+]} & \\log\\mathrm{[Fe^{3+}]} \\end{smallmatrix} \\bigr)$.  There is a $6 \\times 1$ vector of species concentrations as well \\begin{equation}\n",
    "\\mathbf{C}= \\begin{pmatrix}     \n",
    "\\mathrm{[H+]} \\\\ \n",
    "\\mathrm{[Fe^{3+}]} \\\\\n",
    "\\mathrm{[OH-]} \\\\ \n",
    "\\mathrm{[FeOH^{2+}]} \\\\ \n",
    "\\mathrm{[Fe(OH)2+]} \\\\ \n",
    "\\mathrm{[Fe(OH)4-]}\n",
    "\\end{pmatrix}     \\end{equation}\n",
    "\n",
    "Total concentrations are in a $2\\times 1$ vector $\\mathbf{T} = \\bigl( \\begin{smallmatrix}\n",
    "\\text{TOTH} \\\\ \\text{Fe}_{\\mathrm{T}} \\end{smallmatrix} \\bigr)$.  The logK values are summarized in the $ 6\\times 1$ vector\n",
    "\\begin{equation}\n",
    " \\mathbf{K}=\\begin{pmatrix}\n",
    "0\\\\\n",
    "0 \\\\\n",
    "\\log{K_w}\\\\\n",
    "\\log{K_{H1}} \\\\\n",
    " \\log{K_{H2}}\\\\\n",
    "  \\log{K_{H4}}\\\\\n",
    "   \\end{pmatrix} \\end{equation}\n",
    "\n",
    "Finally, we need a $6 \\times 2$ matrix of stoichiometric coefficients\n",
    "\\begin{equation}\n",
    "\\mathbf{A} = \\begin{pmatrix} 1&0\\\\\n",
    "0&1\\\\\n",
    "-1 & 0\\\\\n",
    "-1 & 1\\\\\n",
    "-2 & 1\\\\\n",
    "-4 & 1\\\\ \n",
    " \\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "Now the minimization problem is to determine $\\mathbf{X}$ that minimizes the residuals in the mass balance.  This calculation is performed as follows:\n",
    "\\begin{align}\n",
    " \\text{minimize $\\mathbf{R}$ as a function of $\\mathbf{X}$ where}\\quad\n",
    " \\mathbf{R}&=\\mathbf{A}'\\times(10^\\mathbf{C})-\\mathbf{T}\\\\\n",
    "  \\text{and} \\quad      \\mathbf{C}&=10^{(\\mathbf{K}+\\mathbf{A}\\times\\mathbf{X}')}\n",
    "\\end{align}\n",
    "\n",
    "Minimization can be performed using all element of $\\mathbf{R}$ using Newton-Raphson method for example or using other nonlinear optimization methods on some summation of $\\mathbf{R}$ such as sum of squares or sum of absolute values of residuals (BrassardandBodurtha:2000).  In this case $\\mathbf{R}$ would be a $2 \\times 1$ vector of mass balance residuals.\n",
    "\n",
    "To solve this here first we will add some functions.  One to return the \"tableau\" and one to adjust for fixed pH (for this example we calculate speciation versus pH, so pH will be fixed).  The third function below solves for the minimum in the difference between calculated and actual mass balance.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "edc3e534",
   "metadata": {},
   "outputs": [],
   "source": [
    "function [KSOLUTION,ASOLUTION,SOLUTIONNAMES] = get_equilib_defn \n",
    "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
    "\n",
    "Tableau=[...\n",
    "%H\t  FeIII        logK            species\n",
    " 1\t  0\t            0              {'H'}\n",
    " 0\t  1\t            0              {'Fe'}\n",
    "-1\t  0\t            -14            {'OH'}\n",
    "-1    1             -2.19          {'FeOH'}\n",
    "-2    1             -5.67          {'FeOH2'}\n",
    "-4    1             -21.6          {'FeOH4'}\n",
    "];\n",
    "\n",
    "n=size(Tableau,2);\n",
    "ASOLUTION=cell2mat(Tableau(:,1:n-2));\n",
    "KSOLUTION=cell2mat(Tableau(:,n-1));\n",
    "SOLUTIONNAMES=strvcat(Tableau(:,n));\n",
    "\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6a779c35-c8cb-46ec-9cde-0a2305714500",
   "metadata": {},
   "outputs": [],
   "source": [
    "% ----------- for fixed pH ----------------\n",
    "\n",
    "function [Ksolution,Asolution]=get_equilib_fixed_pH(KSOLUTION,ASOLUTION,pH)\n",
    "\n",
    "    [N,M]=size(ASOLUTION);\n",
    "    Ksolution=KSOLUTION-ASOLUTION(:,1)*pH;\n",
    "    Asolution=[ASOLUTION(:,2:M)];\n",
    "\n",
    "end\n",
    "\n",
    "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "030ce273-23dd-4e1b-a46e-50cc511c7dfd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created file '/run/media/ssmith/T7/ssmith/scottsmith/teaching/CH480_numerical_methods_chemistry/OCTAVE/CH480CH680/nl_massbalancerrnosolid_NR.m'.\n"
     ]
    }
   ],
   "source": [
    "%%file nl_massbalancerrnosolid_NR.m\n",
    "\n",
    "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
    "function [X,F,J,C] = nl_massbalancerrnosolid_NR(X,Asolution,Ksolution,T)\n",
    "\n",
    "[Nc,Nx]=size(Asolution); %Xsolution=X(1:Nx);\n",
    "criteria=1e-15;\n",
    "\n",
    "for i=1:1000\n",
    "\n",
    "logC=(Ksolution)+Asolution*log10(X); C=10.^(logC); % calc species\n",
    "R=Asolution'*C-T;\n",
    "\n",
    "% Evaluate the Jacobian \n",
    "   z=zeros(Nx,Nx); \n",
    "for j=1:Nx \n",
    "\tfor k=1:Nx \n",
    "\t\tfor i=1:Nc; z(j,k)=z(j,k)+Asolution(i,j)*Asolution(i,k)*C(i)/X(k); end\n",
    "   \tend\n",
    "end\n",
    "\n",
    "J = z;\n",
    "\n",
    "deltaX=z\\(-1*R);\n",
    "one_over_del=max([1, -1*deltaX'./(0.5*X')]);\n",
    "del=1/one_over_del; X=X+del*deltaX;\n",
    "    \n",
    "tst=sum(abs(R));\n",
    "if tst<=criteria; break; end\n",
    "\n",
    "end\n",
    "\n",
    "F=[R]; \n",
    "\n",
    "end\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0155c95b-fd0c-4f6c-ac3d-754a56b29f3d",
   "metadata": {},
   "source": [
    "### make the graph\n",
    "\n",
    "So for each pH we will solve the equilibrium problem.  First load the tableau, then modify for the pH of that iternation. Then solve.  Then store the results and loop.\n",
    "\n",
    "At the end make a plot.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0c3d756a-6757-4a60-9e45-a19ac0d3342a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600px\" height=\"600px\" viewBox=\"0 0 560 420\">\n",
       "\n",
       "<title>Gnuplot</title>\n",
       "<desc>Produced by GNUPLOT 5.4 patchlevel 3 </desc>\n",
       "\n",
       "<g id=\"gnuplot_canvas\">\n",
       "\n",
       "<rect x=\"0\" y=\"0\" width=\"560\" height=\"420\" fill=\"none\"/>\n",
       "<defs>\n",
       "\n",
       "\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\" stroke=\"currentColor\"/>\n",
       "\t<path id=\"gpPt0\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M-1,0 h2 M0,-1 v2\"/>\n",
       "\t<path id=\"gpPt1\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M-1,-1 L1,1 M1,-1 L-1,1\"/>\n",
       "\t<path id=\"gpPt2\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\"/>\n",
       "\t<rect id=\"gpPt3\" stroke-width=\"0.222\" stroke=\"currentColor\" x=\"-1\" y=\"-1\" width=\"2\" height=\"2\"/>\n",
       "\t<rect id=\"gpPt4\" stroke-width=\"0.222\" stroke=\"currentColor\" fill=\"currentColor\" x=\"-1\" y=\"-1\" width=\"2\" height=\"2\"/>\n",
       "\t<circle id=\"gpPt5\" stroke-width=\"0.222\" stroke=\"currentColor\" cx=\"0\" cy=\"0\" r=\"1\"/>\n",
       "\t<use xlink:href=\"#gpPt5\" id=\"gpPt6\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<path id=\"gpPt7\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\"/>\n",
       "\t<use xlink:href=\"#gpPt7\" id=\"gpPt8\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<use xlink:href=\"#gpPt7\" id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\"/>\n",
       "\t<use xlink:href=\"#gpPt9\" id=\"gpPt10\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<use xlink:href=\"#gpPt3\" id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\"/>\n",
       "\t<use xlink:href=\"#gpPt11\" id=\"gpPt12\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<path id=\"gpPt13\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\"/>\n",
       "\t<use xlink:href=\"#gpPt13\" id=\"gpPt14\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<filter id=\"textbox\" filterUnits=\"objectBoundingBox\" x=\"0\" y=\"0\" height=\"1\" width=\"1\">\n",
       "\t  <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n",
       "\t  <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n",
       "\t</filter>\n",
       "\t<filter id=\"greybox\" filterUnits=\"objectBoundingBox\" x=\"0\" y=\"0\" height=\"1\" width=\"1\">\n",
       "\t  <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n",
       "\t  <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n",
       "\t</filter>\n",
       "</defs>\n",
       "<g fill=\"none\" color=\"white\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g stroke=\"none\" shape-rendering=\"crispEdges\">\n",
       "\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"64.14,362.40 534.81,362.40 534.81,18.02 64.14,18.02 \"/>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,362.40 L76.74,362.40 M534.82,362.40 L522.22,362.40  \"/>\t<g transform=\"translate(55.75,366.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">0</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,293.52 L76.74,293.52 M534.82,293.52 L522.22,293.52  \"/>\t<g transform=\"translate(55.75,297.42)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">0.2</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,224.64 L76.74,224.64 M534.82,224.64 L522.22,224.64  \"/>\t<g transform=\"translate(55.75,228.54)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">0.4</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,155.77 L76.74,155.77 M534.82,155.77 L522.22,155.77  \"/>\t<g transform=\"translate(55.75,159.67)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">0.6</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,86.89 L76.74,86.89 M534.82,86.89 L522.22,86.89  \"/>\t<g transform=\"translate(55.75,90.79)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">0.8</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,18.01 L76.74,18.01 M534.82,18.01 L522.22,18.01  \"/>\t<g transform=\"translate(55.75,21.91)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">1</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M64.14,362.40 L64.14,349.80 M64.14,18.01 L64.14,30.61  \"/>\t<g transform=\"translate(64.14,384.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">2</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M158.28,362.40 L158.28,349.80 M158.28,18.01 L158.28,30.61  \"/>\t<g transform=\"translate(158.28,384.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">4</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M252.41,362.40 L252.41,349.80 M252.41,18.01 L252.41,30.61  \"/>\t<g transform=\"translate(252.41,384.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">6</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M346.55,362.40 L346.55,349.80 M346.55,18.01 L346.55,30.61  \"/>\t<g transform=\"translate(346.55,384.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">8</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M440.68,362.40 L440.68,349.80 M440.68,18.01 L440.68,30.61  \"/>\t<g transform=\"translate(440.68,384.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">10</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M534.82,362.40 L534.82,349.80 M534.82,18.01 L534.82,30.61  \"/>\t<g transform=\"translate(534.82,384.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">12</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g transform=\"translate(19.18,190.21) rotate(-90)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">fraction of total Fe</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g transform=\"translate(299.48,411.30)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"12.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">pH</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"black\" d=\"M434.27,114.49 L434.27,27.01 L526.43,27.01 L526.43,114.49 L434.27,114.49 Z  \"/></g>\n",
       "\t<g id=\"gnuplot_plot_1a\"><title>Fe</title>\n",
       "<g fill=\"none\" color=\"white\" stroke=\"black\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g transform=\"translate(488.63,41.46)\" stroke=\"none\" fill=\"rgb(0,0,0)\" font-family=\"Arial\" font-size=\"10.80\" text-anchor=\"start\">\n",
       "\t\t<text><tspan font-family=\"Arial\">Fe</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0, 114, 189)\" d=\"M441.83,37.95 L481.07,37.95 M64.14,155.81 L68.85,175.91 L73.55,196.59 L78.26,217.29 L82.97,237.47 L87.67,256.61   L92.38,274.30 L97.09,290.25 L101.79,304.27 L106.50,316.34 L111.21,326.51 L115.91,334.89 L120.62,341.66 L125.33,347.02   L130.04,351.18 L134.74,354.35 L139.45,356.72 L144.16,358.44 L148.86,359.69 L153.57,360.56 L158.28,361.17 L162.98,361.58   L167.69,361.86 L172.40,362.05 L177.10,362.17 L181.81,362.25 L186.52,362.31 L191.22,362.34 L195.93,362.36 L200.64,362.38   L205.34,362.38 L210.05,362.39 L214.76,362.39 L219.46,362.40 L224.17,362.40 L228.88,362.40 L233.58,362.40 L238.29,362.40   L243.00,362.40 L247.71,362.40 L252.41,362.40 L257.12,362.40 L261.83,362.40 L266.53,362.40 L271.24,362.40 L275.95,362.40   L280.65,362.40 L285.36,362.40 L290.07,362.40 L294.77,362.40 L299.48,362.40 L304.19,362.40 L308.89,362.40 L313.60,362.40   L318.31,362.40 L323.01,362.40 L327.72,362.40 L332.43,362.40 L337.13,362.40 L341.84,362.40 L346.55,362.40 L351.25,362.40   L355.96,362.40 L360.67,362.40 L365.38,362.40 L370.08,362.40 L374.79,362.40 L379.50,362.40 L384.20,362.40 L388.91,362.40   L393.62,362.40 L398.32,362.40 L403.03,362.40 L407.74,362.40 L412.44,362.40 L417.15,362.40 L421.86,362.40 L426.56,362.40   L431.27,362.40 L435.98,362.40 L440.68,362.40 L445.39,362.40 L450.10,362.40 L454.80,362.40 L459.51,362.40 L464.22,362.40   L468.92,362.40 L473.63,362.40 L478.34,362.40 L483.05,362.40 L487.75,362.40 L492.46,362.40 L497.17,362.40 L501.87,362.40   L506.58,362.40 L511.29,362.40 L515.99,362.40 L520.70,362.40 L525.41,362.40 L530.11,362.40 L534.82,362.40  \"/>\t<use xlink:href=\"#gpDot\" x=\"461.45\" y=\"37.95\" color=\"rgb(  0, 114, 189)\"/>\n",
       "\n",
       "</g>\n",
       "\t</g>\n",
       "\t<g id=\"gnuplot_plot_2a\"><title>FeOH</title>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g transform=\"translate(488.63,63.33)\" stroke=\"none\" fill=\"rgb(0,0,0)\" font-family=\"Arial\" font-size=\"10.80\" text-anchor=\"start\">\n",
       "\t\t<text><tspan font-family=\"Arial\">FeOH</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(217,  83,  25)\" d=\"M441.83,59.82 L481.07,59.82 M64.14,229.02 L68.85,210.82 L73.55,192.73 L78.26,175.47 L82.97,159.79 L87.67,146.41   L92.38,135.96 L97.09,128.91 L101.79,125.61 L106.50,126.19 L111.21,130.65 L115.91,138.76 L120.62,150.15 L125.33,164.28   L130.04,180.51 L134.74,198.10 L139.45,216.30 L144.16,234.40 L148.86,251.80 L153.57,268.04 L158.28,282.81 L162.98,295.93   L167.69,307.37 L172.40,317.17 L177.10,325.46 L181.81,332.39 L186.52,338.12 L191.22,342.83 L195.93,346.67 L200.64,349.79   L205.34,352.30 L210.05,354.33 L214.76,355.96 L219.46,357.27 L224.17,358.31 L228.88,359.14 L233.58,359.81 L238.29,360.34   L243.00,360.76 L247.71,361.10 L252.41,361.36 L257.12,361.58 L261.83,361.75 L266.53,361.88 L271.24,361.99 L275.95,362.07   L280.65,362.14 L285.36,362.19 L290.07,362.24 L294.77,362.27 L299.48,362.30 L304.19,362.32 L308.89,362.34 L313.60,362.35   L318.31,362.36 L323.01,362.37 L327.72,362.38 L332.43,362.38 L337.13,362.39 L341.84,362.39 L346.55,362.40 L351.25,362.40   L355.96,362.40 L360.67,362.40 L365.38,362.40 L370.08,362.40 L374.79,362.40 L379.50,362.40 L384.20,362.40 L388.91,362.40   L393.62,362.40 L398.32,362.40 L403.03,362.40 L407.74,362.40 L412.44,362.40 L417.15,362.40 L421.86,362.40 L426.56,362.40   L431.27,362.40 L435.98,362.40 L440.68,362.40 L445.39,362.40 L450.10,362.40 L454.80,362.40 L459.51,362.40 L464.22,362.40   L468.92,362.40 L473.63,362.40 L478.34,362.40 L483.05,362.40 L487.75,362.40 L492.46,362.40 L497.17,362.40 L501.87,362.40   L506.58,362.40 L511.29,362.40 L515.99,362.40 L520.70,362.40 L525.41,362.40 L530.11,362.40 L534.82,362.40  \"/>\t<use xlink:href=\"#gpDot\" x=\"461.45\" y=\"59.82\" color=\"rgb(217,  83,  25)\"/>\n",
       "\n",
       "</g>\n",
       "\t</g>\n",
       "\t<g id=\"gnuplot_plot_3a\"><title>FeOH2</title>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g transform=\"translate(488.63,85.20)\" stroke=\"none\" fill=\"rgb(0,0,0)\" font-family=\"Arial\" font-size=\"10.80\" text-anchor=\"start\">\n",
       "\t\t<text><tspan font-family=\"Arial\">FeOH2</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(237, 177,  32)\" d=\"M441.83,81.69 L481.07,81.69 M64.14,357.98 L68.85,356.08 L73.55,353.50 L78.26,350.05 L82.97,345.55 L87.67,339.78   L92.38,332.55 L97.09,323.65 L101.79,312.93 L106.50,300.27 L111.21,285.66 L115.91,269.17 L120.62,251.01 L125.33,231.51   L130.04,211.11 L134.74,190.36 L139.45,169.80 L144.16,149.97 L148.86,131.32 L153.57,114.21 L158.28,98.84 L162.98,85.30   L167.69,73.58 L172.40,63.59 L177.10,55.17 L181.81,48.17 L186.52,42.38 L191.22,37.64 L195.93,33.78 L200.64,30.65   L205.34,28.12 L210.05,26.09 L214.76,24.46 L219.46,23.15 L224.17,22.11 L228.88,21.27 L233.58,20.61 L238.29,20.08   L243.00,19.67 L247.71,19.34 L252.41,19.09 L257.12,18.90 L261.83,18.77 L266.53,18.69 L271.24,18.68 L275.95,18.74   L280.65,18.91 L285.36,19.23 L290.07,19.78 L294.77,20.67 L299.48,22.11 L304.19,24.39 L308.89,27.94 L313.60,33.45   L318.31,41.81 L323.01,54.24 L327.72,72.09 L332.43,96.50 L337.13,127.77 L341.84,164.63 L346.55,204.05 L351.25,242.07   L355.96,275.24 L360.67,301.74 L365.38,321.47 L370.08,335.39 L374.79,344.85 L379.50,351.11 L384.20,355.19 L388.91,357.82   L393.62,359.49 L398.32,360.56 L403.03,361.24 L407.74,361.67 L412.44,361.94 L417.15,362.11 L421.86,362.22 L426.56,362.28   L431.27,362.33 L435.98,362.35 L440.68,362.37 L445.39,362.38 L450.10,362.39 L454.80,362.39 L459.51,362.40 L464.22,362.40   L468.92,362.40 L473.63,362.40 L478.34,362.40 L483.05,362.40 L487.75,362.40 L492.46,362.40 L497.17,362.40 L501.87,362.40   L506.58,362.40 L511.29,362.40 L515.99,362.40 L520.70,362.40 L525.41,362.40 L530.11,362.40 L534.82,362.40  \"/>\t<use xlink:href=\"#gpDot\" x=\"461.45\" y=\"81.69\" color=\"rgb(237, 177,  32)\"/>\n",
       "\n",
       "</g>\n",
       "\t</g>\n",
       "\t<g id=\"gnuplot_plot_4a\"><title>FeOH4</title>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g transform=\"translate(488.63,107.07)\" stroke=\"none\" fill=\"rgb(0,0,0)\" font-family=\"Arial\" font-size=\"10.80\" text-anchor=\"start\">\n",
       "\t\t<text><tspan font-family=\"Arial\">FeOH4</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(126,  47, 142)\" d=\"M441.83,103.56 L481.07,103.56 M64.14,362.40 L68.85,362.40 L73.55,362.40 L78.26,362.40 L82.97,362.40 L87.67,362.40   L92.38,362.40 L97.09,362.40 L101.79,362.40 L106.50,362.40 L111.21,362.40 L115.91,362.40 L120.62,362.40 L125.33,362.40   L130.04,362.40 L134.74,362.40 L139.45,362.40 L144.16,362.40 L148.86,362.40 L153.57,362.40 L158.28,362.40 L162.98,362.40   L167.69,362.40 L172.40,362.40 L177.10,362.40 L181.81,362.40 L186.52,362.40 L191.22,362.40 L195.93,362.40 L200.64,362.40   L205.34,362.40 L210.05,362.40 L214.76,362.40 L219.46,362.40 L224.17,362.40 L228.88,362.40 L233.58,362.39 L238.29,362.39   L243.00,362.38 L247.71,362.37 L252.41,362.36 L257.12,362.34 L261.83,362.30 L266.53,362.24 L271.24,362.15 L275.95,362.00   L280.65,361.76 L285.36,361.39 L290.07,360.80 L294.77,359.87 L299.48,358.40 L304.19,356.11 L308.89,352.53 L313.60,347.01   L318.31,338.63 L323.01,326.19 L327.72,308.34 L332.43,283.93 L337.13,252.65 L341.84,215.79 L346.55,176.36 L351.25,138.34   L355.96,105.17 L360.67,78.67 L365.38,58.95 L370.08,45.02 L374.79,35.56 L379.50,29.30 L384.20,25.22 L388.91,22.59   L393.62,20.92 L398.32,19.85 L403.03,19.17 L407.74,18.74 L412.44,18.47 L417.15,18.30 L421.86,18.19 L426.56,18.13   L431.27,18.08 L435.98,18.06 L440.68,18.04 L445.39,18.03 L450.10,18.02 L454.80,18.02 L459.51,18.01 L464.22,18.01   L468.92,18.01 L473.63,18.01 L478.34,18.01 L483.05,18.01 L487.75,18.01 L492.46,18.01 L497.17,18.01 L501.87,18.01   L506.58,18.01 L511.29,18.01 L515.99,18.01 L520.70,18.01 L525.41,18.01 L530.11,18.01 L534.82,18.01  \"/>\t<use xlink:href=\"#gpDot\" x=\"461.45\" y=\"103.56\" color=\"rgb(126,  47, 142)\"/>\n",
       "\n",
       "</g>\n",
       "\t</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb( 38,  38,  38)\" d=\"M64.15,362.39 L534.83,362.39 M64.15,18.00 L534.83,18.00 M64.15,362.39 L64.15,18.00 M534.83,362.39 L534.83,18.00    \"/></g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"rgb( 38,  38,  38)\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "</g>\n",
       "</svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%plot -s 600,600 -f 'svg'\n",
    "\n",
    "%define the equilibrium sytem (The tableau)\n",
    "[KSOLUTION,ASOLUTION,SOLUTIONNAMES] = get_equilib_defn;\n",
    "% initial guess (just Fe becasue pH is fixed)\n",
    "Feguess=[-5.5]; guess=[10.^Feguess];\n",
    "\n",
    "%set the pH range and total\n",
    "pH=2:0.1:12; FeT=1e-4; T=[FeT];\n",
    "\n",
    "%now for each pH solve\n",
    "for i=1:length(pH)\n",
    "% adjust for fixed pH\n",
    "[Ksolution,Asolution]=get_equilib_fixed_pH(KSOLUTION,ASOLUTION,pH(i));\n",
    "% calculate species using NR\n",
    "[X,F,J,C] = nl_massbalancerrnosolid_NR(guess,Asolution,Ksolution,T);\n",
    "species_summary(:,i)=C; err(i)=F;\n",
    "end\n",
    "\n",
    "for i=1:size(species_summary,1)\n",
    "txt=[SOLUTIONNAMES(i,:),\"=species_summary(i,:);\"]; eval(txt);\n",
    "end\n",
    "\n",
    "h=plot(pH,Fe/T,pH,FeOH/T,pH,FeOH2/T,pH,FeOH4/T);\n",
    "set(gca,\"fontsize\",12); set(h,'linewidth',2); set(gca,'linewidth',2);\n",
    "h=xlabel('pH'); set(h,'fontsize',12); h=ylabel('fraction of total Fe'); set(h,'fontsize',12);\n",
    "legend('Fe','FeOH','FeOH2','FeOH4','Location','northeast','Orientation','vertical');\n",
    "axis([min(pH) max(pH) 0 1])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "7d6b9f15-db43-4267-8d31-2e8eec29ea9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600px\" height=\"600px\" viewBox=\"0 0 560 420\">\n",
       "\n",
       "<title>Gnuplot</title>\n",
       "<desc>Produced by GNUPLOT 5.4 patchlevel 3 </desc>\n",
       "\n",
       "<g id=\"gnuplot_canvas\">\n",
       "\n",
       "<rect x=\"0\" y=\"0\" width=\"560\" height=\"420\" fill=\"none\"/>\n",
       "<defs>\n",
       "\n",
       "\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\" stroke=\"currentColor\"/>\n",
       "\t<path id=\"gpPt0\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M-1,0 h2 M0,-1 v2\"/>\n",
       "\t<path id=\"gpPt1\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M-1,-1 L1,1 M1,-1 L-1,1\"/>\n",
       "\t<path id=\"gpPt2\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\"/>\n",
       "\t<rect id=\"gpPt3\" stroke-width=\"0.222\" stroke=\"currentColor\" x=\"-1\" y=\"-1\" width=\"2\" height=\"2\"/>\n",
       "\t<rect id=\"gpPt4\" stroke-width=\"0.222\" stroke=\"currentColor\" fill=\"currentColor\" x=\"-1\" y=\"-1\" width=\"2\" height=\"2\"/>\n",
       "\t<circle id=\"gpPt5\" stroke-width=\"0.222\" stroke=\"currentColor\" cx=\"0\" cy=\"0\" r=\"1\"/>\n",
       "\t<use xlink:href=\"#gpPt5\" id=\"gpPt6\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<path id=\"gpPt7\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\"/>\n",
       "\t<use xlink:href=\"#gpPt7\" id=\"gpPt8\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<use xlink:href=\"#gpPt7\" id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\"/>\n",
       "\t<use xlink:href=\"#gpPt9\" id=\"gpPt10\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<use xlink:href=\"#gpPt3\" id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\"/>\n",
       "\t<use xlink:href=\"#gpPt11\" id=\"gpPt12\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<path id=\"gpPt13\" stroke-width=\"0.222\" stroke=\"currentColor\" d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\"/>\n",
       "\t<use xlink:href=\"#gpPt13\" id=\"gpPt14\" fill=\"currentColor\" stroke=\"none\"/>\n",
       "\t<filter id=\"textbox\" filterUnits=\"objectBoundingBox\" x=\"0\" y=\"0\" height=\"1\" width=\"1\">\n",
       "\t  <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n",
       "\t  <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n",
       "\t</filter>\n",
       "\t<filter id=\"greybox\" filterUnits=\"objectBoundingBox\" x=\"0\" y=\"0\" height=\"1\" width=\"1\">\n",
       "\t  <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n",
       "\t  <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n",
       "\t</filter>\n",
       "</defs>\n",
       "<g fill=\"none\" color=\"white\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<g stroke=\"none\" shape-rendering=\"crispEdges\">\n",
       "\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"71.31,384.00 534.81,384.00 534.81,18.02 71.31,18.02 \"/>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M71.31,384.00 L83.91,384.00 M534.82,384.00 L522.22,384.00  \"/>\t<g transform=\"translate(62.92,387.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">-1e-18</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M71.31,292.50 L83.91,292.50 M534.82,292.50 L522.22,292.50  \"/>\t<g transform=\"translate(62.92,295.75)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">-5e-19</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M71.31,201.00 L83.91,201.00 M534.82,201.00 L522.22,201.00  \"/>\t<g transform=\"translate(62.92,204.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">0</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M71.31,109.51 L83.91,109.51 M534.82,109.51 L522.22,109.51  \"/>\t<g transform=\"translate(62.92,112.76)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">5e-19</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M71.31,18.01 L83.91,18.01 M534.82,18.01 L522.22,18.01  \"/>\t<g transform=\"translate(62.92,21.26)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"end\">\n",
       "\t\t<text><tspan font-family=\"Arial\">1e-18</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M71.31,384.00 L71.31,371.40 M71.31,18.01 L71.31,30.61  \"/>\t<g transform=\"translate(71.31,405.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">2</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M164.01,384.00 L164.01,371.40 M164.01,18.01 L164.01,30.61  \"/>\t<g transform=\"translate(164.01,405.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">4</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M256.71,384.00 L256.71,371.40 M256.71,18.01 L256.71,30.61  \"/>\t<g transform=\"translate(256.71,405.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">6</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M349.42,384.00 L349.42,371.40 M349.42,18.01 L349.42,30.61  \"/>\t<g transform=\"translate(349.42,405.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">8</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M442.12,384.00 L442.12,371.40 M442.12,18.01 L442.12,30.61  \"/>\t<g transform=\"translate(442.12,405.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">10</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0,   0,   0)\" d=\"M534.82,384.00 L534.82,371.40 M534.82,18.01 L534.82,30.61  \"/>\t<g transform=\"translate(534.82,405.25)\" stroke=\"none\" fill=\"rgb(38,38,38)\" font-family=\"Arial\" font-size=\"10.00\" text-anchor=\"middle\">\n",
       "\t\t<text><tspan font-family=\"Arial\">12</tspan></text>\n",
       "\t</g>\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "\t<g id=\"gnuplot_plot_1a\"><title>gnuplot_plot_1a</title>\n",
       "<g fill=\"none\" color=\"white\" stroke=\"rgb(  0,   0,   0)\" stroke-width=\"0.50\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"0.50\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb(  0, 114, 189)\" d=\"M71.31,220.85 L75.95,205.97 L80.58,166.28 L85.22,178.68 L89.85,188.60 L94.49,208.45 L99.12,205.97 L103.76,186.12   L108.39,181.16 L113.03,188.60 L117.66,223.33 L122.30,210.93 L126.93,191.08 L131.57,168.76 L136.20,225.81 L140.84,213.41   L145.47,198.52 L150.11,201.00 L154.74,198.52 L159.38,193.56 L164.01,178.68 L168.65,193.56 L173.28,208.45 L177.92,201.00   L182.55,210.93 L187.19,181.16 L191.82,181.16 L196.46,205.97 L201.09,171.24 L205.73,193.56 L210.36,196.04 L215.00,188.60   L219.63,161.32 L224.27,153.88 L228.90,143.96 L233.54,156.36 L238.17,205.97 L242.81,210.93 L247.44,163.80 L252.08,213.41   L256.71,255.57 L261.35,196.04 L265.98,163.80 L270.62,198.52 L275.25,250.61 L279.89,243.17 L284.52,253.09 L289.16,193.56   L293.79,245.65 L298.43,255.57 L303.07,196.04 L307.70,270.45 L312.34,253.09 L316.97,210.93 L321.61,218.37 L326.24,188.60   L330.88,201.00 L335.51,210.93 L340.15,215.89 L344.78,223.33 L349.42,223.33 L354.05,176.20 L358.69,235.73 L363.32,228.29   L367.96,94.36 L372.59,139.00 L377.23,235.73 L381.86,329.97 L386.50,245.65 L391.13,250.61 L395.77,317.57 L400.40,250.61   L405.04,196.04 L409.67,265.49 L414.31,297.73 L418.94,258.05 L423.58,240.69 L428.21,267.97 L432.85,292.77 L437.48,253.09   L442.12,245.65 L446.75,196.04 L451.39,270.45 L456.02,208.45 L460.66,163.80 L465.29,213.41 L469.93,265.49 L474.56,223.33   L479.20,168.76 L483.83,193.56 L488.47,215.89 L493.10,275.41 L497.74,225.81 L502.37,230.77 L507.01,203.49 L511.64,72.04   L516.28,168.76 L520.91,119.16 L525.55,111.72 L530.18,188.60 L534.82,213.41  \"/></g>\n",
       "\t</g>\n",
       "<g fill=\"none\" color=\"white\" stroke=\"rgb(  0, 114, 189)\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"2.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"black\" stroke-width=\"0.50\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "\t<path stroke=\"rgb( 38,  38,  38)\" d=\"M71.32,383.99 L534.83,383.99 M71.32,18.00 L534.83,18.00 M71.32,383.99 L71.32,18.00 M534.83,383.99 L534.83,18.00    \"/></g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"rgb( 38,  38,  38)\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "<g fill=\"none\" color=\"black\" stroke=\"currentColor\" stroke-width=\"1.00\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\">\n",
       "</g>\n",
       "</g>\n",
       "</svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(pH,err)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a48c35-e3d7-4992-a688-22decf996578",
   "metadata": {},
   "source": [
    "### exercises 7\n",
    "\n",
    "recalculate the speciation above but at 4e-5 M added salicylic acid, H2L (There is already 1e-4 M Fe(III).\n",
    "\n",
    "Fe + L = FeL logK=17.55\n",
    "\n",
    "pKa1=2.97 pKa2=13.7\n",
    "\n",
    "plot species versus pH, and also a plot of the two mass balance errors versus pH.\n",
    "\n",
    "for matrix dimensions to match, you will need to put totals in like this T=[FeT; LT];\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f311d8d-a84c-45cb-81e2-ce26f35aa399",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "6.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}