function [MSL, Dstd] = trimBins(rawdata,trimFile,ADCP_offset,MSL,Dstd)

%function [MSL, Dstd] = trimBins(rawdata,trimFile,ADCP_offset,MSL,Dstd)
%Use RDsurface.m, a Matlab function that runs the RDI surface program,
%to create new data file where the bins all fall below mean sea level plus 
%the tidal range.
%
%INPUTS:
%	rawdata =  RDI ADCP binary output file (*.000)
%	trimFile = the next file create by trimming the bad ensembles
%	ADCP_offset = is the an attribute of the depth variable, 
%		called xducer_offset_from_bottom in the Netcdf files
%		If not given, the program requests a Netcdf file to find it.
%	MSL = same as output
%	Dstd = same as output
%
%OUTPUTS:
%	MSL = mean sea level based on RDI surface output
%	Dstd = standard of deviation to give an approximate tidal variation

%sub-functions:
%	Netcdf toolbox
%	RDsurface.m
%	DOS version of surface.exe

%tell us what function is running
Mname=mfilename;
disp(' ')
disp([ Mname ' is currently running']);


if nargin < 1, rawdata = '';, end
if nargin < 2, trimFile = '';, end
if nargin < 3, ADCP_offset = '';, end

if isempty(rawdata), rawdata = '*'; end
if isempty(trimFile), trimFile = '*'; end
%ADCP_offset will be taken care of in the RDsurface program if isempty.

%Get the files in case the inputs are empty
% Get ADCP raw data filename.
if any(rawdata == '*')
	[theFile, thePath] = uigetfile(rawdata, 'Select Binary ADCP File:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	rawdata = [thePath theFile];
end

% Get trimmed ADCP data file.
if any(trimFile == '*')
	[theFile, thePath] = uigetfile(trimFile, 'Select Ensemble trimmed ADCP File:');
	if ~any(theFile), return, end
	if thePath(end) ~= filesep, thePath(end+1) = filesep; end
	trimFile = [thePath theFile];
end


%get the mean sea level and tidal range
%[MSL,Dstd] = RDsurface(rawdata,ADCP_offset,progPath,DepthFile);
if nargin < 4
[MSL,Dstd] = RDsurface(rawdata,ADCP_offset);
end

%open the netcdf file to trim and get soem info
f=netcdf(trimFile,'write');
if isempty(f), return, end;
depth=f{'D'};
theVars=var(f);
bad_num = fillval(depth);

B = f('bin'); 
theBinNum = B(:);


if isempty(theBinNum)
    disp(' ## the number of bins not found.')
    close(f)
    return
end

%for files with _FillValue
if ~isempty(bad_num)
   D = autonan(depth,1);
   depth=D(:);
   
else
bad_num=f.VAR_FILL
   
   if ischar(bad_num)==1;
   		bad_num=str2num(bad_num);
   	end
      
   depth=depth(:);
   idgood=find(depth < bad_num & depth > 0);
   depth=depth(idgood);
   clear idgood  
   
end


%find the bins1that fall below the given sea level including tidal fluctuation   
disp(' ')
disp('Finding the good bins');
goodBins=find(depth < MSL+Dstd);

if isempty(goodBins)  
   disp('## No bins were found within the specified depth range');
   close(f)
   return
elseif isequal(goodBins,length(depth));
   disp('No bad bins were found, file is unmodified');
   f{'D'}.water_depth(:) =MSL;
   close(f)
   return
 
end

%Check to make sure that everything is on order
%The goodbins should run from 1(closets to the ADCP) to the depth.
goodBins=sort(goodBins);
%make sure indices are unique
Bd = find(diff(goodBins) == 0);
if any(Bd) 
   goodBins(Bd) = []; 
	disp('Warning: Empty bins were found')   
end  


% compute bin locations
offset=f{'D'}.xducer_offset_from_bottom(:);
bin1=f{'D'}.center_first_bin(:);
binsize=f{'D'}.bin_size(:);
bincnt=length(goodBins);
EndD = (((bincnt-1)*binsize)+bin1) + offset;

%resize the dimension bins if all ok
if depth(goodBins(end)) - EndD > binsize
   disp(' ')
   disp('**Depths are not incremental.')
   disp('**Bad bins were not trimmed.')
   close(f);
   return
else 
   disp('Redefining the "Bin" dimension')
	B=resize(B,length(goodBins));
end
f{'D'}.water_depth(:) =MSL;
close(f)
    
%Done
disp(' ')
disp(['File ' trimFile ' has been modified'])
disp(['##' num2str(theBinNum-bincnt) ' bins were removed from the top of the water column'])
thecomment='Bins were trimmed by trimBins.m based on RDI surface output.';
history(trimFile,thecomment);

starbeam(trimFile);

