function whisk_pcm(theSite)
%
%   whisk_pcm(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 = {'fsi1468-a.nc'; 'vmcm0014-a.nc'; 'vmcm0009-a.nc'; ...
        '5371v-a_d1.nc'};
    listDepths = {'20 mbs'; '50 mbs'; '75 mbs'; '1 mab'};
elseif strcmp(lower(theSite), 'c')
    theFiles = {'5352adc-a_wave.nc'; '5344-a.cdf'};
    listDepths = {'20 mbs'; '45 mbs'};
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];

%Low-passed wind stress
windFile = 'lt_wind.nc';
ncWind = netcdf(fullfile(dataDir, windFile), 'nowrite');
windTime = singlejd(ncWind{'time'}(:), ncWind{'time2'}(:));
windU = ncWind{'WU_422'}(:);
windV = ncWind{'WV_423'}(:);
theBads = find(windU > 1e3 | windU < -1e3 | windV > 1e3 | windV < -1e3);
windU(theBads) = NaN;
windV(theBads) = NaN;

%Set up the new time base & interpolate the data
windTime_new = min(windTime):1/24:max(windTime);
[windU_i] = smart_interp(windTime, windU, windTime_new, 5);
[windV_i] = smart_interp(windTime, windV, windTime_new, 5);
[windU_lp, windTime_lp] = plfilt(windU_i, windTime_new, 6);
[windV_lp] = plfilt(windV_i, windTime_new, 6);

%Estimate the 10 m wind speed (buoy 44013 is at 5 m)
[tauX, tauY, windU_10, windV_10] = wstress(windU_lp, windV_lp, 5);
[tauU, tauV] = wstress(windU_10, windV_10);

figure, orient landscape

subplot(length(theFiles)+1,1,1)
thePos = get(gca, 'position');
[theAngles, theLengths] = cart2pol(tauU, tauV);
theAngles = theAngles*180/pi;
sticksafe(julian2datenum(windTime_lp), ...
    zeros(size(theLengths)), theLengths, theAngles)
xlim([min(theTimeFrame) max(theTimeFrame)])
ylim([-2 2])
set(gca, 'xtick', theTicks)
datetick('x', 6, 'keepticks', 'keeplimits')
ylabel('dynes/cm^2')
title('Low-passed Wind Stress')
box('on')
set(gca, 'position', thePos);
sticksafe

for indexFile = 1:length(theFiles)
    subplot(length(theFiles)+1, 1, indexFile+1)
    thePos = get(gca, 'position');
    
    ncID = netcdf(fullfile(dataDir, theFiles{indexFile}), 'nowrite');
    theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    if strcmp(lower(theSite), 'c') & indexFile == 1;
        theDepth = ncID{'depth'}(:);
        targDepth = value2Index(theDepth, 20);
        theU = ncID{'u_1205'}(:,targDepth);
        theV = ncID{'v_1206'}(:,targDepth);
    else
        theU = ncID{'u_1205'}(:);
        theV = ncID{'v_1206'}(:);
        if isempty(theU)
            theU = ncID{'east'}(:);
            theV = ncID{'north'}(:);
        end
    end
    
    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(['Low-passed Current, ' listDepths{indexFile}])
    box('on')
    set(gca, 'position', thePos);
    sticksafe
end

suptitle({['Site ' upper(theSite)]; ...
    'Low-passed Wind Stress and Current Observations'})