% 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., 2012, Carolinas Coastal Change Processes Project Data 
%  Report for Observations near Diamond Shoals, NC from February 2010: 
%  U.S. Geological Survey Open-File Report 2012-xxxx.
%
% Program written in Matlab 2012b v8.0.0 
% Program ran on PC with Windows 7 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% B.Armstrong 9/7/2012, Pressure stored to data files as the variables
%                       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

addpath('D:\BNA\Matlab\m_cmg\branches\barmstrong\BNA\t_tide_v1.3beta','-end');

% Define data and output directories
dataDirectory ='\\IGSAGIEGGS-FILE\Stg-proc\data\HatterasNearshore10\BBV_reviewed';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\PLOTS';

% Define adcp data filenames.
plot_location=2;%1=north,2=south,3=offshore
if plot_location==1
    %current file names
    dataFiles={'8771aqd-a.nc','8781aqd-a.nc','8791aqd-a.nc',...
        '8801wh-a.nc','8641wh-a.nc','8651aqd-cal.nc'};
    
    %wave file names
    waveFiles={'8771aqwp-cal.cdf','8781aqdwp-cal.nc','8791aqwp-cal.cdf',...
        '8801adwp-cal.cdf','8641adwp-cal.cdf','8651aqwp-cal.cdf'};
    
    %names of stations
    names={'N1','N2','N3','N4','N5','N6'};
    
    %location
    location='North Sites';
elseif plot_location==2
    %current file names
    dataFiles={'8751aqd-a.nc','8711wh-a.nc','8661aqd-cal.nc',...
        '8671Aaqd-cal.nc','8671Baqd-cal.nc'};
    
    %wave file names
    waveFiles={'8751aqwp-cal.cdf','8711whp-cal.cdf','8661aqwp-cal.cdf',...
        '8671Aaqwp-cal.cdf','8671Baqwp-cal.cdf'};
    
    %names of stations
    names={'N8','N9','N11','N12A','N12B'};
    
    %location
    location='South Sites';
else %plot location ==3
    %current file names
    dataFiles={'8721wh-a.nc','8731wh-a.nc','8741wh-a.nc'};
    
    %wave file names
    waveFiles={'8721adwp-cal.nc','8731adwp-cal.cdf','8741Aadwp-cal.cdf','8741Badwp-cal.cdf'};
    
    %names of stations
    names={'O1','O2','O3'};
    
    %location
    location='Offshore Sites';
end
filt='yes';
nFiles=length(dataFiles);
%% define parameters for plots
if plot_location==3
    startTime=julian(2010,2,5,00);
    stopTime=julian(2010,3,24,00);
else
    startTime=julian(2010,2,5,00);
    stopTime=julian(2010,02,23,00);
end
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

timediff=stopTime-startTime;
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)];
xt=[];
for yy=2010:2010
    for mm=2:2
        for dd=1:2:28
            tic=julian(yy,mm,dd,00);
            if tic>=xlims(1) & tic<=xlims(2)
                xt=[xt; tic];
            end
        end
    end
end
% larger than 'maxGap' hours are not interpolated over
delta_t = 900/(3600*24); %for 15 min. data
jd_int = [startTime:delta_t:stopTime]';
maxGap = 1; %hr

% Get the color scheme to distinguish between data sites
col=get(0,'defaultaxescolororder');

%% Initialize variables
hsig_int=zeros(length(jd_int),nFiles);

% Initialize pressure variables
ncVars={'P_1294','P_1','P_4023'};
vars={'pascals','dbars','mbars'};
for var_name=1:length(vars);
    eval([vars{var_name} '=nan(length(jd_int),nFiles);'])
    eval(['dmean_' vars{var_name} '=nan(length(jd_int),nFiles);'])
end
%% The following loop loads the data and creates plots
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'));
    theFillValue=ncreadatt(ncwaves,'wh_4061','_FillValue');
    hsig(hsig >= theFillValue)=nan;
    hsig(hsig <= 0)=nan;
    hsig_int(:,f) = smart_interp([jd_waves], [hsig],...
        jd_int, 3);
end


for f=1:length(names)
    ncobj=[];  
    nc=(fullfile(dataDirectory,dataFiles{f}));
    % Load the data time=double(ncread(nc,'time'));
    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}));
            if ~isempty(ncobj)
                theFillValue=ncreadatt(nc,ncVars{v},'_FillValue');
                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
                display([nc ' does contain ' ncVars{v} ' in ' vars{v}]);
                eval([vars{v},'(:,f) = smart_interp([jd_press], [data],jd_int, 3);']);
                eval(['dmean_',vars{v},'(:,f) =  smart_interp([jd_press], [dmean_data],jd_int, 3);']);
            else
                eval([vars{v},'(:,f) = nan(size(jd_int));'])
                eval(['dmean_',vars{v},'(:,f) = nan(size(jd_int));'])
            end
        catch
            display([nc ' does not contain ' ncVars{v}]);
        end
    end
end
time_dn=datenum(gregorian(jd_int));
%%
if plot_location==3
corder=get(gca,'ColorOrder');
corder2=zeros(11,3);
corder2(1,:)=corder(1,:);
corder2(2,:)=corder(2,:);
corder2(3,:)=corder(3,:);
corder2(4,:)=corder(3,:);
end

clf
set(gcf,'position',[100 255 800 700])
ha = tight_subplot(3,1,[.04 .3],[.1 .1],[.1 .1]) ;

axes(ha(3));plot(time_dn,dmean_dbars,'LineWidth',1.2);
datetick('x',2,'keeplimits','keepticks');
ylabel(ylab{3});

axes(ha(1));
if plot_location==3
    set(gca,'NextPlot','replacechildren');
    set(ha(:),'ColorOrder',corder2);
end
plot(time_dn,hsig_int,'LineWidth',1.2);
if plot_location==3
set(gca,'YTickLabel',[0 1 2 3 4]);
end
plot(time_dn,hsig_int,'LineWidth',1.2);
legend(names,'FontSize',7);
set(gca,'xlim',xlim,'xtick',get(ha(3),'xtick'),'XTickLabel','');
ylabel(ylab{1});

axes(ha(2));plot(time_dn,dmean_pascals,'LineWidth',1.2);
set(gca,'xlim',xlim,'xtick',get(ha(3),'xtick'),'XTickLabel','');
ylabel(ylab{2});
tit1e=get(ha(1),'Title');
set(tit1e,'String',['Mean Pressure Observations ',location])
set(gcf,'PaperPositionMode','Auto','PaperOrientation','landscape')

switch filt
    case 'yes'
        print(fullfile(plotDirectory,[location,'_filt_press_mn']),'-dpdf');
    otherwise
        print(fullfile(plotDirectory,[location,'_press_mn']),'-dpdf');
end
%%

clf
set(gcf,'position',[100 255 800 700])
ha = tight_subplot(3,1,[.04 .3],[.1 .1],[.1 .1]) ;

axes(ha(3));plot(time_dn,dbars,'LineWidth',1.2);
datetick('x',2,'keeplimits','keepticks');
ylabel(ylab{3});

axes(ha(1));
if plot_location==3
    set(gca,'NextPlot','replacechildren');
    set(ha(:),'ColorOrder',corder2);
end
plot(time_dn,hsig_int,'LineWidth',1.2);
if plot_location==3
set(gca,'YTickLabel',[0 1 2 3 4]);
end
plot(time_dn,hsig_int,'LineWidth',1.2);
legend(names,'FontSize',7);
set(gca,'xlim',xlim,'xtick',get(ha(3),'xtick'),'XTickLabel','');
ylabel(ylab{1});


axes(ha(2));plot(time_dn,pascals,'LineWidth',1.2);
set(gca,'xlim',xlim,'xtick',get(ha(3),'xtick'),'XTickLabel','');
ylabel(ylab{2});

tit1e=get(ha(1),'Title');
set(tit1e,'String',['Pressure Observations ',location])
set(gcf,'PaperPositionMode','Auto','PaperOrientation','landscape')
print(fullfile(plotDirectory,[location,'_press']),'-dpdf');