% adv_ts.m  An mfile to plot the (optionally filtered)
%           timeseries of current velocity data from South
%           Carolina ADV 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., 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% C. Sullivan 03/18/05, Raw ADV 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.
%
% Dependencies:
%       the netcdf toolbox
%       smart_interp.m
%       PL66TN.m
%       timeplt.m
%       dolandscape.m
%       suptitle.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={'8545advs-cal.nc','8575advs-cal.nc','8576advs-cal.nc'};
waveFiles={'8542p-cal.nc','8572p-cal.nc','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={ {'Northward','Velocity STD','(cm/s)'},...          %plot ylabels for stick plots
    {'Eastward','Velocity STD','(cm/s)'},...
    {'Northward','Velocity','(cm/s)'},...
    {'Eastward','Velocity','(cm/s)'},...
    {'Sensor Range','to Boundary','(m)'},...
    {'Significant Wave','Height (m)'} };
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
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
ylims=[0 50; 0 50; -80 80; -80 80; 0 .6; 0 5];      %u,v plot ylimits

% '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]';
maxGap = 3;

% initialize variables
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);
varList={'u','v','ustd','vstd','brange','pressure'};
nVars=length(varList);

for f=1:length(dataFiles)
    nc=netcdf(['X:\hatteras09\candidate_final_files\',waveFiles{f},], '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;
    badInd=find(hsig{1} <= 0);
    hsig{1}(badInd)=nan;
    nc=close(nc);
    hsig_all_comb = smart_interp([jd_waves{1}], [hsig{1}],...
        jd_int, 3);

    if ~isempty(dataFiles{f})
        nc=netcdf(fullfile(dataDirectory,dataFiles{f}), 'nowrite');
        jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
        u{f}=nc{'u_1205'}(:);                          %eastward velocity, cm/s
        ustd{f}=nc{'USTD_4097'}(:);                    %eastward velocity standard deviation, cm/s
        v{f}=nc{'v_1206'}(:);                          %northward velocity, cm/s
        vstd{f}=nc{'VSTD_4098'}(:);                    %northward velocity standard deviation, cm/s
        brange{f}=nc{'brange'}(:);                     %sensor range to boundary, m
        pressure{f}=(nc{'P_4022'}(:)-1000)/100;        %pressure, mbar
        depth{f}=nc{'depth'}(:);                       %bin depths below surface, mbs
        theFillValue=nc{'u_1205'}.FillValue_(:);       %netcdf fill value
        delta_t = gmean(diff(jd{f}));                  %time step, julian days
        site=dataFiles{f}(1:4);                         %site #
        nc=close(nc);
    end

    badInd = find( u{f}+v{f} >= theFillValue );
    u{f}(badInd) = nan;
    v{f}(badInd) = nan;

    badInd = find( ustd{f}+vstd{f} >= theFillValue );
    ustd{f}(badInd) = nan;
    vstd{f}(badInd) = nan;
    
    badInd = find(brange{f} >= theFillValue );
    brange{f}(badInd) = nan;

    jd_all=[jd{f}];
    for var=1:nVars
        eval([varList{var},'_all=[',varList{var},'{f}];'])
    end

    % interpolate data onto 'jd_int'
    for var=1:nVars
        eval([varList{var},'_all_int=smart_interp(jd_all, ',...
            varList{var},'_all, jd_int, maxGap);'])
    end


    % (optionally filter the velocities)
    switch filt
        case 'yes'
            disp('Low-pass filtering u,v current velocities');
            for var=1:nVars
                eval([varList{var},'_all_int=PL66TN(',...
                    varList{var},'_all_int, delta_t*24, 33);'])
            end
        otherwise
            disp('Not filtering current velocities');
    end

    % combine data from all sites in an array for plotting

    for var=1:nVars
        eval([varList{var},'_all_comb(:,f)=',...
            varList{var},'_all_int;'])
    end
    %%%%%%%%%%%%%%%%%%%%%% U,V TIME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % call timeplt.m
    xdata=jd_int;
    ydata=[vstd_all_comb(:,f),...
        ustd_all_comb(:,f),...
        v_all_comb(:,f),...
        u_all_comb(:,f),...
        brange_all_comb(:,f),...
        hsig_all_comb];
    panelOrd=[1:min(size(ydata))];
    figure
    han=timeplt2(xdata,ydata,panelOrd,ylims);
    clear global timeplt*
    set(gcf,'ResizeFcn',[])
    set(gcf,'Tag','')
    figCol=get(gcf,'color');
    dolandscape
    fixpaper
    pos=get(han(6),'position');

    % reformat xlims, ylabels, lines
    for h=1:length(han)
        axes(han(h))
        set(gca,'xlim',xlims,'xtick',xt)
        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)])
        if h~=length(han)
            linehan=findobj(gca,'Type','line','color','b');
            set(linehan,'color',col(1,:));
            linehan=findobj(gca,'Type','line','color',[0 .5 0]);
            set(linehan,'color',col(2,:));
        end
    end

    % re-do xlabels
    axes(han(1))
    set(gca,'XTickLabel',{''})
    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
    switch filt
        case 'yes'
            set(gcf,'PaperPositionMode','Auto')
            suptitle([{'Low-Passed ADV Bottom Water Flow Observations',...
                ['Site ',site]}]);
            print(fullfile(plotDirectory,['adv_site',site,'_lp']),'-dpdf');
        otherwise
            set(gcf,'PaperPositionMode','Auto')
            suptitle([{'ADV Bottom Water Flow Observations',...
                ['Site ',site]}]);
            print(fullfile(plotDirectory,['adv_site',site]),'-dpdf');
    end

end


