function soln = PurePlume(k,F00,N);
% PurePlume.m
% A matlab function that solves the equations of motion for a pure plume
% with a point source of buoyancy at z=0.
%
% Inputs: k - the entrainment coefficient (Real number);
%         F00 - the buoyancy flux at z=0 (Real number);
%         N - the buoyancy frequency of the atmosphere (Real number if the atmosphere is statically stable).
%
% Output: soln - a solution structure

% This file (PurePlume.m) is Copyright (C) 2015 Mark Woodhouse.  Permission granted to reproduce for personal and educational use only, provided the author is identified.  Commercial copying, hiring, lending is prohibited.

%% We prefer to work with N2
N2 = N*N;

%% Set the point source boundary conditions
% The point z=0 is a singular point of the differential equations, 
% so we must move away from the point z=0.
Q0 = 6*k/5*(9*k/10)^(1/3)*F00^(1/3);
M0 = (9*k/10)^(2/3)*F00^(2/3);
F1 = -9*N2/20*6*k/5*(9*k/10)^(1/3)*F00^(1/3);

z0 = 1e-8; % Move this tiny distance from the point source.
% The fluxes of mass, momentum and buoyancy at this small distance from the
% point source at z=0 are then given by
Qz0 = Q0.*z0.^(5/3); 
Mz0 = M0.*z0.^(4/3);
Fz0 = F00 + F1.*z0^(8/3);

%% Set a vector containing the initial fluxes of mass, momentum and buoyancy at z0
IC = [Qz0,Mz0,Fz0];

%% Solve the system of equation numerically using Matlab's 4th order Runge-Kutta routine
% Set some solver options:
options = odeset('RelTol',1e-8, ... % Relative error tolerance
    'AbsTol',1e-9.*ones(size(IC)), ... % Absolute error tolerance
    'InitialStep',1e-5, ... % Size of the initial step
    'Nonnegative',[1,2], ... % The mass flux and momentum fluxes are always positive
    'Events',@StopTop); % Set a condition to identify the top and stop the solver, coded in the StopTop function below.
% Solve the system of equations
sol = ode45(@PlumeEqns, ... % Use the ode45 routine, with the differential equations coded in the PlumeEquations function below.
    [z0 100], ... % Integrate the equations between z=z0 and z=100.
    IC, ... % At z=z0 we apply the initial values in the vector IC.
    options); % Use the solver options.

%% Extract some outputs and build our soln structure.

soln.z = sol.x; % The steps taken by the solver.

y = sol.y; % The solution at the steps in soln.z.

% Integral outputs
Q = y(1,:); % The mass flux
M = y(2,:); % The momentum flux
F = y(3,:); % The buoyancy flux

b = Q./sqrt(M); % The plume radius
w = M./Q; % The plume velocity
gr = F./Q; % The reduced gravity of the plume

% Build the solution structure
soln.Q = Q;
soln.M = M;
soln.F = F;
soln.b = b;
soln.w = w;
soln.gr = gr;

return

%% The plume equations
    function dy = PlumeEqns(z,y)
        % This function returns the right hand side of the system of equations,
        % evaluated at z.
        %
        % Inputs: z - the independent variable is the height from z=0 (Real number)
        %         y - the vector of values of the dependent variables (Q,M,F) at z (Real vector)
        %
        % Outputs: dy - a vector of the derivatives (dQ/dz,dM/dz,dF/dz) at z (Real vector)

        dy = zeros(3,1); % dy must be a column vector.

        % Extract Q, M and F from the input y.
        Q = y(1);
        M = y(2);
        F = y(3);
        
        % Write the differential equations:
        % dQ/dz = 2k*sqrt(M)
        dy(1) = 2.*k.*sqrt(M);

        % dM/dz = Q*F/M
        dy(2) = Q.*F./M;

        % dF/dz = -N2*Q
        dy(3) = -N2*Q;
        
        return

    end

    %% Stop integration at the top of the plume %%
    function [VALUE, ISTERMINAL, DIRECTION] = StopTop(z,y)
        % This function is used to stop the integration.
        % We look for the point where the velocity of the plume is zero
        % (but in fact are happy to stop when the velocity gets close to
        % zero).
        %
        % Inputs: z - the independent variable is the height from z=0 [Real number]
        %         y - the vector of values of the dependent variables (Q,M,F) at z [Real vector]
        %
        % Outputs: VALUE - a number that equals zero when we want to stop
        % (so the solver can look for a change in sign of VALUE to stop) [Real number];
        %          ISTERMINAL - a flag to say whether we stop the integration when VALUE changes sign [Real number];
        %          DIRECTION - a flag to say whether we look for +ve to
        %          -ve changes in VALUE, -ve to +ve changes, or any change
        %          across VALUE=0 [Real number].
        
        % Extract Q, M and F from the input y.
        Q = y(1);
        M = y(2);
        F = y(3);
        
        % Extract reduced gravity and velocity from the input y.
        G = F./Q;
        W = M./Q;
        
        % Stop the integration when W becomes less than 5e-3
        VALUE = W-5e-3;
        
        % We want the solver to stop when VALUE cross 0.
        ISTERMINAL = 1;
        
        % Stop if VALUE changes from +ve to -ve or -ve to +ve.
        DIRECTION = 0;
        
        return
        
    end

end
