function [MSL, Dstd] = trimbins(rawcdfFile,trimFile,MSL,Dstd,ADCP_offset,percentwc)

% function [MSL, Dstd] = trimbins(rawcdfFile,trimFile,MSL,Dstd,ADCP_offset,percentwc)
%
% Modifies the data file so the bins all fall below mean sea level(MSL) plus
% the tidal range using a MSL provided by the pressure sensor, findsurface.m
% or user input.
%
% If no imput values are given it will prompt the user for method of
% trimming the bins by using a pressure sensor, findsurface.m, or user input.
% User can fudge the amount trimmed by specifying the percentage of water
% depth to trim to.  The default is 105%
%
% Trimbins also generates the range to boundary information (variable
% brange) for either uplooking or downlooking orientation.
%
%   NOTE:  If user suspects that the pressure sensor did not work properly
%   as to biofouling or other reasons, it is recommended that the user use
%   findsurface.m or input MSL and Dst manually when using trimbins.
%
%   NOTE:  Because data within 6% of the surface is contaminated, the user
%   may wish to keep only 94% of the water column.  However, there is some
%   reason to beleive that useful information can be obtained from bins
%   above the surface.  This is for the user to decide.
%   Reference: Principles of Operation: A Practical Primer (for ADCPs) pg. 38
%
%INPUTS:
%	rawcdfFile =  raw netCDF version of ADCP data, can be the maskFile or trimFile
%   trimFile = trimmed raw netCDF file
%   MSL = mean sea level.  meanRange from findsurface.m, or []
%	Dstd = the standard of deviation of the depth which should
%			 be an approximate tidal fluctuation
%           stdRange from findsurface.m, or []
%	ADCP_offset = is the an attribute of the depth variable,
%		called transducer_offset_from_bottom in the Netcdf files
%		If not given, the rdsurface requests a Netcdf file to find it.
%   percentwc = amount to trim, as a percentage of water column%
%
%OUTPUTS:
%	MSL = mean sea level based on RDI surface output
%	Dstd = standard of deviation to give an approximate tidal variation
%
% Example calls
% Run in non-interacive mode, supplyins a structure of settings:
% [MSL, Dstd] = trimbins(settings.trimbins)
% where settings =
%     settings.method = 'Pressure Sensor';
%     % methods: 'User Input' | 'Pressure Sensor' | 'USGS Surface Program'
%     settings.rawcdfFile = input file; % leave this alone
%     settings.trimFile = output file; % leave this alone
%     % method to remove bins above the surface (see trimbins documentation
%     settings.method = 'Pressure Sensor'; %'RSSI peak detect' | 'User Input'
%     % percent of water column to make sure is preserved when trimming (1 = 100%)
%     settings.percentwc = 1.12; % to capture full range of tide
%     settings.ADCP_offset = settings.rdi2cdf.transducer_offset; % leave this alone
%     % if user input is selected, one may include
%     settings.MSL mean sea level, m
%     settings.Dstd tidal variation, m
%     % or
%     % these control the MATLAB surface detect algorithm findsurface.m
%     settings.findsurface.cdfFile = settings.rawcdf; % leave this alone
%     settings.findsurface.S = 32; % local salinity assumption
%
% examples can be found in runadcp.m

%%% START USGS BOILERPLATE -------------%
% Use of this program is described in:
%
% Acoustic Doppler Current Profiler Data Processing System Manual
% Jessica M. Côté, Frances A. Hotchkiss, Marinna Martini, Charles R. Denham
% Revisions by: Andrée L. Ramsey, Stephen Ruane
% U.S. Geological Survey Open File Report 00-458
% Check for later versions of this Open-File, it is a living document.
%
% Program written in Matlab v7.1.0 SP3
% Program updated in Matlab 7.2.0.232 (R2006a)
% Program ran on PC with Windows XP Professional OS.
%
% "Although this program has been used by the USGS, no warranty,
% expressed or implied, is made by the USGS or the United States
% Government as to the accuracy and functioning of the program
% and related program material nor shall the fact of distribution
% constitute any such warranty, and no responsibility is assumed
% by the USGS in connection therewith."
%
%%% END USGS BOILERPLATE --------------

%sub-functions:
%	Netcdf toolbox
%   pressurecalcs.m
%   findsurface.m

% Written by Jessica M. Cote
% for the U.S. Geological Survey
% Coastal and Marine Geology Program
% Woods Hole, MA
% http://woodshole.er.usgs.gov/
% Please report bugs to jcote@usgs.gov

% 1/13/11 added fix to keep handle following resize() call
% Updated 13-jul-2009 (MM) remove references to goodends.m, surface.exe
% Updated 18-jun-2008 (MM) change to SVN revision info
% updated 12-feb-08 (MM) make sure that if pressure sensor is not present
% and user chose pressure sensor that the range to boundary is used.
% Uplookers only.
% updated 06-feb-08 (MM) change the logic so that while MSL for trimming
% might be computed for the purposes of trimming uplooker data (most
% accurate) the range to boundary information is always computed using
% findsurface.m
% updated 21-dec-07 (MM) replace SURFACE.EXE with findsurface.m
% updated 2-feb-2007 (MM) replace getinfo with inputdlg
% updated 31-jan-2007 (MM) remove batch calls
% updated 30-jan-2007 (MM) make trimming to surface more flexible, improve docs
% updates 25-jan-2007 (MM) abort if ADCP is not uplooking
% updated 01-jan-2007 (MM) - move starbeam call to runadcp
% version 2.0
% 29-dec-2006 (MM) change to allow struct of arguements to be passed
% version 1.0
% updated 20-Dec-2004 (SDR) runs all revisions with batch
% updated 17-Dec-2004 (SDR) can now be run independently and prompt user
%   for method of depth calculation
% updated 29-Nov-2004 (SDR) added comment in history for depth sensor
%   input, if available
% updated 03-Apr-2003 (ALR) clarified note for Height variable
% updated 16-Aug-2001 If using surface.exe, uses 94% of MSL + 1/2 tidal range
%		to trim surface bins (ALR)
% updated 10-Jul-2001 added ability to use one OR two binary raw ADCP files (ALR)
% updated 28-Dec-2000 09:28:32 - added linefeeds to comment/history attribute (ALR)
% updated 12-Jan-2000 11:07:35
%	-runs with batch and has all revisions
% updated 22-Oct-1999 14:40:41
% updated 15-Oct-1999 16:49:45

% get the current SVN version- the value is automatically obtained in svn
% is the file's svn.keywords which is set to "Revision"
rev_info = 'SVN $Revision: 2101 $';
fprintf('%s %s running\n',mfilename,rev_info)

% this management of variables is circuitous, but living with it for
% now for backwards compatability

if nargin==1 && isstruct(rawcdfFile),
    settings = rawcdfFile;
    if isfield(settings,'rawcdfFile'), rawcdfFile = settings.rawcdfFile;
    else rawcdfFile = ''; end
    if isfield(settings,'trimFile'), trimFile = settings.trimFile;
    else trimFile = ''; end
    if isfield(settings,'percentwc'), percentwc = settings.percentwc;
    else percentwc = 1.05; end
    if isfield(settings,'MSL'), MSL = settings.MSL;
    else MSL = ''; end
    if isfield(settings,'Dstd'), Dstd = settings.Dstd;
    else Dstd = ''; end
    if isfield(settings, 'ADCP_offset'), ADCP_offset = settings.ADCP_offset;
    else settings.ADCP_offset = ''; end
    if isfield(settings,'findsurface') && isstruct(settings.findsurface),
        if isfield(settings.findsurface,'S'), S = settings.findsurface.S; end
        if isfield(settings.findsurface,'D'), D = settings.findsurface.D; end
        if isfield(settings.findsurface,'trimrange'),
            trimrange = settings.findsurface.trimrange; end
        if isfield(settings.findsurface,'shiftbin'), shiftbin = settings.findsurface.shiftbin; end
    end
end

if ~exist('rawcdfFile','var'), rawcdfFile = ''; end
if ~exist('ADCP_offset','var'), settings.ADCP_offset = ''; end
if ~exist('MSL','var'), MSL = ''; end
if ~exist('percentwc','var'), percentwc = []; end
if ~exist('S','var'), S = []; end
if ~exist('D','var'), D = []; end
if ~exist('shiftbin','var'), shiftbin = 0; end
if ~exist('trimrange','var'), trimrange = []; end

if isempty(rawcdfFile), rawcdfFile = '*'; end

% Get ADCP data file to trim.
if any(rawcdfFile == '*')
    [theFile, thePath] = uigetfile({'*M.cdf','Pick masked raw netCDF file';...
        '*.cdf','Pick any raw netCDF file'},...
        'Select raw netCDF ADCP File:');
    if ~any(theFile), return, end
    rawcdfFile = fullfile(thePath, theFile);
end

% copy the file 
copyfile(rawcdfFile,trimFile);

%open the netcdf file to trim and get some info
f=netcdf(trimFile,'write');
if isempty(f), return, end;

depth=f{'D'};
bad_num = fillval(depth);

B = f('bin');
theBinNum = B(:);
ensembles = f{'Rec'}(:);
EnsNum = length(ensembles);
fprintf('Pressure sensor available? %s\n',f.depth_sensor(:))

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

% get the range to boundary using findsurface.m
% TODO - allow user to enter these settings interactively in the script
settings.D = D;
settings.S = S;
settings.trimrange = trimrange;
settings.cdfFile = trimFile;
settings.shiftbin = shiftbin;

% if missing ensembles, send just the first chunk
% Rec is the indicator record and had FillValue in place of missing ensembles- 
% the variable called ensembles is a dimension so can't have a
miss_ens=find(f{'Rec'}(:) > 2e9);  % Rec is an int, so has fillVal=2.1475e9
% check as nan in case autonan is on during run
if isempty(miss_ens)
    miss_ens=find(isnan(f{'Rec'}(:)));
end
if isempty(miss_ens)
    settings.ensembles=f{'ensembles'}(:);
else
    settings.ensembles=[1 miss_ens(1)-1];
end


fsdata = findsurface(settings);
if all(isnan(fsdata.Range2Boundary)), noboundary = 1;
    disp('No boundary (surface) was detected')
else noboundary = 0;
end
rnote = 'range to boundary based on findsurface.m output';
methods = {'Pressure Sensor',...
    'USGS Surface Program', ...
    'User Input'}; 
if strcmpi(f.depth_sensor(:),'NO'), methods = methods(2:end); end

if strcmpi(f.orientation(:),'UP'),
    % get the mean sea level and tidal range for the given ensemble numbers
    if isempty(MSL),
        if exist('settings','var') && isfield(settings,'method'),
            buttonname = settings.method;
        else
            [buttonname, ok] = listdlg(...
                'PromptString','What method to trim the bins?',...
                'SelectionMode','single',...
                'ListString',methods);
            if ok, buttonname = methods{buttonname};
            else % cancel was pressed
                buttonname = 'User Input'; % just ask
            end
        end
        % oops, got pressure sensor when there is no pressure sensor...
        if strcmp(buttonname,'Pressure Sensor') && ...
                strcmpi(f.depth_sensor(:),'NO'),
            disp('Pressure sensor method for trimming requested when ADCP has no pressure sensor')
            [buttonname, ok] = listdlg(...
                'PromptString','There is no pressure sensor: what method to trim the bins?',...
                'SelectionMode','single',...
                'ListString',methods);
            if ok, buttonname = methods{buttonname};
            else % cancel was pressed
                buttonname = 'User Input'; % just ask
            end
            fprintf('Using %s instead\n',buttonname)
        end
        switch buttonname
            case 'Pressure Sensor'
                [MSL,Dstd,Dsurf]=pressurecalcs(trimFile);
                MSLnote = 'water depth = MSL from pressure sensor, by trimbins';
                %to append to history later
                thecomment=sprintf('Bins were trimmed by %s %s based on depth sensor input information.\n',...
                    mfilename, rev_info);
                % Dsurf is range from ADCP to surface, not adjust for height                
            case 'USGS Surface Program'
                MSL=fsdata.MSL; % mean range to surface + ADCP offset
                MSLnote = 'water depth = MSL from USGS surface detect algorithm, by trimbins';
                Dstd=fsdata.stdRange; % std of range to surface
                Dsurf=fsdata.Range2Boundary; % Range to surface time series
                % to append history later
                thecomment=sprintf('Bins were trimmed by %s %s using 94%% of the findsurface.m output.\n',...
                    mfilename, rev_info);
                rnote = 'range to boundary based on findsurface.m output';
            case 'User Input'
                prompt  = {'Enter the mean sea level value, m:',...
                    'Enter the half_the_tidal_range, m:'};
                def     = {'0','0'};
                name   = 'User must input the water depth information';
                lineNo  = 1;
                dlgresult  = inputdlg(prompt,name,lineNo,def);
                MSL = str2double(dlgresult{1});
                MSLnote = 'water depth = MSL from user input, by trimbins';
                Dstd = str2double(dlgresult{2});
                Dsurf=[ ];
                %to append history later
                thecomment=sprintf('Bins were trimmed by %s %s based on user input depth information.\n',...
                    mfilename, rev_info);
        end %button switch
    else
        disp('User input water_depth and tidal variation')
        MSLnote = 'water depth = MSL from user input, by trimbins';
        Dsurf=[ ];
        %to append history later
        thecomment=sprintf('Bins were trimmed by %s %s based on user input depth information.\n',...
            mfilename, rev_info);
    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=str2double(bad_num);
        end

        depth=depth(:);
        idgood= depth < bad_num & depth > 0;
        depth=depth(idgood);
        clear idgood

    end


    %find the bins1 that fall below the given sea level including tidal fluctuation
    if isempty(percentwc),
        prompt  = {'Trim percentage of water column:'};
        def     = {'105'};
        title   = sprintf('MSL = %f with a variation of %f m',MSL,Dstd);
        lineNo  = 1;
        dlgresult  = inputdlg(prompt,title,lineNo,def);
        percentwc = str2double(dlgresult{1})./100;
    end

    disp(' ')
    disp('Finding the good bins');
    % goodBins = find(depth <= MSL+Dstd); Only trims bins out of water
    goodBins=find(depth <= (percentwc*(MSL+(Dstd))));
    fprintf('Trimming to %6.3f percent of the mean seal level + Dstd\n',percentwc*100)

    if isempty(goodBins)
        disp('## No bins were found within the specified depth range - cannot continue trimming');
        close(f)
        return
%     elseif length(goodBins) == length(depth);
%         disp('No bad bins were found, file is unmodified');
%         f{'D'}.WATER_DEPTH = ncfloat(MSL);
%         f.WATER_DEPTH = ncfloat(MSL);
%         f{'D'}.WATER_DEPTH_source = ncchar(MSLnote);
%         f.WATER_DEPTH_source = ncchar(MSLnote);
%         close(f)
%         return

    end

    %Check to make sure that everything is in 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'}.transducer_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')
        disp('May take a few minutes...')
        % this effectivle trims the bins.  Mlint would have you not assign the
        % result of resize to B, don't mess with it.
        B=resize(B,length(goodBins));
        f=parent(B);        %needed in 2010b and newer
    end

    %Record some information
    f{'D'}.WATER_DEPTH = ncfloat(MSL);
    f.WATER_DEPTH = ncfloat(MSL);
    f{'D'}.WATER_DEPTH_source = MSLnote;
    f.WATER_DEPTH_source = MSLnote;
    f{'D'}.WATER_DEPTH_datum = 'MSL';
    f.WATER_DEPTH_datum = 'MSL';    
    disp(MSLnote);

    disp(['# in Dsurf ' num2str(length(Dsurf))])
    disp(['# in ensemble ' num2str(EnsNum)])
    disp(' ')
    disp(['File ' trimFile ' has been modified'])
    disp(['## ' num2str(theBinNum-bincnt) ' bins were removed from the top of the water column'])

else
    % we can't do anything in this function with downlooking data...
    disp('trimbins: ADCP is not uplooking, bins were not trimmed.')
    thecomment=sprintf('Uplooking; bins were NOT trimmed by %s %s;\n',...
        mfilename, rev_info);
end

% as per CRS 8/5/09 always write brange
f{'brange'} = ncfloat('ensemble');
f{'brange'}.long_name = ncchar('range to boundary from transducer head');
f{'brange'}.units = ncchar('m');
f{'brange'}.FillValue_ = 1.0e35;
if noboundary, rnote = ['no boundary detected, ' rnote]; end
f{'brange'}.NOTE = ncchar(rnote);
endef(f)

if ~isempty(fsdata.Range2Boundary) && ~noboundary,
    if length(fsdata.Range2Boundary) == EnsNum,
%         % clean up brange a la Kurt
%         % first, see if we have the right functions...
%         if exist('filtfilt','file') && exist('butter','file'),
%             disp('signal processing toolbox functions present, cleaning brange')
%             dT = gmean(diff(f{'TIM'}(:))); % time between samples in julian day
%             fsdata.Range2Boundary = clean_brange(fsdata.Range2Boundary,dT.*24);
%         else
%             disp('signal processing toolbox functions NOT present, brange NOT cleaned')
%         end
        f{'brange'}(1:EnsNum) = fsdata.Range2Boundary;
    elseif length(fsdata.Range2Boundary) < EnsNum
        f{'brange'}(1:length(fsdata.Range2Boundary)) = fsdata.Range2Boundary;
    else   
        f{'brange'}(:) = ones(1,length(f{'Rec'}(:))).*NaN;
        fprintf('%s: EnsNum %d ~= length(fsdata.Range2Boundary) %d\n',...
            mfilename, EnsNum,length(fsdata.Range2Boundary))
        fprintf('%s: brange variable not written to file\n',mfilename)
    end
else
    f{'brange'}(:) = ones(1,length(f{'Rec'}(:))).*NaN;
        fprintf('%s: brange variable not written to file\n',mfilename)
end

fprintf('Data is from bin 1 at %7.2f to last bin at %7.2f meters\n',f{'D'}(1),f{'D'}(end))

% add minimums and maximums
f=add_minmaxvalues(f);

close(f)
if noboundary, thecomment = ['no boundary detected, ' thecomment]; end

%Done
history(trimFile,thecomment);