not finished code

% site_index.m

% this program will calculate maximum productivity for ws1 plots

% we base this on the idea that biomass follows a logistic function. the
% logistic function is defined as B(t) = K/(1+Ce^(-r*t)). We would like to
% fit a logistic function to our biomass curves and then solve for where
% the derivative of biomass (productivity) is maximized (or where the
% second derivative of biomass is 0). To do this we must fit a logistic
% function to the biomass data.

% import biomass data

file = importdata('biomat.csv',',');
data = file;

% years after the first measurement
T =[0;4;8;12;16;22;28];

%transpose data

dataT = data';
prod_rate = zeros(6, length(dataT));

%create a vector "years" that holds the interval length between periods
years = [4 ; 4 ; 4 ; 4 ; 6 ; 6];

%calculate the difference in biomass-- later just import the ANPP matrix
%but we're doing this for now... no MORTALITY included at the moment!

% sample below shows how to do it for just one plot
prod_rate(1,1) = (dataT(2,1)- dataT(1,1))/years(1);
prod_rate(2,1) = (dataT(3,1)- dataT(2,1))/years(2);
prod_rate(3,1) = (dataT(4,1)- dataT(3,1))/years(3);
prod_rate(4,1) = (dataT(5,1)- dataT(4,1))/years(4);
prod_rate(5,1) = (dataT(6,1)- dataT(5,1))/years(5);
prod_rate(6,1) = (dataT(7,1)- dataT(6,1))/years(6);

% loop through all plots (lengthdataT) and for each interval (i) and store
% the results in the array "prod_rate"
for j = 2:length(dataT)
    for i=1:6
        prod_rate(i,j) = (dataT(i+1,j)-dataT(i,j))/years(i);
    end
end

%create a vector telling us what the maximum productivity is:
inflection_value = zeros(length(prod_rate),1);

inflection_value(1,1) = max(prod_rate(:,1));
for j = 2:length(dataT)
    inflection_value(j,1) = max(prod_rate(:,j));
end

%find that maximum productivity and WHEN IT OCCURS by interval time
% these are our initial guesses for the inflection point
indexOfMaxValue = zeros(length(prod_rate),1);

indexOfMaxValue(1,1) = find(prod_rate(:,1) == max(prod_rate(:,1)));
for j = 2:length(prod_rate)
    indexOfMaxValue(j,1) = find(prod_rate(:,j)== max(prod_rate(:,j)));
end

% to find the slope at the inflection point, which will generate the best
% guess value for "r" we need to find

slopeAtInflection = zeros(length(prod_rate),1);

% find the value one position ahead of the max. If the max occurs at the
% end of the interval, then save that value as the "one ahead" value
oneAhead = zeros(length(dataT),1);

oneAhead(1) = dataT((find(prod_rate(:,1) == max(prod_rate(:,1)))+1),1);
for j = 2:length(dataT)
    if indexOfMaxValue(j)~=6
    oneAhead(j) = dataT((find(prod_rate(:,j) == max(prod_rate(:,j)))+1),j);
    else
        oneAhead(j) = dataT((find(prod_rate(:,j) == max(prod_rate(:,j)))),j);
    end
end

% find the value one position behind the max. If the max occurs at the end
% of the interval it's cool.
oneBehind = zeros(length(dataT),1);
oneBehind(1) = prod_rate((find(prod_rate(:,1) == max(prod_rate(:,1)))-1),1);
for j = 2:length(dataT)
    %if indexOfMaxValue(j)~=6
    oneBehind(j) = dataT((find(prod_rate(:,j) == max(prod_rate(:,j)))-1),j);
   % else
      %  oneBehind(j) = prod_rate((find(prod_rate(:,j) == max(prod_rate(:,j)))),j);
   % end
end

% find the slope at the inflection point. This will help us find R.
slopeAtInflection(1,1) = (oneAhead(1)-oneBehind(1))/((indexOfMaxValue(1) + 1) -(indexOfMaxValue(1) -1));

for j = 2:length(dataT)
    if indexOfMaxValue(j) ~=6
        slopeAtInflection(j) = (oneAhead(j)-oneBehind(j))/ ((indexOfMaxValue(j)+1)-(indexOfMaxValue(j)-1));
    else
        slopeAtInflection(j) =(oneAhead(j)-oneBehind(j))/ ((indexOfMaxValue(j))-(indexOfMaxValue(j)-1));
    end
end

% absolute largest biomass on our watershed serves as the resounding
% asymptote for the initial guess

K = max(max(dataT));

% using the power of differential equations we solve to find that r =
% 4*slopeAtInflectionPoint/K

r = zeros(length(slopeAtInflection),1);
r = (4.*slopeAtInflection)/K;

% we will say that "r" is not less than 0.01 nor greater than 0.90, and we
% will graph in 50 equal increments.
% we will say that the inflection point must be between (logically) 1 and 6

linR=linspace(0.01,0.9,50);
t0=linspace(0,28,50);
[linR,t0]=meshgrid(linR,t0);

% for the first plot this is what you do
[m,n] = size(linR);
e = zeros(size(linR));
for i=1:m
    for j=1:n
        e(i,j) = myerror([linR(i,j); t0(i,j)], T, dataT(:,1));
    end
end

% % images of the first plot
% figure(1)
% mesh(linR, t0, e);
% xlabel('r')
% ylabel('inflection point')
% zlabel('e')
% title('plotting error versus parameters for the first plot');
%
% figure(2)
% contour(linR, t0, e,50)
% xlabel('r')
% ylabel('inflection');
% title('contour map');

% find the minimum error of that plot, the first output is the value of "r"
% and the second output is the inflection point when it occurs.
min = fminsearch(@myerror, [0.1, 4],[],T,dataT(:,1));

% pull out the values for the "r" and "inflection point" from your
% minimized error function
rD1 = min(1);
rD2 = min(2);

% now solve for maximum biomass,which is Krep.

H = 1./(1+exp(-rD1*(T-rD2)));
Krep=(H'*dataT(:,1))/(H'*H);

% for the fourth plot this is what you do
[m4,n4] = size(linR);
e = zeros(size(linR));
for i=1:m4
    for j=1:n4
        e(i,j) = myerror([linR(i,j); t0(i,j)], T, dataT(:,4));
    end
end

% find the minimum error of that plot, the first output is the value of "r"
% and the second output is the inflection point when it occurs.
min4 = fminsearch(@myerror, [0.1, 4],[],T,dataT(:,4));

% pull out the values for the "r" and "inflection point" from your
% minimized error function
rD1 = min4(1);
rD2 = min4(2);

% now solve for maximum biomass,which is Krep.

H = 1./(1+exp(-rD1*(T-rD2)));
Krep4=(H'*dataT(:,1))/(H'*H);