% 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, 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 06/20/14, 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:
%       smart_interp.m
%       PL66TN.m
%       julian.m
%       gregorian.m
%       suptitle.m
%       tight_subplot.m
%       gmean.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 ADV data filenames
dataFiles={'9205advs-cal.nc','9206advs-cal.nc','9235advs-cal.nc','9236advs-cal.nc'};
datanames={'Site 1','Site 1','Site 2','Site 2'};
waveFiles={'9241whp-cal-trm.nc','9261awWvs-p-ll-trm.nc',...
    '9281whp-cal-trm.nc','9303sgdataWvs-p-trm.nc','9321awWvs-p-ll-trm.nc'};
wavenameFiles={'Site 2','Site 3','Site 4','Site 5','site 6'};

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)'},...
    {'Eastward','Velocity','(cm/s)'},...
    {'Northward','Velocity','(cm/s)'},...
    {'Eastward','Velocity STD','(cm/s)'},...
    {'Northward','Velocity STD','(cm/s)'},...          %plot ylabels for stick plots
    };
startTime=julian(2012,1,27,00);
stopTime=julian(2012,04,19,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=2012:2012
    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 5;0 1;-30 30;-30 30;0 25; 0 25;];      %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);
% 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
for f=1:length(dataFiles)
       if ~isempty(dataFiles{f})
        nc=(fullfile(dataDirectory,dataFiles{f}));
        jd{f}=double(ncread(nc,'time')) + (double(ncread(nc,'time2'))/3600/1000/24);
        u{f}=squeeze(ncread(nc,'u_1205'));                          %eastward velocity, cm/s
        ustd{f}=squeeze(ncread(nc,'USTD_4097'));                    %eastward velocity standard deviation, cm/s
        v{f}=squeeze(ncread(nc,'v_1206'));                          %northward velocity, cm/s
        vstd{f}=squeeze(ncread(nc,'VSTD_4098'));                    %northward velocity standard deviation, cm/s
        brange{f}=squeeze(ncread(nc,'brange'));                     %sensor range to boundary, m
        try
        pressure{f}=squeeze((ncread(nc,'P_4023')-1000)/100);        %pressure, mbar
        catch
        pressure{f}=nan; 
        end
        depth{f}=squeeze(ncread(nc,'depth'));                       %bin depths below surface, mbs
        try
            theFillValue=ncreadatt(nc,'u_1205','_FillValue');
        catch
            theFillValue=1.0E35
        end
        delta_t = gmean(diff(jd{f}));                  %time step, julian days
        site=dataFiles{f}(1:4);                         %site #
    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
end
    %%%%%%%%%%%%%%%%%%%%%% U,V TIME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%     ydata=[vstd_all_comb(:,f),...
%         ustd_all_comb(:,f),...
%         v_all_comb(:,f),...
%         u_all_comb(:,f),...
%         brange_all_comb(:,f),...
%         hsig_int(:,1)];
    figure;set(gcf,'position',[100 25 850 1100]);
    ha = tight_subplot(6,1,[.04 .3],[.1 .1],[.1 .1]) ;
    
    axes(ha(1));
plot(jd_int,hsig_int,'LineWidth',1.2);ylim(ylims(1,:));
legend(wavenameFiles,'FontSize',7,'Location','NorthWest','Orientation','Horizontal');

axes(ha(2));
plot(jd_int,brange_all_comb,'LineWidth',1.2);ylim(ylims(2,:));
legend(datanames,'FontSize',7,'Location','NorthWest','Orientation','Horizontal');

axes(ha(3));
plot(jd_int,u_all_comb,'LineWidth',1.2);ylim(ylims(3,:));
legend(datanames,'FontSize',7,'Location','SouthWest','Orientation','Horizontal');

axes(ha(4));
plot(jd_int,v_all_comb,'LineWidth',1.2);ylim(ylims(4,:));
legend(datanames,'FontSize',7,'Location','SouthWest','Orientation','Horizontal');

axes(ha(5));
plot(jd_int,ustd_all_comb,'LineWidth',1.2);ylim(ylims(5,:));
legend(datanames,'FontSize',7,'Location','NorthWest','Orientation','Horizontal');

axes(ha(6));
plot(jd_int,vstd_all_comb,'LineWidth',1.2);ylim(ylims(6,:));
legend(datanames,'FontSize',7,'Location','NorthWest','Orientation','Horizontal');


   
    % reformat xlims, ylabels, lines
    for h=1:length(ha)
        axes(ha(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([2012,1,21,0,0,0]) pos(2) pos(3)])
        if h~=length(ha)
            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
    for h=1:length(ha)
    axes(ha(h))
    set(gca,'XTickLabel',{''})
    end
    axes(ha(6))
    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',...
                ['Sites 1 & 2 ']}]);
            print(fullfile(plotDirectory,['adv_sites_lp']),'-dpdf');
        otherwise
            set(gcf,'PaperPositionMode','Auto')
            suptitle([{'ADV Bottom Water Flow Observations',...
                ['Sites 1& 2 ']}]);
            print(fullfile(plotDirectory,['adv_sites']),'-dpdf');
    end




