%BUG report for COCO
% COCO finds a BP and SN point too much (the equilibrium goes
% assymptotically to zero but does not cross).
%
%  STEP      TIME        ||U||  LABEL  TYPE            E0
%     0  00:00:00  1.0983e+001      4  EP     7.0000e+000
%    10  00:00:00  1.1360e+001      5         7.2772e+000
%    16  00:00:01  1.1381e+001      6  FP     7.2888e+000
%    16  00:00:01  1.1381e+001      7  SN     7.2888e+000
%    20  00:00:01  1.1377e+001      8         7.2808e+000
%    30  00:00:01  1.1161e+001      9         7.0876e+000
%    40  00:00:01  9.2802e+000     10         5.1834e+000
%    45  00:00:02  9.3449e+000     11  FP     5.1534e+000
%    45  00:00:02  9.3449e+000     12  SN     5.1534e+000
%    50  00:00:02  9.6267e+000     13         5.2341e+000
%    59  00:00:02  1.3839e+001     14  SN     7.6242e+000<<<< NOT CORRECT
%    59  00:00:02  1.3840e+001     15  BP     7.5910e+000<<<< NOT CORRECT
%    60  00:00:02  1.7627e+001     16         9.8429e+000
%    61  00:00:03  1.7890e+001     17  EP     1.0000e+001
%
%Vegetation turbidity model
%#Ref: Scheffer, M. (1998) "Ecology of Shallow Lakes" Chapman and Hall. 357 pp.
%
%Eeq = E0.*hV./(hV+V);
%Veq = hE.^p./(hE.^p+E.^p);
%V' = rV.*V.*(1-V./Veq)+inpV;
%E' = rE.*E.*(1-E./Eeq)+inpE;
%
function [bd1] = veg
mod=get_model;
% coco settings
prob = coco_prob();
prob = coco_set(prob, 'ode', 'vectorized','on');
prob = coco_set(prob,'EP', 'NSA',true,'BTP',true);
prob = coco_set(prob,'cont','h0',0.009,'h_min',0.009,'h_max',9,'h_fac_min',0.5,'h_fac_max',2,'PtMX',150);
prob = coco_set(prob,'corr','ItMX',10,'SubItMX',7,'TOL',1e-010,'ResTOL',1e-010);

%run coco 1D
bd1 = coco(prob,'1', 'ode', 'isol', 'ep', ...
       mod.coco_ode, mod.coco_jac, mod.coco_jacp, [0.816;1.378], mod.parnames, mod.p0, ... 
       'E0', [1 10]);

function g_l_res=get_model
%parameters:
E0   = 7;             %Vertical light attenuation without vegetation,m-1
hE   = 2;             %Light attenuation at which the lake is half covered with vegetation,m-1
hV   = 0.2;           %Vegetation cover where the light attenuation is reduced by half,-
inpE = 0;             %,m-1 d-1
inpV = 0;             %,d-1
p    = 4;             %parameter in Hill function,-
rE   = 0.1;           %Growth rate of vertical light attenuation (algae),d-1
rV   = 0.05;          %Growth rate of vegetation,d-1
g_l_res.odehandle=@(t,g_X1)curr_ode0(t,g_X1,E0,hE,hV,inpE,inpV,p,rE,rV);
g_l_res.coco_ode=@(g_X1,g_par)curr_ode0(0,g_X1,g_par(1,:),g_par(2,:),g_par(3,:),g_par(4,:),g_par(5,:),g_par(6,:),g_par(7,:),g_par(8,:));
g_l_res.coco_jac=@(g_X1,g_par)reshape(i_vcat(-g_par(8,:).*(g_X1(1,:).*g_par(2,:).^(-g_par(6,:)).*(g_X1(2,:).^g_par(6,:)+g_par(2,:).^g_par(6,:))-1.0)-g_X1(1,:).*g_par(2,:).^(-g_par(6,:)).*g_par(8,:).*(g_X1(2,:).^g_par(6,:)+g_par(2,:).^g_par(6,:)),-(g_X1(2,:).^2.*g_par(7,:))./(g_par(1,:).*g_par(3,:)),-g_X1(2,:).^(g_par(6,:)-1.0).*g_X1(1,:).^2.*g_par(2,:).^(-g_par(6,:)).*g_par(6,:).*g_par(8,:),-g_par(7,:).*((g_X1(2,:).*(g_X1(1,:)+g_par(3,:)))./(g_par(1,:).*g_par(3,:))-1.0)-(g_X1(2,:).*g_par(7,:).*(g_X1(1,:)+g_par(3,:)))./(g_par(1,:).*g_par(3,:))),[2,2,max(size(g_X1,2),size(g_par,2))]);
g_l_res.coco_jacp=@(g_X1,g_par)reshape(i_vcat(0.0,(1.0./g_par(1,:).^2.*g_X1(2,:).^2.*g_par(7,:).*(g_X1(1,:)+g_par(3,:)))./g_par(3,:),-g_X1(1,:).*g_par(8,:).*(g_X1(1,:).*g_par(2,:).^(-g_par(6,:)).*g_par(2,:).^(g_par(6,:)-1.0).*g_par(6,:)-g_X1(1,:).*g_par(2,:).^(-g_par(6,:)-1.0).*g_par(6,:).*(g_X1(2,:).^g_par(6,:)+g_par(2,:).^g_par(6,:))),0.0,0.0,-g_X1(2,:).*g_par(7,:).*(g_X1(2,:)./(g_par(1,:).*g_par(3,:))-(g_X1(2,:).*1.0./g_par(3,:).^2.*(g_X1(1,:)+g_par(3,:)))./g_par(1,:)),0.0,1.0,1.0,0.0,-g_X1(1,:).*g_par(8,:).*(g_X1(1,:).*g_par(2,:).^(-g_par(6,:)).*(g_X1(2,:).^g_par(6,:).*log(g_X1(2,:))+g_par(2,:).^g_par(6,:).*log(g_par(2,:)))-g_X1(1,:).*g_par(2,:).^(-g_par(6,:)).*log(g_par(2,:)).*(g_X1(2,:).^g_par(6,:)+g_par(2,:).^g_par(6,:))),0.0,0.0,-g_X1(2,:).*((g_X1(2,:).*(g_X1(1,:)+g_par(3,:)))./(g_par(1,:).*g_par(3,:))-1.0),-g_X1(1,:).*(g_X1(1,:).*g_par(2,:).^(-g_par(6,:)).*(g_X1(2,:).^g_par(6,:)+g_par(2,:).^g_par(6,:))-1.0),0.0),[2,8,max(size(g_X1,2),size(g_par,2))]);
g_l_res.p0=[E0;hE;hV;inpE;inpV;p;rE;rV];

g_l_res.parnames={'E0','hE','hV','inpE','inpV','p','rE','rV'};

%function created by GRIND
function g_X2=curr_ode0(t,g_X1,E0,hE,hV,inpE,inpV,p,rE,rV)
Eeq=E0.*hV./(hV+g_X1(1,:));
Veq=hE.^p./(hE.^p+g_X1(2,:).^p);
g_X2(1,:)=rV.*g_X1(1,:).*(1-g_X1(1,:)./Veq)+inpV;
g_X2(2,:)=rE.*g_X1(2,:).*(1-g_X1(2,:)./Eeq)+inpE;

function m=i_vcat(varargin)
len=cellfun('size',varargin,2);
len1=len==1;
if any(len1)&&any(~len1)
    flen1=find(len1);
    maxlen=max(len);
    for i=1:length(flen1)
        varargin{flen1(i)}=varargin{flen1(i)}(ones(1,maxlen));
    end;
end;
m=vertcat(varargin{:});
