%Magma Viscosity Calculation 
%Author: Tom D Pering - University of Sheffield
%You are free to use and alter with acknowledgement
%The calculations are based on the method of Shaw

%Shaw(1972).
%Viscosities of magmatic silicate liquids: an empirical method of prediction.
%American Journal of Science 272, pp. 870-893

%I have also used the following Excel program as a guide
%John D. Winter - Available at: http://www.whitman.edu/geology/winter/
%"Magma Viscosity spreadsheet" - Click Tech Support

%Inputs
%Input the wt% of each element listed and the wt% of water (Wh2o) and the magma
%temperature in celsius (mTc)

%Outputs
%Viscosity of the magma in Pa s-1 with and without gas (H2O) and total
%weight of the magma

function [viscosity_gas,viscosity_no_gas,totalweight]=viscosity(Wsio2,Wtio2,Wal2o3,Wfe2o3,Wfeo,Wmno,Wmgo,Wcao,Wna2o,Wk2o,Wp2o5,Wh2o,mTc)

totalElements=(Wsio2+Wtio2+Wal2o3++Wfe2o3+Wfeo+Wmno+Wmgo+Wcao+Wna2o+Wk2o+Wp2o5);
totalGas=Wh2o;
totalweight=(totalElements+totalGas);

%Constants
    %1 Molecular weights 
Msio2=60.085; Mtio2=79.899; Mal2o3=101.961; Mfe2o3=159.69; Mfeo=71.846; Mmno=70.9374; 
Mmgo=40.304; Mcao=56.079; Mna2o=61.979; Mk2o=94.203; Mp2o5=283.89; Mh2o=18.01528; 

    %2 Mol. Prop. 
a=Wsio2/Msio2;
b=Wtio2/Mtio2;
c=Wal2o3/Mal2o3;
d=Wfe2o3/Mfe2o3;
e=Wfeo/Mfeo;
f=Wmno/Mmno;
g=Wmgo/Mmgo;
h=Wcao/Mcao;
i=Wna2o/Mna2o;
j=Wk2o/Mk2o;
k=Wp2o5/Mp2o5;
l=Wh2o/Mh2o;

total_m_prop=a+b+c+d+e+f+g+h+i+j+k+l;
total_m_prop_no_gas=a+b+c+d+e+f+g+h+i+j+k;

    %3 Mol. Fraction
a2=a/total_m_prop; a2alt=a/total_m_prop_no_gas;
b2=b/total_m_prop; b2alt=b/total_m_prop_no_gas;
c2=c/total_m_prop; c2alt=c/total_m_prop_no_gas;
d2=d/total_m_prop; d2alt=d/total_m_prop_no_gas;
e2=e/total_m_prop; e2alt=e/total_m_prop_no_gas;
f2=f/total_m_prop; f2alt=f/total_m_prop_no_gas;
g2=g/total_m_prop; g2alt=g/total_m_prop_no_gas;
h2=h/total_m_prop; h2alt=h/total_m_prop_no_gas;
i2=i/total_m_prop; i2alt=i/total_m_prop_no_gas;
j2=j/total_m_prop; j2alt=j/total_m_prop_no_gas;
k2=k/total_m_prop; k2alt=k/total_m_prop_no_gas;
l2=l/total_m_prop;

total_m_frac=a2+b2+c2+d2+e2+f2+g2+h2+i2+j2+k2+l2;
total_m_frac_no_gas=a2alt+b2alt+c2alt+d2alt+e2alt+f2alt+g2alt+h2alt+i2alt+j2alt+k2alt;

%Si * Xsio2 * Xi
b3=(4.5*a2)*b2; b3alt=(4.5*a2alt)*b2alt;
c3=(6.7*a2)*c2; c3alt=(6.7*a2alt)*c2alt;
d3=(3.4*a2)*d2; d3alt=(3.4*a2alt)*d2alt;
e3=(3.4*a2)*e2; e3alt=(3.4*a2alt)*e2alt;
f3=(3.4*a2)*f2; f3alt=(3.4*a2alt)*f2alt;
g3=(3.4*a2)*g2; g3alt=(3.4*a2alt)*g2alt;
h3=(4.5*a2)*h2; h3alt=(4.5*a2alt)*h2alt;
i3=(2.8*a2)*i2; i3alt=(2.8*a2alt)*i2alt;
j3=(2.8*a2)*j2; j3alt=(2.8*a2alt)*j2alt;
l3=(2.0*a2)*l2;

total_si=b3+c3+d3+e3+f3+g3+h3+i3+j3+l3;
total_si_alt=b3alt+c3alt+d3alt+e3alt+f3alt+g3alt+h3alt+i3alt+j3alt;
slope=total_si/(1-a2); 
slope_alt=total_si_alt/(1-a2alt);
viscosity_gas=exp(slope*(10000/(mTc+273.15))-1.5*slope-6.4)/10;
viscosity_no_gas=exp(slope_alt*(10000/(mTc+273.15))-1.5*slope_alt-6.4)/10;

end