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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Armstrong, Brandy N., Warner, John C., List, Jeff H., Martini, 
% Marinna A., Montgomery, Ellyn T., Voulgaris, George 
% 2014, Coastal Change Processes Project data report for oceanographic 
% observations near Fire Island, New York, January to April  
% 2012 : U.S. Geological Survey Open-File Report 2014-xxxx, available 
% online at http://pubs.usgs.gov/of/2014/xxxx/.
%
% Program written in Matlab v8.0.0.755
% Program ran on PC with Windows 7 64bit
%
% "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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%B. Armstrong 6/20/2014 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:
%       julian.m
%       gregorian.m
%       PL66TN.m
%       smart_interp.m
%       tight_subplot.m
%       draw.m
%       gmean.m
%       replace.m

clear
more off

% define data and output directories
dataDirectory = '\\IGSAGIEGGS-FILE\Stg-proc\data\FireIsland_2012\candidate_final_files';
plotDirectory = 'I:\FireIsland\FI_12_DataDVD\PLOTS';

% define data file names
dataFiles={'9211wh-trm.nc','9241wh-trmbin.nc',...
            '9261aw-cal-ll-trm-single.nc','9281wh-trm.nc','9321aw-cal-ll-trm-single.nc',...
            '9203pcs-cal.nc','9205advs-cal.nc',...
            '9233pcs-cal.nc','9235advs-cal.nc','9303sgWvs-tide-trm.nc'...
            };
        names={'Site 1m','Site 2m','Site 3m','Site 4m','Site 6m','Site 1f','Site 1f','Site 2f','Site 2f','Site 5m'};
wavenameFiles={'Site 2m','Site 3m','Site 4m','Site 5m','site 6m'};
nFiles=length(dataFiles);
waveFiles={'9241whp-cal-trm.nc','9261awWvs-p-ll-trm.nc',...
            '9281whp-cal-trm.nc','9303sgWvs-p-trm.nc','9321awWvs-p-ll-trm.nc'};
% define parameters for plots
filt = 'no';
startTime=julian(2012,1,27,00);                   
stopTime=julian(2012,04,19,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=[datenum(gregorian(startTime)) (datenum(gregorian(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.5;...   
             -0.2e4  0.2e4;...
               -2.2    2.2;...
             -220  220;];
end
xt=[];
for yy=2012:2012
    for mm=1:12
        tic=julian(yy,mm,01,00);
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; datenum(gregorian(tic))];   % place xticks on the first of every month
        end
    end
end
delta_t = 900/(3600*24); %for 15 min. data
jd_int = [startTime:delta_t:stopTime]';
% Initialize pressure variables
ncVars={'P_1294','P_1','P_4023'};
vars={'pascals','dbars','mbars'};
pascals=zeros(length(jd_int),nFiles); % pressure from ADCP pressure sensor
dmean_pascals=zeros(length(jd_int),nFiles);
dbars=zeros(length(jd_int),nFiles); % pressure from ext. pressure sensor
dmean_dbars=zeros(length(jd_int),nFiles);
mbars=zeros(length(jd_int),nFiles); % pressure from ext. pressure sensor connected to ADV or PCADP
dmean_mbars=zeros(length(jd_int),nFiles);

%%
% Load site wave data
for f=1:length(waveFiles)
    ncwaves=(fullfile(dataDirectory,waveFiles{f}));
    time=double(ncread(ncwaves,'time'));
    time2=double(ncread(ncwaves,'time2'));
    jd_waves = time + (time2/3600/1000/24);
    waveDelta_t=gmean(diff(jd_waves));
    hsig=squeeze(ncread(ncwaves,'wh_4061'));
    try
        theFillValue=ncreadatt(ncwaves,'wh_4061','_FillValue');
    catch
        theFillValue=1.0E35;
    end
    hsig(hsig >= theFillValue)=nan;
    hsig(hsig <= 0)=nan;
    hsig_int(:,f) = smart_interp([jd_waves], [hsig],...
        jd_int,3);
end

%% The following loop loads the data and creates
% single site plots
sites_1=[];sites_2=[];sites_3=[];
for f=1:nFiles
    f
    clear ncobj;
    % Load the data
    nc=(fullfile(dataDirectory,dataFiles{f}));
    time=double(ncread(nc,'time'));
    time2=double(ncread(nc,'time2'));
    jd_press=time + (time2/3600/1000/24);
    delta_t = mean(diff(jd_press));
   
    for v=1:length(ncVars)
        try
        ncobj=squeeze(ncread(nc,ncVars{v}));
        data=ncobj(:);
        try
            theFillValue=ncreadatt(nc,ncVars{v},'_FillValue');
        catch
            theFillValue=1.0E35;
        end
            data=ncobj;
            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
            data_int = smart_interp([jd_press],[data],jd_int,3);       
            eval([vars{v},'(:,f) = data_int;'])
            eval(['dmean_',vars{v},'(:,f) = smart_interp([jd_press],[dmean_data],jd_int,3);']) 
            eval(['sites_',num2str(v),'=[sites_',num2str(v),';',num2str(f),'];']);   
        catch
            display ([ncVars{v} ' is not in file ' nc]);
        end
    end   
end

%% plots
clf
set(gcf,'position',[100 25 850 1100])
H = tight_subplot(4,1,[.04 .3],[.1 .1],[.1 .1]) ;

    H(1)=subplot(4,1,1,'align');
    set(H(1),'nextplot','add')
    set(H(1),'box','on')
    plot(H(1), datenum(gregorian(jd_int)), hsig_int);legend(wavenameFiles,'Location','NorthWest','Orientation','Horizontal');
    
    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
        plot(H(2), datenum(gregorian(jd_int)), dmean_pascals(:,sites_1));legend(H(2),names(sites_1),'Location','NorthWest','Orientation','Horizontal');
        plot(H(3), datenum(gregorian(jd_int)), dmean_dbars(:,sites_2));legend(H(3),names(sites_2),'Location','NorthWest','Orientation','Horizontal');
        plot(H(4), datenum(gregorian(jd_int)), dmean_mbars(:,sites_3));legend(H(4),names(sites_3),'Location','NorthWest','Orientation','Horizontal');
        %reformat xlims, ylims, and ylabels on single site plots
        for h=1:length(H)
            set(H(h),'xlim',xlims,'xtick',xt)
            set(H(h),'xticklabel',[])
            set(H(h),'ylim',ylims(h,:))
            ylabhan=get(H(h),'Ylabel');
            set(ylabhan,'String',ylab{h})
            pos=get(ylabhan,'position');
            
            % re-do xlabels
            set(H(h),'xlim',xlims,'xtick',xt)
            set(H(h),'xticklabel',{[datestr(xt,'mmm dd yyyy')]})
            
            % give it a title and export to .pdf
            if h==length(H)
                set(gcf,'PaperPositionMode','Auto')
                switch filt
                    case 'yes'
                        tit1e=get(H(1),'Title');
                        set(tit1e,'String','Low-Passed Mean-Removed Pressure Observations Fire Island 2012'); set(gcf,'PaperPositionMode','Auto')
                        print(fullfile(plotDirectory,['FI12_press_site_lp']),'-dpdf');
                    otherwise
                        tit1e=get(H(1),'Title');
                        set(tit1e,'String','Mean-Removed Pressure Observations Fire Island 2012');set(gcf,'PaperPositionMode','Auto')
                        print(fullfile(plotDirectory,['FI12_press_site']),'-dpdf')
                end
            end
        end
