function v = del2(f,varargin)
%MATLAB Code Generation Library Function

%   Copyright 1984-2014 The MathWorks, Inc.
%#codegen

eml_invariant(nargin > 0, ...
    'MATLAB:minrhs');
eml_invariant(isa(f,'float'), ...
    'Coder:toolbox:del2_2', ...
    'IfNotConst','Fail');
ONE = ones(coder.internal.indexIntClass);
TWO = cast(2,coder.internal.indexIntClass);
onescalarh = nargin == 2 && ...
    eml_is_const(isscalar(varargin{1})) && isscalar(varargin{1});
eml_invariant(nargin == 1 || onescalarh || nargin-1 == ndims(f) || isvector(f), ...
    'Coder:MATLAB:del2_InvalidInputs');
if isvector(f)
    % Streamlined code for vector case. The other reason for treating this
    % case separately is to allow a single h vector input, which is
    % otherwise not supported.
    if size(f,1) >= 3
        if nargin == 1
            v = eml_del2(f,ONE);
        else
            v = eml_del2(f,ONE,varargin{1});
        end
    elseif size(f,2) >= 3
        if nargin == 1
            v = eml_del2(f,TWO);
        else
            v = eml_del2(f,TWO,varargin{1});
        end
    else
        v = eml_expand(coder.internal.scalarEg(f),size(f));
    end
elseif nargin == 1
    v = eml_del2(f,ONE);
    for k = coder.unroll(TWO:coder.internal.ndims(f))
        if size(f,k) >= 3
            v = v + eml_del2(f,k);
        end
    end
elseif onescalarh
    % Special case for two input arguments where the second is a scalar.
    % The second input is the scalar h value for all dimensions.
    v = eml_del2(f,ONE,varargin{1});
    for k = coder.unroll(TWO:coder.internal.ndims(f))
        if size(f,k) >= 3
            v = v + eml_del2(f,k,varargin{1});
        end
    end
else % nargin-1 == ndims(f)
    v = eml_del2(f,ONE,varargin{2});
    if size(f,2) >= 3
        v = v + eml_del2(f,TWO,varargin{1});
    end
    for k = coder.unroll(cast(3,coder.internal.indexIntClass):nargin-1)
        if size(f,k) >= 3
            v = v + eml_del2(f,k,varargin{k});
        end
    end
end
v = v/ndims(f);

%--------------------------------------------------------------------------

function g = eml_del2(f,dim,h)
coder.internal.prefer_const(dim);
vlen = cast(size(f,dim),coder.internal.indexIntClass);
vlenm1 = coder.internal.indexMinus(vlen,1);
vlenm2 = coder.internal.indexMinus(vlen,2);
g = coder.nullcopy(f);
scalarh = nargin == 3 && eml_is_const(isscalar(h)) && isscalar(h);
eml_invariant(nargin < 3 || scalarh || numel(h) == vlen, ...
    'Coder:MATLAB:del2_InvalidInputs');
if vlen < 3
    g(:) = 0;
    return
end
if nargin == 3 && ~scalarh
    for k = 1:vlenm1
        h(k) = h(coder.internal.indexPlus(k,1)) - h(k);
    end
end
vstride = coder.internal.prodsize(f,'below',dim);
vspread = coder.internal.indexTimes(coder.internal.indexMinus(vlen,1),vstride);
npages  = coder.internal.prodsize(f,'above',dim);
i2 = zeros(coder.internal.indexIntClass);
for i = 1:npages
    i1 = i2;
    i2 = coder.internal.indexPlus(i2,vspread);
    for j = 1:vstride
        i1 = coder.internal.indexPlus(i1,1);
        i2 = coder.internal.indexPlus(i2,1);
        % Vector is f(i1:vstride:i2)
        if nargin < 3 % H not supplied for this dimension.
            % Take centered second differences on interior points to compute g/2
            for k = 1:coder.internal.indexMinus(vlen,2)
                % If f and g were simple vectors, this would look like:
                % g(k+1) = ((f(k+2)-f(k+1)) - (f(k+1)-f(k)))/2;
                g(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))) = ( ...
                    (f(coder.internal.indexPlus(i1,coder.internal.indexTimes(coder.internal.indexPlus(k,1),vstride))) - ...
                    f(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride)))) -  ...
                    (f(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))) - ...
                    f(coder.internal.indexPlus(i1,coder.internal.indexTimes(coder.internal.indexMinus(k,1),vstride))))) / 2;
            end
            % Linearly extrapolate second differences from interior
            if vlen > 3
                g(i1) = g(coder.internal.indexPlus(i1,vstride))*2 -  ...
                    g(coder.internal.indexPlus(i1,coder.internal.indexTimes(vstride,2)));
                g(i2) = -g(coder.internal.indexMinus(i2,coder.internal.indexTimes(vstride,2))) +  ...
                    g(coder.internal.indexMinus(i2,vstride))*2;
            else % vlen == 3
                g(i1) = g(coder.internal.indexPlus(i1,vstride));
                g(i2) = g(coder.internal.indexPlus(i1,vstride));
            end
        elseif eml_is_const(isscalar(h)) && isscalar(h) % H is a scalar.
            % Take centered second differences on interior points to compute g/2
            for k = 1:coder.internal.indexMinus(vlen,2)
                % If f and g were simple vectors, this would look like:
                % g(k+1) = ((f(k+2)-f(k+1))/h - (f(k+1)-f(k))/h)/(2*h);
                g(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))) = ( ...
                    (f(coder.internal.indexPlus(i1,coder.internal.indexTimes(coder.internal.indexPlus(k,1),vstride))) - ...
                    f(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))))/h -  ...
                    (f(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))) - ...
                    f(coder.internal.indexPlus(i1,coder.internal.indexTimes(coder.internal.indexMinus(k,1),vstride))))/h)/(2*h);
            end
            % Linearly extrapolate second differences from interior
            if vlen > 3
                g(i1) = g(coder.internal.indexPlus(i1,vstride))*2 -  ...
                    g(coder.internal.indexPlus(i1,coder.internal.indexTimes(vstride,2)));
                g(i2) = -g(coder.internal.indexMinus(i2,coder.internal.indexTimes(vstride,2))) +  ...
                    g(coder.internal.indexMinus(i2,vstride))*2;
            else % vlen == 3
                g(i1) = g(coder.internal.indexPlus(i1,vstride));
                g(i2) = g(coder.internal.indexPlus(i1,vstride));
            end
        else % H is an abscissa vector.
            % Take centered second differences on interior points to compute g/2
            for k = 1:coder.internal.indexMinus(vlen,2)
                % If f and g were simple vectors, this would look like:
                % g(k+1) = ((f(k+2)-f(k+1))/h(k) - (f(k+1)-f(k))/h(k-1))/(h(k+1) + h(k));
                g(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))) = ( ...
                    (f(coder.internal.indexPlus(i1,coder.internal.indexTimes(coder.internal.indexPlus(k,1),vstride))) - ...
                    f(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))))/h(coder.internal.indexPlus(k,1)) -  ...
                    (f(coder.internal.indexPlus(i1,coder.internal.indexTimes(k,vstride))) - ...
                    f(coder.internal.indexPlus(i1,coder.internal.indexTimes(coder.internal.indexMinus(k,1),vstride))))/h(k) ...
                    )/(h(coder.internal.indexPlus(k,1)) + h(k));
            end
            % Linearly extrapolate second differences from interior
            if vlen > 3
                g(i1) = g(coder.internal.indexPlus(i1,vstride))*(h(1)+h(2))/h(2) -  ...
                    g(coder.internal.indexPlus(i1,coder.internal.indexTimes(vstride,2)))*(h(1))/h(2);
                g(i2) = -g(coder.internal.indexMinus(i2,coder.internal.indexTimes(vstride,2)))*(h(vlenm1))/h(vlenm2) +  ...
                    g(coder.internal.indexMinus(i2,vstride))*(h(vlenm1)+h(vlenm2))/h(vlenm2);
            else % vlen == 3
                g(i1) = g(coder.internal.indexPlus(i1,vstride));
                g(i2) = g(coder.internal.indexPlus(i1,vstride));
            end
        end
    end
end

%--------------------------------------------------------------------------
