function mbiwe_pres_lp
%
%   Function to create the plot of low-passed bottom pressure from
%   Sites B and C of the MBIWE, and MBLT-A.

%Soupy Alexander, 7/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.

%Full Record
theTimeFrame = [datenum(1998, 07, 29, 00, 00, 00) datenum(1998, 09, 03, 00, 00, 00)];
theTicks = min(theTimeFrame):5:max(theTimeFrame);
theCode = 6;
dataDir = 'C:\SCHTUFF\MASS_BAY\MBAY_IWAVE\REPORT\DATAFILES';

%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(4,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

siteList = {'MBLT-A'; 'B'; 'C'};

presFiles = {'5301Atcp-a1h.nc'; '5371tcp-a1h.nc'; '5351tcp-a1h.nc'};

for indexFile = 1:length(presFiles)
    ncID = netcdf(fullfile(dataDir, presFiles{indexFile}), 'nowrite');
    theTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    thePres = ncID{'P_4023'}(:);
    
    thePres(find(thePres > 1e4 | thePres < -1e4)) = NaN;
    
    thePres = thePres - nanmean(thePres);
    
    [thePresLP, theTimeLP] = plfilt(thePres, theTime, 6);
    
    subplot(4,1,indexFile+1)
    thePos = get(gca, 'position'); 
    plot(julian2datenum(theTimeLP), thePresLP)
    
    ylim([-25 25])
    xlim(theTimeFrame)
    set(gca, 'xtick', theTicks)
    datetick('x', theCode, 'keepticks', 'keeplimits')
    box('on')
    ylabel('mbar')
    
    title(['Site ' siteList{indexFile}])
    
    set(gca, 'position', thePos)
end

suptitle('Low-passed Pressure Data')