% BubbleMagma.m

params.T = 1273;
params.R = 185;

L = 1000;
params.conduit_r = 5;

params.rho_m = 2700;

params.bubble_diam = 1e-3;

params.g = 9.81;

params.Cd = 1;
params.C = 0.01;

params.P0 = params.rho_m*params.g*L;

params.A = params.rho_m*params.R*params.T/params.P0;
params.B = params.rho_m*params.g*params.bubble_diam/params.P0;
params.sigma = params.conduit_r/params.bubble_diam;

phi0 = 0.05;
rho_g0 = params.P0/params.R/params.T/params.rho_m;
u_m0 = 0.1/sqrt(params.P0/params.rho_m);
u_g0 = u_m0 + sqrt(params.B*rho_g0/params.Cd);

f0 = [phi0,rho_g0,u_g0,u_m0]';

opts = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',0.1,'InitialStep',1e-6);

[z,f] = ode45(@(z,f) BubbleMagmaEqns(z,f,params),[0,1],f0,opts);

plot(f(:,1),z);

function dfdz = BubbleMagmaEqns(z,f,params)
%% The governing equations of motion for rise of a bubbly magma
% Inputs:
% z -- the independent variable, the distance along the conduit with 
%   z = 0 denoting the top of the magma chamber
% f -- the vector of dependent variables with
%   f(1) -- the volume fraction of bubbles
%   f(2) -- the density of bubbles
%   f(3) -- the velocity of bubbles
%   f(4) -- the velocity of magma
% params -- a parameters structure containing:
%   g -- gravitational acceleration
%   rho_m -- magma density
%   Cd -- drag coefficient
%   C -- wall friction coefficient
%   r -- conduit diameter
%   bubble_diam -- bubble diameter
% Outputs:
%   dfdz -- the vector of derivatives df/dz

    phi = f(1);   % volume fraction of bubbles
    rho_g = f(2); % density of bubbles
    u_g = f(3);   % velocity of gas phase
    u_m = f(4);   % velocity of magma
    
    A = zeros(4,4); % A is a 4x4 matrix specifying the differential equations
    b = zeros(4,1); % b is a 4x1 vector specifying the right-hand side of the differential equations
    % In matrix form the equations are A*dfdz = b
    A(1,1) = rho_g*u_g;
    A(1,2) = phi*u_g;
    A(1,3) = rho_g*phi;
    A(2,1) = -u_m;
    A(2,4) = (1-phi);
    A(3,2) = params.A;
    A(3,3) = rho_g*u_g;
    A(4,2) = params.A;
    A(4,4) = u_m;
    
    b(3) = -rho_g - params.Cd/params.B*(u_m-u_g)*abs(u_m-u_g);
    b(4) = -1 - params.C/params.B/params.sigma*u_m*u_m;% + params.Cd/params.B*(u_m-u_g)*abs(u_m-u_g)*phi;
    
    dfdz = A\b;
    
    return
    
end