'''
MagmaForge: Fractional Crystallization
=======================================

Demonstrates how to run an equilibrium crystallization routine using MagmaForge.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load-- please be patient!):

.. image:: https://mybinder.org/badge_logo.svg
  :target: https://mybinder.org/v2/gl/swmatthews-research%2FThermoEngineLite/main?urlpath=%2Fdoc%2Ftree%2F.%2Fdoc%2Fsource%2Fauto_examples%2F_2_magmaforge%2Fplot_magmaforge_crystallization_frac.ipynb
'''

# %% 
# Initialization
# --------------
#
# Import necessary packages:

import thermoengine
from thermoengine import magmaforge

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# %%
# Define a bulk composition:

morb_oxides = pd.Series({ 
            'SiO2':  48.68, 
            'TiO2':   1.01, 
            'Al2O3': 17.64, 
            'Fe2O3':  0.89,
            'Cr2O3':  0.0425, 
            'FeO':    7.59, 
            'MgO':    9.10, 
            'CaO':   12.45, 
            'Na2O':   2.65, 
            'K2O':    0.03, 
            'P2O5':   0.08, 
            'H2O':    0.2},)

# %% 
# System and Calculations
# -----------------------
#
# Define system using magmaforge:

sys = magmaforge.System(comp=morb_oxides, 
                        T_liquidus=True,
                        P_bar=1000.0, 
                        logfO2 = ('QFM', -1),
                        database='MELTS_v1_0',
                        )

# %%
# Crystallize system:

sys.crystallize(method='frac', 
                T_step=4,
                fix_fO2=True, 
                T_final_C = 1000
                )

# %% 
# Output
# ------
#
# Magma Evolution plot

magmaforge.plot.magma_evolution(sys.history)

# %%
# If you want to specify the color for each phase, you can assign them via a dictionary:
phase_colors = {'Spinel': 'gray',
                'Feldspar': 'tan',
                'Olivine': 'springgreen',
                'Clinopyroxene' : 'green',
                'Orthopyroxene': 'darkolivegreen'}

magmaforge.plot.magma_evolution(sys.history, phase_colors=phase_colors)

# %%
# Return useful output variables
# ------------------------------
#
# The following commands demonstrate examples of how to return properties for all states in system's history.
# Other properties can be accessed as well-- see https://thermoenginelite.readthedocs.io/en/latest/api/magmaforge_system.html for more information.

T = sys.history.get_temperatures(unit='C')
P = sys.history.get_pressures(unit='bar')
logfO2 = sys.history.get_logfO2s(buffer='QFM')
liq_comp = sys.history.get_phase_comps(phase_name='Liquid', unit='wt_oxides')
phase_frac_table = sys.history.get_phase_frac_table(phases='present')

print(phase_frac_table) # replace with whichever variable you would like to access

# %%
# If you wish to save the full state history, you can do so with the .save_full_state_history_table() method.
# Soon, you will be able to reload this table into MagmaForge to access its properties.
sys.history.save_full_state_history_table(filename='magmaforge_output.csv', P_unit='bar', T_unit='K', fO2_buffer='QFM')


# %%