'''
Calculator: Activity of Magnetite in Spinel
===========================================

Calculates a_Fe3O4 for spinel oxybarometry applications.

Written by Suzanne Birner (birners@berea.edu) in 2017. 
Last updated June 2025.

Open this code in an executable MyBinder instance:

.. image:: https://mybinder.org/badge_logo.svg
  :target: https://mybinder.org/v2/gl/swmatthews-research%2FThermoEngineLite/main?urlpath=%2Fdoc%2Ftree%2F.%2Fdoc%2Fsource%2Fmanual_examples%2Factivity_of_magnetite%2Factivity_of_magnetite_calculator.ipynb

Uses spinel models from Sack and Ghiorso (1991ab):
* Sack, R. O. & Ghiorso, M. S. (1991). Chromian spinels as petrogenetic indicators: Thermodynamics and petrological applications. American Mineralogist 76, 827–847.
* Sack, R. O. & Ghiorso, M. S. (1991). An internally consistent model for the thermodynamic properties of Fe-Mg-titanomagnetite-aluminate spinels. Contributions to Mineralogy and Petrology 106, 474–505.

User Instructions:

1. Either open this code as an executable MyBinder instance (link above) or download the source code. If you download the source code, you will need to have the thermoengine package installed on your machine.

2. Open the input_data folder to download the example input data OR download a blank input template at https://gitlab.com/skbirner/Spinel_Oxybarometry/blob/master/aFe3O4%20calculator%20input.xlsx. 
 
3. Fill in the input data template with spinel endmember data (columns AA-AG in the "Spinel endmembers" tab of the "fO2 calculation template.xlsx" workbook, located at https://gitlab.com/skbirner/Spinel_Oxybarometry/blob/master/fO2%20calculation%20template.xlsx.
 
4. Upload your data file to the input_data folder. Enter the name of the input data file below.
 
5. Select "Run All Cells" from the Jupyter menu.
 
6. Your data will be output in a pandas dataframe under the section "Output Activity Data". This script automatically calculates activity at T+25°C and T-25°C degrees as well, which is used in the error calculation for fO2.
 
7. This data can then be copy-pasted into columns I-K in the "fO2" tab within the "fO2 calculation template.xlsx" workbook. (Use "Paste Special --> Text" in Excel to get it to format mostly properly).
'''

# %%
# Input Data:
# -----------
#
# Enter the name of your input file below. 
# (By default, the file read is that of the test data. When replacing the filename, be sure to keep the "input_data/" prefix at the beginning to indicate that the data is stored in this subdirectory. Also be sure to include the file extension (e.g., '.xlsx')

Excel_File = 'input_files/aFe3O4_calculator_input_testdata.xlsx'

# %%
# Functions and Code:
# -------------------
#
# This section contains the code that calculates magnetite activity in spinel. You shouldn't need to change anything here.

### Import functions
import thermoengine
import numpy as np
import pandas as pd

### Helper functions

### Main function
def activity_magnetite_in_spinel(spl_endmember_mols, TK, Pbar, sample_ID, ox_or_end, warn=True):
    # Calculates the activity of magnetite in spinel, from endmember moles, P (in bars) and T (in K)
    # Input 'ox' for oxide input and 'end' for endmember input
    # Written by Suzanne Birner, May 2017, last updated June 2025

    # Input endmember order: Chromite, Hercynite, Magnetite, Spinel (MgAl), Ulvospinel

    R = 8.314 #J/(mol*K)

    # Endmember order: Hercynite, Magnesioaluminate, Chromite, Ulvospinel, Magnetite
    reordered_endmember_mols = np.array([spl_endmember_mols[1], spl_endmember_mols[3], spl_endmember_mols[0], spl_endmember_mols[4], spl_endmember_mols[2]])

    # Import spinel phase objects (sample and pure magnetite comparison)
    phs_calc = thermoengine.model.get_phase_calculator('MELTS_Spinel_SG91')
    phs = thermoengine.phases.SolutionPhase('Spl', phs_calc)
    phs_0 = thermoengine.phases.SolutionPhase('Spl', phs_calc)
    
    # Calculate chemical potentials (sample and pure magnetite comparison)
    u = phs.chem_potential(TK, Pbar, mol=reordered_endmember_mols, endmember=4)
    u_0 = phs_0.chem_potential(TK, Pbar, mol=np.array([0,0,0,0,1.0]), endmember=4)

    ln_a = (u - u_0)/(R*TK)
    a = np.exp(ln_a)
    return a[0]

### Read Data

# Load file 
xl = pd.ExcelFile(Excel_File)
input_data = xl.parse('Sheet1').values

# Define variables
sample_ID = input_data[:,0]
sp_endmembers = input_data[:,2:7]
TK = input_data[:,7] + 273.15
Pbar = input_data[:,8]

### Perform calculations

# Pre-allocate output vector
activity = np.zeros(len(TK))
activity_Tpos = np.zeros(len(TK))
activity_Tneg = np.zeros(len(TK))

# Iterate over rows 
for idx, val in enumerate(TK):
    activity[idx] = activity_magnetite_in_spinel(sp_endmembers[idx,:], TK[idx], Pbar[idx], 
                                                      sample_ID[idx], 'end')
    activity_Tpos[idx] = activity_magnetite_in_spinel(sp_endmembers[idx,:], TK[idx]+25, Pbar[idx], 
                                                      sample_ID[idx], 'end', warn=False)
    activity_Tneg[idx] = activity_magnetite_in_spinel(sp_endmembers[idx,:], TK[idx]-25, Pbar[idx], 
                                                      sample_ID[idx], 'end', warn=False)
    
# %%
# Results
# -------

df = pd.DataFrame(
    {'a_Fe3O4': activity,
     'a_Fe3O4, T+': activity_Tpos,
     'a_Fe3O4, T-': activity_Tneg
    })

pd.set_option('display.max_rows', None)
df