function [cbin,dat,n,sd] = pbin02(p,D,wbin,mm,chk)
%
% function [cbin,dat,n,sd] = pbin02(p,D,wbin,mm,chk)
% 24 August 2004 -by- David S. Ciochetto (Oceanograpy at Dalhousie University)
% 
% Modified form of pbin.m. A routine to take a pressure vector (p) and corresponding data (D)
% and a bin width (wbin) and take the average of the data in the given bin sizes. It returns
% the bin centers (cbin) and the binned data (dat). I will use NaNmean such that data may include
% NaN values. p should be from the CTD so it is typically in [db] which for coastal applications
% is roughly [m]. No depth corrections are made.
%
% It will calculate the number of data used to make each bin mean (n) and return that if desired.
% It will also use nanstd to return the standard deviation of the non-nan data in each bin.
%	
% mm is a flag to select the mean = 0 or the median = 1. The default is set to the mean. The
% median is used when there are a lot of outliers makeing the probability distribution non-Gaussian.
%
% Alternately, if wbin is a vector of values (not a single value) corresponding to p, then these
% values will be used as the ends of each bin. Note that one less dat (average) value will be
% returned than length(wbin).
%
%	Note: Given wbin, it is set up to use equal size pressure bins. Data should be in columns.
%				D can be a matrix and it will return a row of values. p can be a vector or a matrix of the
%				same size as D with corresponding columns.
% 
%				23 June 2005 - I had trouble sending it a matrix of values for D. It worked when the
%				data were sent as individual vectors. No time to fix this now. The call is in ADVat10m.m
% USES: 

if ~exist('mm')
	mm=0;
end

% Note that this was set up to work with p as a (+) number...
% So check if p values are positive or negative and make positive.
neg = 'n';
idx = p<0;
if sum(idx)>0.1*length(p)								% If over 10% of the values are (-) then invert it.
	neg = 'y';
	P = -p;
else
	P = p;
end

% Check if the bin widths (wbin) are supplied or they must be calculated.
% If calculated, include all of the data and set the bins to be equal
if length(wbin) == 1
	% (set to include all of the data in integer bins)
	mp = floor(max(P));											% max of p rounded DOWN to an integer value
	mxp = mp + ceil((max(P)-mp)/wbin)*wbin;	% adjust mp to the proper next p bin above the integer value 
	clear mp;
	mp = floor(min(P));		% Set the minimum value then check to make sure it is not less than 0
	if rem((mxp-mp),wbin)~=0
		mnp = mp - wbin + rem((mxp-mp),wbin);
	else
		mnp = mp;
	end
	if mnp < 0						% if it is less than zero, tack the number that it is on to the max p to have equal bins
		mxp = mxp - mnp; 
		mnp = 0;
	end
	pin = [mnp:wbin:mxp]';					% pressure limit index (ends not centers)
	clear mp mxp mnp

	cbin = pin(1:end-1)+wbin/2;			% Location of the bin centers
	if neg == 'y'
		cbin = -cbin;
	end
else	% bin edges are supplied.
	if sum(wbin) < 0	% Make sure that the wbin values are positive to match the pressure we changed to positive
		wbin = -wbin;
	end
	pin = wbin;

	cbin = pin(1:end-1)+diff(wbin)./2;			% Location of the bin centers
	if neg == 'y'
		cbin = -cbin;
	end
end


dat = []; n=[]; sd=[];
if size(P,2) == 1	% If the pressure vector is a vector not a matrix
	for pp = 1:length(pin)-1
	 	% There is a problem here for all three statistics. If there is only one data point in the bin, then it
 		% takes the stat for all the columns. There is no way to limit its dimension in the function call.
 		if size(D(P>=pin(pp) & P<pin(pp+1),:),1) == 1
 			if mm == 0
				dat=[dat; mean(D(P>=pin(pp) & P<pin(pp+1),:),1)];% Works but without NaN consideration, can force columns
			elseif mm == 1
				dat = [dat; median(D(P>=pin(pp) & P<pin(pp+1),:),1)];
			else
				disp('Choose mm = 1 for median and 0 (default) for mean')
			end
 			n   = [n; sum(~isnan(D(P>=pin(pp) & P<pin(pp+1),:)),1)];
 			sd  = [sd; std(D(P>=pin(pp) & P<pin(pp+1),:),0,1)];
		else
 			if mm == 0
				dat = [dat; nanmean(D(P>=pin(pp) & P<pin(pp+1),:))];		% nanmean works on columns by default
 			elseif mm == 1
				dat = [dat; nanmedian(D(P>=pin(pp) & P<pin(pp+1),:))];
			else
				disp('Choose mm = 1 for median and 0 (default) for mean')
			end
			n=[n; sum(~isnan(D(P>=pin(pp) & P<pin(pp+1),:)))];	% number of data used for the mean, ok for matricies
 			sd  = [sd; nanstd(D(P>=pin(pp) & P<pin(pp+1),:))];			% std for the samples in each bin
		end
	end
else
	% The above works when P is simply one vector... there are indexing issues when it is not
	% It does not like the indexing for when P is a matrix too
	for pp = 1:length(pin)-1
		IND = P>=pin(pp) & P<pin(pp+1); IND=+IND;	% Set up the index matrix and make it non-logical
		D1=D; D1(IND==0)=NaN;		% Set up a dummy matrix to ditch unwanted data
		if size(D1,1) == 1
			if mm == 0
				dat = [dat; mean(D1,1)];
			elseif mm == 1
				dat = [dat; median(D1,1)];
			else
				disp('Choose mm = 1 for median and 0 (default) for mean')
			end
			n   = [n; sum(~isnan(D1),1)];
			sd  = [sd; std(D1,0,1)];
		else
			if mm == 0
				dat = [dat; nanmean(D1)];
			elseif mm == 1
				dat = [dat; nanmedian(D1)];
			else
				disp('Choose mm = 1 for median and 0 (default) for mean')
			end
			n   = [n; sum(~isnan(D1))];
			sd  = [sd; nanstd(D1)];
		end
		clear D1 IND
	end
end

if exist('chk') & lower(chk)=='y'
	% Check it
	figure
	plot(D,p,'-k.',dat,cbin,'r-o'); grid on; xlabel('data'); ylabel('p [m]');
	title(['Testing pbin on data (blk) for ' num2str(wbin) ' m bins (red)'])
end

% TEST data (vector P not matrix P)
% data = repmat([1:3 1:5 1:10]',1,8);
% bins = [1 2 4 8 12]';
% z = [0.5:0.5:9]'; % Problem when I have a matrix of z values
% OK - then make some of the values NaN with the array editor and check - OK tested!!
