%-------------------------------------------------------------------------%
%         This software is licensed by Creative Commons BY-NC-SA:         %
%      http://creativecommons.org/licenses/by-nc-sa/3.0/it/legalcode      %
%-------------------------------------------------------------------------%
%
% The file is an example concerning the determination of the basins of 
% attraction of two equilibrium points
%
%-------------------------------------------------------------------------%
clear all
close all
clc 
ngrid = 12; % Number of equispaced points on the edges of the cube
tol = 0.0001; % Tolerance used to stop bisection
ed = 4; % Edge of the cube 
ep = 0.1; % Shape parameter
npu = 4; % Number or PU subdomains
npi = 2; % Number of equilibria partitioning the phase state (tipically
% origin and coexistence equilibrium point)
% Define the system parameters
p = 1; u = 1; a = 1; b = 2; q = 0.5; v = 0.2; c = 0.3; e = 1; r = 2; w = 1;
f = 3; g = 2;
par = [p,u,a,b,q,v,c,e,r,w,f,g]; % Model parameters
% Define the matrix of equilibria
% 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: (repeat) a stable equilibrium point
%    from row npi+4 to row M: remaining equilibria
E = [0 0 par(10); ...
    0 0 0;...
    -par(2)*(par(1)*par(9)*par(5)-par(12)*par(10)*par(1)*par(6)*par(8)-...
    par(3)*par(6)*par(9)*par(5)+par(3)*par(6)*par(9)*par(10)*par(8)+...
    par(4)*par(10)*par(5)*par(12)*par(6)-par(4)*par(10)*par(9)*par(5))...
    /(-par(1)*par(9)*par(5)+par(9)*par(7)*par(2)*par(6)*par(3)+par(11)...
    *par(2)*par(10)*par(4)*par(5)-par(11)*par(2)*par(10)*par(3)*par(6)...
    *par(8)-par(12)*par(10)*par(6)*par(7)*par(2)*par(4)+par(12)*par(10)...
    *par(1)*par(6)*par(8)) -par(6)*(par(1)*par(9)*par(5)-par(1)*par(9)...
    *par(7)*par(2)-par(11)*par(2)*par(10)*par(4)*par(5)+par(9)*par(10)...
    *par(7)*par(2)*par(4)-par(9)*par(10)*par(1)*par(8)+par(11)*par(2)...
    *par(10)*par(1)*par(8))/(-par(1)*par(9)*par(5)+par(9)*par(7)*par(2)...
    *par(6)*par(3)+par(11)*par(2)*par(10)*par(4)*par(5)-par(11)*par(2)...
    *par(10)*par(3)*par(6)*par(8)-par(12)*par(10)*par(6)*par(7)*par(2)*...
    par(4)+par(12)*par(10)*par(1)*par(6)*par(8))  (-par(7)*par(2)*par(1)...
    *par(12)*par(6)+par(9)*par(7)*par(2)*par(6)*par(3)+par(2)*par(1)*...
    par(11)*par(5)-par(2)*par(6)*par(3)*par(11)*par(5)+par(1)*par(5)*...
    par(12)*par(6)-par(1)*par(9)*par(5))*par(10)/(-par(1)*par(9)*par(5)+...
    par(9)*par(7)*par(2)*par(6)*par(3)+par(11)*par(2)*par(10)*par(4)*...
    par(5)-par(11)*par(2)*par(10)*par(3)*par(6)*par(8)-par(12)*par(10)*...
    par(6)*par(7)*par(2)*par(4)+par(12)*par(10)*par(1)*par(6)*par(8));...
    par(2)*par(5)*(-par(1)+par(3)*par(6))/(-par(1)*par(5)+par(7)*par(2)...
    *par(6)*par(3)) par(1)*par(6)*(-par(5)+par(7)*par(2))/(-par(1)*...
    par(5)+par(7)*par(2)*par(6)*par(3)) 0;...
    par(2)*par(5)*(-par(1)+par(3)*par(6))/(-par(1)*par(5)+par(7)*...
    par(2)*par(6)*par(3)) par(1)*par(6)*(-par(5)+par(7)*par(2))/...
    (-par(1)*par(5)+par(7)*par(2)*par(6)*par(3)) 0;...
    par(2)*par(9)*(-par(1)+par(4)*par(10))/(-par(1)*par(9)+par(11)*...
    par(2)*par(10)*par(4)) 0 par(10)*par(1)*(-par(9)+par(11)*par(2))...
    /(-par(1)*par(9)+par(11)*par(2)*par(10)*par(4));...
    0  par(6)*par(9)*(-par(5)+par(10)*par(8))/(-par(5)*par(9)+par(12)...
    *par(6)*par(10)*par(8)) par(10)*par(5)*(-par(9)+par(12)*par(6))...
    /(-par(5)*par(9)+par(12)*par(6)*par(10)*par(8));...
    par(2) 0 0;...
    0 par(6) 0];
tmax = 90; % Integration time
k = 8; % Used to estimate normals
nsp = 2; % Number of stable equilibrium points
% Reconstruct the attraction basins
[A1,A2,A3] = SeparatrixSurfaces(ngrid,tol,ed,ep,npu,npi,par,E,tmax,k,nsp);