This page was generated from docs/Examples/S6_S2_Corrections/CalcS6ST_Muth.ipynb. Interactive online version: Binder badge.

Python Notebook Download

Comparing calculated S6/ST to measured (Muth and Wallace, 2021)

[4]:

#!pip install PySulfSat --upgrade
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import PySulfSat as ss
pd.options.display.max_columns = None
ss.__version__

[4]:

'1.0.5'

Load data (measured MI comps)

[5]:

df_out=ss.import_data('Muth_data_Merged.xlsx', sheet_name='Sheet1')

df_out_trim=df_out.iloc[:, 0:12]
df_out_trim['Sample_ID_Liq']=df_out['MI_Name']
df_out_trim.head()

We have replaced all missing liquid oxides and strings with zeros.

[5]:
SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Sample_ID_Liq
0 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.235010 0.158692 BBL-5-32
1 49.217861 1.288750 16.749303 9.614724 0.171826 6.033436 9.214457 3.699514 0.716508 0.332001 1.232055 0.213707 BBL-5-33
2 51.666197 1.123500 15.420533 9.144514 0.157054 6.318138 10.347291 3.389791 0.470203 0.214040 1.209455 0.181930 BBL-5-34
3 49.585428 1.246884 16.657728 9.493343 0.186833 7.117724 8.805293 3.562566 0.792121 0.280157 0.968013 0.173842 BBL-5-43
4 48.499094 0.981204 18.784365 7.927122 0.143673 7.262274 11.025627 2.775515 0.278404 0.114633 1.024796 0.142214 BBL-5-44 
[6]:

# Loading errors from Muth supplement
df_err=ss.import_data_noise('Muth_data_Merged.xlsx', sheet_name='Sheet1')
df_err_trim=df_err.iloc[:, 0:12]
df_err_trim['Sample_ID_Liq']=df_err['MI_Name']
df_err_trim.head()

We have replaced all missing liquid oxides and strings with zeros.

[6]:
SiO2_Liq_Err TiO2_Liq_Err Al2O3_Liq_Err FeOt_Liq_Err MnO_Liq_Err MgO_Liq_Err CaO_Liq_Err Na2O_Liq_Err K2O_Liq_Err P2O5_Liq_Err H2O_Liq_Err Fe3Fet_Liq_Err Sample_ID_Liq
0 0.678665 0.096405 0.404524 0.303896 0.019579 0.520599 0.047518 0.264134 0.030958 0.023520 0.112563 0.018260 BBL-5-32
1 0.825791 0.026556 0.115552 0.420127 0.007783 0.152587 0.172061 0.212874 0.096533 0.030083 0.113527 0.014959 BBL-5-33
2 1.876183 0.018804 0.083276 0.194047 0.018462 0.149600 0.136053 0.187045 0.020508 0.010884 0.025529 0.017848 BBL-5-34
3 0.314276 0.036322 0.313526 0.424967 0.004634 0.043281 0.065173 0.272982 0.036137 0.017175 0.001549 0.009052 BBL-5-43
4 0.386098 0.036339 0.161871 0.387210 0.001353 0.057336 0.127017 0.097441 0.003957 0.005722 0.079326 0.010618 BBL-5-44

Lets convert with Berry et al. (2018) following O’Neill and Mavrogenes (2022)

[Fe3+/ΣFe]Berry = ([Fe3+/ΣFe]Zhang - 0.048] / (1-0.048)

[7]:

df_out_trim['Fe3Fet_Berry']=(df_out_trim['Fe3Fet_Liq']-0.048)/(1-0.048)
df_err_trim['Fe3Fet_Berry_Err']=df_err_trim['Fe3Fet_Liq_Err']
df_out_trim['Fe3Fet_Berry'].head()

[7]:

0    0.116273
1    0.174062
2    0.140683
3    0.132187
4    0.098964
Name: Fe3Fet_Berry, dtype: float64

Lets calculate a temperature

  • Using their H2O contents.

[9]:

import Thermobar as pt
Temp_3=pt.calculate_liq_only_temp(liq_comps=df_out, equationT='T_Put2008_eq22_BeattDMg', P=3, H2O_Liq=df_out['H2O_Liq'])
Temp_3.head()

[9]:

0    1407.340081
1    1427.006146
2    1426.318219
3    1459.082125
4    1440.880509
dtype: float64

Lets perform calcs straight up first

[10]:

calc_GivenFe3_Rep=ss.calculate_OM2022_S6St(df=df_out_trim, T_K=Temp_3,
                    Fe3Fet_Liq=df_out_trim['Fe3Fet_Liq'])

calc_GivenFe3_Rep.head()

[10]:
S6St_Liq LnCS2_calc LnCS6_calc LnKSO2S2 LnS6S2 deltaQFM_calc SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Sample_ID_Liq Fe3Fet_Berry SiO2_Liq_mol_frac MgO_Liq_mol_frac MnO_Liq_mol_frac FeOt_Liq_mol_frac CaO_Liq_mol_frac Al2O3_Liq_mol_frac Na2O_Liq_mol_frac K2O_Liq_mol_frac TiO2_Liq_mol_frac P2O5_Liq_mol_frac Si_Liq_cat_frac Mg_Liq_cat_frac Mn_Liq_cat_frac Fet_Liq_cat_frac Ca_Liq_cat_frac Al_Liq_cat_frac Na_Liq_cat_frac K_Liq_cat_frac Ti_Liq_cat_frac P_Liq_cat_frac Mg_Number_Liq_NoFe3 Mg_Number_Liq_Fe3 logfo2_calc Fe2_Liq_cat_frac
0 0.213088 -3.316331 14.165511 -19.356183 -1.306411 0.938474 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.235010 0.158692 BBL-5-32 0.116273 0.546628 0.089540 0.001584 0.084517 0.115094 0.108101 0.038423 0.004335 0.010465 0.001312 0.474433 0.077714 0.001375 0.073355 0.099894 0.187648 0.066697 0.007525 0.009083 0.002277 0.514420 0.557370 -8.282959 0.061714
1 0.856566 -3.054445 13.864131 -18.817435 1.787057 1.604158 49.217861 1.288750 16.749303 9.614724 0.171826 6.033436 9.214457 3.699514 0.716508 0.332001 1.232055 0.213707 BBL-5-33 0.174062 0.539109 0.098521 0.001594 0.088074 0.108143 0.108113 0.039284 0.005006 0.010618 0.001539 0.467189 0.085377 0.001381 0.076324 0.093716 0.187380 0.068087 0.008677 0.009202 0.002668 0.527985 0.587220 -7.371922 0.060013
2 0.551477 -3.163752 13.828337 -18.836028 0.206641 1.249443 51.666197 1.123500 15.420533 9.144514 0.157054 6.318138 10.347291 3.389791 0.470203 0.214040 1.209455 0.181930 BBL-5-34 0.140683 0.552219 0.100671 0.001422 0.081738 0.118496 0.097125 0.035123 0.003206 0.009032 0.000968 0.485928 0.088586 0.001251 0.071925 0.104271 0.170931 0.061814 0.005642 0.007948 0.001704 0.551889 0.600875 -7.735105 0.058840
3 0.555898 -2.805222 13.380534 -17.970140 0.224530 1.221990 49.585428 1.246884 16.657728 9.493343 0.186833 7.117724 8.805293 3.562566 0.792121 0.280157 0.968013 0.173842 BBL-5-43 0.132187 0.535713 0.114638 0.001710 0.085773 0.101928 0.106052 0.037313 0.005459 0.010133 0.001281 0.465794 0.099676 0.001487 0.074579 0.088625 0.184422 0.064886 0.009493 0.008810 0.002228 0.572004 0.617985 -7.368103 0.061614
4 0.162507 -3.215526 13.566557 -18.446271 -1.639692 0.801218 48.499094 0.981204 18.784365 7.927122 0.143673 7.262274 11.025627 2.775515 0.278404 0.114633 1.024796 0.142214 BBL-5-44 0.098964 0.523668 0.116897 0.001314 0.071580 0.127555 0.119522 0.029052 0.001917 0.007969 0.000524 0.454962 0.101560 0.001142 0.062189 0.110820 0.207680 0.050482 0.003332 0.006924 0.000910 0.620210 0.655621 -8.005795 0.053345 
[11]:

calc_GivenFe3_Berry=ss.calculate_OM2022_S6St(df=df_out_trim, T_K=Temp_3,
                    Fe3Fet_Liq=df_out_trim['Fe3Fet_Berry'])

calc_GivenFe3_Berry.head()

[11]:
S6St_Liq LnCS2_calc LnCS6_calc LnKSO2S2 LnS6S2 deltaQFM_calc SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Sample_ID_Liq Fe3Fet_Berry SiO2_Liq_mol_frac MgO_Liq_mol_frac MnO_Liq_mol_frac FeOt_Liq_mol_frac CaO_Liq_mol_frac Al2O3_Liq_mol_frac Na2O_Liq_mol_frac K2O_Liq_mol_frac TiO2_Liq_mol_frac P2O5_Liq_mol_frac Si_Liq_cat_frac Mg_Liq_cat_frac Mn_Liq_cat_frac Fet_Liq_cat_frac Ca_Liq_cat_frac Al_Liq_cat_frac Na_Liq_cat_frac K_Liq_cat_frac Ti_Liq_cat_frac P_Liq_cat_frac Mg_Number_Liq_NoFe3 Mg_Number_Liq_Fe3 logfo2_calc Fe2_Liq_cat_frac
0 0.015495 -3.316331 14.201996 -19.356183 -4.151629 0.312720 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.235010 0.158692 BBL-5-32 0.116273 0.546628 0.089540 0.001584 0.084517 0.115094 0.108101 0.038423 0.004335 0.010465 0.001312 0.474433 0.077714 0.001375 0.073355 0.099894 0.187648 0.066697 0.007525 0.009083 0.002277 0.514420 0.557370 -8.908713 0.064825
1 0.446942 -3.054445 13.899122 -18.817435 -0.213034 1.162246 49.217861 1.288750 16.749303 9.614724 0.171826 6.033436 9.214457 3.699514 0.716508 0.332001 1.232055 0.213707 BBL-5-33 0.174062 0.539109 0.098521 0.001594 0.088074 0.108143 0.108113 0.039284 0.005006 0.010618 0.001539 0.467189 0.085377 0.001381 0.076324 0.093716 0.187380 0.068087 0.008677 0.009202 0.002668 0.527985 0.587220 -7.813834 0.063039
2 0.098903 -3.163752 13.862660 -18.836028 -2.209475 0.717337 51.666197 1.123500 15.420533 9.144514 0.157054 6.318138 10.347291 3.389791 0.470203 0.214040 1.209455 0.181930 BBL-5-34 0.140683 0.552219 0.100671 0.001422 0.081738 0.118496 0.097125 0.035123 0.003206 0.009032 0.000968 0.485928 0.088586 0.001251 0.071925 0.104271 0.170931 0.061814 0.005642 0.007948 0.001704 0.551889 0.600875 -8.267211 0.061807
3 0.089050 -2.805222 13.415669 -17.970140 -2.325286 0.660675 49.585428 1.246884 16.657728 9.493343 0.186833 7.117724 8.805293 3.562566 0.792121 0.280157 0.968013 0.173842 BBL-5-43 0.132187 0.535713 0.114638 0.001710 0.085773 0.101928 0.106052 0.037313 0.005459 0.010133 0.001281 0.465794 0.099676 0.001487 0.074579 0.088625 0.184422 0.064886 0.009493 0.008810 0.002228 0.572004 0.617985 -7.929417 0.064720
4 0.007370 -3.215526 13.597361 -18.446271 -4.903005 0.085910 48.499094 0.981204 18.784365 7.927122 0.143673 7.262274 11.025627 2.775515 0.278404 0.114633 1.024796 0.142214 BBL-5-44 0.098964 0.523668 0.116897 0.001314 0.071580 0.127555 0.119522 0.029052 0.001917 0.007969 0.000524 0.454962 0.101560 0.001142 0.062189 0.110820 0.207680 0.050482 0.003332 0.006924 0.000910 0.620210 0.655621 -8.721103 0.056034 
[15]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10,5))
ax1.plot([0, 1], [0, 1], '-r')
ax2.plot([0, 1], [0, 1], '-r')

ax1.plot(df_out['S6+/∑S'],
         calc_GivenFe3_Rep['S6St_Liq'],
         'dk', mfc='c', label='OM2022')

ax2.plot(df_out['S6+/∑S'],
         calc_GivenFe3_Berry['S6St_Liq'],
         'dk', mfc='r', label='OM2022 Berry correction')


ax1.legend()
ax2.legend()
ax1.set_xlabel('Meas S6/ST')
ax1.set_ylabel('Calc S6/St')
ax2.set_xlabel('Meas S6/ST')
ax2.set_ylabel('Calc S6/St')

[15]:

Text(0, 0.5, 'Calc S6/St')
../../_images/Examples_S6_S2_Corrections_CalcS6ST_Muth_12_1.png 
[16]:

pt.calculate_R2(df_out['S6+/∑S'],
         calc_GivenFe3_Berry['S6St_Liq'], xy=False)

[16]:

{'R2': '0.46',
 'RMSE': '0.42',
 'RMSE_num': 0.41869,
 'P_val': '0.000',
 'Median': '-0.30',
 'Mean': '-0.32',
 'Int': array([-0.31040758]),
 'Grad': array([0.98574948])}

[17]:

pt.calculate_R2(df_out['S6+/∑S'],
         calc_GivenFe3_Rep['S6St_Liq'], xy=False)

[17]:

{'R2': '0.51',
 'RMSE': '0.19',
 'RMSE_num': 0.18975,
 'P_val': '0.000',
 'Median': '-0.03',
 'Mean': '-0.03',
 'Int': array([0.2016873]),
 'Grad': array([0.70240837])}

Now lets account for error

[19]:

# Dataframe of errors for input composition.
df_err_trim.head()

[19]:
SiO2_Liq_Err TiO2_Liq_Err Al2O3_Liq_Err FeOt_Liq_Err MnO_Liq_Err MgO_Liq_Err CaO_Liq_Err Na2O_Liq_Err K2O_Liq_Err P2O5_Liq_Err H2O_Liq_Err Fe3Fet_Liq_Err Sample_ID_Liq Fe3Fet_Berry_Err
0 0.678665 0.096405 0.404524 0.303896 0.019579 0.520599 0.047518 0.264134 0.030958 0.023520 0.112563 0.018260 BBL-5-32 0.018260
1 0.825791 0.026556 0.115552 0.420127 0.007783 0.152587 0.172061 0.212874 0.096533 0.030083 0.113527 0.014959 BBL-5-33 0.014959
2 1.876183 0.018804 0.083276 0.194047 0.018462 0.149600 0.136053 0.187045 0.020508 0.010884 0.025529 0.017848 BBL-5-34 0.017848
3 0.314276 0.036322 0.313526 0.424967 0.004634 0.043281 0.065173 0.272982 0.036137 0.017175 0.001549 0.009052 BBL-5-43 0.009052
4 0.386098 0.036339 0.161871 0.387210 0.001353 0.057336 0.127017 0.097441 0.003957 0.005722 0.079326 0.010618 BBL-5-44 0.010618 
[ ]:

N_dups=5000
# Lets set temp error at +-50 K
df_out_trim['T_K_Liq']=Temp_3 # Set at calc temp
df_err_trim['T_K_Liq_Err']=50 # Add error

df_noisy_abs=ss.add_noise_2_dataframes(df_values=df_out_trim, df_err=df_err_trim,
        error_type="Abs", error_dist="normal", N_dups=N_dups, sample_name_col='Sample_ID_Liq')
df_noisy_abs.head()

Yay. columns match in 2 dataframes
SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Fe3Fet_Berry T_K_Liq Sample_ID_Liq
0 49.464891 1.405919 15.872579 9.592507 0.144889 5.201385 9.795767 3.864573 0.649167 0.283707 1.132997 0.166138 0.096248 1462.385232 BBL-5-32
1 49.433187 1.292106 16.499790 9.491953 0.151015 4.379767 9.794181 3.509948 0.619213 0.261549 1.312720 0.164751 0.110122 1426.471496 BBL-5-32
2 50.739490 1.424251 16.882461 8.785630 0.210249 5.648638 9.797237 3.370011 0.615548 0.326386 1.118686 0.161864 0.103601 1430.304368 BBL-5-32
3 49.781261 1.173928 16.377589 9.863109 0.190366 5.474283 9.838843 3.853985 0.597666 0.267870 1.221709 0.150209 0.131069 1347.358580 BBL-5-32
4 48.178742 1.259443 17.053891 9.337702 0.175699 6.027648 9.811798 3.737781 0.579071 0.275786 1.188427 0.136557 0.122865 1480.237478 BBL-5-32

Now put all this synthetic simulated data into the S6 calculation

[ ]:

noisy_ONeill_Fe_S6St_Rep=ss.calculate_OM2022_S6St(df=df_noisy_abs, Fe3Fet_Liq=df_noisy_abs['Fe3Fet_Liq'],
                                          T_K=df_noisy_abs['T_K_Liq'])
noisy_ONeill_Fe_S6St_Rep.head()

noisy_ONeill_Fe_S6St_Berry=ss.calculate_OM2022_S6St(df=df_noisy_abs, Fe3Fet_Liq=df_noisy_abs['Fe3Fet_Berry'],
                                      T_K=df_noisy_abs['T_K_Liq'])
noisy_ONeill_Fe_S6St_Berry.head()
S6St_Liq LnCS2_calc LnCS6_calc LnKSO2S2 LnS6S2 deltaQFM_calc SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Fe3Fet_Berry T_K_Liq Sample_ID_Liq SiO2_Liq_mol_frac MgO_Liq_mol_frac MnO_Liq_mol_frac FeOt_Liq_mol_frac CaO_Liq_mol_frac Al2O3_Liq_mol_frac Na2O_Liq_mol_frac K2O_Liq_mol_frac TiO2_Liq_mol_frac P2O5_Liq_mol_frac Si_Liq_cat_frac Mg_Liq_cat_frac Mn_Liq_cat_frac Fet_Liq_cat_frac Ca_Liq_cat_frac Al_Liq_cat_frac Na_Liq_cat_frac K_Liq_cat_frac Ti_Liq_cat_frac P_Liq_cat_frac Mg_Number_Liq_NoFe3 Mg_Number_Liq_Fe3 logfo2_calc Fe2_Liq_cat_frac
0 0.003156 -2.707814 13.542211 -17.885017 -5.755152 -0.110739 49.464891 1.405919 15.872579 9.592507 0.144889 5.201385 9.795767 3.864573 0.649167 0.283707 1.132997 0.166138 0.096248 1462.385232 BBL-5-32 0.546265 0.085632 0.001355 0.088592 0.115909 0.103295 0.041374 0.004573 0.011679 0.001326 0.474778 0.074425 0.001178 0.076998 0.100741 0.179555 0.071919 0.007949 0.010150 0.002305 0.491496 0.536850 -8.662046 0.069588
1 0.008395 -3.092411 13.813302 -18.831883 -4.771689 0.186181 49.433187 1.292106 16.499790 9.491953 0.151015 4.379767 9.794181 3.509948 0.619213 0.261549 1.312720 0.164751 0.110122 1426.471496 BBL-5-32 0.554646 0.073258 0.001435 0.089065 0.117744 0.109094 0.038178 0.004432 0.010905 0.001242 0.481068 0.063540 0.001245 0.077250 0.102124 0.189245 0.066227 0.007688 0.009458 0.002155 0.451302 0.496153 -8.796480 0.068743
2 0.007376 -3.225447 13.675510 -18.728544 -4.902084 0.134270 50.739490 1.424251 16.882461 8.785630 0.210249 5.648638 9.797237 3.370011 0.615548 0.326386 1.118686 0.161864 0.103601 1430.304368 BBL-5-32 0.551512 0.091530 0.001936 0.079862 0.114100 0.108136 0.035511 0.004268 0.011645 0.001502 0.479819 0.079631 0.001684 0.069480 0.099268 0.188159 0.061789 0.007426 0.010131 0.002613 0.534031 0.577595 -8.801322 0.062282
3 0.028592 -3.756132 15.333922 -21.097335 -3.525604 0.521847 49.781261 1.173928 16.377589 9.863109 0.190366 5.474283 9.838843 3.853985 0.597666 0.267870 1.221709 0.150209 0.131069 1347.358580 BBL-5-32 0.543116 0.089035 0.001759 0.089991 0.115012 0.105294 0.040762 0.004159 0.009634 0.001237 0.471680 0.077324 0.001528 0.078154 0.099885 0.182889 0.070801 0.007224 0.008367 0.002149 0.497324 0.537942 -9.492156 0.067910
4 0.036951 -2.575433 13.289339 -17.431583 -3.260501 0.406455 48.178742 1.259443 17.053891 9.337702 0.175699 6.027648 9.811798 3.737781 0.579071 0.275786 1.188427 0.136557 0.122865 1480.237478 BBL-5-32 0.530943 0.099026 0.001640 0.086057 0.115855 0.110750 0.039932 0.004071 0.010440 0.001286 0.459278 0.085660 0.001419 0.074442 0.100217 0.191602 0.069085 0.007042 0.009031 0.002226 0.535026 0.571301 -7.938221 0.065295

Now calculate the error bar for each row

[ ]:

Stats_Fe_S6_Rep=pt.av_noise_samples_series(calc=noisy_ONeill_Fe_S6St_Rep['S6St_Liq'], sampleID=df_noisy_abs['Sample_ID_Liq'])
Stats_Fe_S6_Rep.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 0.250082 0.209301 0.175566 0.893046 0.001502
1 BBL-5-33 5000 0.817118 0.854546 0.132496 0.991158 0.147848
2 BBL-5-34 5000 0.533757 0.547817 0.209747 0.970988 0.008911
3 BBL-5-43 5000 0.539196 0.557356 0.214825 0.965811 0.011851
4 BBL-5-44 5000 0.204286 0.155821 0.167652 0.880138 0.000666 
[ ]:

Stats_Fe_S6_Berry=pt.av_noise_samples_series(calc=noisy_ONeill_Fe_S6St_Berry['S6St_Liq'], sampleID=df_noisy_abs['Sample_ID_Liq'])
Stats_Fe_S6_Berry.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 0.030180 0.015406 0.042787 0.525115 0.000012
1 BBL-5-33 5000 0.443754 0.434528 0.212082 0.980997 0.015766
2 BBL-5-34 5000 0.138767 0.099460 0.127927 0.802016 0.000249
3 BBL-5-43 5000 0.130735 0.088121 0.127792 0.837637 0.000504
4 BBL-5-44 5000 0.018730 0.007451 0.032474 0.445296 0.000001

Lets do the same using Nash

  • If you just want to vary 1 thing, you can do it this way

  • First, take your dataframe, and duplicate it N times, all the columns will be the same, but will now be Sample1-Sample1-Sample1, Sample2-Sample2-Sample2

[ ]:

Dupdf=ss.duplicate_dataframe(df=df_out_trim, N_dup=N_dups)
Dupdf.head()
SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Sample_ID_Liq Fe3Fet_Berry T_K_Liq
0 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.23501 0.158692 BBL-5-32 0.116273 1407.340081
1 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.23501 0.158692 BBL-5-32 0.116273 1407.340081
2 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.23501 0.158692 BBL-5-32 0.116273 1407.340081
3 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.23501 0.158692 BBL-5-32 0.116273 1407.340081
4 49.928525 1.270821 16.755591 9.230927 0.170783 5.486115 9.811559 3.620198 0.620793 0.283074 1.23501 0.158692 BBL-5-32 0.116273 1407.340081

Now you can make the errors for one column at a time if you prefer

  • E.g. adding a temp error of +- 1sigma = 50 K

[ ]:

Temp_Err=ss.add_noise_series(df_out_trim['T_K_Liq'], error_var=50,
error_type="Abs", error_dist="normal", N_dup=N_dups)
# Then add this to the dataframe
Dupdf['T_K_MC']=Temp_Err

  • And adding a Fe3Fet_Liq error from their spreadsheet

[ ]:

Fe3_Err=ss.add_noise_series(df_out_trim['Fe3Fet_Liq'], error_var=df_err_trim['Fe3Fet_Liq_Err'],
error_type="Abs", error_dist="normal", N_dup=N_dups)
##
Dupdf['Fe3Fet_Liq_MC']=Fe3_Err

Now put these uncertainties into Nash

[ ]:

noisy_Nash_S6St=ss.calculate_S6St_Nash2019(Fe3Fet_Liq=Dupdf['Fe3Fet_Liq_MC'],
                                          T_K=Dupdf['T_K_MC'])
noisy_Nash_S6St.head()

0    0.396509
1    0.397516
2    0.030916
3    0.046192
4    0.194181
dtype: float64

[ ]:

Stats_Nash_S6=pt.av_noise_samples_series(calc=noisy_Nash_S6St, sampleID=Dupdf['Sample_ID_Liq'])
Stats_Nash_S6.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 0.195022 0.134288 0.181301 0.929295 0.000372
1 BBL-5-33 5000 0.750935 0.809510 0.196740 0.995808 0.039207
2 BBL-5-34 5000 0.458549 0.443234 0.249592 0.994307 0.004288
3 BBL-5-43 5000 0.485901 0.492004 0.251244 0.981179 0.002699
4 BBL-5-44 5000 0.152452 0.092361 0.161316 0.871681 0.000131

And Jugo

  • First, need to calculate QFM relative to frost from the Fe3Fet ratio

[ ]:

Buffer=pt.convert_fe_partition_to_fo2(liq_comps=df_noisy_abs,  T_K=Dupdf['T_K_MC'], P_kbar=5,
        model="Kress1991", Fe3Fet_Liq=Dupdf['Fe3Fet_Liq_MC'],
 renorm=False)
Buffer.head()

overwriting Fe3Fet_Liq to that specified in the function input
(150000,)
DeltaQFM_Frost1991 DeltaNNO_Frost1991 fo2_calc SiO2_Liq TiO2_Liq Al2O3_Liq FeOt_Liq MnO_Liq MgO_Liq CaO_Liq Na2O_Liq K2O_Liq P2O5_Liq H2O_Liq Fe3Fet_Liq Fe3Fet_Berry T_K_Liq Sample_ID_Liq FeO_Liq Fe2O3_Liq
0 0.060093 -0.462598 2.812926e-08 49.464891 1.405919 15.872579 9.592507 0.144889 5.201385 9.795767 3.864573 0.649167 0.283707 1.132997 0.155280 0.096248 1462.385232 BBL-5-32 8.102986 1.655022
1 0.369503 -0.148512 1.026760e-08 49.433187 1.292106 16.499790 9.491953 0.151015 4.379767 9.794181 3.509948 0.619213 0.261549 1.312720 0.174378 0.110122 1426.471496 BBL-5-32 7.836766 1.839095
2 -0.180464 -0.692690 3.440496e-10 50.739490 1.424251 16.882461 8.785630 0.210249 5.648638 9.797237 3.370011 0.615548 0.326386 1.118686 0.145696 0.103601 1430.304368 BBL-5-32 7.505600 1.422254
3 -0.346166 -0.863513 1.546122e-09 49.781261 1.173928 16.377589 9.863109 0.190366 5.474283 9.838843 3.853985 0.597666 0.267870 1.221709 0.134465 0.131069 1347.358580 BBL-5-32 8.536862 1.473606
4 0.316464 -0.196554 1.446050e-09 48.178742 1.259443 17.053891 9.337702 0.175699 6.027648 9.811798 3.737781 0.579071 0.275786 1.188427 0.176803 0.122865 1480.237478 BBL-5-32 7.686766 1.834370 
[ ]:

noisy_Jugo=ss.calculate_S6St_Jugo2010_eq10(deltaQFM=Buffer['DeltaQFM_Frost1991'])
noisy_Jugo.head()

0    0.010367
1    0.041734
2    0.003448
3    0.001611
4    0.032988
Name: DeltaQFM_Frost1991, dtype: float64

[ ]:

Stats_Jugo_S6=pt.av_noise_samples_series(calc=noisy_Jugo, sampleID=Dupdf['Sample_ID_Liq'])
Stats_Jugo_S6.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 0.023367 0.010717 0.037736 0.650903 0.000033
1 BBL-5-33 5000 0.405935 0.386462 0.223679 0.954401 0.006194
2 BBL-5-34 5000 0.090598 0.054602 0.103943 0.845817 0.000320
3 BBL-5-43 5000 0.100409 0.063135 0.110625 0.775579 0.000495
4 BBL-5-44 5000 0.012257 0.004779 0.022856 0.356391 0.000006 
[ ]:

S_types=ss.convert_S_types(S_wt=df_out['S'])
S_types.head()
S_wt S_ppm SO2_wt SO2_ppm SO3_wt SO3_ppm SO4_wt SO4_ppm
0 0.133403 1334.034269 0.266528 2665.281066 0.333090 3330.904464 0.399653 3996.527863
1 0.150144 1501.439419 0.299974 2999.741572 0.374889 3748.892648 0.449804 4498.043724
2 0.124034 1240.340506 0.247809 2478.089313 0.309696 3096.963717 0.371584 3715.838121
3 0.108941 1089.410000 0.217654 2176.543672 0.272011 2720.110507 0.326368 3263.677343
4 0.093533 935.330000 0.186871 1868.705623 0.233539 2335.393434 0.280208 2802.081245 
[ ]:

fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, figsize=(12,6))
sam='BBL-5-32'
bins=50
ax1.hist(df_noisy_abs['T_K_Liq'].loc[Dupdf['Sample_ID_Liq']==sam],  fc='black', alpha=0.3, bins=bins);
ax2.hist(df_noisy_abs['Fe3Fet_Liq'].loc[Dupdf['Sample_ID_Liq']==sam], fc='black', alpha=0.3, bins=bins);
ax3.hist(Buffer['DeltaQFM_Frost1991'].loc[Buffer['Sample_ID_Liq']==sam],  fc='black', alpha=0.3, bins=bins);
ax4.hist(df_noisy_abs['FeOt_Liq'].loc[df_noisy_abs['Sample_ID_Liq']==sam], fc='black', alpha=0.3, bins=bins);
ax5.hist(df_noisy_abs['MgO_Liq'].loc[df_noisy_abs['Sample_ID_Liq']==sam],  fc='black', alpha=0.3, bins=bins);
ax6.hist(df_noisy_abs['Na2O_Liq'].loc[df_noisy_abs['Sample_ID_Liq']==sam], fc='black', alpha=0.3, bins=bins);
ax1.set_ylabel('# of simulations')
ax4.set_ylabel('# of simulations')
ax1.annotate("T (K)", xy=(0.87, 0.93), xycoords="axes fraction", fontsize=10)
ax2.annotate("Fe$^{3+}$/Fe$_{T}$", xy=(0.75, 0.93), xycoords="axes fraction", fontsize=10)
ax3.annotate("ΔQFM", xy=(0.8, 0.93), xycoords="axes fraction", fontsize=10)
ax4.annotate("FeOt Liq", xy=(0.75, 0.93), xycoords="axes fraction", fontsize=10)
ax5.annotate("MgO Liq", xy=(0.75, 0.93), xycoords="axes fraction", fontsize=10)
ax6.annotate("Na$_2$O Liq", xy=(0.75, 0.93), xycoords="axes fraction", fontsize=10)

ax1.annotate("A", xy=(0.02, 0.93), xycoords="axes fraction", fontsize=12)
ax2.annotate("B", xy=(0.02, 0.93), xycoords="axes fraction", fontsize=12)
ax3.annotate("C", xy=(0.02, 0.93), xycoords="axes fraction", fontsize=12)
ax4.annotate("D", xy=(0.02, 0.93), xycoords="axes fraction", fontsize=12)
ax5.annotate("E", xy=(0.02, 0.93), xycoords="axes fraction", fontsize=12)
ax6.annotate("F", xy=(0.02, 0.93), xycoords="axes fraction", fontsize=12)

fig.savefig('MonteCarlo_Muth_input.png', dpi=300)
../../_images/Examples_S6_S2_Corrections_CalcS6ST_Muth_38_0.png

Figure for text showing S6/St with errorbars

  • Using reported Fe3FeT

[ ]:


fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10,3.5))
ms1=5




ax1.errorbar(df_out['S6+/∑S'],
         Stats_Fe_S6_Rep['Mean_calc'],
            xerr=0.01, yerr=Stats_Fe_S6_Rep['St_dev_calc'],
               fmt='s', ecolor='grey', elinewidth=0.8, mfc='cyan', ms=ms1, mec='k', capsize=3, label='usingFe')


ax2.errorbar(df_out['S6+/∑S'], Stats_Nash_S6['Mean_calc'],
            xerr=0.01, yerr=Stats_Nash_S6['St_dev_calc'],
               fmt='s', ecolor='grey', elinewidth=0.8, mfc='blue', ms=ms1, mec='k', capsize=3, label='usingFe')

ax3.errorbar(df_out['S6+/∑S'], Stats_Jugo_S6['Mean_calc'],
            xerr=0.01, yerr=Stats_Jugo_S6['St_dev_calc'],
               fmt='s', ecolor='grey', elinewidth=0.8, mfc='red', ms=ms1, mec='k', capsize=3, label='usingFe')

# s=ax1.scatter(df_out['S6+/∑S'],
#          Stats_Fe_S6_Rep['Mean_calc'], marker='s', s=50, edgecolors='k', linewidths=0.5,
#             c=S_types['S_ppm'], cmap='hot', zorder=100)
#fig.colorbar(s, ax=ax3)

ax1.plot([0, 1], [0, 1], '-k')
ax2.plot([0, 1], [0, 1], '-k')
ax3.plot([0, 1], [0, 1], '-k')
ax1.set_xlabel('Measured S$^{6+}$/S$_T$')
ax1.set_ylabel('Calc S$^{6+}$/S$_T$ (OM2022)')
ax2.set_xlabel('Measured S$^{6+}$/S$_T$')
ax2.set_ylabel('Calc S$^{6+}$/S$_T$ (Nash2019)')
ax3.set_xlabel('Measured S$^{6+}$/S$_T$')
ax3.set_ylabel('Calc S$^{6+}$/S$_T$ (Jugo2010)')

# Print stats on fig
O22_stats=pt.calculate_R2(df_out['S6+/∑S'], Stats_Fe_S6_Rep['Mean_calc'], xy=False)
ax1.annotate('G', xy=(0.02, 0.94), xycoords="axes fraction", fontsize=12)

ax1.annotate('R$^{2}$='+str(O22_stats['R2']), xy=(0.02, 0.88), xycoords="axes fraction", fontsize=10)
ax1.annotate('RMSE='+str(O22_stats['RMSE']), xy=(0.02, 0.83), xycoords="axes fraction", fontsize=10)

Nash_stats=pt.calculate_R2(df_out['S6+/∑S'], Stats_Nash_S6['Mean_calc'], xy=False)
ax2.annotate('H', xy=(0.02, 0.94), xycoords="axes fraction", fontsize=12)

ax2.annotate('R$^{2}$='+str(Nash_stats['R2']), xy=(0.02, 0.88), xycoords="axes fraction", fontsize=10)
ax2.annotate('RMSE='+str(Nash_stats['RMSE']), xy=(0.02, 0.83), xycoords="axes fraction", fontsize=10)


Jugo_stats=pt.calculate_R2(df_out['S6+/∑S'], Stats_Jugo_S6['Mean_calc'], xy=False)
ax3.annotate('I', xy=(0.02, 0.94), xycoords="axes fraction", fontsize=12)
ax3.annotate('R$^{2}$='+str(Jugo_stats['R2']), xy=(0.02, 0.88), xycoords="axes fraction", fontsize=10)
ax3.annotate('RMSE='+str(Jugo_stats['RMSE']), xy=(0.02, 0.83), xycoords="axes fraction", fontsize=10)

fig.tight_layout()
fig.savefig('Muth_outputs_TheirFe.png', dpi=300, Transparent=True)

C:\Users\penny\AppData\Local\Temp\ipykernel_29520\3752010702.py:56: MatplotlibDeprecationWarning: savefig() got unexpected keyword argument "Transparent" which is no longer supported as of 3.3 and will become an error in 3.6
  fig.savefig('Muth_outputs_TheirFe.png', dpi=300, Transparent=True)
../../_images/Examples_S6_S2_Corrections_CalcS6ST_Muth_40_1.png

Berry Corrected Fe3

[ ]:


fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10,3.5))
ms1=5




ax1.errorbar(df_out['S6+/∑S'],
         Stats_Fe_S6_Berry['Mean_calc'],
            xerr=0.01, yerr=Stats_Fe_S6_Berry['St_dev_calc'],
               fmt='s', ecolor='grey', elinewidth=0.8, mfc='cyan', ms=ms1, mec='k', capsize=3, label='usingFe')


ax2.errorbar(df_out['S6+/∑S'], Stats_Nash_S6['Mean_calc'],
            xerr=0.01, yerr=Stats_Nash_S6['St_dev_calc'],
               fmt='s', ecolor='grey', elinewidth=0.8, mfc='blue', ms=ms1, mec='k', capsize=3, label='usingFe')

ax3.errorbar(df_out['S6+/∑S'], Stats_Jugo_S6['Mean_calc'],
            xerr=0.01, yerr=Stats_Jugo_S6['St_dev_calc'],
               fmt='s', ecolor='grey', elinewidth=0.8, mfc='red', ms=ms1, mec='k', capsize=3, label='usingFe')

# s=ax1.scatter(df_out['S6+/∑S'],
#          Stats_Fe_S6_Berry['Mean_calc'], marker='s', s=50, edgecolors='k', linewidths=0.5,
#             c=S_types['S_ppm'], cmap='hot', zorder=100)
#fig.colorbar(s, ax=ax3)

ax1.plot([0, 1], [0, 1], '-k')
ax2.plot([0, 1], [0, 1], '-k')
ax3.plot([0, 1], [0, 1], '-k')
ax1.set_xlabel('Measured S$^{6+}$/S$_T$')
ax1.set_ylabel('Calc S$^{6+}$/S$_T$ (OM2022)')
ax2.set_xlabel('Measured S$^{6+}$/S$_T$')
ax2.set_ylabel('Calc S$^{6+}$/S$_T$ (Nash2019)')
ax3.set_xlabel('Measured S$^{6+}$/S$_T$')
ax3.set_ylabel('Calc S$^{6+}$/S$_T$ (Jugo2010)')

# Print stats on fig
O22_stats=pt.calculate_R2(df_out['S6+/∑S'], Stats_Fe_S6_Berry['Mean_calc'], xy=False)
ax1.annotate('f) Berry Corrected Fe', xy=(0.02, 0.95), xycoords="axes fraction", fontsize=10)

ax1.annotate('R$^{2}$='+str(O22_stats['R2']), xy=(0.02, 0.88), xycoords="axes fraction", fontsize=10)
ax1.annotate('RMSE='+str(O22_stats['RMSE']), xy=(0.02, 0.83), xycoords="axes fraction", fontsize=10)

Nash_stats=pt.calculate_R2(df_out['S6+/∑S'], Stats_Nash_S6['Mean_calc'], xy=False)
ax2.annotate('g)', xy=(0.02, 0.95), xycoords="axes fraction", fontsize=10)

ax2.annotate('R$^{2}$='+str(Nash_stats['R2']), xy=(0.02, 0.88), xycoords="axes fraction", fontsize=10)
ax2.annotate('RMSE='+str(Nash_stats['RMSE']), xy=(0.02, 0.83), xycoords="axes fraction", fontsize=10)


Jugo_stats=pt.calculate_R2(df_out['S6+/∑S'], Stats_Jugo_S6['Mean_calc'], xy=False)
ax3.annotate('h)', xy=(0.02, 0.95), xycoords="axes fraction", fontsize=10)
ax3.annotate('R$^{2}$='+str(Jugo_stats['R2']), xy=(0.02, 0.88), xycoords="axes fraction", fontsize=10)
ax3.annotate('RMSE='+str(Jugo_stats['RMSE']), xy=(0.02, 0.83), xycoords="axes fraction", fontsize=10)

fig.tight_layout()
fig.savefig('Muth_outputs_Berry.png', dpi=200)
../../_images/Examples_S6_S2_Corrections_CalcS6ST_Muth_42_0.png

We can also do Monte Carlo simulations for the SCSS (and all other calculations)

Lets make a monte carlo simulation for the sulfide composition, say error is +-0.05 Fe/Fe+Ni+Cu

[ ]:

df_out['FeFeNiCu']=0.634
FeFeNiCu_Err=ss.add_noise_series(df_out['FeFeNiCu'], error_var=0.05,
error_type="Abs", error_dist="normal", N_dup=N_dups)
##
df_noisy_abs['FeFeNiCu_MC']=FeFeNiCu_Err

[ ]:

SCSS_S2017_MC=ss.calculate_S2017_SCSS(df=df_noisy_abs, T_K=df_noisy_abs['T_K_Liq'], Fe_FeNiCu_Sulf=df_noisy_abs['FeFeNiCu_MC'],
                                   P_kbar=5)

Using inputted Fe_FeNiCu_Sulf ratio for calculations.
Index(['SiO2_Liq', 'TiO2_Liq', 'Al2O3_Liq', 'FeOt_Liq', 'MnO_Liq', 'MgO_Liq',
       'CaO_Liq', 'Na2O_Liq', 'K2O_Liq', 'P2O5_Liq', 'H2O_Liq', 'Fe3Fet_Liq',
       'Fe3Fet_Berry', 'T_K_Liq', 'Sample_ID_Liq', 'FeFeNiCu_MC',
       'Fe_FeNiCu_Sulf_calc'],
      dtype='object')
no non ideal SCSS as no Cu/CuFeNiCu

Now lets average these per sample

[ ]:

Stats_SCSS=pt.av_noise_samples_series(calc=SCSS_S2017_MC['SCSS2_ppm_ideal_Smythe2017'], sampleID=df_noisy_abs['Sample_ID_Liq'])
Stats_SCSS.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 632.551388 627.403814 100.880616 1101.012247 334.989681
1 BBL-5-33 5000 686.024489 676.780913 110.286012 1158.645523 371.995340
2 BBL-5-34 5000 677.220586 669.393024 104.987236 1086.462394 367.535899
3 BBL-5-43 5000 720.617802 715.509319 110.811735 1143.578685 390.024109
4 BBL-5-44 5000 724.757508 717.712174 113.706819 1249.729467 382.031792

Lets also calculate the SCAS

[ ]:

SCAS_ZT_MC=ss.calculate_ZT2022_SCAS(df=df_noisy_abs, T_K=df_noisy_abs['T_K_Liq'],
                                   P_kbar=5)

g:\my drive\berkeley_new\pysulfsat\pysulfsat_structure\src\PySulfSat\scas_calc.py:133: UserWarning: you entered a P_kbar, just be aware this function isnt actually pressure sensitive
  w.warn('you entered a P_kbar, just be aware this function isnt actually pressure sensitive')

[ ]:

STot_MC=ss.calculate_S_Total_SCSS_SCAS(SCSS=SCSS_S2017_MC['SCSS2_ppm_ideal_Smythe2017'],
                    SCAS=SCAS_ZT_MC['SCAS6_ppm'], S6St_Liq=noisy_ONeill_Fe_S6St_Rep['S6St_Liq'])
STot_MC.head()

C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2120: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2120: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2120: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
Total_S_ppm S2_Tot_ppm S6_Tot_ppm deltaQFM S6St_Liq SCSS_2_ppm SCAS_6_ppm SCSS_Tot SCAS_Tot S6 in SCSS_Tot S2 in SCAS_Tot
0 1149.291586 794.178826 355.112760 None 0.308984 794.178826 5669.644288 1149.291586 18349.311020 355.112760 12679.666731
1 880.624032 659.182431 221.441601 None 0.251460 659.182431 4820.661515 880.624032 19170.699437 221.441601 14350.037921
2 923.236593 645.139125 278.097468 None 0.301220 645.139125 4715.631793 923.236593 15655.100554 278.097468 10939.468760
3 592.958742 537.704168 55.254574 None 0.093185 537.704168 3312.416955 592.958742 35546.859584 55.254574 32234.442628
4 960.926667 873.421809 87.504858 None 0.091063 873.421809 6235.937726 960.926667 68479.385157 87.504858 62243.447431

Lets average this per sample

[ ]:

Stats_St=pt.av_noise_samples_series(calc=STot_MC['Total_S_ppm'], sampleID=df_noisy_abs['Sample_ID_Liq'])
Stats_St.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 929.867356 811.565310 449.091181 6779.050566 343.782828
1 BBL-5-33 5000 4144.582047 4160.318969 1796.855724 10518.839850 559.080996
2 BBL-5-34 5000 1923.345374 1489.772388 1307.591223 10757.676809 464.636175
3 BBL-5-43 5000 2116.846271 1619.711620 1497.499821 11077.653585 528.072061
4 BBL-5-44 5000 990.963299 871.981210 465.739670 7406.081985 410.488689 
[ ]:

Stats_S6=pt.av_noise_samples_series(calc=STot_MC['SCAS_6_ppm'], sampleID=df_noisy_abs['Sample_ID_Liq'])
Stats_S6.head()
Stats_S2=pt.av_noise_samples_series(calc=STot_MC['SCSS_2_ppm'], sampleID=df_noisy_abs['Sample_ID_Liq'])
Stats_S2.head()
Sample # averaged Mean_calc Median_calc St_dev_calc Max_calc Min_calc
0 BBL-5-32 5000 632.551388 627.403814 100.880616 1101.012247 334.989681
1 BBL-5-33 5000 686.024489 676.780913 110.286012 1158.645523 371.995340
2 BBL-5-34 5000 677.220586 669.393024 104.987236 1086.462394 367.535899
3 BBL-5-43 5000 720.617802 715.509319 110.811735 1143.578685 390.024109
4 BBL-5-44 5000 724.757508 717.712174 113.706819 1249.729467 382.031792

Lets calculate the STot using the measured S6/ST amount

[ ]:

df_out['S6St']=df_out['S6+/∑S']
df_out.loc[df_out['S6+/∑S']>1]=1
SCSS_S2017=ss.calculate_S2017_SCSS(df=df_out, T_K=Temp_3, Fe_FeNiCu_Sulf=0.634,
                                   P_kbar=5)


SCAS_ZT=ss.calculate_ZT2022_SCAS(df=df_out, T_K=Temp_3,
                                   P_kbar=5)

STot=ss.calculate_S_Total_SCSS_SCAS(SCSS=SCSS_S2017['SCSS2_ppm_ideal_Smythe2017'],
                    SCAS=SCAS_ZT['SCAS6_ppm'], S6St_Liq=df_out['S6+/∑S'])
STot.head()

Using inputted Fe_FeNiCu_Sulf ratio for calculations.
Index(['SiO2_Liq', 'TiO2_Liq', 'Al2O3_Liq', 'FeOt_Liq', 'MnO_Liq', 'MgO_Liq',
       'CaO_Liq', 'Na2O_Liq', 'K2O_Liq', 'P2O5_Liq', 'H2O_Liq', 'Fe3Fet_Liq',
       'Ni_Liq_ppm', 'Cu_Liq_ppm', 'Unnamed: 0.1', 'Unnamed: 0',
       'Tephra Sample', 'Unnamed: 2', 'MI_Name', 'IS6+', 'IS4+', 'IS2-',
       'S6+/∑S', '% Correction', 'Centroid (eV)', 'Centroid s.e.',
       'Fe3Fet_Liq_Err', 'log (fO2)', 'Δ QFM', 'CinderCone', 'TephraSample',
       'Inclusion', 'MinimumInclusionWidth(μm)', 'MaxmimumInclusionWidth(μm)',
       'Vaporbubble?', 'VaporbubbleDiameter*(μm)', 'OxideinMI?',
       'SulfideinMI?', 'Unnamed: 10', 'DateAnalyzedEPMA', 'Unnamed: 12',
       'SiO2_Liq_Err', 'TiO2_Liq_Err', 'Al2O3_Liq_Err', 'FeOt_Liq_Err',
       'MnO_Liq_Err', 'MgO_Liq_Err', 'CaO_Liq_Err', 'Na2O_Liq_Err',
       'K2O_Liq_Err', 'P2O5_Liq_Err', 'S', 'S_Err', 'Cl', 'Cl_Err',
       'Unnamed: 37', 'H2O_Liq_Err', 'Unnamed: 40', 'OlivineHostFo%',
       'Unnamed: 42', 'Li', 's.e.', 'B', 's.e..1', 'Sc', 's.e..2', 'V',
       's.e..3', 'Rb', 's.e..4', 'Sr', 's.e..5', 'Y', 's.e..6', 'Zr', 's.e..7',
       'Nb', 's.e..8', 'Ba', 's.e..9', 'La', 's.e..10', 'Ce', 's.e..11', 'Nd',
       's.e..12', 'Sm', 's.e..13', 'Dy', 's.e..14', 'Yb', 's.e..15', 'Pb',
       's.e..16', 'Sample_ID_Liq', 'FeFeNiCu', 'S6St', 'Fe_FeNiCu_Sulf_calc'],
      dtype='object')
no non ideal SCSS as no Cu/CuFeNiCu

g:\my drive\berkeley_new\pysulfsat\pysulfsat_structure\src\PySulfSat\scas_calc.py:133: UserWarning: you entered a P_kbar, just be aware this function isnt actually pressure sensitive
  w.warn('you entered a P_kbar, just be aware this function isnt actually pressure sensitive')
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2120: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2120: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
C:\Users\penny\anaconda3\lib\site-packages\pandas\core\indexing.py:2120: FutureWarning: In a future version, the Index constructor will not infer numeric dtypes when passed object-dtype sequences (matching Series behavior)
  new_ix = Index(new_ix)
Total_S_ppm S2_Tot_ppm S6_Tot_ppm deltaQFM S6St_Liq SCSS_2_ppm SCAS_6_ppm SCSS_Tot SCAS_Tot S6 in SCSS_Tot S2 in SCAS_Tot
0 2045.786001 625.153851 1420.632150 None 0.694419 625.153851 4364.741895 2045.786001 6285.460924 1420.632150 1920.719029
1 4840.565722 30.709456 4809.856266 None 0.993656 685.060485 4809.856266 107982.384252 4840.565722 107297.323767 30.709456
2 1076.254500 673.318103 402.936397 None 0.374388 673.318103 4977.354145 1076.254500 13294.653541 402.936397 8317.299396
3 955.162067 716.391457 238.770610 None 0.249979 716.391457 5478.824513 955.162067 21917.125172 238.770610 16438.300659
4 1093.450261 719.152091 374.298170 None 0.342309 719.152091 5359.471170 1093.450261 15656.809530 374.298170 10297.338360 
[ ]:

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12,5))

ax1.errorbar(S_types['S_ppm'],
         Stats_St['Mean_calc'], xerr=0, yerr=Stats_St['St_dev_calc'],
             fmt='d', ecolor='k', elinewidth=0.8, mfc='cyan', ms=10, mec='k', capsize=5)
ax1.plot([500, 3000], [500, 3000], '-k')
ax2.plot([500, 3000], [500, 3000], '-k')
ax1.set_xlabel('Measured S (ppm)')
ax1.set_ylabel('S Total (SCSS+SCAS)')




ax2.errorbar((df_out['S6+/∑S'])*S_types['S_ppm'],
         Stats_S6['Mean_calc'], xerr=0, yerr=Stats_S6['St_dev_calc'],
             fmt='d', ecolor='k', elinewidth=0.8, mfc='cyan', ms=10, mec='k', capsize=5)


ax2.plot([0, 6000], [0, 6000], '-r')
ax2.set_xlabel('S6+ measured in melt')
ax2.set_ylabel('SCAS ZT22')

ax3.errorbar((1-(df_out['S6+/∑S']))*S_types['S_ppm'],
         Stats_S2['Mean_calc'], xerr=0, yerr=Stats_S2['St_dev_calc'],
             fmt='d', ecolor='k', elinewidth=0.8, mfc='cyan', ms=10, mec='k', capsize=5)


ax3.plot([0, 800], [0, 800], '-r')
ax3.set_xlabel('S2- measured in melt')
ax3.set_ylabel('SCSS S2017')
ax1.annotate("S-saturated", xy=(0.5, 0.01), xycoords="axes fraction", fontsize=10)
ax1.annotate("S undersaturated", xy=(0.5, 0.95), xycoords="axes fraction", fontsize=10)


ax2.annotate("Supersaturated in S6", xy=(0.3, 0.05), xycoords="axes fraction", fontsize=10)
ax2.annotate("Undersaturated in S6", xy=(0.3, 0.9), xycoords="axes fraction", fontsize=10)

ax3.annotate("Supersaturated in S2", xy=(0.3, 0.05), xycoords="axes fraction", fontsize=10)
ax3.annotate("Undersaturated in S2", xy=(0.3, 0.9), xycoords="axes fraction", fontsize=10)

fig.tight_layout()
../../_images/Examples_S6_S2_Corrections_CalcS6ST_Muth_57_0.png 
[ ]:



[ ]:



[ ]: