% pcadp_ts.m  An mfile to plot the (optionally filtered)
%             timeseries of current velocity data from South
%             Carolina PCADP data.

% C. Sullivan 03/14/06, Raw PCADP data is processed and saved to netCDF files.
%                       It is assumed that the netCDF files are stored in a
%                       directory defined by the variable 'dataDirectory'.
%                       To plot low-pass filtered current velocities set the
%                       switch 'filt' to 'yes', otherwise (filt='no')
%                       unfiltered current velocities are plotted.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Armstrong, B.N., Warner, J.C., Voulgaris, G., List, J.H., Thieler, E.R., 
% Martini, M.M., Montgomery, E.T., 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-1156.
%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Dependencies:
%       the netcdf toolbox
%       smart_interp.m
%       PL66TN.m
%       timeplt.m
%       dolandscape.m
%       suptitle.m
%       fixpaper.m

clear
more off

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

% define ADV data filenames
dataFiles={'8544pcs-cal.nc','8574pcs-cal.nc'};
waveFiles={'8542p-cal.nc','8572p-cal.nc'};
% dataFiles={'8574pcs-cal.nc'};
% waveFiles={'8572p-cal.nc'};
nFiles=length(dataFiles);

% define parameters for plots
filt='no';                                         %filtering switch, 'yes' or 'no'
col=get(0,'defaultaxescolororder');                 %plot color order
ylab={ {'Significant Wave','Height (m)'},...
    {'Sensor Range','to Boundary','(m)'},...
    {'Northward','Velocity','(cm/s)'},...
    {'Eastward','Velocity','(cm/s)'},...
    {'Northward','Velocity STD','(cm/s)'},...
    {'Eastward','Velocity STD','(cm/s)'} };
startTime=julian(2009,1,8,00);
stopTime=julian(2009,05,15,00);
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %beyond start time and stop time
ylims=[0 4; 0.5 1.2; -25 55; -40 40; 0 45; 0 45;];      %u,v plot ylimits
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];
        end
    end
end

% 'jd_int' is the time array onto which all
% data is interpolated. data gaps larger than
% 'maxGap' hours are not interpolated over
delta_t = 1/24;
jd_int = (startTime:delta_t:stopTime)';
nTime = length(jd_int);
maxGap = 3;

% filter the waves (optional)
if strcmp(filt,'yes')
    hsig_all = PL66TN(hsig_all,waveDelta_t*24,33);
end

% initialize variables for PCADP data
jd=cell(nFiles,1);
depth=cell(nFiles,1);
u=cell(nFiles,1);
v=cell(nFiles,1);
ustd=cell(nFiles,1);
vstd=cell(nFiles,1);
brange=cell(nFiles,1);
press=cell(nFiles,1);

% define PCADP variables
ncVars={'u_1205','v_1206','USTD_4097','VSTD_4098','brange','P_4023'};
vars={'u','v','ustd','vstd','brange','press'};
nVars=length(vars);

% load the PCADP data
for f=1:nFiles
    site=dataFiles{f}(1:4);
    nc=netcdf(['X:\hatteras09\candidate_final_files\',waveFiles{f},], 'nowrite');
    jd_waves_all{f} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
    waveDelta_t=gmean(diff(jd_waves_all{f}));
    hsig{f}=nc{'wh_4061'}(:);
    theFillValue=nc{'wh_4061'}.FillValue_(:);
    badInd=find(hsig{f} >= theFillValue);
    hsig{f}(badInd)=nan;
    badInd=find(hsig{f} <= 0);
    hsig{f}(badInd)=nan;
    nc=close(nc);
    hsig_all = smart_interp([jd_waves_all{f}], [hsig{f}],...
        jd_int, 3);

    nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
    jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
    delta_t = gmean(diff(jd{f}));
    depth{f}=nc{'depth'}(:);
    for n=1:nVars
        ncobj=nc{ncVars{n}};
        theFillValue=ncobj.FillValue_(:);
        data=ncobj(:);
        badInd=find(data>=theFillValue);
        data(badInd)=nan;
        if strcmp(vars{n},'press')
            % de-mean the pressure signal
            data=(data-1000)/100;
            data=data-gmean(data);
        end
        % filter the data (optional)
        switch filt
            case 'yes'
                [nRows,nCols]=size(data);
                for col=1:nCols
                    data(:,col)=PL66TN(data(:,col), delta_t*24, 33);
                end
            otherwise
                disp('Not filtering current velocities');
        end
        eval([vars{n},'{f}=data;'])
    end
    nc=close(nc);



    % combine data from deployments 1 and 2
    jd_all=[jd{f}];
    depth_all=[depth{f}];
    u_all=[u{f}];
    v_all=[v{f}];
    ustd_all=[ustd{f}];
    vstd_all=[vstd{f}];
    brange_all=[brange{f}];
    press_all=[press{f}];

    %%%%%%%%%%%%%%%%%%%%%%%%%% Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure
    dolandscape
    fixpaper

    Han(1)=subplot(6,1,1,'align');
    plot(jd_int, hsig_all)
    hold on
    Han(2)=subplot(6,1,2,'align');
    plot(jd_all,brange_all(:,1))
    hold on
    plot(jd_all,brange_all(:,2),'r')
    plot(jd_all,brange_all(:,3),'k')
    LEG2=legend('Beam 1','Beam 2','Beam 3');
    set(LEG2,'orientation','horizontal')
    set(LEG2,'location','NorthEast')
    set(LEG2,'box','off')

    % to contour the velocity bins ...
    x=repmat(jd_all,1,length(depth_all));
    z=repmat(depth_all',length(jd_all),1);
    x2=x';
    z2=z';
    [z2_sort,sortInd]=sortrows(z2,1);
    x2_sort=x2(sortInd,:);
    x3=x2_sort';
    z3=z2_sort';
    y=nan(size(x));
    for n=1:4
        eval(['data1 = ',vars{n},'{1};'])
        y(1:length(jd{1}),1:length(depth{1})) = data1;
        y2=y';
        y2_sort=y2(sortInd,:);
        y3=y2_sort';
        Han(n+2)=subplot(6,1,n+2,'align');
        for ii=1:length(depth_all)
            [Ct,Ht]=contour(x3, y3, z3, [z3(1,ii) z3(1,ii)]);
            hold on
            C(ii)=Ct(1);
            H(ii)=Ht(1);
        end
    end

    % reformat xlims, ylabels, lines
    for h=1:length(Han)
        axes(Han(h))
        axpos=get(gca,'position');
        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([2008,12,30,0,0,0]) pos(2) pos(3)])
    end

    % make a legend and reformat xlabels
    LEG=legend(H,num2str(round(C'*100)/100,5));
    figCol=get(gcf,'color');
    set(LEG,'color',figCol)
    set(LEG,'edgecolor',figCol)
    legpos=get(LEG,'position');
    axpos=get(gca,'position');
    set(LEG,'position',[axpos(1)+axpos(3)+.003 axpos(2) legpos(3) legpos(4)])
    set(LEG,'XAxisLocation','top')
    LEGLBL=get(LEG,'Xlabel');
    set(LEGLBL,'String','Depth (mbs)')
    set(LEGLBL,'Color',[0 0 0],'fontweight','bold')
    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)};
        xtlpos=ylims(1,1)-1.5;
        text(xt(tic),xtlpos,xtl,'horizontalalignment','center',...
            'verticalalignment','top')
    end

    % give it a title and export to .pdf
    set(gcf,'paperpositionmode','auto')
    axes(Han(1))
    switch filt
        case 'yes'
            suptitle(['PCADP Low-PAss Filtered Bottom Water Flow Observations Site ',site])%,'horizontalalignment','center','fontsize',14);
            print(fullfile(plotDirectory,['pcadp_site',site,'_lp']),'-dpdf');
        otherwise
            suptitle(['PCADP Bottom Water Flow Observations Site ',site])%,'horizontalalignment','center','fontsize',14);
            print(fullfile(plotDirectory,['pcadp_site',site]),'-dpdf');
    end
end
