function whisk_adcp(theSite)
%
%   whisk_adcp(theSite)
%
%   Function to create whiskers plots of current at Sites B and C, along
%   with low-passed wind stress

%Soupy Alexander, 6/6/04

%This Matlab m-file was used to create portions of U.S. Geological Survey 
%Data Series 85. 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.

%Citation:  Butman, Bradford, Alexander, P. Soupy, Anderson, S.P., 
%Lightsom, F.L., Scotti, Alberto, and Beardsley, R.C., 2004, The 
%Massachusetts Bay Internal Wave Experiment, August 1998:  Data Report:  
%U.S. Geological Survey Data Series 85, 1 DVD-ROM.

if strcmp(lower(theSite), 'b')
    theFiles = {'5372adc-a_wave.nc'};
    listDepths = [10 30 50 70 100];
elseif strcmp(lower(theSite), 'c')
    theFiles = {'5352adc-a_wave.nc'};
    listDepths = [10 20 30 40];
end

%Full Record
theTimeFrame = [datenum(1998, 07, 29, 00, 00, 00) ...
    datenum(1998, 09, 03, 00, 00, 00)];
timeVector = datenum2julian(min(theTimeFrame)):1/24:...
    datenum2julian(max(theTimeFrame));
theTicks = min(theTimeFrame):5:max(theTimeFrame);
theCode = 6;
dataDir = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES';
yLimits = [-15 15];

figure, orient landscape

for indexDepth = 1:length(listDepths)
    subplot(length(listDepths), 1, indexDepth)
    thePos = get(gca, 'position');
    
    ncID = netcdf(fullfile(dataDir, theFiles{1}), 'nowrite');
    theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    theDepths = ncID{'depth'}(:);
    targDepth = value2index(theDepths, listDepths(indexDepth), 50);
    thisDepth = theDepths(targDepth);
    
    theU = ncID{'u_1205'}(:,targDepth);
    theV = ncID{'v_1206'}(:,targDepth);

    theBads = find(((theU + theV)) > 1e3 | ((theU + theV) < -1e3));
    theU(theBads) = NaN;
    theV(theBads) = NaN;
    
    
    for indexTime = 1:length(timeVector)-1;
        inTime = find((theTime >= timeVector(indexTime)) & ...
            (theTime < timeVector(indexTime+1)));
        oldU(indexTime) = nanmean(theU(inTime,:));
        oldV(indexTime) = nanmean(theV(inTime,:));
    end
    
    theBads = find(isnan(oldU + oldV));
    oldTime = timeVector(1:end-1);
    oldTime(theBads) = [];
    oldU(theBads) = [];
    oldV(theBads) = [];
    
    newU = interp1(oldTime, oldU, timeVector(1:end-1));
    newV = interp1(oldTime, oldV, timeVector(1:end-1));
    
    [lpU, lpTime] = plfilt(newU, timeVector, 6);
    [lpV, junk] = plfilt(newV, timeVector, 6);

    [theAngles, theLengths] = cart2pol(lpU, lpV);
    theAngles = theAngles*180/pi;
    sticksafe(julian2datenum(lpTime), ...
        zeros(size(theLengths)), theLengths, theAngles)
    xlim([min(theTimeFrame) max(theTimeFrame)])
    ylim(yLimits)
    set(gca, 'xtick', theTicks)
    datetick('x', 6, 'keepticks', 'keeplimits')
    ylabel('cm/s')
    title([num2str(thisDepth, '%.0f') ' mbs'])
    box('on')
    set(gca, 'position', thePos);
    sticksafe
end

suptitle({['Site ' upper(theSite)]; ...
    'Low-passed Current Observations, ADCP'})