function [m, s, n, med] = mstdgap(x, varargin);
% function [m, s, n, med] = 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.
% If a 4th output argument is requested, it will be the median.
%
%    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.)
% (00/03/26, added median output.)
% (2001/10/30, changed so m = NaN + i * NaN if x is complex and n is zero.)
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);
if any(~isreal(x(:))),
   m = m + i * NaN;
end

% Reduce to 2 dimensions by combining any dimensions higher than 2.
m = m(1,:);
x = x(:,:);
s = m;
if nargout > 3, med = m; end

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 nargout > 3,
   %xz(ibad) = 1e40 * cos(pi * (1:length(ibad)));
   %med = median(xz);
   %ii = find(med > 1e39 | med < -1e39);
   %aa(ii) = NaN;
   % The commented-out statements are an attempt to keep
   % the calculation completely vectorized.  Unfortunately,
   % it yields the same result as the loop below only if the
   % original row dimension and the number of non-NaN rows
   % in each column are both odd.
   % I think it could still be vectorized by using sort(xz),
   % then calculating the middle 1 or two (two indices in
   % either case) and averaging the values at those indices.
   [nr, nc] = size(x);
   if 0 %% columnwise
      for icol = 1:nc,
         if n(icol) == 0,
            med(icol) = NaN;
         else
            xc = x(:, icol);
            ii = find(~isnan(xc));
            med(icol) = median(xc(ii));
         end
      end
   else %% second attempt to vectorize
      xs = sort(x);
      ic_offset = 1 + (0:(nc-1))*nr;
      i0 = floor((n-1)/2) + ic_offset;
      i1 = floor(n/2) + ic_offset;
      ii = find(n >= 1);
      med(ii) = 0.5 * (xs(i0(ii)) + xs(i1(ii))).';
   end
end

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

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

   if Nshift,
      Nback = ndims_orig - Nshift;
      m = shiftdim(m, Nback);
      s = shiftdim(s, Nback);
      n = shiftdim(n, Nback);
      if nargout > 3,
         med = shiftdim(med, Nback);
      end
   end



   if squeeze_flag == 1,
      m = squeeze(m);
      s = squeeze(s);
      n = squeeze(n);
      if nargout > 3,
         med = squeeze(med);
      end
   end

end

