# Import the reaktoro Python package
from reaktoro import *

# Initialize a thermodynamic database
db = Database("supcrt98.xml")

# Define the chemical system
editor = ChemicalEditor(db)
editor.addAqueousPhaseWithElements(
    "H O Na Cl C"
)  # add aqueous phase by the all possible combination of provided elements
editor.addGaseousPhase(["CO2(g)"])  # add one gaseous species

# Construct the chemical system
system = ChemicalSystem(editor)

print(system)

import numpy as np

print(f"List with {system.numSpecies()} species:")
for species in system.species():
    print(species.name())

print(f"List with {system.numPhases()} phases:")
for phase in system.phases():
    print(phase.name())

print(f"List with {(system.numElements())} elements:")
for element in system.elements():
    print(element.name())

print("List of phases with number of species in each phase:")
for phase in system.phases():
    print(f" * Phase {phase.name()} contains {phase.numSpecies()} species")

matrix = system.formulaMatrix()
print(f"Formula matrix of the size {matrix.shape[0]} x {matrix.shape[1]}:\n", matrix)

print("Index of the element H: ", system.indexElement("H"))
print("Index of the phase Aqueous: ", system.indexPhase("Aqueous"))
print("Index of the species Cl-: ", system.indexSpecies("Cl-"))

species = [
    "Cl-",
    "ClO-",
    "ClO2-",
    "ClO3-",
    "ClO4-",
    "HCl(aq)",
    "HClO(aq)",
    "HClO2(aq)",
    "NaCl(aq)",
]
print("Indices of species with Cl: ", system.indicesSpecies(species))

n = np.ones(system.numSpecies())
elements_amount = system.elementAmounts(n)
hydrogen_amount = system.elementAmount(2, n)
print(
    "Element amounts (in mol) provided 1 molal for all species: ",
    elements_amount,
)
print(
    "Hydrogen amounts (in mol) provided 1 molal for all species: ",
    hydrogen_amount,
)

T = 60
P = 100
thermo_properties = system.properties(T, P)

print("List of standard partial molar Gibbs energies of the species:")
for energies, species in zip(
    thermo_properties.standardPartialMolarGibbsEnergies().val,
    system.species()
):
    print(f"\u03B4G ({species.name():>10}) = {energies}")

print("List of standard partial molar enthalpies of the species:")
for enthalpies, species in zip(
    thermo_properties.standardPartialMolarEnthalpies().val,
    system.species()
):
    print(f"\u03B4H ({species.name():>10}) = {enthalpies}")

chemical_properties = system.properties(T, P, n)

print("Chemical potentials of the species:")
for potential, species, index in zip(
    chemical_properties.chemicalPotentials().val,
    system.species(),
    list(range(1, system.numSpecies()+1))
):
    print(f"\u03BC_{index} ({species.name():>10}) = {potential}")

print("Logarithms of activities of the species:")
for activity, species, index in zip(
    chemical_properties.lnActivities().val,
    system.species(),
    list(range(1, system.numSpecies()+1))
):
    print(f"ln (a_{index}) ({species.name():>10}) = {activity}")