% aqdp_curr_ts.m  An mfile to plot the (optionally filtered) timeseries
%                 of current velocities from South Carolina Aquadopp 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 12/01/04, It is assumed that the Aquadopp data files are
%                       stored in a directory defined by the variable
%                       'dataDirectory'.  It assumes the user has run
%                       spit_wad.m by G. Voulgaris in the directory to split
%                       the .wad data file into individual burst .mat files
%                       with pressure and velocity data.  To plot low-pass
%                       filtered current velocities set the switch 'filt' to
%                       'yes', otherwise (filt='no') unfiltered current
%                       velocities are plotted.
%             01/27/05, Aquadopp data has now been processed by running
%                       READ_aquapro_beam.m which loads the raw velocites
%                       and converts them from beam to ENU coordinates.
%                       Spikes and bad data are nan'd by aqdp_clean.m.
%                       Processed data is written to netcdf files by
%                       aqdpCurr2nc.m.  Timeseries plots of current velocites
%                       are created using velocites at the bottom, 6 mbs, 
%                       4 mbs, and the surface.
%             02/24/05, Reformatting of labels and ticks, add additional
%                       plots
%             05/24/06, Update paths to files and file names.  Update to
%                       show Hs in the top panel of the plot. 
%             05/26/06, Add switch to show velocities in either east and
%                       north coordinates or cross-shore and longshore
%                       coordinates
%
% Dependencies:  
%       the netcdf toolbox
%       rot.m
%       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\ADCP';

% define Aquadopp data filenames.  
dataFiles={'7311aqc-a.nc','7511aqc-a.nc'};
nFiles=length(dataFiles);

% define switches
filt='yes';   %filtering switch, 'yes' to filter data, or 'no'
coord='xl';  %coordinates switch,'xl' for cross-shore, longshore; 'en' for
             %east,north
sticks='no';%sticks switch, 'yes' to create magnitude,direction plots and
             %u,v plots, 'no' to create just u,v plots                                 

% set x-limits, ticks; x-limits extend 5% beyond startTime and stopTime
startTime=julian(2003,10,25,00);                   
stopTime=julian(2004,04,25,00);                    
timediff=stopTime-startTime;                               
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; 
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]; 
        end
    end
end

% set y-limits, labels
ylims=[-30 30; -30 30; -30 30; -30 30; 0 2]; 
switch coord
    case 'en'
        yl_bott6={ {'Northward','Velocity','at the Bottom','(cm/s)'},...
             {'Eastward','Velocity','at the Bottom','(cm/s)'},...
             {'Northward','Velocity','at 6 mbs','(cm/s)'},...
             {'Eastward','Velocity','at 6 mbs','(cm/s)'},...
             {'Site 2','Significant Wave','Height (m)'} };
        yl_4surf={ {'Northward','Velocity','at the 4 mbs','(cm/s)'},... 
             {'Eastward','Velocity','at the 4 mbs','(cm/s)'},...
             {'Northward','Velocity','at the Surface','(cm/s)'},...
             {'Eastward','Velocity','at the Surface','(cm/s)'},...
             {'Site 2','Significant Wave','Height (m)'} };
    case 'xl'
        yl_bott6={ {'Cross-shore','Velocity','at the Bottom','(cm/s)'},... 
             {'Longshore','Velocity','at the Bottom','(cm/s)'},...
             {'Cross-shore','Velocity','at 6 mbs','(cm/s)'},...
             {'Longshore','Velocity','at 6 mbs','(cm/s)'},...
             {'Site 2','Significant Wave','Height (m)'} };
        yl_4surf={ {'Cross-shore','Velocity','at the 4 mbs','(cm/s)'},... 
             {'Longshore','Velocity','at the 4 mbs','(cm/s)'},...
             {'Cross-shore','Velocity','at the Surface','(cm/s)'},...
             {'Longshore','Velocity','at the Surface','(cm/s)'},...
             {'Site 2','Significant Wave','Height (m)'} };
    otherwise
        error('Not a valid coordinates switch!')
end
switch sticks
    case 'yes'
        yl_sticks={ {'Velocity','at the Bottom','(cm/s)'},...
             {'Velocity','at 6 mbs','(cm/s)'},...
             {'Velocity','at 4 mbs','(cm/s)'},...
             {'Velocity','at the Surface','(cm/s)'},...
             {'Site 2','Significant Wave','Height (m)'} };
end

% Interpolate all data onto the array 'jd_int' for timeplt.m. Data gaps
% 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');

% 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_bottom=cell(nFiles,1);
u_6mbs=cell(nFiles,1);
u_4mbs=cell(nFiles,1);
u_surface=cell(nFiles,1);
v_bottom=cell(nFiles,1);
v_6mbs=cell(nFiles,1);
v_4mbs=cell(nFiles,1);
v_surface=cell(nFiles,1);
varList={'u','v'};
ncVarList={'u_1205','v_1206'};
nVars=length(ncVarList);

for f=1:length(dataFiles)

        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'}(:);  
        site=nc.DESCRIPT(6);
        for v=1:length(ncVarList)
            ncobj=nc{ncVarList{v}};
            theFillValue=ncobj.FillValue_(:);
            data=ncobj(:);
            bads=find(data>=theFillValue);
            data(bads)=nan;
            questBins=ncobj.bins_questioned(:);
            %velocity in bottom bin
            eval([varList{v},'_bottom{f}=data(:,1);']) 
            %velocity in surface bin where surface is upward-most good bin
            if isempty(questBins)
                eval([varList{v},'_surface{f}=data(:,end);'])
            else
                eval([varList{v},'_surface{f}=data(:,questBins(1)-1);'])
                %uppermost good bin at site 8 is at 5.6mbs for the second
                %deployment so don't show this on the plots
                if f==2
                   u_surface{f}=nan(size(data,1),1);
                   v_surface{f}=nan(size(data,1),1);
                end  
            end
            %velocity in bin 4 mbs
            [binInd,dist2bin]=value2Index(depth{f},4,0.5);
            eval([varList{v},'_4mbs{f}=data(:,binInd);'])
            %velocity in bin 6 mbs
            [binInd,dist2bin]=value2Index(depth{f},6,0.5);
            eval([varList{v},'_6mbs{f}=data(:,binInd);'])
            if v==length(ncVarList)
                % (optionally convert velocities from east, north coordinates
                % to longshore, cross-shore coordinates using rot.m. Assume
                % a shoreline orientation of 45 deg.
                switch coord
                    case 'xl'
                        disp(['Converting from north,east coords '...
                                           'to cross,long coords'])
                        [u_bottom{f},v_bottom{f}]=rot(u_bottom{f},v_bottom{f},-45);
                        [u_6mbs{f},v_6mbs{f}]=rot(u_6mbs{f},v_6mbs{f},-45);
                        [u_4mbs{f},v_4mbs{f}]=rot(u_4mbs{f},v_4mbs{f},-45);
                        [u_surface{f},v_surface{f}]=rot(u_surface{f},v_surface{f},-45);
                    case 'en'
                        disp('Not converting coordinate systems')
                    otherwise
                        disp('Not a valid coordinate system!')
                end
            end
        end
        nc=close(nc);
        
        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},'_bottom_all=[',...
                      varList{var},'_bottom{f-1}; ',varList{var},'_bottom{f}];'])
                eval([varList{var},'_6mbs_all=[',...
                      varList{var},'_6mbs{f-1}; ',varList{var},'_6mbs{f}];'])
                eval([varList{var},'_4mbs_all=[',...
                      varList{var},'_4mbs{f-1}; ',varList{var},'_4mbs{f}];'])
                eval([varList{var},'_surface_all=[',...
                      varList{var},'_surface{f-1}; ',varList{var},'_surface{f}];'])
            end
                    
            % interpolate data onto 'jd_int' 
            for var=1:nVars
                eval([varList{var},'_bottom_all_int=smart_interp(jd_all, ',...
                      varList{var},'_bottom_all, jd_int, maxGap);'])
                eval([varList{var},'_6mbs_all_int=smart_interp(jd_all, ',...
                      varList{var},'_6mbs_all, jd_int, maxGap);'])
                eval([varList{var},'_4mbs_all_int=smart_interp(jd_all, ',...
                      varList{var},'_4mbs_all, jd_int, maxGap);'])
                eval([varList{var},'_surface_all_int=smart_interp(jd_all, ',...
                      varList{var},'_surface_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},'_bottom_all_int=PL66TN(',...
                              varList{var},'_bottom_all_int, delta_t*24, 33);'])
                        eval([varList{var},'_6mbs_all_int=PL66TN(',...
                              varList{var},'_6mbs_all_int, delta_t*24, 33);'])
                        eval([varList{var},'_4mbs_all_int=PL66TN(',...
                              varList{var},'_4mbs_all_int, delta_t*24, 33);'])
                        eval([varList{var},'_surface_all_int=PL66TN(',...
                              varList{var},'_surface_all_int, delta_t*24, 33);'])
                    end
                otherwise
                    disp('Not filtering current velocities')
            end
            
            %%%%%%%%%%%%%%%%%%%%%% STICK PLOT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % call timeplt.m to create individual STICK plots
            switch sticks
                case 'yes'
                    % Make vectors
                    clear i
                    uv_bottom_all_int = u_bottom_all_int + i.*v_bottom_all_int;
                    uv_6mbs_all_int = u_6mbs_all_int + i.*v_6mbs_all_int;
                    uv_4mbs_all_int = u_4mbs_all_int + i.*v_4mbs_all_int;
                    uv_surface_all_int = u_surface_all_int + i.*v_surface_all_int;
                    % (clearing global timeplt*, the ResizeFcn,
                    % and the Tag are necessary to release
                    % timeplt.m's control over your figure!)
                    xdata=jd_int;
                    ydata=[uv_bottom_all_int uv_6mbs_all_int,...
                        uv_4mbs_all_int uv_surface_all_int,...
                        hsig_all_comb];
                    xdata=xdata(1:8:end); %show a stick every 2 hours
                    ydata=ydata(1:8:end,:);
                    panelOrd=[1:min(size(ydata))];
                    figure
                    han=timeplt(xdata,ydata,panelOrd,ylims);
                    clear global timeplt*
                    set(gcf,'ResizeFcn',[])
                    set(gcf,'Tag','')
                    dolandscape
                    fixpaper

                    % reformat xlims,ylabels, and line color
                    for h=1:length(han)
                        axes(han(h))
                        set(gca,'xlim',xlims,'xtick',xt)
                        ylabhan=get(gca,'Ylabel');
                        set(ylabhan,'String',yl_sticks{h})
                        pos=get(ylabhan,'position');
                        set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) %align ylabels nicely
                        if h~=length(han)
                            obj=findobj(gca,'type','line','color','b');
                            set(obj,'color',col(7,:))
                        end
                        if h==4
                            xtxt=julian(2004,3,15);
                            ytxt=15;
                            text(xtxt,ytxt,'No Data','horizontalalignment','center')
                        end
                    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)}; %2 line label
                        xtlpos=ylims(1,1)-2;
                        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 Water Flow Observations',['Site ',site]});
                            print(fullfile(plotDirectory,['s',site,'_spdir_s46b_lp']),'-dpdf');
                        otherwise
                            set(gcf,'PaperPositionMode','Auto')
                            suptitle({'Water Flow Observations',['Site ',site]});
                            print(fullfile(plotDirectory,['s',site,'_spdir_s46b']),'-dpdf');
                    end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            %%%%%%%%%%%%%%%%%%%%%% U,V PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % call timeplt.m to create individual U, V plots at the BOTTOM
            % and 6 MBS
            xdata=jd_int;
            ydata=[v_bottom_all_int u_bottom_all_int,...
                  v_6mbs_all_int u_6mbs_all_int,...
                  hsig_all_comb];
            panelOrd=[1:min(size(ydata))]; 
            figure
            han=timeplt(xdata,ydata,panelOrd,ylims);
            clear global timeplt*
            set(gcf,'ResizeFcn',[])
            set(gcf,'Tag','')
            dolandscape
            fixpaper

            % reformat xlims, ylabels, and line color
            for h=1:length(han)
                axes(han(h))
                set(gca,'xlim',xlims,'xtick',xt)
                ylabhan=get(gca,'Ylabel');
                set(ylabhan,'String',yl_bott6{h})
                pos=get(ylabhan,'position');
                set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) 
                if h~=length(han)
                    obj=findobj(gca,'type','line','color','b');
                    set(obj,'color',col(7,:))
                end
            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 Water Flow Observations',['Site ',site]});
                    switch coord
                        case 'en'
                            print(fullfile(plotDirectory,['s',site,'_en_6b_lp']),'-dpdf');
                        case 'xl'
                            print(fullfile(plotDirectory,['s',site,'_xl_6b_lp']),'-dpdf');
                        otherwise
                            disp('Not a valid coordinate system!')
                    end
                otherwise
                    set(gcf,'PaperPositionMode','Auto')
                    suptitle({'Water Flow Observations',['Site ',site]});
                    switch coord
                        case 'en'
                            print(fullfile(plotDirectory,['s',site,'_en_6b']),'-dpdf');
                        case 'xl'
                            print(fullfile(plotDirectory,['s',site,'_xl_6b']),'-dpdf');
                        otherwise
                            disp('Not a valid coordinate system!')
                    end
            end
            
            % call timeplt.m to create individual U, V plots at the 4 MBS
            % and the SURFACE
            xdata=jd_int;
            ydata=[v_4mbs_all_int u_4mbs_all_int,...
                   v_surface_all_int u_surface_all_int,...
                   hsig_all_comb];
            panelOrd=[1:min(size(ydata))]; 
            figure
            han=timeplt(xdata,ydata,panelOrd,ylims);
            clear global timeplt*
            set(gcf,'ResizeFcn',[])
            set(gcf,'Tag','')
            dolandscape
            fixpaper

            % reformat xlims, ylabels, and line color
            for h=1:length(han)
                axes(han(h))
                set(gca,'xlim',xlims,'xtick',xt)
                ylabhan=get(gca,'Ylabel');
                set(ylabhan,'String',yl_4surf{h})
                pos=get(ylabhan,'position');
                set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)]) 
                if h~=length(han)
                    obj=findobj(gca,'type','line','color','b');
                    set(obj,'color',col(7,:))
                end
                if h==3 || h==4
                    xtxt=julian(2004,3,15);
                    ytxt=15;
                    text(xtxt,ytxt,'No Data','horizontalalignment','center')
                end
            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 Water Flow Observations',['Site ',site]});
                    switch coord
                        case 'en'
                            print(fullfile(plotDirectory,['s',site,'_en_s4_lp']),'-dpdf');
                        case 'xl'
                            print(fullfile(plotDirectory,['s',site,'_xl_s4_lp']),'-dpdf');
                        otherwise
                            disp('Not a valid coordinate system!')
                    end
                otherwise
                    set(gcf,'PaperPositionMode','Auto')
                    suptitle({'Water Flow Observations',['Site ',site]});
                    switch coord
                        case 'en'
                            print(fullfile(plotDirectory,['s',site,'_en_s4']),'-dpdf');
                        case 'xl'
                            print(fullfile(plotDirectory,['s',site,'_xl_s4']),'-dpdf');
                        otherwise
                            disp('Not a valid coordinate system!')
                    end
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            
        end
        
end