Worked Examples

This section provides the same set of worked examples that is contained in the ChemApp documentation provided by GTT-Technologies, but using the friendly interface instead of the basic one. When you compare the basic and friendly worked examples, the benefits of using the friendly interface should be clear.

Equilibrium Composition of Mixture Phases

"""
Equilibrium composition of mixture phases.
"""

# Import the ChemApp for Python modules that we will use in our calculations.
from chemapp.core import (TemperatureUnit,
                            PressureUnit, VolumeUnit, AmountUnit,
                            EnergyUnit)
from chemapp.friendly import Units as cau
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import StreamCalculation as casc

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# Initialise ChemApp for Python.
# Open, read and close the thermochemical data-file cosi.dat (system
# C-O-Si).
cats.load(str(data_dir / 'cosi.dat'))

# Change default units to desired units.
cau.set(T = TemperatureUnit.K,
        P = PressureUnit.bar,
        V = VolumeUnit.dm3,
        A = AmountUnit.mol,
        E = EnergyUnit.J)

# For the present example, the input composition consists of 2 mol of
# CO/GAS/ and 3 mol of quartz.
ico = cats.get_index_pc(0, 'CO')  # index of gas phase is always 0.
caec.set_IA_pc(0, ico, 2.0)

isio2 = cats.get_index_ph('SiO2(quartz)')
caec.set_IA_pc(isio2, 0, 3.0)

# The temperature at which the equilibrium is to be calculated is set to
# 1500 K.
caec.set_eq_T(1500.0)

# Calculate the equilibrium and print a ChemSage result table.
casc.calculate_eq(print_results=True)

# First, get the equilibrium amount of the gas phase, which is the only
# mixture phase in the system. Note that since we changed the units, the
# amount unit (mol) is used.
agas = caec.get_eq_A_ph(0)

print(f'\nEquilibrium amount of gas phase in mol: {agas:.4f}')

# To determine the mole fraction of the gas phase, relative to the whole
# system, we first need to get the amount contained in the system.
asys = caec.get_eq_A()

print(f'\nMole fraction of the gas phase: {agas / asys:.5f}')

# The equilibrium amounts of all phase constituents of the gas phase can
# be obtained by using the get_eq_A_pcs_in_ph function and the phase
# index.
values = caec.get_eq_A_pcs_in_ph(0)

# Print the equilibrium amount of each phase constituent.
npcon = cats.get_count_pcs_in_ph(0)
print('\nEquilibrium amounts of gas phase constituents:')

names = cats.get_name_pcs_in_ph(0)
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5E} mol')

# The mole fractions of the phase constituents can be calculated using
# the values just retrieved.

# Print the mole fractions of each phase constituent.
print(f'\nMole fractions of gas phase constituents:')
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20}  {values[i] / agas:.5E}')

# Information relating to the system components (in the present case
# C, O, and Si) to the phases can also be obtained.

# Get the mole fraction of the system components in the gas phase.
nscom = cats.get_count_scs()
print(f'\nMole fractions of system components in the gas phase:')
names = cats.get_name_scs()
values = caec.get_eq_X_scs_in_ph(0)
for i in range(nscom):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5}')

# To get the fraction of a system component that reports to the gas
# phase, we have to make use of the information contained in the
# stoichiometry matrix.
isi = cats.get_index_sc('Si')
fracsi = 0.0

# Get the equilibrium amount of each gas phase constituent, then use the
# stoichiometry matrix to determine how much Si is contained in the
# constituent, and sum it all up.
for i in range(npcon):
    wmass = cats.get_mm_pc(0, i)
    stoi = cats.get_stoich_pc(0, i)
    value = caec.get_eq_A_pc(0, i)
    fracsi += stoi[isi] * value

# Get the total amount of Si in the system.
sitot = caec.get_eq_A_sc(isi)

print(f'\nMole fraction of system component Si that reports to the gas '
      f'phase: {fracsi / sitot:.5E}')

# To demonstrate how ChemApp for Python can be used to check the relationship
# between partial and integral values of extensive properties, a
# different data-file for a system that contains non-ideal mixture
# phases is used.

# Open a new thermochemical data-file (pbsn.dat, system Pb-Sn) for
# reading.
cats.load(str(data_dir / 'pbsn.dat'))

# Find out the system component index numbers for Pb and Sn, as well as
# the phase index number for the liquid phase:
ipb = cats.get_index_sc('Pb')
isn = cats.get_index_sc('Sn')
iliq = cats.get_index_ph('LIQUID')

# Using the index numbers just determined, we can enter our incoming
# amounts, in this case an alloy consisting of 0.7 mol Sn and 0.3 mol
# Pb:
xsn = 0.7
caec.set_IA_sc(isn, xsn)
caec.set_IA_sc(ipb, 1.0 - xsn)

# Set the temperature to 600 K. Under these conditions, the liquid phase
# is the only stable phase in this system (see the Pb-Sn phase diagram
# provided in worked example 3).
caec.set_eq_T(600)

# Calculate the equilibrium.
caec.calculate_eq()

# Get the integral enthalpy of the liquid phase.
hint1 = caec.get_eq_H_ph(iliq)

# Get the partial molar enthalpies of Sn/LIQUID/ and Pb/LIQUID/.
hsn = caec.get_eq_HM_pc(iliq, isn)
hpb = caec.get_eq_HM_pc(iliq, ipb)

# Calculate the integral enthalpy of the liquid phase from the two
# partial values.
hint2 = xsn * hsn + (1.0 - xsn) * hpb

# Print both values.
print(f'\nhint1 = {hint1:.3f} J')
print(f'hint2 = {hint2:.3f} J')

# As expected, the two values are the same.

Equilibrium Composition of Mixture Phases

"""
Equilibrium composition of mixture phases.
"""

# Import the ChemApp for Python modules that we will use in our calculations.
from chemapp.core import (TemperatureUnit,
                            PressureUnit, VolumeUnit, AmountUnit,
                            EnergyUnit)
from chemapp.friendly import Units as cau
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import StreamCalculation as casc

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# Initialise ChemApp for Python.
# Open, read and close the thermochemical data-file cosi.dat (system
# C-O-Si).
cats.load(str(data_dir / 'cosi.dat'))

# Change default units to desired units.
cau.set(T = TemperatureUnit.K,
        P = PressureUnit.bar,
        V = VolumeUnit.dm3,
        A = AmountUnit.mol,
        E = EnergyUnit.J)

# For the present example, the input composition consists of 2 mol of
# CO/GAS/ and 3 mol of quartz.
ico = cats.get_index_pc(0, 'CO')  # index of gas phase is always 0.
caec.set_IA_pc(0, ico, 2.0)

isio2 = cats.get_index_ph('SiO2(quartz)')
caec.set_IA_pc(isio2, 0, 3.0)

# The temperature at which the equilibrium is to be calculated is set to
# 1500 K.
caec.set_eq_T(1500.0)

# Calculate the equilibrium and print a ChemSage result table.
casc.calculate_eq(print_results=True)

# First, get the equilibrium amount of the gas phase, which is the only
# mixture phase in the system. Note that since we changed the units, the
# amount unit (mol) is used.
agas = caec.get_eq_A_ph(0)

print(f'\nEquilibrium amount of gas phase in mol: {agas:.4f}')

# To determine the mole fraction of the gas phase, relative to the whole
# system, we first need to get the amount contained in the system.
asys = caec.get_eq_A()

print(f'\nMole fraction of the gas phase: {agas / asys:.5f}')

# The equilibrium amounts of all phase constituents of the gas phase can
# be obtained by using the get_eq_A_pcs_in_ph function and the phase
# index.
values = caec.get_eq_A_pcs_in_ph(0)

# Print the equilibrium amount of each phase constituent.
npcon = cats.get_count_pcs_in_ph(0)
print('\nEquilibrium amounts of gas phase constituents:')

names = cats.get_name_pcs_in_ph(0)
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5E} mol')

# The mole fractions of the phase constituents can be calculated using
# the values just retrieved.

# Print the mole fractions of each phase constituent.
print(f'\nMole fractions of gas phase constituents:')
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20}  {values[i] / agas:.5E}')

# Information relating to the system components (in the present case
# C, O, and Si) to the phases can also be obtained.

# Get the mole fraction of the system components in the gas phase.
nscom = cats.get_count_scs()
print(f'\nMole fractions of system components in the gas phase:')
names = cats.get_name_scs()
values = caec.get_eq_X_scs_in_ph(0)
for i in range(nscom):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5}')

# To get the fraction of a system component that reports to the gas
# phase, we have to make use of the information contained in the
# stoichiometry matrix.
isi = cats.get_index_sc('Si')
fracsi = 0.0

# Get the equilibrium amount of each gas phase constituent, then use the
# stoichiometry matrix to determine how much Si is contained in the
# constituent, and sum it all up.
for i in range(npcon):
    wmass = cats.get_mm_pc(0, i)
    stoi = cats.get_stoich_pc(0, i)
    value = caec.get_eq_A_pc(0, i)
    fracsi += stoi[isi] * value

# Get the total amount of Si in the system.
sitot = caec.get_eq_A_sc(isi)

print(f'\nMole fraction of system component Si that reports to the gas '
      f'phase: {fracsi / sitot:.5E}')

# To demonstrate how ChemApp for Python can be used to check the relationship
# between partial and integral values of extensive properties, a
# different data-file for a system that contains non-ideal mixture
# phases is used.

# Open a new thermochemical data-file (pbsn.dat, system Pb-Sn) for
# reading.
cats.load(str(data_dir / 'pbsn.dat'))

# Find out the system component index numbers for Pb and Sn, as well as
# the phase index number for the liquid phase:
ipb = cats.get_index_sc('Pb')
isn = cats.get_index_sc('Sn')
iliq = cats.get_index_ph('LIQUID')

# Using the index numbers just determined, we can enter our incoming
# amounts, in this case an alloy consisting of 0.7 mol Sn and 0.3 mol
# Pb:
xsn = 0.7
caec.set_IA_sc(isn, xsn)
caec.set_IA_sc(ipb, 1.0 - xsn)

# Set the temperature to 600 K. Under these conditions, the liquid phase
# is the only stable phase in this system (see the Pb-Sn phase diagram
# provided in worked example 3).
caec.set_eq_T(600)

# Calculate the equilibrium.
caec.calculate_eq()

# Get the integral enthalpy of the liquid phase.
hint1 = caec.get_eq_H_ph(iliq)

# Get the partial molar enthalpies of Sn/LIQUID/ and Pb/LIQUID/.
hsn = caec.get_eq_HM_pc(iliq, isn)
hpb = caec.get_eq_HM_pc(iliq, ipb)

# Calculate the integral enthalpy of the liquid phase from the two
# partial values.
hint2 = xsn * hsn + (1.0 - xsn) * hpb

# Print both values.
print(f'\nhint1 = {hint1:.3f} J')
print(f'hint2 = {hint2:.3f} J')

# As expected, the two values are the same.

Simple Process Using Streams

"""
A simple process using streams.
"""
# Note that this program contains phase target calculations, which
# cannot be performed with the 'light' version of ChemApp for Python.

from chemapp.core import(ResultVariable, ConditionVariable,
                           StreamVariable, Quantity,
                           AmountUnit, TemperatureUnit, PressureUnit,
                           VolumeUnit, EnergyUnit, TargetVariable)
from chemapp.friendly import Info as cai
from chemapp.friendly import Units as cau
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import StreamCalculation as casc

# import packages for plotting graphs.
import matplotlib.pyplot as plt
import pandas as pd

# import package for 3D-plot.
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.mlab import griddata
import scipy.interpolate

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# The following function calculates the sum of all equilibrium amounts
# of the condensed phases. It starts with phase #1, since the gas phase
# is #0. The index number of the last phase to be considered is passed
# as an argument to the function and is retrieved in the main program.

def getacp(lastph):
    getacp = 0.0
    # Add up the equilibrium amounts of all condensed phases.
    for i in range(lastph):
        amount = caec.get_eq_A_ph(i)
        getacp += amount
    return amount

# Main program.

# Initialise ChemApp for Python.
# Open data-file for reading. In this case, a data-file containing the
# thermochemical data for the system carbon-nitrogen-oxygen is selected.
# Read data-file.
# Close data-file.

cats.load(str(Path(data_dir) / 'cno.dat'))

# Since we will perform an extensive property target calculation later,
# we will use streams. This way we can assign an initial state regarding
# temperature and pressure to the incoming substances, which will be
# taken into account when ChemApp for Python calculates the extensive property
# (e.g. heat) balance.

# Since we are going to use degrees Celsius as temperature unit
# throughout the calculations, we will change the unit now.
cau.set(T = TemperatureUnit.C,
        V = VolumeUnit.dm3,
        P = PressureUnit.bar,
        E = EnergyUnit.J,
        A = AmountUnit.mol)

# According to the process scheme shown above, define the first stream.
# It is the incoming stream of mixed CO2 and CO at 200 C and 1 bar,
# which we will call 'CO2CO'. The amount of CO2 and CO in this stream
# is set to 0.5 and 2.5 mol, respectively.
casc.create_st('CO2CO', 200.0, 1.0)

# Set the constituent amounts of CO2 and CO. To do this, get the index
# number of these constituents of the gas phase, which is phase number
# 0.
ico2 = cats.get_index_pc(0, 'CO2')
ico = cats.get_index_pc(0, 'CO')

# Use the index numbers just obtained to set the constituent amount of
# CO2 in the stream 'CO2CO' to 0.5 mol.
casc.set_IA_pc('CO2CO', 0, ico2, 0.5)

# Similarly, set the amount of CO in the same stream to 2.5 mol.
casc.set_IA_pc('CO2CO', 0, ico, 2.5)

# Define the second stream ('O2N2'), which is air supplied at room
# temperature and 1 bar for the combustion. The air is assumed to
# consist of O2 and N2 at a ratio of 1:4. The amount of O2 is
# initially set to 0.5 mol (thus N2 = 2.0 mol), and will later
# be varied.
casc.create_st('O2N2', 25.0, 1.0)

# Get the index numbers of O2 and N2.
in2 = cats.get_index_pc(0, 'N2')
io2 = cats.get_index_pc(0, 'O2')

# Set the amount of O2.
casc.set_IA_pc('O2N2', 0, io2, 0.5)

# Set the amount of N2.
casc.set_IA_pc('O2N2', 0, in2, 2.0)

# Set the temperature of the reactor at which the equilibrium reaction
# (i.e. the combustion) will take place. It is initially set to 600 C
# at first, but will be varied later.
casc.set_eq_T(600.0)

# In the result table to be printed out subsequently, the amount unit
# for the condensed phases formed should be 'gram', while the input
# was given in 'mol'.
cau.set_A_unit(AmountUnit.gram)

# Calculate the equilibrium and print out the result table for a
# quick overview on the equilibrium state.
casc.calculate_eq(print_results = True)

# As can be seen from the above table, both the gas phase and carbon
# are stable under the conditions chosen. It can also be seen that
# N2, CO2, and CO are the only majority gas phase species. All
# other gas phase species are produced in negligible amounts only.

# The values for the amounts of condensed phases and the gas phase
# can be retrieved using get_count_phs(). Since the gas phase, if
# present in the data-file, has always the phase index number 0,
# this means that we have to add up all the amounts of the phases
# from phase #2 to the last one in order to retrieve the total
# amount of condensed phases.

# Find out how many phases we have in the data-file, and use this
# number in all subsequent calls to getacp.
nphase = cats.get_count_phs()

# Get the total equilibrium amount of the condensed phases, using
# the function getacp declared above.
acond = getacp(nphase)

print(f'\nAmount of condensed material: {acond:.4f} gram.')


# Retrieve the amount of gas phase. Note: The current amount unit
# is still 'gram'. If the amount of the gas phase is desired in
# 'mol' instead, tqcsu would need to be called again.
agas = casc.get_eq_A()

print(f'\nAmount of gas phase: {agas:.2f} gram.')

# Retrieve values for some extensive properties.  Since we used
# streams for this example, these extensive property values are
# calculated relative to the incoming streams, when the whole
# system is selected.

# Get the enthalpy change associated with this reaction.
result = casc.get_eq_H()

print(f'\nDelta H: {result:.4E} J')

# The negative sign of this value indicates that energy is released
# by the system.

# Get the change in volume of the process, which is the volume of
# the equilibrium gas phase minus the sum of the volumes of the
# incoming streams.
result =  casc.get_eq_V()
print(f'Delta V: {result:.2f} dm^3')

# Note: using the option 'V' with tqgetr normally would retrieve
# the _total_ volume of the system. There is one exception, in
# which case it retrieves the _change_ in volume, and that is the
# case when streams are used _and_ the whole system is selected.
# To retrieve the _total_ volume instead in the present case, use
# the option 'VT':
result =  casc.get_eq_VT()
print(f'Total V: {result:.2f} dm^3')

# Compare the two values for the volume with the corresponding
# ones in the output table produced earlier. The volume printed
# at the top of such an output table is always the _total_ volume,
# whereas the line specifying the extensive property balance toward
# the end of the table shows the _change_ in volume.

# tqstxp can be used to retrieve extensive property values for the
# incoming streams. Get the volume of the incoming CO2/CO stream:
result =  casc.get_st_V('CO2CO')
print(f'\nV(CO2CO): {result:.2f} dm^3')

# Get the value of the incoming air stream:
result =  casc.get_st_V('O2N2')
print(f'V(O2N2): {result:.3f} dm^3')

# For users of ChemApp versions prior to 3.2.0: Please note that
# the values returned by tqstxp always er to the _last_ calculated
# equilibrium state. Even if one uses tqstxp to retrieve properties
# of _incoming_ streams, one thus needs to perform an equilibrium
# calculation before the correct values are obtained.

# As can be seen from the above results, this process, under the
# conditions we have chosen, results in solids to be formed in the
# reactor. We will now further investigate this systems by looking
# at how the amount of condensed material changes when
# a) the reaction temperature is varied, and
# b) the amount of air supplied is varied.

# *********************************************************************
# Varying the reaction temperature.
# *********************************************************************

# Since several individual equilibria will be calculated, the results
# (temperature and amount of condensed phases) will be stored to a
# file for further processing.
file_11 = open('caf103_1.rst', 'w')

# Vary the temperature between 200 C and 1000 C in steps of 5 C.
for j in range(200, 1001, 5):

    # Use the variable temp to hold the current temperature.
    temp = j
    casc.set_eq_T(temp)

    # Calculate the equilibrium.
    casc.calculate_eq()


    # Get the amount of condensed phases.
    acond = getacp(nphase)

    # Write the temperature and the amount of condensed phases to the
    # result file.
    file_11.write(f'{temp},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure.

# Read the created caf103_1.rst file, and name the first column
# 'temp' and name the second 'acond'.
data_1 = pd.read_csv('caf103_1.rst', names = ['temp', 'acond'])
plt.axis([200, 1000, 0, 10])
plt.ylabel('Amount/gram')
plt.xlabel('Temperature/C')
plt.plot(data_1['temp'], data_1['acond'], 'r-')
plt.show()

# *********************************************************************
# Varying the amount of air supplied.
# *********************************************************************

# Instead of the temperature, the amount of incoming air will be varied
# now. The temperature will be kept constant at the initial value of
# 600 C. To plot the results later on, we will again store them in a
# separate file.
file_11 = open('caf103_2.rst', 'w')

# Set the temperature for the equilibrium calculation (600 C).
casc.set_eq_T(600.0)

# Calculate 201 equilibria...
for j in range(201):

    # ... with the O2 amount varying between 0 and 1 mol...
    iao2 = j * 0.005
    # ... and the N2 amount between 0 and 4 mol.
    ian2 = iao2 * 4.0

    # To make sure ChemApp for Python interprets these values as 'mol', set the
    # amount unit accordingly.
    cau.set_A_unit(AmountUnit.mol)

    # Set the amount of O2 in the air stream.
    casc.set_IA_pc('O2N2', 0, io2, iao2)

    # Set the amount of N2 in the air stream.
    casc.set_IA_pc('O2N2', 0, in2, ian2)

    # Calculate the equilibrium.
    casc.calculate_eq()

    # Set the amount unit to 'gram'.
    cau.set_A_unit(AmountUnit.gram)

    # Get the amount of condensed phases.
    acond = getacp(nphase)

    # Write the amounts of O2 and condensed phases to the result file.
    file_11.write(f'{iao2},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure:
data_2 = pd.read_csv('caf103_2.rst', names = ['iao2', 'acond'])
plt.axis([0, 1, 0, 12])
plt.ylabel('Amount/gram')
plt.xlabel('n(O2)/mol')
plt.plot(data_2['iao2'], data_2['acond'], 'r-')
plt.show()

# *********************************************************************
# Varying the amount of air supplied _and_ the temperature.
# *********************************************************************

# To investigate the combined influence of temperature and amount of air
# supplied on the resulting amount of condensed phases, both parameters
# are varied in nested loops. The results are again stored in a file for
# further processing.
file_11 = open('caf103_3.rst', 'w')

# The outer loop will vary the temperature between 200 C and 1000 C.
for j in range(200, 1001, 20):

    temp = j
    casc.set_eq_T(temp)

    # The inner loop varies the amount of incoming air as above.
    for k in range(0, 201, 5):
        iao2 = k * 0.005
        ian2 = iao2 * 4.0

        cau.set_A_unit(AmountUnit.mol)
        # Set the incoming amount of O2.
        casc.set_IA_pc('O2N2', 0, io2, iao2)

        # Set the incoming amount of N2.
        casc.set_IA_pc('O2N2', 0, in2, ian2)

        casc.calculate_eq()
        cau.set_A_unit(AmountUnit.gram)

        # Get the amount of condensed phases
        acond = getacp(nphase)

        # Write the temperature and the amounts of O2 and condensed
        # phases to the result file.
        file_11.write(f'{temp},{iao2},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure:

data_3 = pd.read_csv('caf103_3.rst', names = ['temp', 'iao2', 'acond'])
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')
ax.set_xlim3d(200, 1000)
ax.set_ylim3d(0, 1)
ax.set_zlim3d(-10, 20)
ax.set_xlabel('Temperature/C')
ax.set_ylabel('n(O2)/mol')
ax.set_zlabel('Amount/gram')

x = data_3['temp']
y = data_3['iao2']
z = data_3['acond']
ax.scatter(x, y, z, color = 'red')
plt.show()

# Check if we are working with the 'light' version. If we do, omit
# the following phase target calculation(s).
islite = cai.is_light()
if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n '
          f'*** since they are not possible with the\n '
          f'*** ''light''version of ChemApp for Python.')

else:

    # *****************************************************************
    # Varying the amount of air supplied _and_ the temperature,
    # searching for the explicit condition where the amount of condensed
    # phase is zero.
    # *****************************************************************

    # Assuming that one is interested in the above process _not_
    # producing any condensed phases (i.e. complete combustion), it would
    # be helpful to have a diagram which corresponds to the bottom plane
    # of the above figure: A graph which, as a function of the incoming
    # amount of air, shows the temperature at which the amount of
    # condensed phases formed is zero.

    # For the calculation procedure with ChemApp for Python, this results in the
    # following steps:
    # 1) We are looking for the temperature at which condensed phases
    # _just_ start to form. In terms of ChemApp for Python, we thus have to define
    # a precipitation target with respect to the gas phase.
    # 2) For each value of the amount of air supplied, define the 'O2N2'
    # stream with the proper amounts of O2 and N2.
    # 3) ChemApp for Python is then supposed to vary the temperature until the
    # target (a phase precipitates from the gas phase) is met.
    # Temperature is thus the target variable.
    # 4) The equilibrium is calculated for each of the amounts of air
    # supplied, then the equilibrium temperature is retrieved and stored
    # in the result file.
    file_11 = open('caf103_4.rst', 'w')

    # Reset the amount unit to 'mol'.
    cau.set_A_unit(AmountUnit.mol)

    # Define the target as a precipitation target for the gas
    # phase (phase index 0).
    casc.set_target_precipitation_from_ph(0)

    # The initial estimate for the temperature is passed via the first
    # element of vals.
    vals = 600.0

    # Vary the amount of supplied air.
    for j in range(201):
        iao2 = j * 0.005
        ian2 = iao2 * 4.0

        # Set the incoming amount of O2.
        casc.set_IA_pc('O2N2', 0, io2, iao2)

        # Set the incoming amount of N2.
        casc.set_IA_pc('O2N2', 0, in2, ian2)

        # Define the target variable (temperature), with the initial
        # estimate in the first element of vals, and calculate the
        # equilibrium.
        casc.calculate_eq_T(vals)

        # Retrieve the calculated equilibrium temperature.
        temp = casc.get_eq_T()

        # Store the amount of O2 supplied and the calculated
        # equilibrium temperature to the result file.
        file_11.write(f'{iao2},{temp}\n')
    # Close data-file.
    file_11.close()

    # When plotted, the resulting graph looks like the following Figure.
    data_4 = pd.read_csv('caf103_4.rst', names = ['iao2', 'temp'])
    plt.axis([0, 1, 500, 800])
    plt.ylabel('Temperature/C')
    plt.xlabel('n(O2)/mol')
    plt.plot(data_4['iao2'], data_4['temp'], 'r-')
    plt.show()

One-dimensional Phase Mapping

"""
One-dimensional phase mapping calculations.

Note that this program contains phase target calculations, which cannot be
performed with the 'light' version of ChemApp for Python.
"""

from chemapp.core import (Status,
    ConditionVariable, TargetVariable, ResultVariable, PhaseMapVariable,
    AmountUnit)
from chemapp.friendly import Units as cau
from chemapp.friendly import Info as cai
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import PhaseMapCalculation as capmc

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# import packages for plotting graphs.
import matplotlib.pyplot as plt

# The function psp prints the value of the 'search variable' (temperature,
# pressure, or composition), the names of the stable phases, and their amounts.
# It is used after every call to calculate_ph_txs.
# The variable_str is set to 'Temperature/K' for temperature calls,
# 'Pressure/bar' for pressure etc. This function also prints a header for the result
# table that is populated by psp.
def psp(variable_str, variable, results):
    print(f'\n{variable_str:<26}{"Stable phases":<18}{"Amount/mol"}')

    eps = 1E-8

    for r in results:
        print(f'{getattr(r, variable):.2f}')

        for ph in r.phs.values():
            if ph.A > eps:
                print(f'{"":<26}{ph.name:<24}{ph.A:.2f}')


# Since the calculate_ph_txs_IA is different, a psp_ia function was written and
# it works on the same principles as psp.
def psp_ia(results):
    print(f'\n{"X(Fe2SiO4)":<26}{"X(Mg2SiO4)":<26}{"Stable phases":<18}'
          f'{"Amount/mol"}')

    eps = 1E-8

    for r in results:
        print(f'{r.IAs[0].A:<26.3f}{r.IAs[1].A:<26.3f}')

        for ph in r.phs.values():
            if ph.A > eps:
                print(f'{"":<52}{ph.name:<24}{ph.A:.2f}')


asize = 100
eps = 1E-3
nused = 0

# Initialise ChemApp for Python, open, read and close the thermochemical data-file
# pbsn.dat (system Pb-Sn).
cats.load(str(data_dir / 'pbsn.dat'))


# Check if we are working with the 'light' version.
# If we are, omit the following phase target calculation(s).
islite =  cai.is_light()

if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n'
          f'*** since they are not possible with the\n'
          f'*** ''light'' version of ChemApp for Python.')

    exit(0)


# The data-file for system Pb-Sn contains a second copy of the BCT_A5 phase,
# because of the inherent metastable miscibility gap. Since this miscibility
# gap will stay metastable under the conditions used in this example, we
# will eliminate the second copy of that phase, which means it will be
# disregarded completely in the following calculations.
bct2 = cats.get_index_ph('BCT_A5#2')
cats.set_status_ph(bct2, Status.ELIMINATED)

# Determine the index numbers for the system components Pb and Sn.
ipb = cats.get_index_sc('Pb')
isn = cats.get_index_sc('Sn')

# As indicated in the Pb-Sn phase diagram, label "1", the
# first phase map should be done with the temperature as the search
# variable, for X(Sn)=0.2.

# Set the incoming amounts for the Sn and Pb.
capmc.set_IA_sc(isn, 0.2)
capmc.set_IA_sc(ipb, 0.8)

tarray = list(range(410, 600, 10))

# The purpose of the array tarray is to hold the temperature values at which
# the equilibria will be calculated later.
# This set of temperature values consists of,
# 1) a number of evenly spaced temperature intervals,
# 2) and the temperatures at which calculate_ph_txs_T finds phase boundaries.

# One call is made to ChemApp for Python friendly, and all of the transitions are
# calculated and sorted.
results = capmc.calculate_ph_txs_T(400.0, 600.0)
psp('Temperature/K', 'T', results)

# We are now ready to calculate the equilibria for all the temperatures in
# tarray. These equilibrium calculations now being performed are _not_
# target calculations of any kind any more, that's why we reset _all_
# previous conditions including all targets and target variables.
caec.remove_eq_conditions_all()

# The following equilibria are to be calculated for the same alloy
#  concentration of course, so we use the same incoming amounts again.
caec.set_IA_sc(isn, 0.2)
caec.set_IA_sc(ipb, 0.8)

# The graphical representation of the one-dimensional phase map will be
# done using two x/y graphs: one displaying the phase amounts vs. the
# temperature, the other one the phase activities vs. the
# temperature.
# For every temperature in tarray, calculate the equilibrium:
for t in tarray:
    caec.set_eq_T(t)

    # Note again that these are regular equilibrium calculations, no
    # targets are involved:
    caec.calculate_eq()
    results.append(caec.get_result_object())

# Sort the result list by temperature. This is required because results are appended
# to the list.
results.sort(key=lambda r: r.T)

# When plotted, the two graphs, displaying the equilibrium phase
# amounts and phase activities vs. temperature look like Figures
# 3 and 4.

T_array = [r.T for r in results]
A_liq_array = [r.phs['LIQUID'].A for r in results]
A_bct_array = [r.phs['BCT_A5#1'].A for r in results]
A_fcc_array = [r.phs['FCC_A1'].A for r in results]

# Print the sorted values in T_array in ascending order.
for i in range(0, len(T_array)):
    print(f'TARRAY({i:>2}) = {T_array[i]:>8.5f}')

# Phase amounts vs. temperature.
plt.axis([400, 600, 0, 1])
plt.ylabel('Amount/mol')
plt.xlabel('Temperature/K')
plt.plot(T_array, A_liq_array, 'r-+',
         linewidth=0.5, markersize=5, label='LIQUID')
plt.plot(T_array, A_bct_array, 'g-x',
         linewidth=0.5, markersize=5, label='BCT_A5#1')
plt.plot(T_array, A_fcc_array, 'm-s',
         linewidth=0.5, markersize=5, markerfacecolor='none',
         label='FCC_A1')
plt.legend(loc='best')
# plt.show()

AC_liq_array = [r.phs['LIQUID'].AC for r in results]
AC_bct_array = [r.phs['BCT_A5#1'].AC for r in results]
AC_fcc_array = [r.phs['FCC_A1'].AC for r in results]

# Activities vs. temperature.
plt.axis([400, 600, 0.55, 1])
plt.ylabel('Activity')
plt.xlabel('Temperature/K')
plt.plot(T_array, AC_liq_array, 'r-+',
         linewidth=0.5, markersize=5, label='LIQUID')
plt.plot(T_array, AC_bct_array, 'g-x',
         linewidth=0.5, markersize=5, label='BCT_A5#1')
plt.plot(T_array, AC_fcc_array, 'm-s',
         linewidth=0.5, markersize=5, markerfacecolor='none',
         label='FCC_A1')
plt.legend(loc='best')
# plt.show()

# For the next one-dimensional phase mapping calculations we will only
# use tabular output. First, remove all previously set conditions.
capmc.remove_eq_conditions_all()

# We will now calculate another one-dimensional phase map in the same
# system, again with temperature as the search variable. The new
# alloy composition will cause calculate_ph_txs_T to find the eutectic.
capmc.set_IA_sc(isn, 0.5)
capmc.set_IA_sc(ipb, 0.5)

# For this, we call calculate_ph_txs_T to calculate the phase transitions
# according to temperature.
results = capmc.calculate_ph_txs_T(400.0, 600.0)
psp('Temperature/K', 'T', results)

# Note that the above table lists a 3-phase equilibria now, which for the
# current system means we found the eutectic temperature (compare this
# with the results for the same equilibrium determined in worked example
# 3: Phase equilibria and phase target calculations).

# So far we have calculated one-dimensional phase maps for a
# temperature/composition phase diagram, in our case keeping the
# composition (and implicitly the pressure) constant and using the
# temperature as the search variable. The next phase map we will
# calculate at a constant temperature, this time varying the
# pressure. For this purpose, we will use a different thermochemical
# data-file for the system Fe2SiO4-Mg2SiO4 which includes
# pressure-dependent Gibbs-energy data.

# Open, Read and close the thermochemical data-file femgsio4.dat
# (system Fe2SiO4-Mg2SiO4).
cats.load(str(data_dir / 'femgsio4.dat'))

# Determine the index numbers for the stoichiometric compounds Fe2SiO4
# and Mg2SiO4. These two phases will be used to define the incoming
# amounts for the next calculations. Note that they are included in the
# thermochemical data-file only for this purpose (to allow for a
# somewhat easier definition of incoming amounts), since the Gibbs
# energy functions for these two stoichiometric compounds are only
# dummies.

# Select a composition of 2 mol-% Fe2SiO4, the rest is Mg2SiO4:
capmc.set_IA_pc('Fe2SiO4', 0, 0.02)
capmc.set_IA_pc('Mg2SiO4', 0, 0.98)

# Select the range for the search variable, which is the pressure in the
# current case. We will search for any phase boundaries between 120 kbar
# and 150 kbar (see also the line marked '1' in Figure 2.)

# The temperature is set to 1000 K.
capmc.set_eq_T(1000.0)

# Call calculate_ph_txs_P to calculate the transitions pressure.
results = capmc.calculate_ph_txs_P(120000.0, 150000.0)

# Call psp (print stable phases), passing the value for the pressure
# we just retrieved. psp will print a formatted entry for the result
# table.
psp('Pressure/bar', 'P', results)

# The result table is complete now. Compare the phase boundaries found
# with the phase diagram in Figure 2. Note that the total
# amount in the system is always calculated to be 2 mol. This is due to
# the fact that the phase constituents of all mixture phases in this
# data-file (Spinell, Olivine, and \DF-phase) have been scaled with a
# factor of 0.5 in order to obtain the correct ideal entropy between Fe and
# Mg.

# The last one-dimensional phase map will use composition as the search
# variable, keeping both temperature and pressure constant. As with
# other phase target calculations which use composition as the target
# variable, the input in this case is a bit more complicated (see also
# worked example 3: Phase equilibria and phase target calculations).
# This is due to the fact that when one is interested in comparing the
# result for this type of calculation with graphical phase diagram
# information, it is necessary to keep the total substance amount in the
# system constant (i.e. use concentrations rather than absolute
# amounts).

# First, reset all conditions and targets again:
capmc.remove_eq_condition(-1)

# Set the pressure to 135 kbar (see Figure 2):
capmc.set_eq_P(135000.0)

# The composition range for Fe2SiO4 has a minimum of 0 mol and a maximum of 0.2
# mol, compare again to Figure 2. The inputs required by calculate_ph_txs_IA
# are defined in the list "ranges". Here the ph, pc (None for this example),
# range_start and range_end are defined.

# The results are printed by psp_ia.
ranges = [capmc.IA_range('Fe2SiO4', '', 0.0, 0.2),
          capmc.IA_range('Mg2SiO4', '', 1.0, 0.8)]
results = capmc.calculate_ph_txs_IA(ranges)
psp_ia(results)