function data = findsurface(cdfFile, ensembles, S, D, trimrange, verbose, allout)
%    findsurface(settings) also works
% findsurface find range from ADCP head to sea surface using ADCP backscatter
% function data = findsurface(cdfFile, ensembles, S, D, trimrange, verbose, allout)
%
% Inputs
% cdfFile = raw netCDF data file with AGC data
% ensembles = ensemble number range, [1 n] or [] or for all ensembles
% for sound speed correction, use the following
% S = salinity to assume (preferably what was set in the ADCP at deployment)
%       S = [] the nominal value of 32 will be assumed
% D = depth, m to use, this will override a pressure sensor if present.
% trimrange = [minDepth maxDepth] | [];
%       how to trim outliers, surface detech is often noisy and
%       needs a thumb-finger trim for bets results
%       if no trim is wanted, use [] to override
% verbose - turn on the plotting feature and watch live.  This
%       will seriously slow things down and is used for testing.
% allout = 1 - pass out all data used, temperature, ensemble number, etc. for
%       diagnostic purposes or stand alone use of this function
% inputs may be entered as a struct 
%   (how this is called by the ADCP toolbox)
%   settings.cdfFile
%   settings.ensembles
%   settings.S
%   settings.D
%   settings.verbose
%   settings.allout
%   settings.trimrange
%   settings.shiftbin = 0; % or N bins, available only for script invocation with the
%           settings structure.  findsurface starts looking for a peak from the bin
%           furthest from the ADCP head to the middle of the profile.  Sometimes,
%           if there are many extra bins, a second acoustic bounce may make a
%           second surface appear.  Shiftbin allows the user to try to shift the
%           search area closer to the ADCP by N bins.
%
% Outputs
% data structure with the following fields
%   data.Range2Boundary = range from ADCP head to the water surface, time series
%   data.meanRange = mean of the Range2Boundary time series
%   data.stdRange = STD of the Range2Boundary time series
%   data.MSL = mean sea level estimated from range to boundary
%   data.ADCP_offset = offset from ADCP to bed used to compute MSL
% returned if allout == TRUE
%   data.ensembles = vector of ensemble numbers (cdf record)
%   data.peakfit = the peak in AGC determined by parabolic fit for each beam
%   data.WMVpeak = peak in water mass volume as computed using S, D, etc.
%   data.pitch = pitch data
%   data.roll = roll data
%   data.press = pressure data if present in ADCP
%   data.Tx = temperature data
%   data.S = salinity used
%   data.SVelADCP = soubd velocity computed internally by the ADCP

% This program reads selected ensembles, finds the peak in each beam 
% associated with the sea surface, computes the range to the surface for each 
% beam and the average range.  To do this, it corrects for range attenuation, 
% assuming typical salt water absorption, and beam spread. It finds the 
% largest echo intensity then computes a parabolic fit. 
% 
% This method is from :
%   Visbeck, M. and J. Fischer, "Sea surface conditions remotely sensed by 
%      upward-looking ADCPs." J. Atmosph. and Oceanic Technol., 12, 141-149,
%      1995.
% Augmented by ideas gleaned from
%   Shcherbina, A. Y., D. L. Rudnick, et al. 
%      "Ice-Draft Profiling from Bottom-Mounted ADCP Data." 
%      Journal of Atmospheric and Oceanic Technology 22(8): 1249-1266.
%   Deines, Kent, "Backscatter Estimation Using Broadbad Acoustic Doppler 
%      Current Profilers," unpublished RD Instruments technical document, 1999.
% and information from the Teledyme RDI 
%   "Principles of Operation, A Practical Primer"
%      originally written in 1996 for RD Instruments by Lee Gordon.
% and the basics from
%   _Fundamentals of Acoustics_, Kinsler, Frey, Coppens & Sanders, 3rd Ed., 1982.

%%% 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 --------------

% Written by Marinna Martini 21 dec 07
% for the U.S. Geological Survey
% Woods Hole Science Center
% used by trimbins in the ADCP toolbox

% Update 23-jan-08 (MM) apply an ADCP height offset and return real mean
% sealevel
% Update 22-jan-08 (MM) use names that reflect the true meaning of the data
% Update 18-jan-08 (MM) fix a stupid indexing mistake
% Update 16-jan-08 (MM) profile & restructure for better speed
% Update 10-jan-08 (MM) make compatible with ADCP Toolbox
rev_info = 'SVN $Revision: 2101 $';
fprintf('%s %s running\n',mfilename,rev_info)

% translate inputs from struct of settings
if nargin==1 && isstruct(cdfFile),
    settings = cdfFile;
    fields = {'cdfFile', 'ensembles', 'S', 'D', 'trimrange',...
        'verbose', 'allout', 'shiftbin'};
    for ifield = 1:length(fields),
        if isfield(settings,fields{ifield}), 
            eval(sprintf('%s = settings.%s;',fields{ifield},fields{ifield}))
        else eval(sprintf('%s = '''';',fields{ifield}))
        end
    end
end

% test inputs
if (~exist('cdfFile','var') || isempty(cdfFile)),
    [theFile, thePath] = uigetfile('*.cdf', 'Select netCDF raw ADCP File:');
    if ~any(theFile), return, end
    if thePath(end) ~= filesep, thePath(end+1) = filesep; end
    cdfFile = [thePath theFile];
end
if ~isfield(settings,'S'),
    S = 32;
end
if ~isfield(settings,'D'),
    D = [];
end
if ~isfield(settings,'trimrange'),
    trimrange = [];
end
if ~isfield(settings,'verbose') || isempty(verbose),
    verbose = 0;
end
if ~isfield(settings,'allout') || isempty(allout),
    allout = 0;
end
if ~isfield(settings,'shiftbin'),
    shiftbin = 0;
end
if shiftbin,
    fprintf('User has elected to shift the peak search area towards the ADCP by %d bins\n',shiftbin)
end

cdf = netcdf(cdfFile,'r');
ADCP_offset=cdf{'D'}.transducer_offset_from_bottom(:);
if isempty(ADCP_offset) || (ADCP_offset <= 0),
    disp(sprintf('%s: bad or missing transducer_offset_from_bottom in %s',...
        mfilename, cdfFile))
    disp(sprintf('%s: ADCP offset applied to pressure data is %6.2f m',...
        mfilename,ADCP_offset))
end

% get basic information
nbeams = cdf.beams_in_velocity_calculation(:);
if nbeams > 4, 
    disp('More than 4 beams may cause issues with tilt compensation'); 
end

nbins = length(cdf{'D'});
BinSize = cdf{'D'}.bin_size(1);
ADCP_depth = cdf{'D'}.WATER_DEPTH(1)-ADCP_offset;
if isempty('D'), D = ADCP_depth; end
% This field, set by the WT-command, contains the length of the
% transmit pulse. When the Workhorse receives a <BREAK>
% signal, it sets the transmit pulse length as close as possible to
% the depth cell length (WS-command). This means the Workhorse
% uses a WT command of zero. However, the WT field
% contains the actual length of the transmit pulse used.
% Scaling: LSD = 1 centimeter; Range = 0 to 65535 cm (2150
% feet)
xmitPlen = cdf.transmit_pulse_length(:)./100; % cm to m
% LagD / Transmit lag distance
% Lag length displayed by file viewer in winadcp
% This field, determined mainly by the setting of the WMcommand,
% contains the distance between pulse repetitions.
% Scaling: LSD = 1 centimeter; Range = 0 to 65535 centimeters
lagD = cdf.transmit_lag_distance(:)./100; % cm to m
beamAngle = cdf.beam_angle(:);
% we need some basic acoustic information
% From Deines, 1999
%         STEP 5 – DETERMINE TRANSMIT POWER
% To obtain absolute backscatter data, transmit power must
% be estimated. Table 1 shows the expected transmit power
% (dBW) for RDI BBADCP’s. In general this power is proportional
% to the input voltage (dBV). Some low frequency
% Broadband systems have a high power module that removes
% this dependence upon input voltage. For systems running
% from alkaline batteries the input varies 6 dB over the life of
% the battery. This variation is only 3 dB over the middle 80%
% of the battery life. The estimate of power in Table 1 for
% Workhorse systems is at 33 volts, which is in the center of the
% plateau of the battery discharge curve.
% Workhorse:
% System Frequency,kHz, Name C, dB, PdBW, Rayleigh Distance,meters
% 76.8 WH Long Ranger -159.1 24.0 1.30
% 307.2 WH Sentinel -143.5 14.0 0.98
% WH Monitor -143.0 14.0 0.98
% 614.4 RioGrande -139.3 9.0 1.96
% 1228.8 WH Sentinel -129.1 4.8 1.67
% WH Monitor -129.1 4.8 1.67
NominalFreq = cdf.frequency(:);
switch NominalFreq,
    case 300
        %f = 307200; % Hz transmit frequency from .dlg file
        SL = 14.0; % PdBW source level
        %SL1 = 216.3; % dB re micropascal at 1 m
        %SL21 = 180.0; % dB re micropascal at 34.6 m
        %ffm = 34.6; % m
        Ba = 3.7; % degrees beam apeture from .dlg file
        Alpha = mean([0.062 0.082]); % from the primer
    case 600
        %f = 614400; % Hz
        SL = 9.0; % PdBW source level
        %SL1 = 217.1; % dB re micropascal at 1 m
        %SL21 = 180.0; % dB re micropascal at 21.6 m
        %ffm = 21.6; % m
        Ba = 3.7; % degrees beam apeture from .dlg file
        Alpha = mean([0.14 0.20]); % from the primer
    case 1200
        %f = 1228800; % Hz
        SL = 4.8; % PdBW source level
        %SL1 = 214.0; % dB re micropascal at 1 m
        %SL21 = 180.0; % dB re micropascal at 8.7 m
        %ffm = 8.7; % m
        Ba = 3.7; % degrees beam apeture from .dlg file
        Alpha = mean([0.44 0.66]); % from the primer
end
% H = cdf{'D'}.ADCP_depth(1) - cdf{'D'}.transducer_offset_from_bottom(1);
% % HsH = Hs/H = 0.027 for 20 deg beam angle and 4 deg beam aperture
% HsH = (sin(beamAngle*pi/180).*sin(Ba*pi/180));
% HsH = HsH./(cos((beamAngle+Ba)*pi/180).*cos((beamAngle-Ba)*pi/180));

Ball = 0:nbins-1;
CtrBin1 = cdf{'D'}.center_first_bin(1);
% Rv = distance from transducer to depth cell, vertically
% we recompute here from scratch to avoid confusion if cdf{'d'} was
% adjusted (those bin depths are relative to the sea bed, we want Range)
Rv = CtrBin1+ BinSize*Ball;

% find the ensemble numbers requested in the data
ensnums = cdf{'Rec'}(:);
if ~exist('ensembles','var') || isempty(ensembles), % no ensemble numbers indicated
    eidx(1) = 1;
    eidx(2) = length(ensnums);
else
    tempidx = find(ensnums >= ensembles(1),1,'first');
    if isempty(tempidx),
        eidx(1) = 1;
        disp('first ensemble number not found, starting at beginning of cdf file')
    else
        eidx(1) = tempidx;
    end
    tempidx = find(ensnums >= ensembles(2),1,'first');
    if isempty(tempidx),
        eidx(2) = length(ensnums);
        disp('last ensemble number not found going to the end of cdf file')
    else
        eidx(2) = tempidx;
    end
end
disp(sprintf('using ensembles %d to %d',ensnums(eidx(1)),ensnums(eidx(2))))
clear ensnums

% pull all data into memory for speed
data.ensembles = cdf{'Rec'}(eidx(1):eidx(2));
agc = cell(nbeams,1);
for ibeam = 1:nbeams,
    varname = sprintf('AGC%d',ibeam);
    agc{ibeam} = cdf{varname}(eidx(1):eidx(2),:).*cdf{varname}.norm_factor(1); % converts to dB
end
Tx = cdf{'Tx'}(eidx(1):eidx(2));
data.pitch = cdf{'Ptch'}(eidx(1):eidx(2));
data.roll = cdf{'Roll'}(eidx(1):eidx(2));
data.SVelADCP = cdf{'sv'}(eidx(1):eidx(2)); % SVelADCP = Sound Velocity at ADCP head
if ~isempty(cdf{'Pressure'}),
    %convert from pascals to decibars: db = pascal/10000
    press = cdf{'Pressure'}(eidx(1):eidx(2)).*1e-4;
elseif ~isempty(D), % user provided an override for depth
    press = ones(length(data.ensembles),1).*D;
else
    press = []; % we have no depth to work with.
end

peakfit = ones(length(data.ensembles),nbeams).*NaN;
WMVpeak = peakfit;
%press = ones(nrec,1).*NaN;
Rs = ones(nbeams, nbins).*NaN;
n = ceil(nbins/2); % number of points to look back for the peak.


tic
for irec = 1:length(data.ensembles),
    if verbose,
        clf;
        hold on;
    end
    for ibeam = 1:nbeams,
        agcpall = agc{ibeam}(irec,:); % pull out hte current profile
        pitch = data.pitch(irec);
        roll = data.roll(irec);
        SVelADCP = data.SVelADCP(irec); % SVelADCP = Sound Velocity at ADCP head

        % correct the beam angle for the effective beam angle for this
        % deployment, inferred from Shcherbina:  "because the echo form
        % this angle arrives earlier than the one from the nominal beam
        % angle, th distance to the surface is underestimated."
        % 
        % EffBeamAngle = beamAngle + Ba/2;
        EffBeamAngle = beamAngle;
        % correct the vertical distance from transducer to depth cell, 
        % Rv, for sound speed variation if pressure and salinity are known
        % compute sound speed based on ADCP variables
        if ~isempty(press),
            SVel = 1449.2+4.6.*Tx(irec)-0.055.*(Tx(irec)^2)+0.00029.*(Tx(irec)^3)+(1.34-0.01*Tx(irec))*(S-35)+0.016*press(irec);
        end
        % assumption of sound speed rather than a sounds speed profile
        % Visbeck uses a harmonic mean soundspeed for his projection, we
        % use actual soundspeed based on temperature at the ADCP
        % This will be a little inaccurate in highly stratified conditions.
        if ~isempty(S) && ~isempty(press),
            % comment the next line to overide sound speed correction
            Rv = CtrBin1.*(SVelADCP/SVel)+ BinSize.*(SVelADCP/SVel).*Ball;
        end

        % map the nominal Rv onto the beam for the orientation of *this*
        % ensemble = a slant range for each beam
        % we need this as an estimate of range to bin for other corrections
        switch ibeam
            case 1,
                Rs = Rv./(cos((roll+EffBeamAngle)*pi/180).*cos(pitch*pi/180));
            case 2,
                Rs = Rv./(cos((roll-EffBeamAngle)*pi/180).*cos(pitch*pi/180));
            case 3,
                Rs = Rv./(cos((pitch-EffBeamAngle)*pi/180).*cos(roll*pi/180));
            case 4,
                Rs = Rv./(cos((pitch+EffBeamAngle)*pi/180).*cos(roll*pi/180));
        end
        % overide tilt correction
        %Rs = Rv;

        % calculate target strength, correct for absorption
        % sonar equation from Kinsler: EL = SL-2TL+TS
        % Echo Level returned = Source Level - two way Transmission Loss + Target Strength, db
        % from the primer
        % EI = SL + SV + constant - 20log(R) - 2aR
        % EI = echo intensity (dB) (AGC or RSSI) 
        % SL = source level(dB) (transmit power)
        % SV = water mass volume backscattering strength (dB) (or, TS term)
        % TL = two way transmission loss term ( - 20log(R) - 2aR)
        % i.e. the target strength of a mass of particles in the sampling volume
        % Alpha = absorption Coeff (dB/m)
        % constant - we will omit this, use relative measurements for surface detection
        % select transducer output specifics per RDI memo 18Nov2005
        % SV = EI - SL + 20log(R) + 2aR
        % convert agcp from counts to dB 0.45 dB per Workhorse count.
        % 
        TL = 20.*log(Rs)+2.*Alpha.*Rs; % transmission losses per bin
        WMV = agcpall - SL + TL;
        % everything improved when I removed this!
        %WMV = agcpall; % override absorption correction
        
        % select just the end of the profile to prevent false positives
        % have not yet figured out how to deal with multible bounces
        % user has option to shift this search area if htere is a double
        % peak by the value shiftbin.  So shiftbins in from the end to
        % shiftbin + n bins from the end will be the search area
        if n+shiftbin > nbins,
            disp(sprintf('shiftbin value of %d was too large, omitting', shiftbin))
        end
        WMV = WMV(end-n-shiftbin:end-shiftbin); % target strength, one beam, some bins
        Rs = Rs(end-n-shiftbin:end-shiftbin); % slant ranges

        % this method uses standard MATLAB methods
        dRs = 0.1; %(Rs(end)-Rs(1))/ifactor;
        Rsi = Rs(1):dRs:Rs(end);
        WMVi = interp1(Rs,WMV,Rsi,'spline');
        [m, idx] = max(WMVi);
        if irec == 1, %>8999, % diagnostic stuff
            fprintf('irec %d, size WMV [%d %d]\n',irec,size(WMV,1),size(WMV,2))
        end
        WMVpeak(irec) = WMVi(idx);
        %WMVpeak(irec,ibeam) = WMVi(idx);
        Rspeak = Rsi(idx);
        
        % this method depends on the signal processing toolbox interp
        % % find the peak
        % ifactor = 5; % granularity of interpolation
        % WMVi = interp(WMV,ifactor);
        % [m, idx] = max(WMVi);
        % WMVpeak(irec) = WMVi(idx);
        % % interpolate slant range, find the range to match the peak bin
        % Rsi = interp(Rs,ifactor);
        % Rspeak = Rsi(idx);
        
        % Adjust for the filtering of the RSSI output as described in Shcherbina 
        % dH = 0.5 * tauBL + BBlag
        % Rspeak = Rspeak + 0.5.*BinSize + xmitPlen; % this is Shcherbina
        % translate the slant range result back to vertical
        % I'm not sure that translation from vertical to slant range and
        % back is actually doing much except provide a more exact range for 
        % absorption corrections, but it doesn't hurt...
        switch ibeam
            case 1,
                Rspeak = Rspeak.*cos((roll+EffBeamAngle)*pi/180).*cos(pitch*pi/180);
            case 2,
                Rspeak = Rspeak.*cos((roll-EffBeamAngle)*pi/180).*cos(pitch*pi/180);
            case 3,
                Rspeak = Rspeak.*cos((pitch-EffBeamAngle)*pi/180).*cos(roll*pi/180);
            case 4,
                Rspeak = Rspeak.*cos((pitch+EffBeamAngle)*pi/180).*cos(roll*pi/180);
        end
        peakfit(irec,ibeam) = Rspeak;
        
        if verbose,
            h(1) = plot(Rs,WMV,'o');
            hold on
            text(0.1,0.75,sprintf('Pitch %4.2f, Roll %4.2f',pitch,roll),'units','normalized');
            text(0.1,0.8,sprintf('Ensemble %d',data.ensembles(irec)),'units','normalized');
            xlabel('Range, m')
            ylabel('volume backscatter')
            h(2) = plot(Rsi,WMVi,'-'); 
            text(Rsi(end),WMVi(end),sprintf('%1d',ibeam),'units','data','color','b');
            h(3) = plot(Rspeak,WMVpeak(irec),'rx');
        end
    end % end for ibeam =
    if verbose, 
        ylims = get(gca,'ylim');
        if ~isempty(press), h(4) = plot(press(irec), median(ylims), 'gp'); end
        h(5) = plot(gmean(peakfit(irec,:)), median(ylims), 'rx');
        if ~isempty(press),
            text(0.1,0.85,sprintf('press - fit = %5.2f',press(irec)-gmean(peakfit(irec,:))),'units','normalized');
            if ~isnan(press(irec)-gmean(peakfit(irec,:))),
                legend(h,{'WMV','WMV interp','WMV interp max','pressure','mean fit'})
            end
            text(0.1,0.9,sprintf('SVelADCP/SVel = %5.1f/%5.1f = %3.2f',SVelADCP,SVel,SVelADCP/SVel),'units','normalized');
        end
        text(0.1,0.50,sprintf('xmit lag disance %4.2f',lagD),'units','normalized');        
        text(0.1,0.45,sprintf('lag length %4.2f',xmitPlen),'units','normalized');        
        text(0.1,0.7,sprintf('Bin size %4.2f',BinSize),'units','normalized');        
        text(0.1,0.65,sprintf('ADCP Depth %5.2f m',ADCP_depth),'units','normalized');        
        text(0.1,0.6,sprintf('ADCP height %5.2f m',ADCP_offset),'units','normalized');        
        text(0.1,0.55,sprintf('Salinity %5.2f m',S),'units','normalized');        
        pause(0.25)
    end
    if (irec<1000 && ~rem(irec,100)) || (irec>1000 && ~rem(irec,1000)),
        fprintf('%d ensembles read, at %d in %d min\n',irec-1,data.ensembles(irec),toc/60),
    end
end

close(cdf)

%%

% remove the aggregious outliers before returning the data
% also do some statistics
% was using Dm, now use Rm for mean RANGE 22-jan-08 (MM)
Rm = gmean(peakfit')'; % mean for the four beams
Rmts=gmean(Rm); % scalar mean for time series
Rmin=gmin(Rm);
Rmax=gmax(Rm);
disp(['findsurface found the following range statistics in meters: min = '...
   num2str(Rmin) ' max = ' num2str(Rmax) ' Mean = ' num2str(Rmts)])
%are we really getting only the good data?
p=figure; 
plot(data.ensembles,Rm,'r.')
xlabel('ensemble number')
ylabel('meters')
hold

% check with the user
if isempty(trimrange),
    prompt={'Minimum depth in meters:',...
        'Maximum depth in meters:'};
    title='Check trim values for surface data';
    DefAns={num2str(Rmin),num2str(Rmax)};
    dlgresult=inputdlg(prompt,title,1,DefAns);
    if ~isempty(dlgresult),
        Rmin=str2double(dlgresult{1});
        Rmax=str2double(dlgresult{2});
    end
else
    Rmin = trimrange(1);
    Rmax = trimrange(2);
end
% replace the wild points with NaN
% keep the length the same
idbad=find(Rm < Rmin | Rm > Rmax);
Rm(idbad) = ones(length(idbad),1).*NaN;
Rmts = gmean(Rm);
Rstd = gstd(Rm);
fprintf('User modified depth constraints in meters: min = %f, max = %f, mean = %f\n',...
    Rmin, Rmax, Rmts);
clf, plot(data.ensembles, Rm,'r.');
xlabel('ensemble number')
ylabel('meters')
hold;
        

fprintf('ADCP measured %f m mean range from ADCP to surface\n',Rmts);
disp(['The tidal range is approximately ' num2str(Rstd) ' m']);
disp(' ')

pause(3)
close(p)

%%
data.Range2Boundary = Rm;  % was data.Dsurf = Dm;
data.meanRange = Rmts; % was data.MSL;
data.stdRange = Rstd;
data.MSL = Rmts+ADCP_offset;
data.ADCP_offset = ADCP_offset;
if allout,
    data.peakfit = peakfit;
    data.WMVpeak = WMVpeak;
    data.press = press;
    data.Tx = Tx;
    data.S = S;
    data.xmitPlen = xmitPlen;
    data.lagD = lagD;
else 
    data = rmfield(data,{'pitch','roll','SVelADCP','ensembles'});
end