'''
IN PROGRESS: Calculator: Solid Phase Properties
===============================================

This notebook calculates the activity of endmembers given the composition of a solid phase.
It replaces the functionality previously implemented by the MELTS Supplemental Calculator (https://melts.ofm-research.org/CalcForms/)

This calculator is a work in progress -- please let us know if you find any errors!

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_4_unorganized%2Fplot_supplemental_calculator_solids.ipynb

'''

# %% [markdown]
# Information
# -----------
#
# The following solution phases are implemented in this calculator:
#
# - 'Ol'   (olivine)
# - 'Fsp'  (feldspar)
# - 'Opx'  (orthopyroxene)
# - 'Cpx'  (clinopyroxene)
# - 'SplS' (spinel solid solution)
# - 'Gt'   (garnet)
#
# **This notebook is experimental and has not yet been extensively tested! Proceed with caution!**

# %%
# User Input
# ----------
#
# 1. Define phase composition in wt% oxides, temperature, and pressure below. Components not included in the phase model will be ignored.
# 2. Once you have input your data, click "Run" --> "Run All Cells" at the top of the notebook. This will run the entire notebook and produce an output table at the bottom of the file.
#
phase = 'SplS'

oxide_wt_dict = { 
        'SiO2':     0, 
        'TiO2':     0.118399, 
        'Al2O3':    31.4741, 
        'Fe2O3':    3.623535369,
        'Cr2O3':    33.6659, 
        'FeO':      15.77978837, 
        'MnO':      0.248043,
        'MgO':      13.0798, 
        'NiO':      0.140058,
        'CoO':      0.000,
        'CaO':      0.019049, 
        'Na2O':     0.0039, 
        'K2O':      0.000, 
        'P2O5':     0.000, 
 }

TK = 794.12 + 273.15
Pbar = 6000

# Advanced options
database='MELTS_v1_2'


# %%
# Calculations
# ------------
#
# The user does not need to modify this code-- it will run automatically when "Run All Cells" is triggered.

###-----------------------------------------------------------------------------------------
### 0. Import modules
###-----------------------------------------------------------------------------------------

import numpy as np
import pandas as pd
from thermoengine import model, rockychem

###-----------------------------------------------------------------------------------------
### 1. Create phase instance
###-----------------------------------------------------------------------------------------

phs_obj = model.Database(database_name="MELTS_v1_2").get_phase(phase)
phs_endmems = rockychem.convert_phase_comp(phase, database=database, wt_oxides=oxide_wt_dict, to='mol_endmems')

if phase == 'SplS':
    # Custom method for calculating endmembers; more accurate for magnetite activity than default endmember calculation
    mol_oxides = rockychem.convert_bulk_comp(wt_oxides=oxide_wt_dict, to='mol_oxides')

    mol_cat = pd.Series([mol_oxides['TiO2'], 
                        2*mol_oxides['Al2O3'],
                        2*mol_oxides['Cr2O3'],
                        2*mol_oxides['Fe2O3'],
                        mol_oxides['FeO'],
                        mol_oxides['MgO']],
                        ['Ti', 'Al', 'Cr', 'Fe3', 'Fe2', 'Mg'])

    cat_O = pd.Series([2, 3/2, 3/2, 3/2, 1, 1],
        ['Ti', 'Al', 'Cr', 'Fe3', 'Fe2', 'Mg'])
    mol_O = sum(mol_cat*cat_O)
    norm_cat = 4*mol_cat/mol_O

    magnesioaluminate = norm_cat['Mg']
    chromite = norm_cat['Cr']/2
    ulvospinel = norm_cat['Ti']
    magnetite = norm_cat['Fe3']/2
    hercynite = norm_cat['Fe2'] - (chromite + 2*ulvospinel + magnetite)

    phs_endmems = np.array([hercynite, magnesioaluminate, chromite, ulvospinel, magnetite])

###-----------------------------------------------------------------------------------------
### 2. Calculate Gibbs Free Energies of Formation (Composition-Independent)
###-----------------------------------------------------------------------------------------

end_array = np.identity(len(phs_endmems))
G_array = np.zeros_like(phs_endmems)
for i in range(len(phs_endmems)):
    G = phs_obj.gibbs_energy(TK,Pbar,mol=end_array[i,:])
    G_array[i] = G

###-----------------------------------------------------------------------------------------
### 3. Calculate Activities (Composition-Dependent)
###-----------------------------------------------------------------------------------------

activities = phs_obj.activity(TK, Pbar, mol=phs_endmems)

###-----------------------------------------------------------------------------------------
### 4. Combine and Display Data
###-----------------------------------------------------------------------------------------

data = {'Endmember': phs_obj.endmember_names, 
        'Formula': phs_obj.endmember_formulas,
        'Gibbs Energy (kJ)': G_array/1000, 
        'Activity': activities}
endmember_df = pd.DataFrame(data)
pd.set_option('display.max_columns', None)
pd.set_option('display.max_rows', None)
endmember_df

# %%