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

# # Calculating Viscosity from liquid compositions
# - This notebook shows how to calculate Viscosity using Giordano et al. (2008)
# - You can download the Excel spreadsheet from: https://github.com/PennyWieser/Thermobar/blob/main/docs/Examples/Other_features/Viscoity_Giordano.xlsx

# In[1]:


# If you haven't done so, pip install Thermobar by removing the # symbol
#!pip install Thermobar


# In[2]:


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


# ## Lets load in some melt compositions from a MELTS model published in Wieser et al. (2022)

# In[3]:


Liqs_import2=pt.import_excel('Viscoity_Giordano.xlsx', sheet_name='MELTSTest', suffix="_Liq")
Liqs2=Liqs_import2['Liqs']
Liqs_input2=Liqs_import2['my_input']


# ## Inspect the liquid data you have loaded in to make sure it makes sense

# In[4]:


Liqs2.head()


# ## Lets calculate viscosity at the temperature stored in the column "Temp HT1987_K"
# - Here, we had already calculated temperature using Helz and Thornber, which was stored in the input spreadsheet in a column named 'Temp HT1987_K'
# - The dataframe Liqs_input2 contains all input columns, so we can access the values stored in this column using Liqs_input2['Temp HT1987_K']
# - This temperature needs to be in Kelvin!

# In[5]:


Calc_ExcelT=pt.calculate_viscosity_giordano_2008(liq_comps=Liqs2, 
                                                  T=Liqs_input2['Temp HT1987_K'])


# In[10]:


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.plot(Calc_ExcelT['MgO_Liq'], Calc_ExcelT['n_melt'], '-k')
ax2.plot(Calc_ExcelT['SiO2_Liq'], Calc_ExcelT['n_melt'], '-k')

ax1.set_ylabel('Viscosity (PaS)')
ax1.set_xlabel('MgO content Liq (Wt%)')
ax2.set_xlabel('SiO2 content Liq (Wt%)')


# ## Using a different thermometer for temperature
# - You can get a list of all thermometers in Thermobar using the help function

# In[11]:


help(pt.calculate_liq_only_temp)


# ### Lets use "T_Put2008_eq13"

# In[19]:


CalcT_eq13=pt.calculate_liq_only_temp(liq_comps=Liqs2, equationT="T_Put2008_eq13")
Calc_Puteq13=pt.calculate_viscosity_giordano_2008(liq_comps=Liqs2, 
                                                  T=CalcT_eq13)


# In[20]:


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.plot(Calc_Puteq13['MgO_Liq'], Calc_Puteq13['n_melt'], '-r', label='TPut13')
ax2.plot(Calc_Puteq13['SiO2_Liq'], Calc_Puteq13['n_melt'], '-r')


ax1.plot(Calc_ExcelT['MgO_Liq'], Calc_ExcelT['n_melt'], '-k', label='HT87')
ax2.plot(Calc_ExcelT['SiO2_Liq'], Calc_ExcelT['n_melt'], '-k')
ax1.legend()
ax1.set_ylabel('Viscosity (PaS)')
ax1.set_xlabel('MgO content Liq (Wt%)')
ax2.set_xlabel('SiO2 content Liq (Wt%)')


# ## With different F2O contents
# - By default, we perform calculations with no F, to use the same input structure as the rest of the liquids
# - However, Giordano parameterize in terms of F2O, so you can enter this straight in the function
# - We have 2 functions, allowing you to convert from F2O to F and back

# In[21]:


F2O_calc=pt.convert_F_to_F2O(F_ppm=1000)
F2O_calc


# In[22]:


F_calc=pt.convert_F2O_to_F_ppm(F2O_wt=F2O_calc)
F_calc


# In[23]:


WithF=pt.calculate_viscosity_giordano_2008(liq_comps=Liqs2, 
                                    T=Liqs_input2['Temp HT1987_K'],
                                     F2O_content=F2O_calc)


# In[25]:


plt.plot( Calc_ExcelT['MgO_Liq'], Calc_ExcelT['n_melt'], '-k', label='No F')
plt.plot( WithF['MgO_Liq'], WithF['n_melt'], '-g', label='With 1000 ppm F')
plt.legend()
plt.ylabel('Viscosity')
plt.xlabel('MgO content Liq')