close all, 
clear all
%% Community parameters
    S = 15; % number of species in the community
    x_f = 0; % fraction of habitat in reserve
    c = 5; % propagule production rate
    r = 1; % recovery rate
    h_scroll = 0:.05:5; % harvest rates 
    w = 0.05; % cost of effort
    facil_choice = 2; % choose which method to use to calculate; method 1 is encouraged for protected diversity...
        % versus diversity gain as the accuracy is higher; otherwise, method
        % 2 is much faster to run
        
%% Other parameters
t_end = 200; % number of time steps to run through to find equilibrium solution

piF = zeros(3,length(h_scroll)); % pre-allocation
alphaF = zeros(1,length(h_scroll)); 
piN = zeros(1,length(h_scroll)); 
alphaN = zeros(1,length(h_scroll)); 

%% Null model
x_phi = h_scroll*(1-x_f)./(r+h_scroll);
        
a = -c; % terms in quadratic
b = (-h_scroll - c.*x_f + c*(1-x_phi-x_f));
d = c*x_f.*(1-x_phi-x_f);

x_else = (-b-sqrt(b.^2-(4*a*d)))./(2*a); % the + sign

alpha_sum = zeros(1,length(h_scroll)); 
     
            for l = 1:S
                temp1 = factorial(S)/(factorial(l)*factorial(S-l))...
                    *(x_else./(1-x_phi-x_f)).^l.*(1-(x_else)./(1-x_phi-x_f)).^(S-l).*(1-x_phi-x_f)*l; % what proportion of patches will have l species * l species
                    alpha_sum = alpha_sum + temp1;                     
            end

        l = 0; % calculates how many patches are empty, but habitable
        temp2 = factorial(S)/(factorial(l)*factorial(S-l))...
                    *(x_else./(1-x_phi-x_f)).^l.*(1-(x_else)./(1-x_phi-x_f)).^(S-l).*(1-x_phi-x_f); 
        alphaN = alpha_sum + x_f*S; % calculates alpha diversity
        piN = h_scroll.*alpha_sum - h_scroll*w; % calculates profit, for each species contributing 1 to a patch's value

%% Facilitation model
for j = 1:length(h_scroll)
    h = h_scroll(j);     
    % Facilitation Model: Equilibrium Method 1 (Simulation)
    if facil_choice == 1
        x0 = zeros(1,S+2)+(1-x_f)/(S+2); % creates an initial condition

        n = length(x0)-2;
        A = [zeros(n,2), triu(ones(n,n))];
        res = 1-sum(x0);

        [~,x] = ode45(@facil,[0 t_end], x0);

        x = x(end,:); % gives the proportion of patches in [x_phi, x_0, x_1,...x_S]
    end

    % Facilation Model: Equilibrium Method 2 (Quadratic)
    if facil_choice == 2
        x_phi = h*(1-x_f) / (h+r); % defines x_phi; % added rp (reserve proportion?) 1/10/12

        % solve for x_0 
        a = c; % terms in the quadratic
        b = -h - c + c*x_phi;
        e = r*x_phi; 

        x_0 = (-b - sqrt(b^2 - 4*a*e))/(2*a); % uses the negative solution only to solve

        % solve for everything else
        x = zeros(1,S); % pre-allocates space for x

        for m = 1:S-1 % cycles through each m                          
            y_i = sum(x(1:m-1)); % sums all below the current index
            a = c;
            b = -h - c*(1-x_phi-x_0) + c*y_i; 

            if m>1
                  e = c*x(m-1)*(1-x_phi-x_0-y_i);
            else
                  e = c*x_0*(1-x_phi-x_0-y_i); 
            end
            x(m) = (-b - sqrt(b^2 - 4*a*e))/(2*a);
        end
        
        if sum(x<0)>0 % if any solutions are negative
            idx = find(x<0); 
            x(idx) = 1 - (x_phi + x_0 + sum(x(1:idx-1))); 
        else
            x(S) = 1 - (x_phi + x_0 + sum(x(1:S-1))); 
        end
        x = real(x); % removes non-real part, which is O(10^-16)      
        x = [x_phi x_0 x]; % gives the proportion of patches in [x_phi, x_0, x_1,...x_S]
        x(end) = x(end) - x_f; % sets x_S to not include fraction in reserves      
    end
    %% Faciliation: Calculation of Diversity and Profit

    % price vectors for patches in each state; each column represents a
    % different value of kappa
    p_i(:,1) = cumsum(1*0.9.^[0:S-1]'.*[zeros(S,1)+1]); % kappa = 0.9
    p_i(:,2)= cumsum(1*1.^[0:S-1]'.*[zeros(S,1)+1]); % kappa = 1
    p_i(:,3) = cumsum(1*1.1.^[0:S-1]'.*[zeros(S,1)+1]); % kappa = 1.1

    for i = 1:3 % runs through each pricing vector
        p = p_i(:,i); 
        piF(i,j) = h*(sum([0; 0; p].*x')-w);
    end

        x(end) = x(end) + x_f; % adds the reserve fraction to the x_S proportion for diversity calculation
        alphaF(j) = [0 0 1:S]*x'; 
        
end
   
subplot(2,1,1), 
plot(h_scroll,piN,'LineWidth',2), hold on, plot(h_scroll,piF,'--','LineWidth',2)
axis([0 max(h_scroll) 0 max(piN)])
xlabel('h'), ylabel('\pi'), title('Profit v. h'), legend('Null','Facilitation, k=0.9','Facilitation, k=1','Facilitation, k=1.1')

subplot(2,1,2),
plot(h_scroll,alphaN,'LineWidth',2), hold on, plot(h_scroll,alphaF,'--','LineWidth',2)
axis([0 max(h_scroll) 0 S])
xlabel('h'), ylabel('\alpha'), title('\alpha diversity v. h'), legend('Null','Facilitation')

