#  # Carbon Deposition in a C-H-O Ternary Diagram
#  
#  This example calculates the carbon deposition boundaries for gas mixtures
#  of the elements carbon, hydrogen and oxygen and shows them in a ternary
#  diagram. On the carbon-rich side of a boundary the formation of solid carbon
#  is thermodynamically favored.

import gaspype as gp
import numpy as np
import matplotlib.pyplot as plt
import ternary


#  Setting temperatures and pressure. The temperatures get an additional
#  dimension, so they are broadcast against the compositions later on:

t_range = np.array([500, 600, 700, 800])  # °C
t = t_range[:, None] + 273.15  # K

p = 1e5  # Pa


#  The search for the boundaries below runs through arbitrary atomic compositions.
#  Oxygen (`O2`) and atomic carbon gas (`C`) are part of the fluid system to
#  make sure that an equilibrium composition exists for each of them, also with
#  an excess of oxygen or carbon:

fs = gp.fluid_system(['H2', 'H2O', 'CO2', 'CO', 'CH4', 'O2', 'C'])


#  The mixtures are defined by their atomic composition. Therefore we create
#  an `elements` object with one mol of atoms for each of the three elements:

carbon = gp.elements({'C': 1}, fs)
hydrogen = gp.elements({'H': 1}, fs)
oxygen = gp.elements({'O': 1}, fs)


#  Each mixture is described by its molar fraction of carbon `xc` and by the
#  hydrogen fraction of the remaining atoms `h_ratio = H / (H + O)`. The following
#  function returns the molar fractions of hydrogen and oxygen:

h_ratio = np.linspace(0.01, 0.99, 50)


def get_h_o_fractions(xc):
    xh = (1 - xc) * h_ratio
    xo = 1 - xc - xh
    return xh, xo


#  `gp.carbon_activity` calculates the activity of solid carbon (graphite) for
#  the equilibrium state of a mixture. The shapes of the mixture, the temperature
#  and the pressure are broadcast against each other, so a single call returns
#  the activities of all compositions for all temperatures:

def get_carbon_activity(xc):
    xh, xo = get_h_o_fractions(xc)
    mix = carbon * xc + hydrogen * xh + oxygen * xo
    return gp.carbon_activity(mix, t, p)


#  At the carbon deposition boundary the carbon activity is 1. For each hydrogen
#  ratio and each temperature the carbon fraction of the boundary is searched by
#  bisection: The carbon fraction is increased if the activity is below 1 and
#  decreased otherwise, while the step size is halved with every iteration.
#  All boundaries are searched at once:

xc = np.full((len(t_range), len(h_ratio)), 0.5)
step_size = 0.5

while step_size > 1e-5:
    activity = get_carbon_activity(xc)

    step_size /= 2
    xc = np.where(activity < 1, xc + step_size, xc - step_size)

activity = get_carbon_activity(xc)
xh, xo = get_h_o_fractions(xc)


#  A search has only found a boundary if its activity ended up close to 1. The
#  other points are discarded:

on_boundary = (activity > 0.9) & (activity < 1.1)


#  Plot the carbon deposition boundaries. The points of the ternary diagram are
#  tuples of the coordinates for the bottom, the right and the left axis in
#  percent:

scale = 100
fig, ax = plt.subplots(figsize=(10, 8))
tax = ternary.TernaryAxesSubplot(ax=ax, scale=scale)

# Draw boundary and gridlines
tax.boundary(linewidth=2.0)
tax.gridlines(color="black", multiple=10)

tax.left_axis_label("H", fontsize=20, offset=0.2)
tax.right_axis_label("C", fontsize=20, offset=0.2)
tax.bottom_axis_label("O", fontsize=20, offset=0.2)

# Set ticks
tax.ticks(ticks=list(range(0, 110, 20)), axis='lbr', linewidth=1, fontsize=15, offset=0.04)

# Remove default matplotlib ticks and box
tax.clear_matplotlib_ticks()
tax.get_axes().axis('off')

for i, tc in enumerate(t_range):
    mask = on_boundary[i]
    plot_data = list(zip(xo[i, mask] * scale, xc[i, mask] * scale, xh[i, mask] * scale))
    tax.plot(plot_data, label=f'{tc} °C')

tax.legend(fontsize=15)

# The axis labels are only drawn by tax.show() or by this call
tax._redraw_labels()