function pcut = addEvalVariableCuts(p)

pcut = p;
if ~isempty(p.evalMap)
   % temp = p.F_struc(1:p.K.f,:);
   % p.F_struc(1:p.K.f,:) = [];
    %pcut = p;
    pcut = emptyNumericalModel;
    for i = 1:length(p.evalMap)
        y = p.evalVariables(i);
        x = p.evalMap{i}.variableIndex;
        xL = p.lb(x);
        xU = p.ub(x);
        
        % Generate a convex hull polytope
        if xL<xU
            if ~isempty(p.evalMap{i}.properties.convexhull)
                % A convex hull generator function is available!
                K = [];
                try
                    if isfield(p.evalMap{i},'oldhull') & isequal(p.evalMap{i}.oldhull.xL,xL) & isequal(p.evalMap{i}.oldhull.xU,xU)
                        Ax = p.evalMap{i}.oldhull.Ax;
                        Ay = p.evalMap{i}.oldhull.Ay;
                        b = p.evalMap{i}.oldhull.b;
                        K = p.evalMap{i}.oldhull.K;
                    else
                        [Ax,Ay,b,K]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
                        p.evalMap{i}.oldhull.xL = xL;
                        p.evalMap{i}.oldhull.xU = xU;
                        p.evalMap{i}.oldhull.Ax = Ax;
                        p.evalMap{i}.oldhull.Ay = Ay;
                        p.evalMap{i}.oldhull.b = b;
                        p.evalMap{i}.oldhull.K = K;
                    end
                catch
                    if isfield(p.evalMap{i},'oldhull') & isequal(p.evalMap{i}.oldhull.xL,xL) & isequal(p.evalMap{i}.oldhull.xU,xU)
                        Ax = p.evalMap{i}.oldhull.Ax;
                        Ay = p.evalMap{i}.oldhull.Ay;
                        b = p.evalMap{i}.oldhull.b;
                        K = p.evalMap{i}.oldhull.K;
                    else
                        [Ax,Ay,b]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
                        p.evalMap{i}.oldhull.xL = xL;
                        p.evalMap{i}.oldhull.xU = xU;
                        p.evalMap{i}.oldhull.Ax = Ax;
                        p.evalMap{i}.oldhull.Ay = Ay;
                        p.evalMap{i}.oldhull.b = b;
                        p.evalMap{i}.oldhull.K = K;
                    end
                end
            else
                if length(xL)>1
                    disp(['The ' p.evalMap{i}.fcn ' operator does not have a convex hull operator'])
                    disp('This is required for multi-input single output operators');
                    disp('Sampling approximation does not work in this case.');
                    error('Missing convex hull operator');
                end
                % sample function
                z = linspace(xL,xU,100);
                
                if isequal(p.evalMap{i}.fcn,'power_internal2')
                    % Special code for automatically converting sigmonial
                    % terms to be solvable with bmibnb
                    fz = feval(p.evalMap{i}.fcn,z,p.evalMap{i}.arg{2});
                    
                else
                    arg = p.evalMap{i}.arg;
                    arg{1} = z;
                    fz = feval(p.evalMap{i}.fcn,arg{1:end-1});
                    % end
                    [minval,minpos] = min(fz);
                    [maxval,maxpos] = max(fz);
                    xtestmin = linspace(z(max([1 minpos-5])),z(min([100 minpos+5])),100);
                    xtestmax = linspace(z(max([1 maxpos-5])),z(min([100 maxpos+5])),100);
                    arg{1} = xtestmin;
                    fz1 = feval(p.evalMap{i}.fcn,arg{1:end-1});
                    arg{1} = xtestmax;
                    fz2 = feval(p.evalMap{i}.fcn,arg{1:end-1});
                    z = [z(:);xtestmin(:);xtestmax(:)];
                    fz = [fz(:);fz1(:);fz2(:)];
                    [z,sorter] = sort(z);
                    fz = fz(sorter);
                    [z,ii,jj]=unique(z);
                    fz = fz(ii);
                end
                
                %                    fz = feval(p.evalMap{i}.fcn,z);
                %   end
                % create 4 bounding planes
                % f(z) < k1*(x-XL) + f(xL)
                % f(z) > k2*(x-XL) + f(xL)
                % f(z) < k3*(x-XU) + f(xU)
                % f(z) > k4*(x-XU) + f(xU)
                k1 = max((fz(2:end)-fz(1))./(z(2:end)-xL))+1e-12;
                k2 = min((fz(2:end)-fz(1))./(z(2:end)-xL))-1e-12;
                k3 = min((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))+1e-12;
                k4 = max((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))-1e-12;
                Ax = [-k1;k2;-k3;k4];
                Ay = [1;-1;1;-1];
                b =  [k1*(-z(1)) + fz(1);-(k2*(-z(1)) + fz(1));k3*(-z(end)) + fz(end);-(k4*(-z(end)) + fz(end))];
                K = [];
            end
            if ~isempty(b)
                if isempty(K)
                    % Compatibility with old code
                    K.f = 0;
                    K.l = length(b);
                end
                F_structemp = zeros(size(b,1),length(p.c)+1);
                F_structemp(:,1+y) = -Ay;
                F_structemp(:,1+x) = -Ax;
                F_structemp(:,1) = b;
                localModel = createNumericalModel(F_structemp,K);
                pcut = mergeNumericalModels(pcut,localModel);
              %  pcut.F_struc = [F_structemp; pcut.F_struc];
              %  pcut.K.l = pcut.K.l + K.l;%length(b);
              %  pcut.K.f = pcut.K.f + K.f;%length(b);
            end
        end
    end
    
    pcut = mergeNumericalModels(p,pcut);
    %pcut.F_struc = [temp;pcut.F_struc];
end

