function N = samppm(p,meth,sf)
% function N = samppm(p)
% 12 May 2004 - David S. Ciochetto
% 
% A function to take a pressure signal (p) and calculate the average number of samples per m (N).
% N is the profile average.
%
% It can apply several methods selected by meth
% 	meth = 1: Calculate the average based on integer spacing [ p = 3, 4, 5, 6, 7, 8, ...]
% 	meth = 2: Calculate the average based on a moving average with a width of one meter
%		meth = 3: Calculate the average simply bases on the max range and the sampling frequency (sf) which
%							must be supplied (default 8 Hz).
%
% The default is method 1.
% 
% NOTES:
%		It can work with profiles with NaN values.
%		It can not work with Matricies of p data since each column may be different.
%		Assumes that p is monotonically increasing S.T. max(p) is the last non NaN value.
%		(Maybe I was trying to say that there should not be gappy data also. I have not checked
%		what it does for data dropouts. Do not need it so keep it in mind for the future use.)
%
%		For CTD data from the Scotian Shelf, the numbers from the two methods match when rounded.
%
%		Case 1 will return a NaN if there is fewer than 1 m of data.
%
% 	If any case fails to return a value, case 3 is used to make a rough estimate.
%		The user must supply the sampling rate and the estimate will assume that the data taken with no gaps,
% 	i.e. consecuitively. This method is less accurate since it is simply # samples/max range. If there are
%		few data and the instrument bobs, this could be a problem, it will not average over changes in vertical
%		velocity of the sensor package.
% 
% USES:
%

% Make p positive values
p = abs(p);

if ~exist('meth')
	meth = 1;
end

if ~exist('sf')
	sf = 8; %[Hz]
end

switch meth
case 1	% Default: Method 1. Calculate the average based on regular integer spacing.
	% Set up a reference grid of unit intervals
	% I want the average per meter so I will chop off data that does not occupy a full m bin at the top and bottom
	grd = [ceil(min(p)): floor(max(p))]';
	
	S = [];	% samples per meter for each meter bin
	for i = 1:max(grd)-min(grd)
		idx = p > grd(i) & p <= grd(i+1);
		idx = idx - 0;
		S = [S; sum(idx)];
	end
	
case 2
	S = []; i = 1;
	while p(i)+1 <= max(p)
		idx  = p > p(i) & p <= p(i)+1;
		idx = idx - 0;
		S = [S; sum(idx)];
		i = i + 1;
	end
	
case 3		% Simple overall average (Look into fixing this better with a runlength and sf and W? (see 2004-05.txt)
	S = sum(~isnan(p))/(max(p)-min(p));
end
	
% Now the average value of samples per meter for the profile
N = mean(S);	% N = nanmean(S);

% In case there are short casts with less than one m data and method 1 is used.... return something if there is data.
if isnan(N) & sum(isnan(p)) < length(p)
	N = sum(~isnan(p))/(max(p)-min(p));
end
