from reaktoro import *

db = Database("supcrt98.xml")

editor = ChemicalEditor(db)

editor.addAqueousPhaseWithElements("H O Na Cl C Ca Mg Si")

editor.addGaseousPhase("H2O(g) CO2(g)")

editor.addMineralPhase("Halite")
editor.addMineralPhase("Calcite")
editor.addMineralPhase("Magnesite")
editor.addMineralPhase("Dolomite")
editor.addMineralPhase("Quartz")

system = ChemicalSystem(editor)

problem = EquilibriumProblem(system)

problem.setTemperature(70, "celsius")
problem.setPressure(100, "bar")
problem.add("H2O", 1.0, "kg")
problem.add("CO2", 2.0, "mol")
problem.add("NaCl", 1.0, "mol")
problem.add("CaCO3", 10.0, "g")
problem.add("MgCO3", 5.0, "g")
problem.add("Quartz", 1.0, "mol")

state = equilibrate(problem)

T = state.temperature()
P = state.pressure()
n = state.speciesAmounts()

print(f"T = {T} K")
print(f"P = {P} Pa")
print(f"n (in mol) = \n{n}")

print("Species names : n (in mol)")
for species in system.species():
    name = species.name()
    amount = state.speciesAmount(name)
    print(f"{name:>13} = {amount}")

state.output("state.txt")

properties = ChemicalProperties(system)

properties.update(T, P, n)

properties = state.properties()

lna = properties.lnActivities().val

a = numpy.exp(lna)
print("Species names : activities")
for i, species in enumerate(system.species()):
    print(f"{species.name():>13} : {a[i]}")

mu = properties.chemicalPotentials().val
print("Species names : potentials (in kJ/mol)")
for i, species in enumerate(system.species()):
    print(f"{species.name():>13} : {mu[i]}")

evaluate_pH = ChemicalProperty.pH(system)
pH = evaluate_pH(properties)

print(f"The pH of the aqueous phase is {pH.val}.")
print(f"Its sensitivity with respect to speciation is:")
print("Species names : ∂(pH)/∂n (in 1/mol)")
for i, species in enumerate(system.species()):
    print(f"{species.name():>13} : {pH.ddn[i]}")