% Plot control burn
% 
% function [pine, hard] = plot_cb(cbf, years, grab_fgrab)
%
% cbf = control burn frequency (how often you would like burns to occur)
% years = total time span you would like to consider
% grab_fgrab = switch to select GRAB vs FGRAB computing method
%   select 0 for GRAB
%   select 1 for FGRAB
%   default is FGRAB
%

function [popsize_pine_new, popsize_hard_new] = plot_cb(cb, total_time, grab_type)

% VARIABLES
%
% popsize_pine = vector containing the population size of the pines for each year
% popsize_pine_new = vector containing the population size of the pines for each year
%   with effects of hardwood crowding
% popsize_hard = vector containing the population size of the hardwoods for each year
% popsize_hard_new = vector containing the population size of the hardwoods for each year
%   with effects of pine tree crowding
%

% Get matrices
cb_matrices;

% Set initial vectors to zero and enter initial pop size 
popsize_pine = 0;
popsize_pine(1) = sum(pine);
popsize_pine_new = 0;
popsize_pine_new(1) = sum(pine);
popsize_hard = 0;
popsize_hard(1) = sum(hard);
popsize_hard_new = 0;
popsize_hard_new(1) = sum(hard);

% Run through how population dynamics change each year over the total years.
for I = 0 : total_time,
    year = I + 0;
    K = I + 1;
    
    % Find size of pine tree population
    [X] = control_burn(year, cb, Gpine, Rpine, Dpine, pine,grab_type);
    popsize_pine(K + 1) = sum(X); 
    
    % Find size of hard wood population
    [Y] = control_burn(year, cb, Ghard, Rhard, Dhard, hard,grab_type);
    popsize_hard(K + 1) = sum(Y) - Y(1);

    % sum(X) = total number of pine trees
    % sum(Y) = total number of hardwoods
    
    % If there are more than three times as many pines than hardwoods
    if sum(X) >= 3.0 * sum(Y)
        Y(1) = 0.00 * Y(1);
        Y(2) = 0.25 * Y(2);
        Y(3) = 0.35 * Y(3);
        popsize_hard_new(K + 1) = sum(Y);
        popsize_pine_new(K + 1) = sum(X);
        
    % If there are between 1.25 and 3 times as many pines as hardwoods    
    elseif sum(X) >= 1.25 * sum(Y) & sum(X) < 3.0 * sum(Y)
        Y(1) = 0.10 * Y(1);
        Y(2) = 0.20 * Y(2);
        Y(3) = 0.30 * Y(3);
        popsize_hard_new(K + 1) = sum(Y);
        popsize_pine_new(K + 1) = sum(X);
        
    % If there are between 1 and 1.25 times as many pines as hardwoods    
    elseif sum(X) >= sum(Y) & sum(X) < 1.25 * sum(Y)
        X = 1.00 .* X;
        popsize_pine_new(K + 1) = sum(X);
    
    % If there are between 1 and 2 times as many hardwoods as pines    
    elseif sum(Y) > sum(X) & sum(Y) < 2.0 * sum(X)
        X = 0.75 .* X;
        popsize_pine_new(K + 1) = sum(X);
    
    % If there are between 2 and 5 times as many hardwoods as pines    
    elseif sum(Y) >= 2.0 * sum(X) & sum(Y) < 5.0 * sum(X)
        X = 0.50 .* X;
        popsize_pine_new(K + 1) = sum(X);
    
    % If there are more than 5 times as may hardwoods as pines    
    else
        X = 0.25 .* X;
        popsize_pine_new(K + 1) = sum(X);

    end
     
    % Carrying capacities
    C_pine = 500;
   
    % Carrying Capacity on pine tree population (somewhat adhoc)
    if sum(X) > C_pine
        X(1) = 0 * X(1);
        X(2) = 0 * X(2);
        X(3) = 0.25 * X(3);
        X(4) = 0.75 * X(4);
        X(5) = 0.90 * X(5);
        popsize_pine_new(K + 1) = sum(X);
    end
   
end

% Determine Lenth of Arrays
pinesize = length(popsize_pine_new);
hardsize = length(popsize_hard);

% Construct an array representing time
years = [1:1:pinesize];

% Check to see that population data was collected for an equal number of years for
%  pines and hardwoods
if pinesize ~= hardsize
    error = 'An error hass occured in calculations'
else
    % Plot populations
    plot(years,popsize_pine_new,'g', years,popsize_hard,'b')
    xlabel('Time in Years')
    ylabel('Number of Trees')
    title('Control Burn Every n Years')
    legend('Pine Trees','Hardwoods');    
end