% pressure_ts.m  An mfile to plot the timeseries of North Carolina pressure
%                data.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
%  Armstrong, B.N., Warner, J.C., Voulgaris, G., List, J.H., Thieler, E.R., 
%  Martini, M.M., 2011, Carolinas Coastal Change Processes Project Data 
%  Report for Observations near Diamond Shoals, NC from January-April 2009: 
%  U.S. Geological Survey Open-File Report 2011-xxxx.
%
% Program written in Matlab v7.1.0 SP3
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% C. Sullivan 05/16/06, Pressure stored to data files as the variables
%                       hght_18 (m), P_1294 (deca-Pascals), P_1 (Db), and
%                       P_4023 (mbar). The variable is determined by the
%                       type of data file in which pressure is stored 
%                       (ie: in adv/pcadp data files, pressure is P_4023
%                       (mbar), in adcp data files, pressure is P_1294
%                       (deca-Pascals)).  
%                     
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       dolandscape.m
%       fixpaper.m

clear
more off

% define data and output directories
dataDirectory = 'Z:\data\hatteras09\BBV_reviewed';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2009DVD\mfiles';

% define data file names
dataFiles={ '8533pt-a.nc','8542wh-a.nc',...
            '8544pcs-cal.nc',...
            '8545advs-cal.nc',...
            '8551wh-a.nc','8563pt-a.nc',...
            '8572wh-a.nc',...
            '8574pcs-cal.nc',...
            '8575advs-cal.nc',...
            '8576advs-cal.nc','8581wh-a.nc'};
nFiles=length(dataFiles);
waveFiles={'8542p-cal.nc','8551p-cal.nc',...
    '8572p-cal.nc','8581p-cal.nc'};
% define parameters for plots
filt = 'no';
startTime=julian(2009,1,8,00);
stopTime=julian(2009,05,15,00);
cols=get(0,'defaultaxescolororder');     %plot color order
ylab={ {'Significant Wave','Height (m)'},...   %plot ylabels
       {'Pressure','(Pa)'},...
       {'Pressure','(dbar)'},...
       {'Pressure','(mbar)'} };
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %past start time and stop time
switch filt
    case 'yes'
        ylims=[ 0     2;...   
             -1e3   1e3;...
              -.1     .1;...
              -50    50;];
    otherwise
        ylims=[ 0     4;...   
             -0.15e4  0.15e4;...
               -2    2;...
             -200  200;];
end
xt=[];
for yy=2009:2009
    for mm=1:12
        tic=julian(yy,mm,01,00);
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic];   % place xticks on the first of every month
        end
    end
end

% Initialize pressure variables
ncVars={'P_1294','P_1','P_4023'};
vars={'pascals','dbars','mbars'};
jd_press=cell(nFiles,1);
pascals=cell(nFiles,1); % pressure from ADCP pressure sensor
dmean_pascals=cell(nFiles,1);
dbars=cell(nFiles,1); % pressure from ext. pressure sensor
dmean_dbars=cell(nFiles,1);
mbars=cell(nFiles,1); % pressure from ext. pressure sensor connected to ADV or PCADP
dmean_mbars=cell(nFiles,1);

% Load site wave data
nc=netcdf(['Z:\data\hatteras09\candidate_final_files\',waveFiles{2},], 'nowrite');
jd_waves{1} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
waveDelta_t=gmean(diff(jd_waves{1}));
hsig{1}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{1} >= theFillValue);
hsig{1}(badInd)=nan;
nc=close(nc);

nc=netcdf(['Z:\data\hatteras09\candidate_final_files\',waveFiles{3},], 'nowrite');
jd_waves{2} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
waveDelta_t=gmean(diff(jd_waves{2}));
hsig{2}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{2} >= theFillValue);
hsig{2}(badInd)=nan;
nc=close(nc);

% The following loop loads the data and creates
% single site plots
for f=1:nFiles
    f
    % Load the data
    nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
    jd_press{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
    delta_t = mean(diff(jd_press{f}));
    for v=1:length(ncVars)
        ncobj=nc{ncVars{v}};
        if ~isempty(ncobj)
            data=ncobj(:);
            theFillValue=ncobj.FillValue_(:);
            bads=find(data >= theFillValue);
            data(bads)=nan;
            dmean_data = data - gmean(data);
            switch filt
                case 'yes'
                    dmean_data=pl66tn(dmean_data,delta_t*24,33);
            end
            eval([vars{v},'{f} = data;'])
            eval(['dmean_',vars{v},'{f} = dmean_data;'])
        else
            eval([vars{v},'{f} = nan(size(jd_press{f}));'])
            eval(['dmean_',vars{v},'{f} = nan(size(jd_press{f}));'])
        end
    end
    nc=close(nc)
    
end

%% Site plots
% Define the site and color
for  st=1:3
    % set up single site plots
    F=figure(f);
    dolandscape
    fixpaper

    H(1)=subplot(4,1,1,'align');
    set(H(1),'nextplot','add')
    set(H(1),'box','on')
    if st==1
    plot(H(1), jd_waves{1}, hsig{1})
    elseif st==2
    plot(H(1), jd_waves{2}, hsig{2})
    else
        hsig_wavestip=zeros(length(jd_waves{1}),1);
        hsig_wavestip(:)=NaN;
    plot(H(1), jd_waves{1}, hsig_wavestip)        
    end
    H(2)=subplot(4,1,2,'align');
    set(H(2),'nextplot','add')
    set(H(2),'box','on')
    H(3)=subplot(4,1,3,'align');
    set(H(3),'nextplot','add')
    set(H(3),'box','on')
    H(4)=subplot(4,1,4,'align');
    set(H(4),'nextplot','add')
    set(H(4),'box','on')

    % plot to single site plots
    col={'r','b','g','c','m','r','b','g','c','m','b','r','g','c','m'};
    if st==1
    for f=1:5
        figure(F)
        plot(H(2), jd_press{f}, dmean_pascals{f}, col{f})
        plot(H(3), jd_press{f}, dmean_dbars{f}, col{f})
        plot(H(4), jd_press{f}, dmean_mbars{f}, col{f})
    end
    end
if st==2
    for f=6:10
        figure(F)
        plot(H(2), jd_press{f}, dmean_pascals{f}, col{f})
        plot(H(3), jd_press{f}, dmean_dbars{f}, col{f})
        plot(H(4), jd_press{f}, dmean_mbars{f}, col{f})
    end
end
    if st==3
    for f=11
        figure(F)
        plot(H(2), jd_press{f}, dmean_pascals{f}, col{f})
        plot(H(3), jd_press{f}, dmean_dbars{f}, col{f})
        plot(H(4), jd_press{f}, dmean_mbars{f}, col{f})
    end
    end
    %reformat xlims, ylims, and ylabels on single site plots
    for h=1:length(H)
        axes(H(h))
        set(gca,'xlim',xlims,'xtick',xt)
        set(gca,'xticklabel',[])
        set(gca,'ylim',ylims(h,:))
        ylabhan=get(gca,'Ylabel');
        set(ylabhan,'String',ylab{h})
        pos=get(ylabhan,'position');
        if h==length(H)
            if st==1
                    site='North Site';
                    text(julian(2009,3,1),230,'ADV and PCADP (Near Bottom)','horizontalalignment','center','Parent',H(4))
                    text(julian(2009,3,1),2.3,'PT (Near Surface)','horizontalalignment','center','Parent',H(3))
                    text(julian(2009,3,1),1900,'ADCP (Near Bottom)','horizontalalignment','center','Parent',H(2))
            elseif st== 2 %South site
                    site='South Site';
                    text(julian(2009,3,1),230,'ADV and PCADP (Near Bottom)','horizontalalignment','center','Parent',H(4))
                    text(julian(2009,3,1),2.3,'PT (Near Surface)','horizontalalignment','center','Parent',H(3))
                    text(julian(2009,3,1),1900,'ADCP (Near Bottom)','horizontalalignment','center','Parent',H(2))
            elseif st== 3 %Tip site
                    site='Tip Site';
                    text(julian(2009,3,1),1900,'ADCP (Near Bottom)','horizontalalignment','center','Parent',H(2))
            end
        end

        % re-do xlabels
        axes(H(end))
        oldhan=findobj(gca,'Tag','xtick labels');
        delete(oldhan);
        for tic=1:length(xt)
            xtg=gregorian(xt(tic));
            xtd=datenum(xtg);
            xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)};
            text(xt(tic),-210,xtl,'horizontalalignment','center',...
                'verticalalignment','top')
        end

        % give it a title and export to .pdf
        if h==length(H)
            set(gcf,'PaperPositionMode','Auto')
            switch filt
                case 'yes'
                    suptitle({'Low-Passed Mean-Removed Pressure Observations ',site})
                    print(fullfile(plotDirectory,['press_site',site,'_lp']),'-dpdf')
                otherwise
                    suptitle({'Mean-Removed Pressure Observations ',site})
                    print(fullfile(plotDirectory,['press_site',site]),'-dpdf')
            end
        end
    end
end
