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

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% 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_4 (Pascals), P_1 (dbar), 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_4
%                       (pascals)).  
%                     
% Dependencies:
%       the netcdf toolbox
%       julian.m
%       dolandscape.m
%       fixpaper.m

clear
more off

% define data and output directories
dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\PRESS';

% define data file names
dataFiles={ '7201adc-a.nc','7401adc-a.nc',...
            '7221adc-a.nc','7421adc-a.nc',...
            '7223pt-a.nc','7423pt-a.nc',...
            '7241adc-a.nc','7441adc-a.nc',...
            '7243pt-a.nc','7443pt-a.nc',...
            '7331adc-a.nc','7531adc-a.nc',...
            'adv7261vp-cal2.nc','adv7461vp-cal2.nc',...
            'adv7262vp-cal2.nc','adv7462vp-cal2.nc',...
            '7263arc-a.nc','7463arc-a.nc',...
            'adv7281vp-cal.nc','adv7481vp-cal.nc',...
            '728pcvp-cal5.nc','748pcvp-cal6.nc',...
            '7321adc-a.nc','7521wh-a.nc',...
            '7311aqc-a.nc','7511aqc-a.nc' };
nFiles=length(dataFiles);

% define parameters for plots
filt = 'yes';
startTime=julian(2003,10,25,00);        %plot start stime
stopTime=julian(2004,04,25,00);         %plot stop time  
cols=get(0,'defaultaxescolororder');     %plot color order
cols(7,:)=[.6 .6 .6];
cols(8,:)=[0 0 0];
ylab={ {'Site 2','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;...   
             -5e3   5e3;...
              -.5     .5;...
              -50    50;];
    otherwise
        ylims=[ 0     2;...   
             -2e4  2e4;...
               -2    2;...
             -200  200;];
end
xt=[];
for yy=2003:2004
    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_4','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 2 wave data
nc=netcdf(fullfile(dataDirectory,'7221adwp-cal.nc'), '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;
switch filt
    case 'yes'
        hsig{1}=pl66tn(hsig{1},1,33);
end
nc=close(nc);

nc=netcdf(fullfile(dataDirectory,'7521adwp-cal.nc'), 'nowrite');
jd_waves{2} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
hsig{2}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{2} >= theFillValue);
hsig{2}(badInd)=nan;
switch filt
    case 'yes'
        hsig{2}=pl66tn(hsig{2},1,33);
end
nc=close(nc);

% Combine wave data from deployments 1 and 2
jd_waves_all=[jd_waves{1}; nan; jd_waves{2}];
hsig_all=[hsig{1}; nan; hsig{2}];

% The following loop loads the data and creates
% single site plots
for f=1:nFiles
    
    % 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);
    
    % plots
    if ismember(f,[1 3 7 11 13 19 23 25])
        % Define the site and color
        switch f
            case 1
                site='1'; col=cols(1,:);
            case 3
                site='2'; col=cols(2,:);
            case 7
                site='3'; col=cols(3,:);
            case 11
                site='4'; col=cols(4,:);
            case 13
                site='6'; col=cols(5,:);
            case 19
                site='7A'; col=cols(6,:);
            case 23
                site='7B'; col=cols(7,:);
            case 25
                site='8'; col=cols(8,:);
        end
        
        % 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')
        plot(H(1), jd_waves{1}, hsig{1})
        plot(H(1), jd_waves{2}, hsig{2})
        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')
        
        if f==1
            % set up multi-site plot
            FH=figure;
            dolandscape
            fixpaper

            HA(1)=subplot(4,1,1,'align');
            set(HA(1),'nextplot','add')
            set(HA(1),'box','on')
            plot(HA(1), jd_waves{1}, hsig{1})
            plot(HA(1), jd_waves{2}, hsig{2})
            HA(2)=subplot(4,1,2,'align');
            set(HA(2),'nextplot','add')
            set(HA(2),'box','on')
            HA(3)=subplot(4,1,3,'align');
            set(HA(3),'nextplot','add')
            set(HA(3),'box','on')
            HA(4)=subplot(4,1,4,'align');
            set(HA(4),'nextplot','add')
            set(HA(4),'box','on')
        end
    end
    
    % plot to single site plots
    figure(F)
    plot(H(2), jd_press{f}, dmean_pascals{f}, 'color', col)
    plot(H(3), jd_press{f}, dmean_dbars{f}, 'color', col)
    plot(H(4), jd_press{f}, dmean_mbars{f}, 'color', col)
    
    % plot to multi-site plot
    figure(FH)
    P(f)=plot(HA(2), jd_press{f}, dmean_pascals{f}, 'color', col);
    P(f)=plot(HA(3), jd_press{f}, dmean_dbars{f}, 'color', col);
    P(f)=plot(HA(4), jd_press{f}, dmean_mbars{f}, 'color', col);
    
    if ismember(f,[2 6 10 12 18 22 24 26])
        %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');
            set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)])
            if h==length(H)
            switch f
                case 2 %site 1
                    text(julian(2004,2,1),0,'No ADV or PCADP','horizontalalignment','center','Parent',H(4))
                    text(julian(2004,2,1),0,'No external pressure sensor','horizontalalignment','center','Parent',H(3))
                    text(julian(2003,12,1),0,'No Data','horizontalalignment','center','Parent',H(2))
                case 6 %site 2
                    text(julian(2004,2,1),0,'No ADV or PCADP','horizontalalignment','center','Parent',H(4))
                case 10 %site 3
                    text(julian(2004,2,1),0,'No ADV or PCADP','horizontalalignment','center','Parent',H(4))
                    text(julian(2003,12,1),0,'No Data','horizontalalignment','center','Parent',H(2))
                case 12 %site 4
                    text(julian(2004,2,1),0,'No ADV or PCADP','horizontalalignment','center','Parent',H(4))
                    text(julian(2004,2,1),0,'No external pressure sensor','horizontalalignment','center','Parent',H(3))
                case 18 %site 6
                    text(julian(2004,2,1),0,'No external pressure sensor','horizontalalignment','center','Parent',H(3))
                case 22 %site 7A
                    text(julian(2004,2,1),0,'No external pressure sensor','horizontalalignment','center','Parent',H(3))
                    text(julian(2004,2,1),0,'No ADCP','horizontalalignment','center','Parent',H(2))
                case 24 %site 7B
                    text(julian(2004,2,1),0,'No ADV or PCADP','horizontalalignment','center','Parent',H(4))
                    text(julian(2004,2,1),0,'No external pressure sensor','horizontalalignment','center','Parent',H(3))
                    text(julian(2003,12,1),0,'No Data','horizontalalignment','center','Parent',H(2))
                case 26
                    text(julian(2004,2,1),0,'No ADV or PCADP','horizontalalignment','center','Parent',H(4))
                    text(julian(2004,2,1),0,'No external pressure sensor','horizontalalignment','center','Parent',H(3))
            end
            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
        set(gcf,'PaperPositionMode','Auto')
        switch filt
            case 'yes'
                suptitle({'Low-Passed Mean-Removed Pressure Observations',['Site ',site]})
                print(fullfile(plotDirectory,['press_site',site,'_lp']),'-dpdf')
            otherwise
                suptitle({'Mean-Removed Pressure Observations',['Site ',site]})
                print(fullfile(plotDirectory,['press_site',site]),'-dpdf')
        end
        
        if f==26
           %reformat xlims, ylims, and ylabels on multi-site plot
           for h=1:length(HA)
               axes(HA(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');
               set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)])
               if h==1
               LEG=legend(P([2 6 10 12 18 22 24 26]),'Site 1','Site 2','Site 3','Site 4',...
                                 'Site 6','Site 7A','Site 7B','Site 8');
               set(LEG,'Location','NorthOutside')
               set(LEG,'Orientation','Horizontal')
               figCol=get(gcf,'color');
               set(LEG,'color',figCol)
               set(LEG,'edgecolor',figCol)
               pos=get(LEG,'position');
               set(LEG,'position',[pos(1) .9 pos(3) pos(4)])

               end
           end
           % re-do xlabels
           axes(HA(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
           set(gcf,'PaperPositionMode','Auto')
           switch filt
            case 'yes'
                suptitle({'Low-Passed Mean-Removed Pressure Observations','All Sites'})
                print(fullfile(plotDirectory,'press_all_lp'),'-dpdf')
            otherwise
               suptitle({'Mean-Removed Pressure Observations','All Sites'})
               print(fullfile(plotDirectory,'press_all'),'-dpdf')
           end
        end
        
    end
        
end