%Magma Density 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 Bottinga & Weill - ref below

%Bottinga & Weill (1970).
%Densities of liquid silicate systems calculated from partial molar volumes
%of oxide components.
%American Journal of Science 269, pp. 169-182

%I have also used the following Excel programs as a guide
%John D. Winter - Available at: http://www.whitman.edu/geology/winter/
%"Magma Density spreadsheet"
%GabbroSoft - http://www.gabbrosoft.org/spreadsheets.html
%"MAGMA-DENSITY"

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

function [density_gas,density_no_gas,totalweight]=density(Wsio2,Wtio2,Wal2o3,Wfeo,Wmgo,Wcao,Wna2o,Wk2o,Wh2o,mTc) 
    %1 wt% of composition
    
totalElements=(Wsio2+Wtio2+Wal2o3+Wfeo+Wmgo+Wcao+Wna2o+Wk2o);

    %2 wt% of gases
totalGas=Wh2o;
    
    %3 Totals (can't exceed 100%)
totalweight=(totalElements+totalGas);

%Constants
    %1 Molecular weights 
Msio2=60085; Mtio2=79899; Mal2o3=101961; Mfeo=71846; Mmgo=40304; Mcao=56079; Mna2o=61979; Mk2o=94203; Mh2o=18015.28; 

    %2 Molecular volume at 1400C
MVsio2=26.75;  MVtio2=22.45;  MVal2o3=37.8; MVfeo=13.94; MVmgo=12.32; MVcao=16.59; MVna2o=29.03; MVk2o=46.3;
MVh2o=17; 

    %3 Coefficient of thermal expansion (*10^5)
CTsio2=0.000001; CTtio2=0.000371; CTal2o3=0.000026; CTfeo=0.000347; CTmgo=0.000122; CTcao=0.000167; CTna2o=0.000259; CTk2o=0.000359; CTh2o=0; 

%Calculations
    %Step 1 Mol Prop
a=Wsio2/Msio2; b=Wtio2/Mtio2; c=Wal2o3/Mal2o3; d=Wfeo/Mfeo; e=Wmgo/Mmgo; f=Wcao/Mcao;
g=Wna2o/Mna2o; h=Wk2o/Mk2o; j=Wh2o/Mh2o; 
new=(a+b+c+d+e+f+g+h+j);
    
    %Step 2 Mol Fraction * Mol Vol
a2=(a/new)*(MVsio2*(1+(CTsio2*(1400-mTc)))); b2=(b/new)*(MVtio2*(1+(CTtio2*(1400-mTc))));
c2=(c/new)*(MVal2o3*(1+(CTal2o3*(1400-mTc)))); d2=(d/new)*(MVfeo*(1+(CTfeo*(1400-mTc))));
e2=(e/new)*(MVmgo*(1+(CTmgo*(1400-mTc)))); f2=(f/new)*(MVcao*(1+(CTcao*(1400-mTc))));
g2=(g/new)*(MVna2o*(1+(CTna2o*(1400-mTc)))); h2=(h/new)*(MVk2o*(1+(CTk2o*(1400-mTc))));
j2=(j/new)*(MVh2o*(1+(CTh2o*(1400-mTc))));
new2=(a2+b2+c2+d2+e2+f2+g2+h2+j2);

    %Step 3 Mol weight* weight 
a3=(a/new)*Msio2; b3=(b/new)*Mtio2; c3=(c/new)*Mal2o3; d3=(d/new)*Mfeo;
e3=(e/new)*Mmgo; f3=(f/new)*Mcao; g3=(g/new)*Mna2o; h3=(h/new)*Mk2o;
j3=(j/new)*Mh2o; 
new3=(a3+b3+c3+d3+e3+f3+g3+h3+j3);

density_gas=(new3/new2);
density_no_gas=(a3+b3+c3+d3+e3+f3+g3+h3)/(a2+b2+c2+d2+e2+f2+g2+h2); %without so2 and co2

clear Wsio2 Wtio2 Wal2o3 Wfeo Wmgo Wcao Wna2o Wk2o Wco2 Wh2o Wso2
clear Msio2 Mtio2 Mal2o3 Mfeo Mmgo Mcao Mna2o Mk2o Mco2 Mh2o Mso2
clear MVsio2 MVtio2 MVal2o3 MVfeo MVmgo MVcao MVna2o MVk2o MVco2 MVh2o MVso2
clear CTsio2 CTtio2 CTal2o3 CTfeo CTmgo CTcao CTna2o CTk2o CTco2 CTh2o CTso2
clear a b c d e f g h i j k new a2 b2 c2 d2 e2 f2 g2 h2 i2 j2 k2 new2
clear a3 b3 c3 d3 e3 f3 g3 h3 i3 j3 k3 new3
end