from ChemPlugin import *

print("ChemPlugin example -- pH titration\n")

# Create ChemPlugin instance
cp = ChemPlugin("stdout")

# Configure the instance
cmds =  ("Ca++ = 1 mmol/kg",
         "Na+ = 1 mmol/kg",
         "Cl- = 3 mmol/kg",
         "HCO3- = 2 mmol/kg",
         "pH = 4")
        
cmds += ("react 3 mmol/kg NaOH", "delxi = 0.1")

for cmd in cmds:
    cp.Config(cmd)

# Initialize the instance
cp.Initialize()
spec = cp.Report("species")
nsp = len(spec)

# Time Marching loop
while True:
    deltat = cp.ReportTimeStep()
    if cp.AdvanceTimeStep(deltat): break
    if cp.AdvanceChemical(): break

    conc = cp.Report("concentration aqueous", "umol/kg")
    for i in range(nsp):
        if conc[i] >= 10.0:
            print("%-10s = %12.2f umol/kg" % (spec[i], conc[i]))
    print("")

input()