%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% File: SeparatrixSurfaces(ngrid,tol,ed,ep,npu,npi,par,...
%       E,tmax,k,nsp)
%
% Goal: script that recontructs the basins of attraction in 3D dynamical
%       systems with two or three stable equilibira
%
% Inputs: ngrid: number of grid points in one direction of the cube
%         tol:   tolerance used during the bisection algorithm
%         ed:    edge of the cube
%         ep:    shape parameter for RBF. Determines size of support of
%                basis function. 
%         npu:   number of PU cells in one direction
%         npi:   number of equilibrium points to be interpolated
%                (tipically origin and coexistence saddle point)
%         par:   vector of model parameters
%         E:     matrix of the equilibria Mx3, where M is the number of
%                equilibria. The matrix is organized as follows:
%                row 1: a stable equilibrium point
%                from row 2 to row npi+1: equilibrium points to be 
%                interpolated
%                row npi+2: a stable equilibrium point
%                row npi+3: a stable equilibrium point
%                from row npi+4 to row M: remaining equilibria
%        tmax:   integration time
%        k:      number of nearest points used for estimate the normals
%        nsp:    number of stable equilibria
%
% Outputs: A1 A2 and A3: matrices of separatrix points and plot of the
%                        attraction basins
%
% Calls on: SeparatrixPoints.m, PointCloud3D_PUCS.m
%
% Remark:   if the integration time is not sufficient or if there are
%           not two or three stable equilibria the algorithm returns
%           an error message
%
%-------------------------------------------------------------------------%
function [A1,A2,A3]=SeparatrixSurfaces(ngrid,tol,ed,ep,npu,npi,par,...
    E,tmax,k,nsp)
% Use symbolic calculation for computing the system Jacobian
syms x y z
J=Jacobian(par);
size_E=size(E,1);
lambda=ones(3,size_E); % Initialize
% For all equilibria find the eigenvalues
for i=1:size_E
    lambda(:,i)=eig(subs(J, [x y z], [E(i,1) E(i,2) E(i,3)]));
end
% Check if the equilibrium points stored in the matrix E satisfy 
% the stability and feasibility conditions
if (E(1,1)<0) || (E(1,2)<0) || (E(1,3)<0)
    disp('First equilibrium not feasible')
    return
end
if (lambda(1,1)>0) || (lambda(2,1)>0) || (lambda(3,1)>0)
    disp('First equilibrium not stable')
    return
end
if (E(npi+2,1)<0) || (E(npi+2,2)<0) || (E(npi+2,3)<0)
    disp('Second equilibrium not feasible')
    return
end
if (lambda(1,npi+2)>0) || (lambda(2,npi+2)>0) || (lambda(3,npi+2)>0)
    disp('Second equilibrium not stable')
    return
end
if (E(npi+3,1)<0) || (E(npi+3,2)<0) || (E(npi+3,3)<0)
    disp('Third equilibrium not feasible')
    return
end
if (lambda(1,npi+3)>0) || (lambda(2,npi+3)>0) || (lambda(3,npi+3)>0)
    disp('Third equilibrium not stable')
    return
end
% Check if the npi equilibrium points to be interpolated are feasible
% but not stable
if npi~=0
    for i=2:npi+1
        if (E(i,1)<0) || (E(i,2)<0) || (E(i,3)<0)
            disp('Point to be interpolated not feasible')
            return
        end
        if (lambda(1,i)<0) && (lambda(2,i)<0) && (lambda(3,i)<0)
            disp('Point to be interpolated stable')
            return
        end
    end
end
% Check in the remaining equilibria if there exist other feasible and
% stable equilibrium points
if (npi+3<size_E)
    for i=npi+4:size_E
        if ((E(i,1)>0) && (E(i,2)>0) && (E(i,3)>0)) && ...
                ((lambda(1,i)<0) && (lambda(2,i)<0) && (lambda(3,i)<0))
            disp('Exist an other stable equilibrium')
            return
        end
    end
end
% Declare equispaced points on the edge of the cube
x_grid = linspace(0,ed,ngrid); y_grid = linspace(0,ed,ngrid);
z_grid = linspace(0,ed,ngrid);
npts1=0; npts2=0; npts3=0; % Initialize
% Calculate separatrix points
for j=1:ngrid
    for i=1:ngrid
        % Declare initial conditions used to find a separatrix point
        p1=[x_grid(i) x_grid(i); 0 ed; z_grid(j) z_grid(j)];
        % Find a separatrix point
        [point1 point2 point3] = SeparatrixPoints(p1,tmax,tol,...
            npi,par,E,nsp);
        % If the two initial conditions evolve towards the same stable
        % equilibrium point then bisection cannot be performed
        if ((point1(1)==0) && (point1(2)==0) && (point1(3)==0)) && ...
                ((point2(1)==0) && (point2(2)==0) && (point2(3)==0)) ...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            disp('Bisection cannot be performed')
        % If the integration time is not sufficient the algorithm 
        % returns an error message
        elseif ((point1(1)==-1) && (point1(2)==-1) && (point1(3)==-1))...
                && ((point2(1)==-1) && (point2(2)==-1) && ...
                (point2(3)==-1)) && ((point3(1)==-1) && ...
                (point3(2)==-1) && (point3(3)==-1))
            disp('Integration is time not sufficient')
            return
            % Store the separatrix point into one of the three matrices
        elseif ((point2(1)==0) && (point2(2)==0) && (point2(3)==0))...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            npts1 = npts1+1;
            A1(npts1,1) = point1(1);
            A1(npts1,2) = point1(2);
            A1(npts1,3) = point1(3);
        elseif ((point1(1)==0) && (point1(2)==0) && (point1(3)==0))...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            npts2 = npts2+1;
            A2(npts2,1) = point2(1);
            A2(npts2,2) = point2(2);
            A2(npts2,3) = point2(3);
        else
            npts3 = npts3+1;
            A3(npts3,1) = point3(1);
            A3(npts3,2) = point3(2);
            A3(npts3,3) = point3(3);
        end
    end
end
% Use the second set of initial conditions to find separatrix points
for j=1:ngrid
    for i=1:ngrid
        % Declare initial conditions used to find a separatrix point
        p1=[0 ed; y_grid(i) y_grid(i); z_grid(j) z_grid(j)];
        % Find a separatrix point
        [point1 point2 point3] = SeparatrixPoints(p1,tmax,tol,...
            npi,par,E,nsp);
        % If the two initial conditions evolve towards the same stable
        % equilibrium point then bisection cannot be performed
        if ((point1(1)==0) && (point1(2)==0) && (point1(3)==0)) && ...
                ((point2(1)==0) && (point2(2)==0) && (point2(3)==0)) ...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            disp('Bisection cannot be performed')
        % If the integration time is not sufficient the algorithm 
        % returns an error message
        elseif ((point1(1)==-1) && (point1(2)==-1) && (point1(3)==-1))...
                && ((point2(1)==-1) && (point2(2)==-1) && ...
                (point2(3)==-1)) && ((point3(1)==-1) && ...
                (point3(2)==-1)  && (point3(3)==-1))
            disp('Integration time is not sufficient')
            return
            % Store the separatrix point into one of the three matrices
        elseif ((point2(1)==0) && (point2(2)==0) && (point2(3)==0)) ...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            npts1 = npts1+1;
            A1(npts1,1) = point1(1);
            A1(npts1,2) = point1(2);
            A1(npts1,3) = point1(3);
        elseif ((point1(1)==0) && (point1(2)==0) && (point1(3)==0)) ...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            npts2 = npts2+1;
            A2(npts2,1) = point2(1);
            A2(npts2,2) = point2(2);
            A2(npts2,3) = point2(3);
        else
            npts3 = npts3+1;
            A3(npts3,1) = point3(1);
            A3(npts3,2) = point3(2);
            A3(npts3,3) = point3(3);
        end
    end
end
% Use the third set of initial conditions to find separatrix points
for j=1:ngrid
    for i=1:ngrid
        % Declare initial conditions used to find a separatrix point
        p1=[x_grid(i) x_grid(i); y_grid(j) y_grid(j); 0 ed];
        % Find a separatrix point
        [point1 point2 point3] = SeparatrixPoints(p1,tmax,tol,...
            npi,par,E,nsp);
        % If the two initial conditions evolve towards the same stable
        % equilibrium point then bisection cannot be performed
        if ((point1(1)==0) && (point1(2)==0) && (point1(3)==0)) && ...
                ((point2(1)==0) && (point2(2)==0) && (point2(3)==0)) && ...
                ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            disp('Bisection cannot be performed')
            % Check if the integration time is sufficient
        elseif ((point1(1)==-1) && (point1(2)==-1) && (point1(3)==-1))...
                && ((point2(1)==-1) && (point2(2)==-1) && ...
                (point2(3)==-1)) && ((point3(1)==-1) && (point3(2)==-1) ...
                && (point3(3)==-1))
            disp('Integration time is not sufficient')
            return
        % If the integration time is not sufficient the algorithm 
        % returns an error message
        elseif ((point2(1)==0) && (point2(2)==0) && (point2(3)==0)) ...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            npts1 = npts1+1;
            A1(npts1,1) = point1(1);
            A1(npts1,2) = point1(2);
            A1(npts1,3) = point1(3);
        elseif ((point1(1)==0) && (point1(2)==0) && (point1(3)==0)) ...
                && ((point3(1)==0) && (point3(2)==0) && (point3(3)==0))
            npts2 = npts2+1;
            A2(npts2,1) = point2(1);
            A2(npts2,2) = point2(2);
            A2(npts2,3) = point2(3);
        else
            npts3 = npts3+1;
            A3(npts3,1) = point3(1);
            A3(npts3,2) = point3(2);
            A3(npts3,3) = point3(3);
        end
    end
end
% If the number of stable equilibria is two then the algorithm returns
% only A1
if nsp == 2
    A2 = 1; A3 = 1;
    % Add the npi equilibrium points to be interpolated
    if npi~=0
        s = 1; B = zeros(npi,3); i = 2;
        while(i<=npi+1)
            B(s,1) = E(i,1); B(s,2) = E(i,2); B(s,3) = E(i,3);
            s = s+1; i = i+1;
        end
        A1 = [B; A1];
    end
    % Plot separatrix points
    figure
    hold on
    plot3(A1(:,1),A1(:,2),A1(:,3),'m*')
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis([0 ed 0 ed 0 ed])
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
    % Interpolation with inmplicit PUM
    figure
    hold on
    PointCloud3D_PUCS(A1,ep,npu,k)
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis([0 ed 0 ed 0 ed])
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
else
    % Define the matrices of separatrix points describing the three basins
    % of attraction
    A1 = sortrows(A1,1); A2 = sortrows(A2,1); A3 = sortrows(A3,1);
    A4 = [A1; A2]; A5 = [A1; A3]; A6 = [A3; A2];
    % Add the npi equilibrium points to be interpolated
    if npi~=0
        s = 1; B = zeros(npi,3); i = 2;
        while(i<=npi+1)
            B(s,1) = E(i,1); B(s,2) = E(i,2); B(s,3) = E(i,3);
            s = s+1; i = i+1;
        end
        A4 = [B; A4]; A5 = [B; A5]; A6 = [B; A6]; 
    end
    figure
    hold on
    plot3(A1(:,1),A1(:,2),A1(:,3),'m*')
    plot3(A2(:,1),A2(:,2),A2(:,3),'b*')
    plot3(A3(:,1),A3(:,2),A3(:,3),'r*')
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis([0 ed 0 ed 0 ed])
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
    plot3(E(npi+3,1),E(npi+3,2),E(npi+3,3),'.b','MarkerSize',20)
    % Interpolation with inmplicit PUM
    figure
    hold on
    PointCloud3D_PUCS(A4,ep(1),npu(1),k(1))
    axis([0 ed 0 ed 0 ed])
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
    plot3(E(npi+3,1),E(npi+3,2),E(npi+3,3),'.b','MarkerSize',20)
    figure
    hold on
    PointCloud3D_PUCS(A5,ep(2),npu(2),k(2))
    axis([0 ed 0 ed 0 ed])
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
    plot3(E(npi+3,1),E(npi+3,2),E(npi+3,3),'.b','MarkerSize',20)    
    figure
    hold on
    PointCloud3D_PUCS(A6,ep(3),npu(3),k(3))
    axis([0 ed 0 ed 0 ed])
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
    plot3(E(npi+3,1),E(npi+3,2),E(npi+3,3),'.b','MarkerSize',20)
    % Plot the three basins of attraction
    figure
    hold on
    PointCloud3D_PUCS(A4,ep(1),npu(1),k(1))
    PointCloud3D_PUCS(A5,ep(2),npu(2),k(2))
    PointCloud3D_PUCS(A6,ep(3),npu(3),k(3))
    axis([0 ed 0 ed 0 ed])
    xlabel('x','fontsize',13)
    ylabel('y','fontsize',13)
    zlabel('z','fontsize',13)
    axis square
    grid on
    plot3(E(1,1),E(1,2),E(1,3),'.b','MarkerSize',20)
    i = 2;
    if npi~=0
        while(i<=npi+1)
            plot3(E(i,1),E(i,2),E(i,3),'.k','MarkerSize',20)
            i = i+1;
        end
    end
    plot3(E(npi+2,1),E(npi+2,2),E(npi+2,3),'.b','MarkerSize',20)
    plot3(E(npi+3,1),E(npi+3,2),E(npi+3,3),'.b','MarkerSize',20)
end
end