function [m,s,n] = mstdgap(x, sq);
% function [m,s,n] = mstdgap(x, 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.
%
%    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.)

squeeze_flag = 0;
if nargin == 2,
   if sq(1) == 's',
      squeeze_flag = 1;
   end
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);
if ~isempty(ibad),
   x(ibad) = zeros(size(ibad));
end
sx = sum(x);
m(isg) = sx(isg) ./ n(isg);

xdm = x - m(ones(dimx(1),1),:);
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 squeeze_flag == 1,
      m = squeeze(m);
      s = squeeze(s);
      n = squeeze(n);
   end

end