% Software for the computation of multi-component factorizations and the set of feasible solutions 
% using analytical approach and inner boundary curve.
%
% Copyright (c) 2021 Tomass Andersons, University of Rostock, Department of Mathematics.
% 
% The software can be used for academic, research and other similar noncommercial uses. The user
% acknowledges that this software is still in the development stage and that it is provided by the
% copyright holders and contributors "as is" and any express or implied warranties, fitness for a
% particular purpose are disclaimed. In no event shall the copyright owner or contributors be liable
% for any direct, indirect, incidental, special, exemplary, or consequential damages.
% 
% The copyright holders provide no reassurances that the source code provided does not infringe any
% patent, copyright, or any other intellectual property rights of third parties. The copyright
% holders disclaim any liability to any recipient for claims brought against recipient by any third
% party for infringement of that parties intellectual property rights.

%% WARNING!
% The algorithm in its current form is not intended for noisy data or
% cases when INNPOL is very close to FIRPOL. The results might be 
% incorrect. Use polygon inflation algorithm in these cases.

%% 
% This GUI accompanies:
% T. Andersons, M. Sawall, K. Neymeyr,
% Analytical enclosure of the set of solutions of the three-species multivariate curve resolution problem
% J. Math. Chem. 60, 1750-1780 (2022).
% DOI: 10.1002/cem.3381

%%
% This is the calculation file of the inner boundary curve algorithm
% This file is required for the GUI in GUI4IBCA.m
% It can also be used for standalone calculations without the GUI.

% After an instance of the class IBCA has been declared (e.g. g = IBCA)
% and the data matrix has been loaded (g.D = D), 
% the following functions can be used:
% g.calcSVD
% g.calcTSVD
% g.calcBoundaries
% g.calcIBCF
% g.calcnIBC

% There are also two different functions for the calculation of the SFS:
% g.calcSFSxor 
% for the simple calculation of SFS (sometimes incorrect) and
% g.calcIntersection
% g.calcSFS
% for the analytical calculation

% Further plotting functions are available
% g.plotBoundary(1)
% g.plotBoundary(2)
% g.plotIBC
% g.plotSFS

% If you wish to load an existing result in to the IBCA object, then use 
% g.loadRes
% (it is the responsibility of the user to load a valid IBCA result file, 
% there are no checks for that)

% It is also possible to return the coordinates of the boundary
% triangle for a specific rotaton with the following function
% g.alpha = 1; % This is the rotation
% g.calcTriangle % This carries out the calculation
% The results are saved in the object itself:
% g.pTriangle % Coordinates of the primal boundary triangle
% g.pTangents % Coordinates of the tangential points of the primal boundary triangle
% g.dTriangle % Coordinates of the dual boundary triangle
% g.dTangents % Coordinates of the tangential points of the dual boundary triangle
% g.iTr % Number of the subfunction for this rotation angle
% g.pTriangle(3,:) % Coordinates of the free vertex on the IBC for the given rotation

%%
% WARNING!
% The algorithm does not function well, if INNPOL is very close to FIRPOL.
% The fast calculation of SFS can sometimes be incorrect. If this is
% suspected, use analytical SFS.
% Results should be validated with the polygon inflation algorithm
% This file is used with Matlab R2019B

%%
classdef IBCA < handle
    
    %% TODO
    % Plot SFS for a line
    % Break if incorrect data
    % Plots unused in GUI for debug
    %% Properties
    properties
        D       % Nonnegative Matrix
        T       % Time vector (optional) containing the time evolution beginning at reaction start
        X       % Channel vector (optional) containing the wavelengths/numbers of the spectra
        U       % Singular value decomposition
        S       % Singular value decomposition
        V       % Singular value decomposition
        F1      % Primal FIRPOL
        kF1     % Primal FIRPOL coefficients
        I1      % Primal INNPOL
        kI1     % Primal INNPOL coefficients
        F2      % Dual FIRPOL
        kF2     % Dual FIRPOL coefficients
        I2      % Dual INNPOL
        kI2     % Dual INNPOL coefficients
        TV      % Transition values and parameters
        IBCF    % Inner boundary curve function coefficients
        nIBC    % Approximation of inner boundary curve
        isec    % Intersection between IBC and FIRPOL
        SFS     % Segments of the set of feasible solutions saved in cells
        
        alpha   % Used to calculate a particular boundary triangle with this angle
        pTriangle   % Primary boundary triangle
        pTangents   % Tangential points of the primary triangle
        dTriangle   % Dual boundary triangle
        dTangents   % Tangential points of the dual triangle
        iTr         % Index of TV and IBCA for the current alpha
    end
    
    properties(GetAccess=private)
        para
        debug
    end
    
    %% Ordinary methods
    methods
        %% Constructor
        function obj = IBCA(para)
            if nargin==0
                para.debug = 0;           % Debug mode
                para.maxSVD = 4;          % Number of singular values
                para.INNFIRtol = 1E-12;   % calcBoundaries INNPOL point outside of FIRPOL tolerance
                para.isIsectol = 1E-12;   % isIntersection tolerance
                para.calcIsectol = 1E-12; % calcIntersection tolerance
                para.rotAngle = 0.005;    % calcIBCAprox distAngle (rotation parameter)
                para.overlaptol = 1E-14;  % calcIBCAprox overlapTol (parameter for overlapping points)
                obj.para = para;
            elseif nargin==1
                obj.para=para;
            end
        end
        %% Auxilarry for para
        function setpara(obj,para)
            obj.para=para;
        end
        function para=getpara(obj)
            para=obj.para;
        end
        %% calcSVD: Calculate SVD
        function calcSVD(obj)
            % In: D
            % Out: U, S, V
            [obj.U,obj.S,obj.V]=svds(obj.D,obj.para.maxSVD);
        end
        %% calcTSVD: Calculate the truncated SVD with correct orientation
        function calcTSVD(obj)
            % In: U, S, V
            % Out: U, S, V
            p = zeros(3);
            for i = 1:3
                if max(obj.V(:,i))>=max(-obj.V(:,i))
                    p(i,i) = 1;
                else
                    p(i,i) = -1;
                end
            end
            obj.U = obj.U(:,1:3)*p;
            obj.S = obj.S(1:3,1:3);
            obj.V = obj.V(:,1:3)*p;
        end
        %% calcBoundaries: Calculate INNPOL and FIRPOL for both factors
        function calcBoundaries(obj)
            % In: D,U,S,V
            % Out: F1,kF1,I1,kI1,F2,kF2,I2,kI2
            sizeD = [size(obj.U,1) size(obj.V,1)];
            US = obj.U*obj.S;
            
            % Primal INNPOL
            O1 = zeros(sizeD(1));
            for i=1:sizeD(1)
                O1(i,i) = 1/US(i,1);
            end
            obj.I1 = O1*obj.D*obj.V(:,2:3);
            obj.kI1 = convhull(obj.I1);
            
            % Dual FIRPOL
            sizeTemp = size(obj.kI1,1)-1;
            obj.F2 = zeros(sizeTemp,2);
            obj.F2(:,1) = (-obj.I1(obj.kI1(1:sizeTemp),2)+obj.I1(obj.kI1(2:sizeTemp+1),2))./...
                (obj.I1(obj.kI1(1:sizeTemp),2).*obj.I1(obj.kI1(2:sizeTemp+1),1)-obj.I1(obj.kI1(2:sizeTemp+1),2).*obj.I1(obj.kI1(1:sizeTemp),1));
            obj.F2(:,2) = (-obj.I1(obj.kI1(1:sizeTemp),1)+obj.I1(obj.kI1(2:sizeTemp+1),1))./...
                (obj.I1(obj.kI1(1:sizeTemp),1).*obj.I1(obj.kI1(2:sizeTemp+1),2)-obj.I1(obj.kI1(2:sizeTemp+1),1).*obj.I1(obj.kI1(1:sizeTemp),2));
            obj.kF2 = convhull(obj.F2);
            
            % Dual INNPOL
            O2 = zeros(sizeD(2));
            for i=1:sizeD(2)
                O2(i,i) = 1/obj.V(i,1);
            end
            obj.I2 = O2*obj.V(:,2:3);
            obj.kI2 = convhull(obj.I2);
            
            % Primal FIRPOL
            sizeTemp = size(obj.kI2,1)-1;
            obj.F1 = zeros(sizeTemp,2);
            obj.F1(:,1) = (-obj.I2(obj.kI2(1:sizeTemp),2)+obj.I2(obj.kI2(2:sizeTemp+1),2))./...
                (obj.I2(obj.kI2(1:sizeTemp),2).*obj.I2(obj.kI2(2:sizeTemp+1),1)-obj.I2(obj.kI2(2:sizeTemp+1),2).*obj.I2(obj.kI2(1:sizeTemp),1));
            obj.F1(:,2) = (-obj.I2(obj.kI2(1:sizeTemp),1)+obj.I2(obj.kI2(2:sizeTemp+1),1))./...
                (obj.I2(obj.kI2(1:sizeTemp),1).*obj.I2(obj.kI2(2:sizeTemp+1),2)-obj.I2(obj.kI2(2:sizeTemp+1),1).*obj.I2(obj.kI2(1:sizeTemp),2));
            obj.kF1 = convhull(obj.F1);
            
            % Check if INNPOL in FIRPOL and move near points
            in = inpolygon(obj.I1(obj.kI1(1:end-1),1),obj.I1(obj.kI1(1:end-1),2),obj.F1(obj.kF1,1),obj.F1(obj.kF1,2));
            changeF2 = 0;
            for i=1:size(in,1)
                if in(i)==0
                    for j=1:size(obj.kF1)-1
                        % Find the distance of the linear projection of the point of INNPOL
                        % on each edge of FIRPOL (extension of each edge)
                        l=dot((obj.F1(obj.kF1(j+1),:)-obj.F1(obj.kF1(j),:)),...
                            (obj.I1(obj.kI1(i),:)-obj.F1(obj.kF1(j),:)))/(norm(obj.F1(obj.kF1(j+1),:)-obj.F1(obj.kF1(j),:)))^2;
                        % Find the distance to line segment
                        I1Dist = ones(size(obj.kF1,1)-1,1);
                        if l<0
                            I1Dist(j) = norm(obj.I1(obj.kI1(i),:)-obj.F1(obj.kF1(j),:));
                        elseif l>1
                            I1Dist(j) = norm(obj.I1(obj.kI1(i),:)-obj.F1(obj.kF1(j+1),:));
                        else
                            I1Dist(j) = norm((obj.I1(obj.kI1(i),:)-obj.F1(obj.kF1(j),:)) - l.*(obj.F1(obj.kF1(j+1),:)-obj.F1(obj.kF1(j),:)));
                        end
                        I1MinDist = min(I1Dist);
                        if I1MinDist < obj.para.INNFIRtol
                            % Primal INNPOL is changed. Therefore, dual FIRPOL must be changed too.
                            changeF2 = 1;
                            % Move the INNPOL point in the direction of 0
                            obj.I1(obj.kI1(i),:) = obj.I1(obj.kI1(i),:)-sign(obj.I1(obj.kI1(i),:)).*[1E-12 1E-12];
                        else
                            % TODO Break!!!
                        end
                    end
                end
            end
            if changeF2 == 1
                % Dual FIRPOL
                sizeTemp = size(obj.kI1,1)-1;
                obj.F2 = zeros(sizeTemp,2);
                obj.F2(:,2) = (-obj.I1(obj.kI1(1:sizeTemp),1)+obj.I1(obj.kI1(2:sizeTemp+1),1))./...
                    (obj.I1(obj.kI1(1:sizeTemp),1).*obj.I1(obj.kI1(2:sizeTemp+1),2)-...
                    obj.I1(obj.kI1(2:sizeTemp+1),1).*obj.I1(obj.kI1(1:sizeTemp),2));
                obj.F2(:,1) = (-1-obj.I1(obj.kI1(1:sizeTemp),2).*obj.F2(:,2))./obj.I1(obj.kI1(1:sizeTemp),1);
                obj.kF2 = convhull(obj.F2);
            end
            
            % Debug output
            if obj.para.debug
                figure(80)
                subplot(1,2,1)
                hold on
                plot(obj.F1(obj.kF1,1),obj.F1(obj.kF1,2)) % Primal FIRPOL
                plot(obj.I1(obj.kI1,1),obj.I1(obj.kI1,2)) % Primal INNPOL
                plot(obj.I1(:,1),obj.I1(:,2),'o') % Primal INNPOL
                
                subplot(1,2,2)
                hold on
                plot(obj.F2(obj.kF2,1),obj.F2(obj.kF2,2)) % Dual FIRPOL
                plot(obj.I2(obj.kI2,1),obj.I2(obj.kI2,2)) % Dual INNPOL
                plot(obj.I2(:,1),obj.I2(:,2),'o') % Dual INNPOL
            end
        end
        
        %% calcIBCF: Find the transition values and the inner boundary function
        function calcIBCF(obj)
            % In: F1,kF1,I1,kI1,F2,kF2,I2,kI2
            % Out: TV, IBCF
            sizeI = size(obj.kF2,1)-1;
            sizeF = size(obj.kF1,1)-1;
            % Transition values and parameters
            obj.TV = zeros(3*sizeI+2*sizeF,12);
            
            % Calculate transition values of the first kind
            obj.TV(1:sizeI,:) = [atan2(obj.F2(obj.kF2(1:sizeI),2),obj.F2(obj.kF2(1:sizeI),1)) ...
                obj.I1(obj.kI1(2:sizeI+1),1) obj.I1(obj.kI1(2:sizeI+1),2) NaN(sizeI,1) NaN(sizeI,1) ...
                NaN(sizeI,1) NaN(sizeI,1) NaN(sizeI,1) NaN(sizeI,1) NaN(sizeI,1) NaN(sizeI,1) ones(sizeI,1)];
            
            % Calculate transition values of the second kind
            for i=1:sizeF
                % Find the convex hull of primal INNPOL and one vertex of primal FIRPOL
                % to find tangents of primal INNPOL.
                [z1l,z2l,z1r,z2r] = IBCA.calcTCH(obj.I1,obj.F1,obj.kF1,i);
                % Left tangent -> for right transition value
                % Right tangent -> for left transition value
                
                % Parameter for edge
                [a1,a2] = IBCA.calcLEQ(obj.F1(obj.kF1(i),:),obj.F1(obj.kF1(i+1),:));
                % Add into the parameter list together with the transition value from
                % tangential points
                obj.TV(sizeI+(2*i),:) = [atan2(z2r,z1r) NaN NaN a1 a2 NaN NaN NaN NaN NaN NaN 2]; % u_l
                obj.TV(sizeI+(2*i-1),:) = [atan2(z2l,z1l) NaN NaN NaN NaN a1 a2 NaN NaN NaN NaN 3]; % u_r
            end
            
            % Calculate transition values of the third kind
            for i=1:sizeI
                % Find the convex hull of dual INNPOL and one vertex of dual FIRPOL
                % to find tangents of dual INNPOL.
                [z1l,z2l,z1r,z2r] = IBCA.calcTCH(obj.I2,obj.F2,obj.kF2,i);
                % Left tangent -> for right transition value
                % Right tangent -> for left transition value
                
                % Find the convex hull of primal INNPOL and dual point of the right tangent
                % to find tangents of primal INNPOL and left transition values.
                [~,~,q1l,q2l] = IBCA.calcTCH(obj.I1,[z1r z2r],[1],1);
                
                % Find the convex hull of primal INNPOL and dual point of the left tangent
                % to find tangents of primal INNPOL and right transition values.
                [q1r,q2r,~,~] = IBCA.calcTCH(obj.I1,[z1l z2l],[1],1);
                
                % Parameter for edge
                [a1,a2] = IBCA.calcLEQ(obj.F2(obj.kF2(i),:),obj.F2(obj.kF2(i+1),:));
                % Add into the parameter list together with the transition value
                obj.TV(sizeI+2*sizeF+(2*i-1),:) = [atan2(q2l, q1l) NaN NaN NaN NaN NaN NaN a1 a2 NaN NaN 4]; % o_l
                obj.TV(sizeI+2*sizeF+(2*i),:) = [atan2(q2r, q1r) NaN NaN NaN NaN NaN NaN NaN NaN a1 a2 5]; % o_r
            end
            
            % Sorting transition values and filling the parameter table
            obj.TV = sortrows(obj.TV);
            % Filling in the NaN values
            for i = 2:size(obj.TV,1)
                for j = 1:size(obj.TV,2)
                    if isnan(obj.TV(i,j))
                        obj.TV(i,j)=obj.TV(i-1,j);
                    end
                end
            end
            for j = 1:size(obj.TV,2)
                if isnan(obj.TV(1,j))
                    obj.TV(1,j)=obj.TV(size(obj.TV,1),j);
                end
            end
            for i = 2:size(obj.TV,1)
                for j = 1:size(obj.TV,2)
                    if isnan(obj.TV(i,j))
                        obj.TV(i,j)=obj.TV(i-1,j);
                    end
                end
            end
            
            % Calculating subfunction
            obj.IBCF = zeros(3*sizeI+2*sizeF,10);
            for i = 1:size(obj.IBCF,1)
                obj.IBCF(i,1) = obj.TV(i,1);
                obj.IBCF(i,2:10) = IBCA.calcExpression(obj.TV(i,2:11));
            end
        end
        
        %% calcnIBC: Calculate the approximation of IBC with equidistant angles
        function calcnIBC(obj)
            % In: IBCF
            % Out: nIBC
            alp=-pi:obj.para.rotAngle:pi;
            obj.nIBC=zeros(numel(alp),2);
            for i=1:numel(alp)
                for j=1:size(obj.IBCF,1)-1
                    k = size(obj.IBCF,1);
                    if alp(i)<=obj.IBCF(1,1)
                        obj.nIBC(i,1) = (obj.IBCF(k,2)+obj.IBCF(k,3)*tan(alp(i))+obj.IBCF(k,4)*(tan(alp(i)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(alp(i))+obj.IBCF(k,10)*(tan(alp(i)))^2);
                        obj.nIBC(i,2) = (obj.IBCF(k,5)+obj.IBCF(k,6)*tan(alp(i))+obj.IBCF(k,7)*(tan(alp(i)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(alp(i))+obj.IBCF(k,10)*(tan(alp(i)))^2);
                    elseif alp(i)<=obj.IBCF(j+1,1) && alp(i)>obj.IBCF(j,1)
                        obj.nIBC(i,1) = (obj.IBCF(j,2)+obj.IBCF(j,3)*tan(alp(i))+obj.IBCF(j,4)*(tan(alp(i)))^2)/...
                            (obj.IBCF(j,8)+obj.IBCF(j,9)*tan(alp(i))+obj.IBCF(j,10)*(tan(alp(i)))^2);
                        obj.nIBC(i,2) = (obj.IBCF(j,5)+obj.IBCF(j,6)*tan(alp(i))+obj.IBCF(j,7)*(tan(alp(i)))^2)/...
                            (obj.IBCF(j,8)+obj.IBCF(j,9)*tan(alp(i))+obj.IBCF(j,10)*(tan(alp(i)))^2);
                    elseif alp(i)>obj.IBCF(k,1)
                        obj.nIBC(i,1) = (obj.IBCF(k,2)+obj.IBCF(k,3)*tan(alp(i))+obj.IBCF(k,4)*(tan(alp(i)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(alp(i))+obj.IBCF(k,10)*(tan(alp(i)))^2);
                        obj.nIBC(i,2) = (obj.IBCF(k,5)+obj.IBCF(k,6)*tan(alp(i))+obj.IBCF(k,7)*(tan(alp(i)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(alp(i))+obj.IBCF(k,10)*(tan(alp(i)))^2);
                    end
                end
            end
            
            % Removing overlapping points, but preserving the order of the rows
            [~,idx] = uniquetol(obj.nIBC,obj.para.overlaptol,'ByRows',true);
            obj.nIBC = obj.nIBC(sort(idx),:);
        end
        %% calcIntersection: Calculate points on inner boundary curve and the SFS
        function calcIntersection(obj)
            % In: IBCF,F1,kF1
            % Out: isec
            F1Edge = zeros(2,1);
            obj.isec = [];
            
            for k = 1:size(obj.IBCF,1)
                if obj.para.debug
                    % Plot ellipse
                    alpha2=-10:0.005:10;
                    alpha_res2=zeros(numel(alpha2),2);
                    for j = 1:numel(alpha2)
                        alpha_res2(j,1) = (obj.IBCF(k,2)+obj.IBCF(k,3)*tan(alpha2(j))+obj.IBCF(k,4)*(tan(alpha2(j)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(alpha2(j))+obj.IBCF(k,10)*(tan(alpha2(j)))^2);
                        alpha_res2(j,2) = (obj.IBCF(k,5)+obj.IBCF(k,6)*tan(alpha2(j))+obj.IBCF(k,7)*(tan(alpha2(j)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(alpha2(j))+obj.IBCF(k,10)*(tan(alpha2(j)))^2);
                    end
                    figure(81)
                    hold on
                    p1 = plot(alpha_res2(:,1),alpha_res2(:,2),'-r');
                end
                for i = 1:size(obj.kF1,1)-1
                    [F1Edge(1),F1Edge(2)] = IBCA.calcLEQ(obj.F1(obj.kF1(i),:),obj.F1(obj.kF1(i+1),:));
                    % Calculate both intersections of the conic section with an edge of FIRPOL
                    is1Talp = -1/2/(F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*F1Edge(1)*obj.IBCF(k,3)-1/2/...
                        (F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*F1Edge(2)*obj.IBCF(k,6)-1/2/...
                        (F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*obj.IBCF(k,9)+1/2/...
                        (F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*(-4*F1Edge(1)^2*obj.IBCF(k,2)*...
                        obj.IBCF(k,4)+F1Edge(1)^2*obj.IBCF(k,3)^2-4*F1Edge(1)*F1Edge(2)*obj.IBCF(k,2)*obj.IBCF(k,7)+2*...
                        F1Edge(1)*F1Edge(2)*obj.IBCF(k,3)*obj.IBCF(k,6)-4*F1Edge(1)*F1Edge(2)*obj.IBCF(k,4)*obj.IBCF(k,5)-...
                        4*F1Edge(2)^2*obj.IBCF(k,5)*obj.IBCF(k,7)+F1Edge(2)^2*obj.IBCF(k,6)^2-4*F1Edge(1)*obj.IBCF(k,8)*...
                        obj.IBCF(k,4)+2*F1Edge(1)*obj.IBCF(k,9)*obj.IBCF(k,3)-4*F1Edge(1)*obj.IBCF(k,10)*obj.IBCF(k,2)-...
                        4*F1Edge(2)*obj.IBCF(k,8)*obj.IBCF(k,7)+2*F1Edge(2)*obj.IBCF(k,9)*obj.IBCF(k,6)-...
                        4*F1Edge(2)*obj.IBCF(k,10)*obj.IBCF(k,5)-4*obj.IBCF(k,8)*obj.IBCF(k,10)+obj.IBCF(k,9)^2)^(1/2);
                    
                    is1x = (obj.IBCF(k,2)+obj.IBCF(k,3)*is1Talp+obj.IBCF(k,4)*(is1Talp)^2)/...
                        (obj.IBCF(k,8)+obj.IBCF(k,9)*is1Talp+obj.IBCF(k,10)*(is1Talp)^2);
                    is1y = (obj.IBCF(k,5)+obj.IBCF(k,6)*is1Talp+obj.IBCF(k,7)*(is1Talp)^2)/...
                        (obj.IBCF(k,8)+obj.IBCF(k,9)*is1Talp+obj.IBCF(k,10)*(is1Talp)^2);
                    
                    is2Talp = -1/2/(F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*F1Edge(1)*obj.IBCF(k,3)-1/2/...
                        (F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*F1Edge(2)*obj.IBCF(k,6)-1/2/...
                        (F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*obj.IBCF(k,9)-1/2/...
                        (F1Edge(1)*obj.IBCF(k,4)+F1Edge(2)*obj.IBCF(k,7)+obj.IBCF(k,10))*(-4*F1Edge(1)^2*obj.IBCF(k,2)*...
                        obj.IBCF(k,4)+F1Edge(1)^2*obj.IBCF(k,3)^2-4*F1Edge(1)*F1Edge(2)*obj.IBCF(k,2)*obj.IBCF(k,7)+2*...
                        F1Edge(1)*F1Edge(2)*obj.IBCF(k,3)*obj.IBCF(k,6)-4*F1Edge(1)*F1Edge(2)*obj.IBCF(k,4)*obj.IBCF(k,5)-...
                        4*F1Edge(2)^2*obj.IBCF(k,5)*obj.IBCF(k,7)+F1Edge(2)^2*obj.IBCF(k,6)^2-4*F1Edge(1)*obj.IBCF(k,8)*...
                        obj.IBCF(k,4)+2*F1Edge(1)*obj.IBCF(k,9)*obj.IBCF(k,3)-4*F1Edge(1)*obj.IBCF(k,10)*obj.IBCF(k,2)-...
                        4*F1Edge(2)*obj.IBCF(k,8)*obj.IBCF(k,7)+2*F1Edge(2)*obj.IBCF(k,9)*obj.IBCF(k,6)-...
                        4*F1Edge(2)*obj.IBCF(k,10)*obj.IBCF(k,5)-4*obj.IBCF(k,8)*obj.IBCF(k,10)+obj.IBCF(k,9)^2)^(1/2);
                    is2x = (obj.IBCF(k,2)+obj.IBCF(k,3)*is2Talp+obj.IBCF(k,4)*(is2Talp)^2)/...
                        (obj.IBCF(k,8)+obj.IBCF(k,9)*is2Talp+obj.IBCF(k,10)*(is2Talp)^2);
                    is2y = (obj.IBCF(k,5)+obj.IBCF(k,6)*is2Talp+obj.IBCF(k,7)*(is2Talp)^2)/...
                        (obj.IBCF(k,8)+obj.IBCF(k,9)*is2Talp+obj.IBCF(k,10)*(is2Talp)^2);
                    
                    if obj.para.debug
                        p4 = plot([-2 2], [(F1Edge(1)*(-2)+1)/(-F1Edge(2)) (F1Edge(1)*(2)+1)/(-F1Edge(2))],'r');
                        p5 = plot([(F1Edge(2)*(-2)+1)/(-F1Edge(1)) (F1Edge(2)*(2)+1)/(-F1Edge(1))],[-2 2], 'r');
                    end
                    
                    % Find valid intersections
                    if isIntersection(obj,is1Talp,is1x,is1y,F1Edge,i,k)
                        if isempty(obj.isec)
                            obj.isec = [obj.isec;[is1Talp,is1x,is1y,i,k]];
                        else
                            if ~any(ismembertol(obj.isec(:,2:3),[is1x,is1y],obj.para.calcIsectol,'ByRows',true))
                                obj.isec = [obj.isec;[is1Talp,is1x,is1y,i,k]];
                            end
                        end
                    end
                    if isIntersection(obj, is2Talp,is2x,is2y,F1Edge,i,k)
                        if isempty(obj.isec)
                            obj.isec = [obj.isec;[is2Talp,is2x,is2y,i,k]];
                        else
                            if ~any(ismembertol(obj.isec(:,2:3),[is2x,is2y],obj.para.calcIsectol,'ByRows',true))
                                obj.isec = [obj.isec;[is2Talp,is2x,is2y,i,k]];
                            end
                        end
                    end
                    if obj.para.debug
                        p2 = plot(is1x,is1y,'*r');
                        p3 = plot(is2x,is2y,'*b');
                        delete(p2)
                        delete(p3)
                        delete(p4)
                        delete(p5)
                    end
                end
                if obj.para.debug
                    delete(p1)
                end
            end
        end
        %% isIntersection: Intersection point corectness tests
        function res = isIntersection(obj,is1Talp,is1x,is1y,F1Edge,i,k)
            % Test whether a given point is a valid intersection of FP and IBC
            % In: is1Talp,is1x,is1y,F1Edge,F1,kF1,i,IBCF,k
            % Out: res
            res = 0; % Means that it is not a valid intersection
            % Test whether the solution is real
            realityTest = (abs(imag(is1x))<=obj.para.isIsectol && abs(imag(is1y))<=obj.para.isIsectol);
            if realityTest
                is1x = real(is1x);
                is1y = real(is1y);
                % Test whether the point is on a line of FP
                onFPTest = (abs(F1Edge(1)*is1x+F1Edge(2)*is1y+1)<obj.para.isIsectol);
                if onFPTest
                    % Test whether the point is on the correct interval of the FP edge
                    % with equations
                    xIntervalTest = (is1x>=obj.F1(obj.kF1(i),1)-obj.para.isIsectol && is1x<=obj.F1(obj.kF1(i+1),1)+obj.para.isIsectol) ||...
                        (is1x>=obj.F1(obj.kF1(i+1),1)-obj.para.isIsectol && is1x<=obj.F1(obj.kF1(i),1)+obj.para.isIsectol);
                    yIntervalTest = (is1y>=obj.F1(obj.kF1(i),2)-obj.para.isIsectol && is1y<=obj.F1(obj.kF1(i+1),2)+obj.para.isIsectol) ||...
                        (is1y>=obj.F1(obj.kF1(i+1),2)-obj.para.isIsectol && is1y<=obj.F1(obj.kF1(i),2)+obj.para.isIsectol);
                    if xIntervalTest && yIntervalTest
                        % Test whether the point is on the correct interval of IBC
                        % (The next test is most likely better for this)
                        isAngle = atan2(is1y,is1x);
                        startx = (obj.IBCF(k,2)+obj.IBCF(k,3)*tan(obj.IBCF(k,1))+obj.IBCF(k,4)*(tan(obj.IBCF(k,1)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(obj.IBCF(k,1))+obj.IBCF(k,10)*(tan(obj.IBCF(k,1)))^2);
                        starty = (obj.IBCF(k,5)+obj.IBCF(k,6)*tan(obj.IBCF(k,1))+obj.IBCF(k,7)*(tan(obj.IBCF(k,1)))^2)/...
                            (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(obj.IBCF(k,1))+obj.IBCF(k,10)*(tan(obj.IBCF(k,1)))^2);
                        startAngle = atan2(starty,startx);
                        if k<size(obj.IBCF,1)
                            endx = (obj.IBCF(k,2)+obj.IBCF(k,3)*tan(obj.IBCF(k+1,1))+obj.IBCF(k,4)*(tan(obj.IBCF(k+1,1)))^2)/...
                                (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(obj.IBCF(k+1,1))+obj.IBCF(k,10)*(tan(obj.IBCF(k+1,1)))^2);
                            endy = (obj.IBCF(k,5)+obj.IBCF(k,6)*tan(obj.IBCF(k+1,1))+obj.IBCF(k,7)*(tan(obj.IBCF(k+1,1)))^2)/...
                                (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(obj.IBCF(k+1,1))+obj.IBCF(k,10)*(tan(obj.IBCF(k+1,1)))^2);
                            endAngle = atan2(endy,endx);
                        else
                            endx = (obj.IBCF(k,2)+obj.IBCF(k,3)*tan(obj.IBCF(1,1))+obj.IBCF(k,4)*(tan(obj.IBCF(1,1)))^2)/...
                                (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(obj.IBCF(1,1))+obj.IBCF(k,10)*(tan(obj.IBCF(1,1)))^2);
                            endy = (obj.IBCF(k,5)+obj.IBCF(k,6)*tan(obj.IBCF(1,1))+obj.IBCF(k,7)*(tan(obj.IBCF(1,1)))^2)/...
                                (obj.IBCF(k,8)+obj.IBCF(k,9)*tan(obj.IBCF(1,1))+obj.IBCF(k,10)*(tan(obj.IBCF(1,1)))^2);
                            endAngle = atan2(endy,endx);
                        end
                        if startAngle<=endAngle
                            tanTest = (isAngle>=startAngle-obj.para.isIsectol && isAngle<=endAngle+obj.para.isIsectol);
                        else
                            tanTest = (isAngle>=startAngle-obj.para.isIsectol || isAngle<=endAngle+obj.para.isIsectol);
                        end
                        if tanTest
                            % Test whether the tangens of alpha is in the correct interval for the IBC
                            starTalp = tan(obj.IBCF(k,1));
                            if k<size(obj.IBCF,1)
                                endTalp = tan(obj.IBCF(k+1,1));
                            else
                                endTalp = tan(obj.IBCF(1,1));
                            end
                            if starTalp<=endTalp
                                isTalpTest = (is1Talp>=starTalp-obj.para.isIsectol && is1Talp<=endTalp+obj.para.isIsectol);
                            else
                                isTalpTest = (is1Talp>=starTalp-obj.para.isIsectol || is1Talp<=endTalp+obj.para.isIsectol);
                            end
                            if isTalpTest
                                res = 1;
                            end
                        end
                    end
                end
            end
        end
        %% calcSFS: Calculate the SFS
        function calcSFS(obj)
            % In: isec,nIBC,F1,kF1
            % Out: SFS
            obj.SFS = {};
            % If intersection points were not found
            if isempty(obj.isec)
                % Check whether there exists a solution (one donut segment)
                if inpolygon(obj.nIBC(1,1),obj.nIBC(1,2),obj.F1(obj.kF1(:),1),obj.F1(obj.kF1(:),2))
                    obj.SFS(end+1) = {obj.F1(obj.kF1(:),:)};
                    obj.SFS(end+1) = {obj.nIBC};
                end
            else
                % Sort intersection points by angles
                isecAngle = atan2(obj.isec(:,3),obj.isec(:,2));
                [isecAngle,idx] = sort(isecAngle);
                obj.isec = obj.isec(idx,:);
                
                % Sort FIRPOL vertices by angles in the opposite order
                FAngle = atan2(obj.F1(obj.kF1(1:size(obj.kF1,1)-1),2),obj.F1(obj.kF1(1:size(obj.kF1,1)-1),1));
                [FAngle,idx] = sort(FAngle,'descend');
                F = obj.F1(obj.kF1(idx),:);
                
                for i = 1:size(obj.isec,1)
                    % Find the closest nIBC index to a given intersection point
                    isecMat = [obj.isec(i,2)*ones(size(obj.nIBC,1),1) obj.isec(i,3)*ones(size(obj.nIBC,1),1)];
                    [~, closestIndex] = min(sqrt(sum((isecMat - obj.nIBC).^2,2)));
                    % Find the index of the point left and right of the closest index
                    if closestIndex == size(obj.nIBC,1)
                        closestIndexPlus1 = 1;
                        closestIndexMinus1 = closestIndex-1;
                    elseif closestIndex == 1
                        closestIndexPlus1 = closestIndex+1;
                        closestIndexMinus1 = size(obj.nIBC,1);
                    else
                        closestIndexPlus1 = closestIndex+1;
                        closestIndexMinus1 = closestIndex-1;
                    end
                    % If the closest nIBC index plus one is in FIRPOL
                    if inpolygon(obj.nIBC(closestIndexPlus1,1),obj.nIBC(closestIndexPlus1,2),obj.F1(obj.kF1(:),1),obj.F1(obj.kF1(:),2))
                        % Find the next isec index
                        if i == size(obj.isec,1)
                            iNext = 1;
                        else
                            iNext = i+1;
                        end
                        % Find the next nIBC index that is closest to the next isec
                        isecMatNext = [obj.isec(iNext,2)*ones(size(obj.nIBC,1),1) obj.isec(iNext,3)*ones(size(obj.nIBC,1),1)];
                        [~, closestIndexNext] = min(sqrt(sum((isecMatNext - obj.nIBC).^2,2)));
                        if not(inpolygon(obj.nIBC(closestIndexNext,1),obj.nIBC(closestIndexNext,2),...
                                obj.F1(obj.kF1(:),1),obj.F1(obj.kF1(:),2)))
                            if closestIndexNext == 1
                                closestIndexNext = size(obj.nIBC,1);
                            else
                                closestIndexNext = closestIndexNext-1;
                            end
                        end
                        % Check whether the nIBC indices loop over the end of the array and
                        % save them
                        if closestIndexPlus1>closestIndexNext
                            nIBCsegIndex = [closestIndexPlus1:size(obj.nIBC,1),1:closestIndexNext];
                        else
                            nIBCsegIndex = [closestIndexPlus1:closestIndexNext];
                        end
                        % Find the possible FIRPOL vertices and check for loop over the end
                        % of the array
                        if isecAngle(i)<isecAngle(iNext)
                            Fidx = and((FAngle > isecAngle(i)),(FAngle < isecAngle(iNext)));
                            segmentSFS = [obj.isec(i,2:3); obj.nIBC(nIBCsegIndex,:); obj.isec(iNext,2:3); F(Fidx,:)];
                        else
                            Fidx = (FAngle < isecAngle(iNext));
                            Fidx2 = (FAngle > isecAngle(i));
                            segmentSFS = [obj.isec(i,2:3); obj.nIBC(nIBCsegIndex,:); obj.isec(iNext,2:3); F(Fidx,:); F(Fidx2,:)];
                        end
                        % Save the segment
                        obj.SFS(end+1) = {segmentSFS};
                        if obj.para.debug
                            figure(81)
                            plot(polyshape(obj.SFS{end}(:,1),obj.SFS{end}(:,2)))
                        end
                        % If the neither the next nor the previous nIBC index is in FIRPOL
                    elseif  not(inpolygon(obj.nIBC(closestIndexMinus1,1),obj.nIBC(closestIndexMinus1,2),...
                            obj.F1(obj.kF1(:),1),obj.F1(obj.kF1(:),2)))
                        % Add as an isolated point
                        obj.SFS(end+1) = {obj.isec(i,2:3)};
                    end
                end
            end
        end
        %% calcSFSxor: Calculate SFS without the use of isec
        function calcSFSxor(obj)
            % In: nIBC,F1,kF1
            % Out: SFS
            if obj.para.debug==0
                warning off
            end
            p1 = polyshape(obj.F1(obj.kF1,1),obj.F1(obj.kF1,2));
            p2 = polyshape(obj.nIBC(:,1),obj.nIBC(:,2));
%             p3 = intersect(p1,p2);
%             p4 = xor(p1,p3);
            p3 = xor(p1,p2);
            p4 = intersect(p1,p3);
            if obj.para.debug==0
                warning on
            end
            fullSFS = p4.Vertices;
            if ~isempty(fullSFS)
                idx = find(isnan(fullSFS(:,1)));
                if isempty(idx) % One segment (should be impossible)
                    obj.SFS{1} = fullSFS;
                else
                    for i = 1:length(idx)+1
                        if i==1
                            obj.SFS{i} = fullSFS(1:idx(i),:);
                        elseif i==length(idx)+1
                            obj.SFS{i} = fullSFS(idx(i-1):end,:);
                        else
                            obj.SFS{i} = fullSFS(idx(i-1):idx(i),:);
                        end
                    end
                end
            end
            % Remove NaN and last row = first row
            for i = 1:size(obj.SFS,2)
                obj.SFS{i} = obj.SFS{i}(all(~isnan(obj.SFS{i}),2),:);
                obj.SFS{i} = [obj.SFS{i};obj.SFS{i}(1,:)];
            end
        end
        %% plotBoundary: Plot INNPOL and FIRPOL
        function [pFP,pIP] = plotBoundary(obj,factor,ax)
            if nargin==1
                figure(81)
                hold on
                ax = gca;
                factor = 1;
            elseif nargin==2
                figure(81)
                hold on
                ax = gca;
            end
            if factor==1
                F = obj.F1;
                kF = obj.kF1;
                I = obj.I1;
                kI = obj.kI1;
                colF = 'b';
                colI = 'r';
            elseif factor==2
                F = obj.F2;
                kF = obj.kF2;
                I = obj.I2;
                kI = obj.kI2;
                colF = 'r';
                colI = 'b';
            end
            pFP = plot(ax,F(kF,1),F(kF,2),'Color',colF,'Visible','off');
            pIP = plot(ax,I(kI,1),I(kI,2),'-o','Color',colI,'Visible','off');
            if nargin<3
                set(pFP,'Visible','on');
                set(pIP,'Visible','on');
            end
        end
        %% plotIBC: Plot inner boundary curve
        function pIBC = plotIBC(obj,ax)
            if nargin==1
                figure(81)
                hold on
                ax = gca;
            end
            pIBC = plot(ax, obj.nIBC(:,1),obj.nIBC(:,2),'-k','Visible','off');
            if nargin==1
                set(pIBC,'Visible','on');
            end
        end
        %% plotSFS: Plot the set of feasible solutions
        function pSFS = plotSFS(obj,ax)
            if nargin==1
                figure(81)
                hold on
                ax = gca;
            end
            if obj.para.debug==0
                warning off
            end
            if size(obj.SFS,2)==2
                poly1 = polyshape(obj.SFS{1}(:,1),obj.SFS{1}(:,2));
                poly2 = polyshape(obj.SFS{2}(:,1),obj.SFS{2}(:,2));
                pSFS{1} = plot(ax,xor(poly1,poly2),'Visible','off');
            elseif size(obj.SFS,2)>2
                for i = 1:size(obj.SFS,2)
                    if size(obj.SFS{i},1) == 1
                        pSFS{i} = plot(ax,obj.SFS{i}(1,1),obj.SFS{i}(1,2),'o','Visible','off');
                        set(pSFS{i}, 'markerfacecolor', get(pSFS{i}, 'color'));
                    elseif  size(obj.SFS{i},1) == 2
                        pSFS{i} = plot(ax,obj.SFS{i}(:,1),obj.SFS{i}(:,2),'Visible','off');
                    else
                        pSFS{i} = plot(ax,polyshape(obj.SFS{i}(:,1),obj.SFS{i}(:,2)),'Visible','off');
                    end
                end
            else
                pSFS{1} = plot([],[]);
            end
            if obj.para.debug==0
                warning on
            end
            if nargin==1
                if size(obj.SFS,2)==2
                    set(pSFS{1},'Visible','on');
                else
                    for i = 1:size(obj.SFS,2)
                        set(pSFS{i},'Visible','on');
                    end
                end
            end
        end
        %% calcTriangle: Calculate the coordinates for the boundary triangles
        function calcTriangle(obj)            
            % In: alpha, TV
            % Out: pTriangle, pTangents, dTriangle, dTangents, iTr
            obj.pTriangle = zeros(3,2);
            obj.pTangents = zeros(3,2); 
            obj.dTriangle = zeros(3,2);
            obj.dTangents = zeros(2,2);
            % Find correct parameters
            for j=1:size(obj.IBCF,1)-1
                k = size(obj.IBCF,1);
                if obj.alpha<=obj.IBCF(1,1) || obj.alpha>obj.IBCF(k,1)
                    obj.iTr = k;
                elseif obj.alpha<=obj.IBCF(j+1,1) && obj.alpha>obj.IBCF(j,1)
                    obj.iTr = j;
                end
            end
            % Find the triangles
            % v
            obj.pTangents(1,:) = [obj.TV(obj.iTr,2) obj.TV(obj.iTr,3)];
            % t
            obj.dTriangle(1,1) = -1/(obj.TV(obj.iTr,2)+obj.TV(obj.iTr,3)*tan(obj.alpha));
            obj.dTriangle(1,2) = -tan(obj.alpha)/(obj.TV(obj.iTr,2)+obj.TV(obj.iTr,3)*tan(obj.alpha));
            % u_l and u_r
            obj.dTangents(1,:) = [obj.TV(obj.iTr,4) obj.TV(obj.iTr,5)];
            obj.dTangents(2,:) = [obj.TV(obj.iTr,6) obj.TV(obj.iTr,7)];
            % s_l and s_r
            [obj.pTriangle(1,1),obj.pTriangle(1,2)] = obj.calcLEQ(obj.dTriangle(1,:),obj.dTangents(1,:));
            [obj.pTriangle(2,1),obj.pTriangle(2,2)] = obj.calcLEQ(obj.dTriangle(1,:),obj.dTangents(2,:));
            % o_l and o_r
            obj.pTangents(2,:) = [obj.TV(obj.iTr,8) obj.TV(obj.iTr,9)];
            obj.pTangents(3,:) = [obj.TV(obj.iTr,10) obj.TV(obj.iTr,11)];
            % q_l and q_r
            [obj.dTriangle(2,1),obj.dTriangle(2,2)] = obj.calcLEQ(obj.pTriangle(1,:),obj.pTangents(2,:));
            [obj.dTriangle(3,1),obj.dTriangle(3,2)] = obj.calcLEQ(obj.pTriangle(2,:),obj.pTangents(3,:));
            % p
            [obj.pTriangle(3,1),obj.pTriangle(3,2)] = obj.calcLEQ(obj.dTriangle(2,:),obj.dTriangle(3,:));
        end
        %% Load result
        function loadRes(obj)
            [fname, fpath] = uigetfile({'*.mat'},...
                'Select *.mat file with an existing IBCA result');
            if max(fname)>0 && max(fpath)>0
                Sf=load([fpath fname]);
                names=fieldnames(Sf);
                if not(isempty(names))
                    if ismember('D',names)
                        obj.D = getfield(Sf,'D');
                    end
                    if ismember('T',names)
                        obj.T = getfield(Sf,'T');
                    end
                    if ismember('X',names)
                        obj.X = getfield(Sf,'X');
                    end
                    if ismember('U',names)
                        obj.U = getfield(Sf,'U');
                    end
                    if ismember('S',names)
                        obj.S = getfield(Sf,'S');
                    end
                    if ismember('V',names)
                        obj.V = getfield(Sf,'V');
                    end
                    if ismember('FIRPOL_1',names)
                        obj.F1 = getfield(Sf,'FIRPOL_1');
                        obj.kF1 = 1:size(obj.F1,1);
                    end
                    if ismember('INNPOL_1',names)
                        obj.I1 = getfield(Sf,'INNPOL_1');
                        obj.kI1 = 1:size(obj.I1,1);
                    end
                    if ismember('FIRPOL_2',names)
                        obj.F2 = getfield(Sf,'FIRPOL_2');
                        obj.kF2 = 1:size(obj.F2,1);
                    end
                    if ismember('INNPOL_2',names)
                        obj.I2 = getfield(Sf,'INNPOL_2');
                        obj.kI2 = 1:size(obj.I2,1);
                    end
                    if ismember('Contact_Change_Val_Param',names)
                        obj.TV = getfield(Sf,'Contact_Change_Val_Param');
                    end
                    if ismember('Inn_Bound_Cur_Fun_Coeffs',names)
                        obj.IBCF = getfield(Sf,'Inn_Bound_Cur_Fun_Coeffs');
                    end
                    if ismember('Inn_Bound_Cur_Approx',names)
                        obj.nIBC = getfield(Sf,'Inn_Bound_Cur_Approx');
                    end
                    if ismember('Intersect_FIRPOL_Inn_Bound_Cur',names)
                        obj.isec = getfield(Sf,'Intersect_FIRPOL_Inn_Bound_Cur');
                    end
                    if ismember('Set_of_Feasib_Sol',names)
                        obj.SFS = getfield(Sf,'Set_of_Feasib_Sol');
                    end
                end
            end
        end 

        %% plotConic: Plot the whole conic section of a section of the IBC
        % This is auxiliary function that plots corresponding conic sections
        % In the if part the correct interval hast to be given
        % Alternatively you can delete the if-part and use breakpoints
        function plotConic(obj,ax)
            % j=16, j=18
            if nargin==1
                figure(81)
                hold on
                ax = gca;
            end
            
            alp=0:0.005:pi;
            for j=1:size(obj.IBCF,1)-1
                if j==17 % use 5 and 17 for 9x8 model problem
                    for i=1:numel(alp)
                        conic(i,1) = (obj.IBCF(j,2)+obj.IBCF(j,3)*tan(alp(i))+obj.IBCF(j,4)*(tan(alp(i)))^2)/...
                            (obj.IBCF(j,8)+obj.IBCF(j,9)*tan(alp(i))+obj.IBCF(j,10)*(tan(alp(i)))^2);
                        conic(i,2) = (obj.IBCF(j,5)+obj.IBCF(j,6)*tan(alp(i))+obj.IBCF(j,7)*(tan(alp(i)))^2)/...
                            (obj.IBCF(j,8)+obj.IBCF(j,9)*tan(alp(i))+obj.IBCF(j,10)*(tan(alp(i)))^2);
                    end
                    alpo=obj.IBCF(j,1):0.005:obj.IBCF(j+1,1);
                    for i=1:numel(alpo)
                        conicoverlap(i,1) = (obj.IBCF(j,2)+obj.IBCF(j,3)*tan(alpo(i))+obj.IBCF(j,4)*(tan(alpo(i)))^2)/...
                            (obj.IBCF(j,8)+obj.IBCF(j,9)*tan(alpo(i))+obj.IBCF(j,10)*(tan(alpo(i)))^2);
                        conicoverlap(i,2) = (obj.IBCF(j,5)+obj.IBCF(j,6)*tan(alpo(i))+obj.IBCF(j,7)*(tan(alpo(i)))^2)/...
                            (obj.IBCF(j,8)+obj.IBCF(j,9)*tan(alpo(i))+obj.IBCF(j,10)*(tan(alpo(i)))^2);
                    end
                    plot(ax, conic(:,1),conic(:,2),'-g')
                    plot(ax, conicoverlap(:,1),conicoverlap(:,2),'color',[0 0.5 0])
                end
            end
        end
    end
    %% Static methods
    methods(Static)
        %% calcTCH: Find the tangent using convex hull
        function [z1l,z2l,z1r,z2r] = calcTCH(I,F,kF,i)
            % Find the convex hull of INNPOL and one vertex of FIRPOL
            % to find the tangents of INNPOL.
            Itg = [I;F(kF(i),:)];
            kItg = convhull(Itg);
            % Index of the FIRPOL vertex
            iF = find(kItg==size(I,1)+1,1);
            
            % If the point is too close to INNPOL, take the INNPOL point
            if size(iF,1)==0
                iF = find(all(ismembertol(Itg(kItg,:),Itg(size(I,1)+1,:),10e-14)==1,2)==1);
            end
            % If the point is near the convex hull edge, take neigbouring
            % points
            if size(iF,1)==0
                for i = 1:size(kItg)-1
                    distToEdge = IBCA.calcMinDistPointLine(Itg(size(I,1)+1,:),Itg(kItg(i),:),Itg(kItg(i+1),:));
                    if distToEdge<1e-3
                        kItg = [kItg(1:i) ;size(I,1)+1; kItg(i+1:end)];
                        break;
                    end
                end
                iF = find(kItg==size(I,1)+1,1);
            end
            iF = iF(1);
                        
            % Index of the left and right tangential points
            if iF==1
                iTl = size(kItg,1)-1;
                iTr = 2;
            else
                iTl = iF-1;
                iTr = iF+1;
            end
            % Left tangent
            [z1l,z2l] = IBCA.calcLEQ(Itg(kItg(iF),:),Itg(kItg(iTl),:));
            % Right tangent
            [z1r,z2r] = IBCA.calcLEQ(Itg(kItg(iF),:),Itg(kItg(iTr),:));
        end
        %% calcMinDistPointLine: Calculate the shortest distance between a point an a line
        function dist = calcMinDistPointLine(p,le1,le2)
            % p point, le1 first line endpoint, le2 second line endpoint
            testBetween = -((le1(1)-p(1))*(le2(1)-le1(1))+(le1(2)-p(2))*(le2(2)-le1(2)))/((le2(1)-le1(1))^2+(le2(2)-le1(2))^2);
            if testBetween>=0 && testBetween <=1
                dist = abs((le2(1)-le1(1))*(le1(2)-p(2))-(le2(2)-le1(2))*(le1(1)-p(1)))/sqrt((le2(1)-le1(1))^2+(le2(2)-le1(2))^2);
            else
                dist = min(sqrt((le1(1)-p(1))^2+(le1(2)-p(2))^2),sqrt((le2(1)-p(1))^2+(le2(2)-p(2))^2));
            end
        end
        %% calcLEQ: Calculate the intersection of 2 lines or the line equation from 2 points
        function [e1, e2] = calcLEQ(g, h)
            e1 = (-g(2)+h(2))/(g(2)*h(1)-g(1)*h(2));
            e2 = (-g(1)+h(1))/(g(1)*h(2)-g(2)*h(1));
        end
        %% calcExpression: Analytical expression for the inner boundary curve
        function res = calcExpression(A)
            % A(:,1)  is  v_1
            % A(:,2)  is  v_2
            % A(:,3)  is  u_l1
            % A(:,4)  is  u_l2
            % A(:,5)  is  u_r1
            % A(:,6)  is  u_r2
            % A(:,7)  is  o_l1
            % A(:,8)  is  o_12
            % A(:,9)  is  o_r1
            % A(:,10) is  o_r2
            
            % res(:,1) is c_x0
            % res(:,2) is c_x1
            % res(:,3) is c_x2
            % res(:,4) is c_y0
            % res(:,5) is c_y1
            % res(:,6) is c_y2
            % res(:,7) is c_d0
            % res(:,8) is c_d1
            % res(:,9) is c_d2
            
            res = zeros(size(A,1),9);
            res(:,1) = -A(:,7).*A(:,9).*A(:,3).*A(:,6).*A(:,1)+A(:,7).*A(:,9).*A(:,4).*A(:,5).*A(:,1)+...
                A(:,7).*A(:,10).*A(:,4).*A(:,6).*A(:,1)+A(:,7).*A(:,3).*A(:,6).*A(:,1).^2-A(:,8).*A(:,9).*A(:,4).*A(:,6).*A(:,1)+...
                A(:,8).*A(:,4).*A(:,6).*A(:,1).^2-A(:,9).*A(:,4).*A(:,5).*A(:,1).^2-A(:,10).*A(:,4).*A(:,6).*A(:,1).^2+...
                A(:,7).*A(:,9).*A(:,4)-A(:,7).*A(:,9).*A(:,6)+A(:,7).*A(:,6).*A(:,1)-A(:,9).*A(:,4).*A(:,1);
            res(:,2) = -A(:,7).*A(:,9).*A(:,3).*A(:,6).*A(:,2)+A(:,7).*A(:,9).*A(:,4).*A(:,5).*A(:,2)-...
                A(:,7).*A(:,10).*A(:,3).*A(:,6).*A(:,1)+A(:,7).*A(:,10).*A(:,4).*A(:,6).*A(:,2)+2.*A(:,7).*A(:,3).*A(:,6).*A(:,1).*A(:,2)+...
                A(:,8).*A(:,9).*A(:,4).*A(:,5).*A(:,1)-A(:,8).*A(:,9).*A(:,4).*A(:,6).*A(:,2)+2.*A(:,8).*A(:,4).*A(:,6).*A(:,1).*A(:,2)-...
                2.*A(:,9).*A(:,4).*A(:,5).*A(:,1).*A(:,2)-2.*A(:,10).*A(:,4).*A(:,6).*A(:,1).*A(:,2)-A(:,7).*A(:,9).*A(:,3)+...
                A(:,7).*A(:,9).*A(:,5)+A(:,7).*A(:,10).*A(:,4)+A(:,7).*A(:,3).*A(:,1)+A(:,7).*A(:,6).*A(:,2)-...
                A(:,8).*A(:,9).*A(:,6)+A(:,8).*A(:,4).*A(:,1)+A(:,8).*A(:,6).*A(:,1)-A(:,9).*A(:,4).*A(:,2)-...
                A(:,9).*A(:,5).*A(:,1)-A(:,10).*A(:,4).*A(:,1)-A(:,10).*A(:,6).*A(:,1)+A(:,7)-A(:,9);
            res(:,3) = -A(:,7).*A(:,10).*A(:,3).*A(:,6).*A(:,2)+A(:,7).*A(:,3).*A(:,6).*A(:,2).^2+A(:,8).*A(:,9).*A(:,4).*A(:,5).*A(:,2)+...
                A(:,8).*A(:,4).*A(:,6).*A(:,2).^2-A(:,9).*A(:,4).*A(:,5).*A(:,2).^2-A(:,10).*A(:,4).*A(:,6).*A(:,2).^2-...
                A(:,7).*A(:,10).*A(:,3)+A(:,7).*A(:,3).*A(:,2)+A(:,8).*A(:,9).*A(:,5)+A(:,8).*A(:,4).*A(:,2)+A(:,8).*A(:,6).*A(:,2)-...
                A(:,9).*A(:,5).*A(:,2)-A(:,10).*A(:,4).*A(:,2)-A(:,10).*A(:,6).*A(:,2)+A(:,8)-A(:,10);
            res(:,4) = -A(:,7).*A(:,10).*A(:,3).*A(:,6).*A(:,1)-A(:,7).*A(:,3).*A(:,5).*A(:,1).^2+A(:,8).*A(:,9).*A(:,4).*A(:,5).*A(:,1)-...
                A(:,8).*A(:,4).*A(:,5).*A(:,1).^2+A(:,9).*A(:,3).*A(:,5).*A(:,1).^2+A(:,10).*A(:,3).*A(:,6).*A(:,1).^2-...
                A(:,7).*A(:,10).*A(:,6)-A(:,7).*A(:,3).*A(:,1)-A(:,7).*A(:,5).*A(:,1)+A(:,8).*A(:,9).*A(:,4)-A(:,8).*A(:,4).*A(:,1)+...
                A(:,9).*A(:,3).*A(:,1)+A(:,9).*A(:,5).*A(:,1)+A(:,10).*A(:,6).*A(:,1)-A(:,7)+A(:,9);
            res(:,5) = A(:,7).*A(:,10).*A(:,3).*A(:,5).*A(:,1)-A(:,7).*A(:,10).*A(:,3).*A(:,6).*A(:,2)-...
                2.*A(:,7).*A(:,3).*A(:,5).*A(:,1).*A(:,2)-A(:,8).*A(:,9).*A(:,3).*A(:,5).*A(:,1)+A(:,8).*A(:,9).*A(:,4).*A(:,5).*A(:,2)-...
                A(:,8).*A(:,10).*A(:,3).*A(:,6).*A(:,1)+A(:,8).*A(:,10).*A(:,4).*A(:,5).*A(:,1)-2.*A(:,8).*A(:,4).*A(:,5).*A(:,1).*A(:,2)+...
                2.*A(:,9).*A(:,3).*A(:,5).*A(:,1).*A(:,2)+2.*A(:,10).*A(:,3).*A(:,6).*A(:,1).*A(:,2)+A(:,7).*A(:,10).*A(:,5)-...
                A(:,7).*A(:,3).*A(:,2)-A(:,7).*A(:,5).*A(:,2)-A(:,8).*A(:,9).*A(:,3)+A(:,8).*A(:,10).*A(:,4)-A(:,8).*A(:,10).*A(:,6)-...
                A(:,8).*A(:,4).*A(:,2)-A(:,8).*A(:,5).*A(:,1)+A(:,9).*A(:,3).*A(:,2)+A(:,9).*A(:,5).*A(:,2)+A(:,10).*A(:,3).*A(:,1)+...
                A(:,10).*A(:,6).*A(:,2)-A(:,8)+A(:,10);
            res(:,6) = A(:,7).*A(:,10).*A(:,3).*A(:,5).*A(:,2)-A(:,7).*A(:,3).*A(:,5).*A(:,2).^2-A(:,8).*A(:,9).*A(:,3).*A(:,5).*A(:,2)-...
                A(:,8).*A(:,10).*A(:,3).*A(:,6).*A(:,2)+A(:,8).*A(:,10).*A(:,4).*A(:,5).*A(:,2)-A(:,8).*A(:,4).*A(:,5).*A(:,2).^2+...
                A(:,9).*A(:,3).*A(:,5).*A(:,2).^2+A(:,10).*A(:,3).*A(:,6).*A(:,2).^2-A(:,8).*A(:,10).*A(:,3)+A(:,8).*A(:,10).*A(:,5)-...
                A(:,8).*A(:,5).*A(:,2)+A(:,10).*A(:,3).*A(:,2);
            res(:,7) = A(:,7).*A(:,10).*A(:,4).*A(:,6)+A(:,7).*A(:,4).*A(:,5).*A(:,1)-A(:,8).*A(:,9).*A(:,4).*A(:,6)+...
                A(:,8).*A(:,4).*A(:,6).*A(:,1)-A(:,9).*A(:,3).*A(:,6).*A(:,1)-A(:,10).*A(:,4).*A(:,6).*A(:,1)+A(:,3).*A(:,6).*A(:,1).^2-...
                A(:,4).*A(:,5).*A(:,1).^2+A(:,7).*A(:,4)-A(:,9).*A(:,6)-A(:,4).*A(:,1)+A(:,6).*A(:,1);
            res(:,8) = -A(:,7).*A(:,10).*A(:,3).*A(:,6)-A(:,7).*A(:,10).*A(:,4).*A(:,5)-A(:,7).*A(:,3).*A(:,5).*A(:,1)+...
                A(:,7).*A(:,4).*A(:,5).*A(:,2)+A(:,8).*A(:,9).*A(:,3).*A(:,6)+A(:,8).*A(:,9).*A(:,4).*A(:,5)-...
                A(:,8).*A(:,3).*A(:,6).*A(:,1)+A(:,8).*A(:,4).*A(:,6).*A(:,2)+A(:,9).*A(:,3).*A(:,5).*A(:,1)-...
                A(:,9).*A(:,3).*A(:,6).*A(:,2)+A(:,10).*A(:,4).*A(:,5).*A(:,1)-A(:,10).*A(:,4).*A(:,6).*A(:,2)+...
                2.*A(:,3).*A(:,6).*A(:,1).*A(:,2)-2.*A(:,4).*A(:,5).*A(:,1).*A(:,2)-A(:,7).*A(:,3)+A(:,8).*A(:,4)+...
                A(:,9).*A(:,5)-A(:,10).*A(:,6)+A(:,3).*A(:,1)-A(:,4).*A(:,2)-A(:,5).*A(:,1)+A(:,6).*A(:,2);
            res(:,9) = A(:,7).*A(:,10).*A(:,3).*A(:,5)-A(:,7).*A(:,3).*A(:,5).*A(:,2)-A(:,8).*A(:,9).*A(:,3).*A(:,5)-...
                A(:,8).*A(:,3).*A(:,6).*A(:,2)+A(:,9).*A(:,3).*A(:,5).*A(:,2)+A(:,10).*A(:,4).*A(:,5).*A(:,2)+...
                A(:,3).*A(:,6).*A(:,2).^2-A(:,4).*A(:,5).*A(:,2).^2-A(:,8).*A(:,3)+A(:,10).*A(:,5)+A(:,3).*A(:,2)-A(:,5).*A(:,2);
        end
    end
end
