function [m,s,n] = mstdgap(x, varargin);
% function [m,s,n] = mstdgap(x, N, squeeze);
% Mean, Standard Deviation, and Number of good points
% in an array with NaN-filled gaps.  Calculations are
% done column-wise (or on the first non-singleton dimension)
% as for the normal mean and std.
%
%    N: optional argument; dimension along which to calculate
%
%    squeeze: optional string argument; if 'squeeze', then
%           the squeeze operator is applied to the outputs.
%
% EF
% (Modified 98/07/15 to handle arrays with more than 2 dimensions.)
% (99/11/30, added optional argument N.)
%

squeeze_flag = 0;
Nshift = 0;
ndims_orig = ndims(x);
if nargin > 1,
   for ii = 1:length(varargin)
      arg = varargin{ii};
      if strcmp(arg, 'squeeze'),
         squeeze_flag = 1;
      elseif arg > 0 & arg <= ndims_orig
         Nshift = arg-1;
      else
         error('invalid optional argument')
      end
   end
end

if Nshift, x = shiftdim(x, Nshift); end

% Strip off leading singleton dimensions.
[x, nshift] = shiftdim(x);


dimx = size(x);
dimout = dimx;
dimout(1) = 1;
ndims = length(dimx);

m = NaN * zeros(dimout);

% Reduce to 2 dimensions by combining any dimensions higher than 2.
m = m(1,:);
x = x(:,:);
s = m;

goodmask = ~isnan(x);
ibad = find(~goodmask);
n  = sum(goodmask);
isg = find(n ~= 0);
xz = x;
if ~isempty(ibad),
   xz(ibad) = zeros(size(ibad));
end
sx = sum(xz);
m(isg) = sx(isg) ./ n(isg);

xdm = x - m(ones(dimx(1),1),:);
if ~isempty(ibad),
   xdm(ibad) = zeros(size(ibad));
end
sxs = sum(xdm .* conj(xdm));

s(isg) = sqrt(sxs (isg) ./ n(isg));

if ndims > 2,
   % Restore the dimensions higher than 2.
   m = reshape(m, dimout);
   s = reshape(s, dimout);
   n = reshape(n, dimout);

   if nshift > 0 & squeeze_flag == 0,
      % Replace the leading singleton dimensions.
      m = shiftdim(m, -nshift);
      s = shiftdim(s, -nshift);
      n = shiftdim(n, -nshift);
   end

   if Nshift,
      Nback = ndims_orig - Nshift;
      m = shiftdim(m, Nback);
      s = shiftdim(s, Nback);
      n = shiftdim(n, Nback);
   end



   if squeeze_flag == 1,
      m = squeeze(m);
      s = squeeze(s);
      n = squeeze(n);
   end

end