function [x_bar, st_dev, n] = smartmean(m,dim);
%
% function [x_bar, st_dev, n] = smartmean(m,dim);
%
% Calculates a smart nan mean/std !
% you'll be able to make nanmeans/std from n-dimensional arrays 
% Inputs:
%			m = input matrix
%        dim = dimension over which you want to average.
%              must be <= number of m_dimensions 
%	        (dim is optional, if not specified the program works
%                over the first dimension).
% Outputs:
%
%        x_bar  = de-nan mean (m_dimensions-1)
%        st_dev = Stand. Deviation (normalized) of x_bar (m_dimensions-1)
%        n = number of non-nan observations used to obtain x_bar (m_dimensions-1)
%
%      *obs: if any of the results show NaN's , it means you did not have enough good
%            data (at least 2) in your original matrix (over the chosen dim).
%
% Mauricio , ADL 02/10/1999

if isempty(m) % Check for empty input.
    error('input matrix is empty !');
end

if (nargin<1 | nargin>2),
   error('Requires at least one input arguments and a max. of two!'); 
end

if nargout>3 , 
   error('Too many output arguments.'); 
end

if nargin==1, 
   dim=1; 
end

if dim > length(size(m)),
   error('You cannot set "dim" > than max number of dimensions in your original matrix , you dummy !');
end
   
  


ind_mat=ind2sub(size(m),~isnan(m)); % N-Dim Mat with zeros where nan were present
n=sum(ind_mat,dim); % Sum of the good index in the n-dim here, the cycles (time)
m(find(isnan(m))) = 0;
sum_m = sum(m,dim);
sum_m_2 = sum_m.^2;
%
m2 = m.^2;
sum_m2 = sum(m2,dim);

%%%% at least 2 valid data to be able to calculate the variance
n_control = find(n == 0 | n == 1);
if(~isempty(n_control)),
   n(n_control) = nan;
end   
%%%%

x_bar = squeeze( sum_m./n ); % mean
st_dev = squeeze( sqrt( (sum_m2 - sum_m_2./n)./(n-1)) ); % std. Dev. normalized to sample variance
n = squeeze(n); % number of valid observations; no nan's !



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% THE END %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






