#!/usr/bin/env python
# coding: utf-8

# # Amphibole Thermobarometry and Chemometry
# - This workbook demonstrates how to use the various amphibole-only, and amphibole-liquid functions
# - You can download the excel spreadsheet with data here:
# https://github.com/PennyWieser/Thermobar/blob/main/docs/Examples/Amphibole/Amphibole_Liquids.xlsx

# ### You need to install Thermobar once on your machine, if you haven't done this yet, uncomment the line below (remove the #)

# In[1]:


#!pip install Thermobar


# ## Load python things

# In[2]:


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import Thermobar as pt
pd.options.display.max_columns = None


# ## Load data

# In[3]:


out=pt.import_excel('Amphibole_Liquids.xlsx', sheet_name="Amp-Liq")
my_input=out['my_input']
Amps=out['Amps']
Liqs=out['Liqs']


# ## Checking data read in correctly (check no columns are full of zeros that you think you inputted)

# In[4]:


display(Amps.head())
display(Liqs.head())


# ## Example 1 - Amphibole-only pressure
# - can use help functions to see equation options

# In[5]:


help(pt.calculate_amp_only_press)


# ### Ridolfi 2021 - Minerals algorithm
# - here, we specify equationP="P_Ridolfi2021", which combines equations 1a-1e from Ridolfi 2012 and chooses the "correct" equation using a number of filters (see Ridolfi 2021, Minerals)
# - The output dataframe shows the calulated pressure, as well as the equation selected
# - The code also prints a warning that some of the rows have oxide sums <90, so a Nan is returned for these following Ridolfi 2021 supplement

# In[6]:


P_Ridolfi2021=pt.calculate_amp_only_press(amp_comps=Amps, equationP="P_Ridolfi2021")
P_Ridolfi2021.head()


# #### Without quality filters
# - If you want to use the pressures which yield "input errors", you can specify "Ridolfi_Filter=False"
# - This is more for looking at what ones it has discarded, we don't necessarily recomend this, we just have it as a feature you can use at your peril. 

# In[7]:


P_Ridolfi2021=pt.calculate_amp_only_press(amp_comps=Amps, equationP="P_Ridolfi2021", 
                                          Ridolfi_Filter=False)
P_Ridolfi2021.head()


# ###  Mutch et al. (2016)

# In[8]:


P_Mutch2016=pt.calculate_amp_only_press(amp_comps=Amps, equationP="P_Mutch2016")
P_Mutch2016.head()


# ## Example 2 - Amphibole-only temperatures
# - Unlike pressures, some thermometers are sensitive to P
# - here, we calculate tempeatures for P=5 kbar

# In[9]:


P_Put2016_5kbar=pt.calculate_amp_only_temp(amp_comps=Amps, equationT="T_Put2016_eq6", P=5)-273.15
P_Put2016_5kbar.head()


# ## Example 3 - Iterating P and T
# - We can also iterate pressure and temperature - most barometers are actually T-insensitive, so it isn't actually iterating, but instead calculates the pressure first, then uses this to calculate the temperature (which does depend on P).

# In[10]:


PT_Mutch=pt.calculate_amp_only_press_temp(amp_comps=Amps, equationT="T_Put2016_eq6", 
                                          equationP="P_Mutch2016")
PT_Mutch.head()


# In[11]:


# Same, but iterating Ridolfi 2021 and Ridolfi 2012
PT_Ridolfi=pt.calculate_amp_only_press_temp(amp_comps=Amps, equationT="T_Ridolfi2012", 
                                            equationP="P_Ridolfi2021")
PT_Ridolfi.head()


# In[12]:


# As above, but without quality filters
PT_Ridolfi=pt.calculate_amp_only_press_temp(amp_comps=Amps, equationT="T_Ridolfi2012", 
                                            equationP="P_Ridolfi2021", Ridolfi_Filter=False)
PT_Ridolfi.head()


# ## Example 4 - Amphibole-liquid thermometry
# - as for Amp-only, except also specify liquid composition
# - help function details your options

# In[13]:


help(pt.calculate_amp_liq_temp)


# ### Using T_Put2016_eq4a_amp_sat. 
# - Note, this calculates the temperature at which a liquid is saturated in amphibole, not the temperature of the liquid

# In[14]:


eq4a_T=pt.calculate_amp_liq_temp(liq_comps=Liqs, amp_comps=Amps, 
                                 equationT="T_Put2016_eq4a_amp_sat")-273.15
eq4a_T.head()


# ### Assessing equilibrium 
# - As for other functions, you can also specify eq_tests=True, to get the Kd Fe-Mg value Putirka suggest may be a useful equilibrium test (+0.28+-0.11)

# In[15]:


eq4a_T_EqTest=pt.calculate_amp_liq_temp(liq_comps=Liqs, amp_comps=Amps, 
                                        equationT="T_Put2016_eq4a_amp_sat", 
                                        eq_tests=True)
eq4a_T_EqTest.head()


# ### Overwriting H2O_Liq.
# - These equations are a little sensitive to melt H2O, because they are parameterized in terms of hydrous fractions, and some contain a separate T term too. You can overwrite the H2O_Liq value in the entered spreadsheet

# In[16]:


eq4a_T_1H2O=pt.calculate_amp_liq_temp(liq_comps=Liqs, amp_comps=Amps,
                        equationT="T_Put2016_eq4a_amp_sat", H2O_Liq=1)-273.15
eq4a_T_1H2O.head()


# ## Example 5 - Amphibole-Liquid barometry
# - Here we are overwriting the input H2O content with 1, if you remove this, reads from spreadsheet

# In[17]:


eq7a_T_1H2O=pt.calculate_amp_liq_press(liq_comps=Liqs, amp_comps=Amps, 
                            equationP="P_Put2016_eq7a", H2O_Liq=1)
eq7a_T_1H2O.head()


# ## Example 6 - Iterating equations for pressure and temperature for amphibole liquid
# - The function prints to tell you actually you selected 2 functions independent of P and T, so only 1 iterative step is being done.

# In[18]:


eq7a_eq4b=pt.calculate_amp_liq_press_temp(liq_comps=Liqs, amp_comps=Amps, 
                                          equationP="P_Put2016_eq7a", 
                                          equationT="T_Put2016_eq4b")
eq7a_eq4b.head()


# ### Both are actually pretty sensitive to water. 

# In[19]:


eq7a_eq4b_1H2O=pt.calculate_amp_liq_press_temp(liq_comps=Liqs, amp_comps=Amps, 
                                               equationP="P_Put2016_eq7a", 
                                               equationT="T_Put2016_eq4b", 
                                               H2O_Liq=1)
eq7a_eq4b_1H2O.head()


# ### can also specify eq_tests=True, and it calculates Kd Fe-Mg, and then classifies within this is within the suggested range of Putirka (2008)

# In[20]:


eq7a_eq4b_1H2O_EqTests=pt.calculate_amp_liq_press_temp(liq_comps=Liqs, amp_comps=Amps, 
                                                       equationP="P_Put2016_eq7a", 
                                                       equationT="T_Put2016_eq4b", 
                                                       H2O_Liq=1, eq_tests=True)
eq7a_eq4b_1H2O_EqTests.head()


# ## Example 8 - Calculating melt compositions from amphibole compositions.
# - You need to specify amphibole compositions.
# - For some equations, e.g., equation 3 and 5 from Zhang, Putirka (2016) SiO2 content you also need to specify temperature. 
# - If you don't, it returns a warning telling you that you won't get results for 3 equations unless you enter a T.

# In[21]:


Melt_Comps=pt.calculate_amp_only_melt_comps(amp_comps=Amps)
Melt_Comps.head()


# - Calculations at a fixed temperature (1300 Kelvin!)

# In[29]:


Melt_Comps2=pt.calculate_amp_only_melt_comps(amp_comps=Amps, T=1300)
Melt_Comps2.head()


# ### Or, we can use an amphibole only thermobarometers to get temperature, to then feed into the chemometry function

# In[30]:


PT_Rid=pt.calculate_amp_only_press_temp(amp_comps=Amps, equationP="P_Ridolfi2021", 
                                        equationT="T_Ridolfi2012", Ridolfi_Filter=False)
Melt_Comps3=pt.calculate_amp_only_melt_comps(amp_comps=Amps, T=PT_Rid['T_K_calc'])
Melt_Comps3.head()


# ## Comparing calculated melts from different equations

# In[31]:


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.plot(Melt_Comps3['SiO2_Eq1_Zhang17'], Melt_Comps3['SiO2_Eq10_Put2016'], 'ok')
ax1.set_xlabel('SiO$_2$ wt% Eq 1 Zhang17')
ax1.set_ylabel('SiO$_2$ wt% Eq 10 Putirka 2016')

ax2.plot(Melt_Comps3['SiO2_Eq1_Zhang17'], Melt_Comps3['TiO2_Eq6_Zhang17'], 'ok')
ax2.set_xlabel('SiO$_2$ wt% Eq 1 Zhang17')
ax2.set_ylabel('TiO$_2$ wt% Eq 6 Zhang17')

ax1.plot([50, 85], [50, 85], '-r')


# In[ ]: