function pres_plot(theYears, dataDirectory, plotDirectory, siteA, siteB)
%
%   pres_plot(theYears, dataDirectory, plotDirectory, siteA, siteB)
%
%   Function to plot the low-passed wind stress and pressure observations
%   from Sites A and B.  Output files are named SITE_pres_YYYY.

%Soupy Alexander, 10/30/03

%Load in the wind stress
[windTime, tauU, tauV] = mblt_windstress(dataDirectory);

presTimeA = [];
presDataA = [];
presTimeA_lp = [];
presDataA_lp = [];

moorNum = {};

%Load in the site A data
for indexFile = 1:length(siteA)
    fileVars = siteA(indexFile).theVars;
    fileName = siteA(indexFile).fileList;
    fileMoor = fileName(1:4);
    if ~strcmp(fileVars, 'P_4023') %| ~isempty(strmatch(fileMoor, moorNum))
        continue
    end
    moorNum = [moorNum; {fileMoor}];
    presFileA = netcdf(fullfile(dataDirectory, siteA(indexFile).fileList), 'nowrite');
    fileTime = singlejd(presFileA{'time'}(:), presFileA{'time2'}(:));
    fileData = presFileA{'P_4023'}(:);  
    bads = find(fileData > 1e5);
    fileTime(bads) = [];
    fileData(bads) = [];
    fileData = fileData - nanmean(fileData);
    [fileData_lp, fileTime_lp] = plfilt(fileData, fileTime, 6);
    
    presTimeA = [presTimeA(:); fileTime(:)];
    presDataA = [presDataA(:); fileData(:)];
    presTimeA_lp = [presTimeA_lp(:); fileTime_lp(:)];
    presDataA_lp = [presDataA_lp(:); fileData_lp(:)];
end

presTimeB = [];
presDataB = [];
presTimeB_lp = [];
presDataB_lp = [];

%Load in the site B data
for indexFile = 1:length(siteB)
    fileVars = siteB(indexFile).theVars;
    if ~strcmp(fileVars, 'P_4023')
        continue
    end
    presFileB = netcdf(fullfile(dataDirectory, siteB(indexFile).fileList), 'nowrite');
    fileTime = singlejd(presFileB{'time'}(:), presFileB{'time2'}(:));
    fileData = presFileB{'P_4023'}(:);
    bads = find(fileData > 1e5);
    fileTime(bads) = [];
    fileData(bads) = [];
    fileData = fileData - nanmean(fileData);
    [fileData_lp, fileTime_lp] = plfilt(fileData, fileTime, 6);
        
    presTimeB = [presTimeB(:); fileTime(:)];
    presDataB = [presDataB(:); fileData(:)]; 
    presTimeB_lp = [presTimeB_lp(:); fileTime_lp(:)];
    presDataB_lp = [presDataB_lp(:); fileData_lp(:)];
end

%Interpolate and low-pass the pressure data
newTimeA = min(presTimeA):1/24:max(presTimeA);
[presDataA_i] = smart_interp_pres(presTimeA, presDataA, newTimeA, 1.5);
%newTimeA_lp = min(presTimeA_lp):6/24:max(presTimeA_lp);
%[presDataA_lp] = smart_interp_pres(presTimeA_lp, presDataA_lp, newTimeA_lp, 6.5);

[presDataA_lp, newTimeA_lp] = plfilt(presDataA_i, newTimeA, 6);

newTimeB = min(presTimeB):1/24:max(presTimeB);
[presDataB_i] = smart_interp_pres(presTimeB, presDataB, newTimeB, 1.5);
newTimeB_lp = min(presTimeB_lp):6/24:max(presTimeB_lp);
[presDataB_lp] = smart_interp_pres(presTimeB_lp, presDataB_lp, newTimeB_lp, 6.5);


%Do the plots
for indexYear = 1:length(theYears)
    targetYear = theYears(indexYear);
    figure, orient landscape
    
    xTick = [];
    for indexMonth = 1:12
        xTick = [xTick(:); datenum(targetYear, indexMonth, 1)];
    end
    
    inYear = find((windTime >= julian(targetYear, 1, 1, 00)) & ...
        windTime <= julian(targetYear, 12, 31, 24));
    windTimeIn = windTime(inYear);
    
    [theAngleRad, theLength] = cart2pol(tauU(inYear), tauV(inYear));
    theAngle = theAngleRad * 180/pi;
    
    goods = find(~isnan(theAngle));
    theAngle = theAngle(goods);
    theLength = theLength(goods);
    
    plotX = julian2datenum(windTimeIn(goods));
    plotY = zeros(size(plotX));
    
    subplot(3,1,1)
    sticksafe(plotX, plotY, theLength, theAngle)
    ylim([-4 4])
    set(gca, 'xtick', xTick, 'ytick', [-4:2:4])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('dynes/cm^2')
    title('Low-passed Wind Stress')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    sticksafe
    
    subplot(3,1,2)
    plot(julian2datenum(newTimeA), presDataA_i)
    ylim([-300 300])
    set(gca, 'xtick', xTick, 'ytick', [-300:100:300])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('mb')
    title('Hourly Averaged Pressure')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    subplot(3,1,3)
    plot(julian2datenum(newTimeA_lp), presDataA_lp)
    ylim([-100 100])
    set(gca, 'xtick', xTick, 'ytick', [-200:25:200])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('mb')
    title('Low-passed Pressure')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    suptitle({'Low-passed Wind Stress and Pressure Observations'; ['Site A, ' ...
                num2str(targetYear, '%.0f')]});
    print(fullfile(plotDirectory, ['a_pres_' num2str(targetYear, '%.0f')]), '-dpdf')
    close
    
    %Site B
    if targetYear ~= 2000
        continue
    end
    figure, orient landscape
    subplot(3,1,1)
    sticksafe(plotX, plotY, theLength, theAngle)
    ylim([-4 4])
    set(gca, 'xtick', xTick, 'ytick', [-4:2:4])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('dynes/cm^2')
    title('Low-passed Wind Stress')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    sticksafe
    
    subplot(3,1,2)
    plot(julian2datenum(newTimeB), presDataB_i)
    ylim([-300 300])
    set(gca, 'xtick', xTick, 'ytick', [-300:50:300])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('mb')
    title('Hourly Averaged Pressure')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    subplot(3,1,3)
    plot(julian2datenum(newTimeB_lp), presDataB_lp)
    ylim([-50 50])
    set(gca, 'xtick', xTick, 'ytick', [-200:25:200])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('mb')
    title('Low-passed Pressure')
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    suptitle({'Low-passed Wind Stress and Pressure Observations'; ['Site B, ' ...
                num2str(targetYear, '%.0f')]});
    print(fullfile(plotDirectory, ['b_pres_' num2str(targetYear, '%.0f')]), '-dpdf')
    close
end