'''
Calculator: Liquid Properties
===========================================

This notebook calculates the activity of endmembers within a liquid.
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%2Fsupplemental_calculator.ipynb

'''

# %%
# User Input
# ----------
#
# 1. Input liquid composition in wt% oxides, temperature, pressure, and log(fO2) below
# 2. Once you have input your data in the cell above, click "Run" --> "Run All Cells" at the top of the notebook. This will run the entire notebook and produce the output table at the bottom of the file.

oxide_wt_dict = {
    'SiO2':    48.4012,	
    'TiO2':     0.9208,	
    'Al2O3':   14.9349, 
    'FeO':      9.0825,
    'MnO':      0.1587, 
    'MgO':      8.1285, 
    'CaO':      14.1368, 
    'Na2O':      1.3124, 
    'K2O':       0.1462, 
    'P2O5':      0.001
    }

TK = 1500 + 273.15 # temperature value in K
Pbar = 15000 # pressure value in bar
logfO2 = None # can replace this with a value for absolute logfO2 (not relative to a buffer). Otherwise, calculates with given Fe2O3 value.

# %%
# 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, equilibrate, chemistry
from thermoengine.core import chem
from thermoengine import redox

###-----------------------------------------------------------------------------------------
### 1. Calculate oxide moles
###-----------------------------------------------------------------------------------------

oxide_moles = chem.format_mol_oxide_comp(oxide_wt_dict, convert_grams_to_moles=True)
oxide_moles_ser = pd.Series(oxide_moles, index=chemistry.OXIDES)
oxide_dict_for_fO2 = {'Liq': oxide_moles}

###-----------------------------------------------------------------------------------------
### 2. Adjust composition to be consistent with log fO2 input
###-----------------------------------------------------------------------------------------

if logfO2 is not None:
    ln_Fe_oxide_ratio = redox.redox_state(TK, Pbar, oxide_comp=oxide_dict_for_fO2, logfO2=logfO2,
                phase_of_interest='Liq', method='Kress91')
    Fe2O3_FeO = np.exp(ln_Fe_oxide_ratio)
    oxide_moles_ser = redox.impose_redox_ratio(Fe2O3_FeO, oxide_moles_ser)

###-----------------------------------------------------------------------------------------
### 3. Create liquid phase instance
###-----------------------------------------------------------------------------------------

liquid = model.Database(database_name="MELTS_v1_2").get_phase('Liq')
liq_end = liquid.calc_endmember_comp(oxide_moles_ser.values, method='least_squares',
                            output_residual=False, normalize=False,
                            decimals=10)

###-----------------------------------------------------------------------------------------
### 4. Calculate Gibbs Free Energies of Formation
###-----------------------------------------------------------------------------------------

end_array = np.identity(liquid.endmember_num)
G_array = np.zeros(liquid.endmember_num)
for i in range(liquid.endmember_num):
    G = liquid.gibbs_energy(TK,Pbar,mol=end_array[i,:])
    G_array[i] = G

###-----------------------------------------------------------------------------------------
### 5. Calculate Activities
###-----------------------------------------------------------------------------------------

activities = liquid.activity(TK, Pbar, mol=liq_end)

###-----------------------------------------------------------------------------------------
### 6. Combine and display data
###-----------------------------------------------------------------------------------------

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

# %%