% Input:
% S: function handle (support function)
% n: vector
% gamma : double
% options: 'fmincon/fminbnd/fminsearch/fminunc' !!! fminsearch/fminunc
% calculates n*x = gamma

% Ouput:
% f: function handle (support function of the intersection of S with n*x <=
% or == gamma)

% Description:
% Computes the intersection of a set S with a Hyperplane defined by
% n*x <= gamma or n*x = gamma for fminsearch und fminunc

function [f] = directintersection(S,n,gamma,option)

if strcmp(option{1},'fmincon')
    f = @(l) intersection(S,n,gamma,l);
elseif strcmp(option{1},'fminbnd')  % golden section search over the intervall 0 pi
    f = @(l) fminbnd_handle(S,n,gamma,l);
elseif strcmp(option{1},'fminsearch')   % searches on entire R ! -> different result then the others !
                                        % => calculates the intersection
                                        % with a hyperplane instead of a
                                        % halfspace
    f = @(l) fminsearch_handle(S,n,gamma,l);
    
elseif strcmp(option{1},'fminunc')
    f = @(l) fminunc_handle(S,n,gamma,l);
elseif strcmp(option{1},'RayAlgo_inequality')
    f = @(l) fRay_handle_inequality(S,n,gamma,l);
elseif strcmp(option{1},'RayAlgo_equality')
    f = @(l) fRay_handle_equality(S,n,gamma,l);
end;

function [fval] = intersection(S,n,gamma,l)
fun = @(x) (S(l-x*n) + x*gamma);
[~,fval] = fmincon(fun,0,-1,0,[],[],[],[],[],optimset('Display','off','Algorithm','active-set'));

function [fval] = fminbnd_handle(S,n,gamma,l)
fun = @(x) (S(l-x*n) + x*gamma);
[~,fval] = fminbnd(fun,0,pi);

function [fval] = fminsearch_handle(S,n,gamma,l)    % defining pi/2 as an estimate start value
fun = @(x) (S(l-x*n) + x*gamma);
[~,fval] = fminsearch(fun,pi/2);

function [fval] = fminunc_handle(S,n,gamma,l)
fun = @(x) (S(l-x*n) + x*gamma);
[~,fval] = fminunc(fun,zeros(size(l,1),1));

function [fval] = fRay_handle_inequality(S,n,gamma,l)
fun = @(x) (S(l-x*n) + x*gamma);
fval = RayAlgo(0,10000000,0.000001,fun);

function [fval] = fRay_handle_equality(S,n,gamma,l)
fun = @(x) (S(l-x*n) + x*gamma);
fval = RayAlgo(-10000000,10000000,0.000001,fun);