% 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:
%
% Sullivan, C.M., Warner, J.C., Martini, M.A., Voulgaris, G., 
% Work, P.A., Haas, K.A., and Hanes, D.H. (2006) 
% South Carolina Coastal Erosion Study Data Report for Observations
% October 2003 - April 2004., USGS Open-File Report 2005-1429.
%
% 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.
%             03/18/05, Create 'all observations' version for SCCES project
%                       meeting on 29Mar05.
%             04/22/05, Save mfile that creates 29Mar05 meeting version as
%                       adv_ts.m, and use for OFR plots.
%             05/23/05, Update paths to files and file names.  Update to
%                       show Hs from site 2 in the top panel of the plot. 
%                       Format title and ylabels like PCADP plots.
%
% Dependencies:  
%       the netcdf toolbox
%       smart_interp.m
%       PL66TN.m
%       timeplt.m
%       dolandscape.m
%       suptitle.m

clear
more off

% define data and output directories
dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\ADV';

% define ADV data filenames 
dataFiles={'adv7281vp-cal.nc','adv7481vp-cal.nc',...
           'adv7282vp-cal2.nc','adv7482vp-cal2.nc',...
           'adv7261vp-cal2.nc','adv7461vp-cal2.nc',...
           'adv7262vp-cal2.nc','adv7462vp-cal2.nc'};
nFiles=length(dataFiles);

% define parameters for plots
filt='yes';                                         %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)'},...
       {'Site 2','Significant Wave','Height (m)'} };
startTime=julian(2003,10,25,00);                            %plot start stime
stopTime=julian(2004,04,25,00);                             %plot stop time  
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=2003:2004
    for mm=1:12
        tic=julian(yy,mm,01,00); 
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic];   % xticks on the first of every month
        end
    end
end
ylims=[0 20; 0 20; -12 7; -12 7; .30 .55; 0 2];      %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;

% load site 2 wave data for the top panel
% of the plot
nc=netcdf(fullfile(dataDirectory,'7221adwp-cal.nc'), '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;
nc=close(nc);

nc=netcdf(fullfile(dataDirectory,'7421adwp-cal.nc'), 'nowrite');
jd_waves{2} = nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
hsig{2}=nc{'wh_4061'}(:);
theFillValue=nc{'wh_4061'}.FillValue_(:);
badInd=find(hsig{2} >= theFillValue);
hsig{2}(badInd)=nan;
nc=close(nc);

hsig_all_comb = smart_interp([jd_waves{1};jd_waves{2}], [hsig{1};hsig{2}],...
                         jd_int, 3);
                     
switch filt
    case 'yes'
        hsig_all_comb=pl66tn(hsig_all_comb, 1, 33);
end

% 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)
 
        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_4023'}(:)-1000)/100;        %pressure, mbar
        if f<3 || f>4
            mean_press{f}=gmean(pressure{f});
            pressure{f}=(pressure{f}-mean_press{f});
        end
        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=nc.DESCRIPT(6);                           %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;
        
        if mod(f,2)==0
            
            % combine data from deployment 1 and
            % deployment 2
            jd_all=[jd{f-1}; jd{f}];
            for var=1:nVars
                eval([varList{var},'_all=[',varList{var},'{f-1}; ',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
            ind=f/2;
            for var=1:nVars
                eval([varList{var},'_all_comb(:,ind)=',...
                      varList{var},'_all_int;'])
            end

        end
        
end

%%%%%%%%%%%%%%%%%%%%%% U,V TIME %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% call timeplt.m 
xdata=jd_int;
ydata=[vstd_all_comb,...
       ustd_all_comb,...
       v_all_comb,...
       u_all_comb,...
       brange_all_comb,...
       hsig_all_comb];
panelOrd=[repmat(1,1,nFiles/2),...
          repmat(2,1,nFiles/2),...
          repmat(3,1,nFiles/2),...
          repmat(4,1,nFiles/2),...
          repmat(5,1,nFiles/2),...
          repmat(6,1,1)];
figure
han=timeplt(xdata,ydata,panelOrd,ylims);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
LEG=legend(han(5),'Site 7A ADV1','Site 7A ADV2','Site 6 ADV1','Site 6 ADV2');
set(LEG,'Location','NorthOutside')   
set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
dolandscape
fixpaper
legpos=get(LEG,'position');
pos=get(han(6),'position');
set(LEG,'position',[pos(1) pos(2)+pos(4) pos(3) legpos(4)])

% reformat xlims, ylabels
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([2003,10,6,0,0,0]) pos(2) pos(3)])
end

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