B=summ(A,varargin)

function B=summ(A,varargin)
    %SUMM Sum of consecutive pairs in an array, analogous to DIFF.
    %   SUMM(X), for a vector X, is [X(2)+X(1)  X(3)+X(2) ... X(n)+X(n-1)].
    %   SUMM(X), for a matrix X, is the matrix of row sums,
    %      [X(2:n,:) + X(1:n-1,:)].
    %   SUMM(X), for an N-D array X, is the sum along the first
    %      non-singleton dimension of X.
    %   SUMM(X,N) is the N-th order sum along the first non-singleton 
    %      dimension (denote it by DIM). If N >= size(X,DIM), SUMM takes 
    %      successive sums along the next non-singleton dimension.
    %   SUMM(X,N,DIM) is the Nth sum function along dimension DIM. 
    %      If N >= size(X,DIM), SUMM returns an empty array.
    %
%   Examples:
%
%      If X = [3 7 5
%              0 9 2]
%      then summ(X,1,1) is [3 16 7], summ(X,1,2) is [10 12
%                                                     9 11],
%      summ(X,2,2) is the 2nd order sum along the dimension 2, and
%      summ(X,3,2) is the empty matrix.
    
    % set default input arguments
    niter=1;
    dims=size(A);
    dim=find(dims>1,1,'first');
    if nargin==1 % no varargs
        if isempty(dim)
            B=[];
            return
        end
    elseif nargin==2 % N is given
        if isempty(dim)
            B=[];
            return
        end
        niter=varargin{1};
    elseif nargin==3 % N and dim are given
        if (isempty(dim))
            newdims=[1 1];
            newdims(varargin{2})=0;
            B=zeros(newdims);
            return
        elseif (varargin{2}>numel(dims))
            B=zeros([dims 0]);
            return
        end
        niter=varargin{1};
        dim=varargin{2};
    end
    
    
    
    for jx=1:niter
        % helpersumm does the actual sum
        temp=helpersumm(A,dims,dim);
        dims=size(temp);
        A=temp;
    end
    
    B=A;
end

function B=helpersumm(A,dims,dim)
    % Use a helper subfunction for actual summ to allow for iteration
    % Example of engine:
    % A=rand(1,10);
    % B=A(1:end-1)+A(2:end);
    %
    % Just need to determine correct index to use
    
    idx1=cell(numel(dims),1);
    idx1(:)={':'};
    idx2=idx1;
    
    L1=true(1, dims(dim));
    L2=L1;
    L1(dims(dim))=false;
    L2(1)=false;
    
    idx1{dim} = L1;
    idx2{dim} = L2;

    B=A(idx1{:})+A(idx2{:});
end