% Control Burn Program
% 
% function [X] = control_burn(years, cbf, growth, replace, death, forest_composition, grab_type)
%
% years = the total number of years you would like to consider
% cbf = control burn frequency (how often you would like to prescribe burns)
% growth = the matrix with growth information for tree population (forest_composition)
% replace = the matrix with reproductive information for tree population
% death = the matrix with death information for the tree population
% grab_type = switch to select GRAB vs FGRAB computing method
%   select 0 for GRAB
%   select 1 for FGRAB
%   default is FGRAB
%

function [X] = control_burn(years, c, growth, replace, death, forest_composition, grab_type);

% Transpose initial forest composition for ease of calculations
x = forest_composition';

% Find sizes of matrices
[M,N] = size(growth);
[L,l] = size(replace);
[m,n] = size(death);

% Make sure the matrices are of the appropriate length
if N ~= l - 1 |n ~= l,
    error = 'control_burn: Matricies are not of the same length.  N = L - 1 and n = L'
    N
    n
    l
    break;
end

% Select either GRAB or FGRAB
if grab_type == 0 
    % Run model using GRAB
    % Run the model over the appropriate number of years
    for I = 0 : years,
        [G,R,B] = grab(I, c, growth, replace, death);
        x = G*x + R*x - B*x;
    
        % Change all negative numbers to zeros
        mask = x < 1;
        x(mask) = 0;
    end
else
    % Run model using FGRAB
    % Run the model over the appropriate number of years
    for I = 0 : years,
        [G,R,B] = fgrab(I, c, growth, replace, death);
        x = G*x + R*x - B*x;
    
        % Change all negative numbers to zeros
        mask = x < 1;
        x(mask) = 0;
    end
end


% Return the forest composition at time years
X = x;