function [Ydm,mYvec]=demean(Y,dim)


if and(size(Y,3)<=1,dim<=2)
    %dims2D=[2,1];
    mY=nanmean(Y,dim);   mYvec=mY(:);
    newdims=ones(1,2);
    newdims(dim)=size(Y,dim);
    
    Ydm=Y-repmat(mY,newdims);
else
    error('Demean function is not able to process objects with more than 2 dimensions')
end