%Stokes Rise of bubbles with a given gas content
%Author: Tom D Pering - University of Sheffield
%You are free to use and alter with acknowledgement
function [u,rs,slug_tt,Re,Pbubble,Vbubble] = stokesF(asl,aSO2,aCO2,aH2O,rho,eta,mTk,length) 
%asl = height above sea level
%aSO2 is the amount of SO2 released per slug in kg
%aCO2 and aH2O in kg
%rho is the density of the magma
%eta is the visocosity of the magma
%mTk is magma temperature in kelvin
%length is the length of the conduit that you want to model bubble
%speed/travel time for

%NOTE - Stokes rise IS NOT suitable for slug flow (Taylor bubbles)

%Calculating Atmopsheric Pressure - input asl to run
slp=101325; %sea level standard pressure
slt=298.15; %sea level standard temperature
tlr=0.0065; %temperature lapse rate
mmDA=0.0289644; %molar mass of dry air
g=9.81; %acceleration due to gravity;
R=8.3144621; %universal gas constant
%also needs asl (height above sea level)
atmP=slp*(1-((tlr*asl)/slt))^((g*mmDA)/(R*tlr));
clear slp slt tlr mmDA

load gas_properties.mat;
mmSlug=(((aSO2*1000)/mmSO2)+((aCO2*1000)/mmCO2)+((aH2O*1000)/mmH2O)); %calculating molar mass of slug
total=(aSO2+aCO2+aH2O); % total molar mass of slug
propSO2=aSO2/total;
propCO2=aCO2/total;
propH2O=aH2O/total;
dSlug=(dSO2*propSO2)+(dCO2*propCO2)+(dH2O*propH2O); %density of the slug

for depth=drange(1:length)
    Pmagma(depth)=(depth*g*rho)+atmP;
    Pbubble(depth)=(depth*g*dSlug)+atmP;    
end

Vbubble=(R*mTk*mmSlug)./Pmagma;
rs_pre=(Vslug./((4/3)*pi));
rs=rs_pre.^(1/3);
u=(2*(rho-dSlug)*g*(rs.^2))/(9*eta);
slug_tt=1./u; %travel time for each bubble at each metre iteration
slug_ttC=cumsum(slug_tt); %cumulative travel time
clear rs_pre

Re=(2*rho.*u.*rs)./eta; %Reynolds number

clear mmCO mmCO2 mmCl mmF mmH2O mmH2S mmSO2 dCO dCO2 dCl dF dH2O dH2S dSO2 total propSO2 propCO2 propH2O g R
end