% adcp_curr_ts.m  An mfile to plot the (optionally filtered) timeseries
%                 of current velocities from South Carolina ADCP 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 12/20/04, Raw ADCP 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
%       rot.m
%       smart_interp.m
%       PL66TN.m
%       timeplt.m
%       dolandscape.m
%       suptitle.m

clear
more off

% define data and output directories
dataDirectory = '\\IGSAGIEGGS-FILE\Stg-proc\data\hatteras09\BBV_reviewed\';
plotDirectory = 'C:\barmstrong\data\Hatteras_2009\';

% define ADCP data filenames 
dataFiles={'8542wh-a.nc','8551wh-a.nc',...
            '8572wh-a.nc','8581wh-a.nc'};
waveFiles={'8542p-cal.nc','8551p-cal.nc',...
            '8572p-cal.nc',''};
        nameFiles={'North FloBee','North Minipod',...
            'South FloBee','Tip Micropod'};

nFiles=length(dataFiles);

% define switches
filt='no';   %filtering switch, 'yes' to filter data, or 'no'
coord='en';  %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(2009,1,8,00);                   
stopTime=julian(2009,05,15,00);                    
timediff=stopTime-startTime;                               
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; 
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

% set y-limits, labels
ylims=[0 360; 0 150; 0 360; 0 150; 0 4];
yl_topbot={ {'Direction','at the Surface','(degrees)'},...
    {'Magnitude','at the Surface','(cm/s)'},...
    {'Direction','at the Bottom','(degrees)'},...
    {'Magnitude','at the Bottom','(cm/s)'},...
    {'Significant','Wave','Height (m)'} };
    
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)'},...
             {'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');


% 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)
    if f~=4 %tip site
        nc=netcdf([dataDirectory,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);
    end
    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(1:11)
    for v=1:length(ncVarList)
        ncobj=nc{ncVarList{v}};
        theFillValue=1.0E35;
        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);'])
        end
        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_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);
    close(nc);


    jd_all=[jd{f}];
    for var=1:nVars
        eval([varList{var},'_bottom_all=[',varList{var},'_bottom{f}];'])
        eval([varList{var},'_surface_all=[',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},'_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},'_surface_all_int=PL66TN(',...
                    varList{var},'_surface_all_int, delta_t*24, 33);'])
            end
        otherwise
            disp('Not filtering current velocities');
    end


    %%%%%%%%%%%%%%%%%%%%%%MAG and DIR Plots%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % interpolate data onto 'jd_int'
    
    angle_r=zeros(size(u_bottom_all_int));
    velb=u_bottom_all_int  + v_bottom_all_int  * sqrt(-1);
    vels=u_surface_all_int  + v_surface_all_int  * sqrt(-1);
    velb=velb .* exp ( sqrt(-1) * angle_r);
    vels=vels .* exp ( sqrt(-1) * angle_r);
    qub=real(velb);
    qus=real(vels);
    qvb=imag(velb);
    qvs=imag(vels);
    dirb=(atan2(qvb,qub)-(pi()/2))*180/pi();
    dirb(dirb<0)=dirb(dirb<0)+360;
    dirs=(atan2(qvs,qus)-(pi()/2))*180/pi();
    dirs(dirs<0)=dirs(dirs<0)+360;
    magb=sqrt(qub.^2+qvb.^2);clear qub qvb velb;
    mags=sqrt(qus.^2+qvs.^2);clear qus qvs vels;

    %%%%%%%%%%%%%%%%%%%%%% Mag,Dir PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % call timeplt.m to create individual U, V plots at the BOTTOM
    % and TOP
    xdata=jd_int;
     if f==4
        hsig_all_comb(:)=NaN;
    end
    ydata=[dirb magb,...
        dirs mags,...
        hsig_all_comb];
    panelOrd=[1:min(size(ydata))];
    figure
    han=timeplt2(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,'LineStyleOrder','.')
        ylabhan=get(gca,'Ylabel');
        set(ylabhan,'String',yl_topbot{h})
        pos=get(ylabhan,'position');
        set(get(gca,'Ylabel'),'position',[julian(2008,12,25) pos(2) pos(3)])
    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'
            tit1e=get(han(5),'Title');
            set(tit1e,'String',['Low-Passed Water Flow Observations ',nameFiles{f}])
            set(gcf,'PaperPositionMode','Auto')
            switch coord
                case 'en'
                    print(fullfile(plotDirectory,['s',dataFiles{f}(1:4),'_en_tb_lp']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,['s',dataFiles{f}(1:4),'_xl_tb_lp']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
        otherwise
            tit1e=get(han(5),'Title');
            set(tit1e,'String',['Water Flow Observations ',nameFiles{f}])
            set(gcf,'PaperPositionMode','Auto')

            switch coord
                case 'en'
                    print(fullfile(plotDirectory,['s',dataFiles{f}(1:4),'_magdir_sb']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,['s',dataFiles{f}(1:4),'_magdirrot_sb']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
    end



end