% adcp_curr_ts.m  An mfile to plot the (optionally filtered) timeseries
%                 of current velocities from North Carolina aquadopp and ADCP 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, 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 6/20/14, 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:  
%       rot.m
%       smart_interp.m
%       PL66TN.m
%       gmean.m
%       draw.m
%       julian.m
%       gregorian.m
%       replace.m
%       tight_subplot

clear

addpath('D:\BNA\Matlab\m_cmg\branches\barmstrong\BNA\t_tide_v1.3beta','-end');
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 ADCP data filenames 
dataFiles={'9211wh-trm.nc','9241wh-trm.nc',...
            '9261aw-cal-ll-trm.nc','9281wh-trm.nc','9321aw-cal-ll-trm.nc'};
waveFiles={'9241whp-cal-trm.nc','9261awWvs-p-ll-trm.nc',...
            '9281whp-cal-trm.nc','9303sgdataWvs-p-trm.nc','9321awWvs-p-ll-trm.nc'};
names={'Site 1','Site 2','Site 3','Site 4','Site 6'};
wavenameFiles={'Site 2','Site 3','Site 4','Site 5','site 6'};


nFiles=length(names);

% 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(2012,1,27,00);                   
stopTime=julian(2012,04,19,00);         
timediff=stopTime-startTime;
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; 
xt=[];
for yy=2012:2012
    for mm=1:5
%         for dd=1:2:28
            tic=julian(yy,mm,1,0);
            if tic>=xlims(1) & tic<=xlims(2)
                xt=[xt; tic];
            end
%         end
    end
end

% set y-limits, labels
ylims=[ 0 4;-100 100; -100 100; -50 50; -50 50]; 

switch coord
    case 'en'
         yl_topbot={ {'Significant Wave','Height (m)'},{'Eastward Velocity','Surface (cm/s)'},...
             {'Northward Velocity','Surface (cm/s)'},...
             {'Eastward Velocity','Bottom (cm/s)'},...
             {'Northward Velocity','Bottom (cm/s)'}};
    case 'xl'
        yl_topbot={ {'Significant Wave','Height (m)'},{'Longshore Velocity','Surface (cm/s)'},... 
             {'Cross-shore','Velocity','at the Surface','(cm/s)'},...
             {'Longshore','Velocity','at the Bottom','(cm/s)'},...
             {'Cross-shore','Velocity','at the Bottom','(cm/s)'}};
    otherwise
        error('Not a valid coordinates switch!')
end
switch sticks
    case 'yes'
        yl_sticks={{'Significant Wave','Height (m)'}, {'Velocity Bottom','(cm/s)'},...
             {'Velocity Surface','(cm/s)'}};
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
u_bottom_int=zeros(length(jd_int),nFiles);
u_surface_int=zeros(length(jd_int),nFiles);
v_bottom_int=zeros(length(jd_int),nFiles);
v_surface_int=zeros(length(jd_int),nFiles);
hsig_int=zeros(length(jd_int),nFiles);
varList={'u','v'};
ncVarList={'u_1205','v_1206'};
nVars=length(ncVarList);
%%
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(names)
    nc=(fullfile(dataDirectory,dataFiles{f}));
    time=double(ncread(nc,'time'));
    time2=double(ncread(nc,'time2'));
    jd= time + (time2/3600/1000/24);
    delta_t=gmean(diff(jd));
    depth=ncread(nc,'depth');
    site=names{f};
    for v=1:length(ncVarList)
        ncobj=squeeze(ncread(nc,ncVarList{v}));
        theFillValue=ncreadatt(nc,ncVarList{v},'_FillValue');
        data=ncobj;
        bads=find(data>=theFillValue);
        data(bads)=nan;
        %velocity in bottom bin
        eval([varList{v},'_bottom=data(1,:);'])
        %velocity in surface bin where surface is upward-most good bin
        eval([varList{v},'_surface=data(end-10,:);'])
        eval(['plot(',varList{v},'_surface);']);
        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,v_bottom]=rot(u_bottom,v_bottom,-45);
                    [u_surface,v_surface]=rot(u_surface,v_surface,-45);
                case 'en'
                    disp('Not converting coordinate systems')
                otherwise
                    disp('Not a valid coordinate system!')
            end
        end
    end


    % interpolate data onto 'jd_int'
    for var=1:nVars
        eval([varList{var},'_bottom_int(:,f)=smart_interp(jd, ',...
            varList{var},'_bottom, jd_int,3);'])
        eval([varList{var},'_surface_int(:,f)=smart_interp(jd, ',...
            varList{var},'_surface, jd_int,3);'])

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

    %% %%%%%%%%%%%%%%%%%%%% U,V PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % create U, V plots at the BOTTOM and TOP
    xdata=jd_int;
    
    clf
    set(gcf,'position',[100 25 850 1100])
    
    ha = tight_subplot(5,1,[.04 .3],[.1 .1],[.1 .1]) ;
    axes(ha(5));plot(time_dn,v_bottom_int,'LineWidth',1.2);
    set(gca,'xtick',datenum(gregorian(xt)));
    datetick('x',2,'keeplimits','keepticks');
    ylabel(yl_topbot{5});
    
    axes(ha(1));
    plot(time_dn,hsig_int,'LineWidth',1.2);
    legend(names,'FontSize',7);
    set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
    set(gca,'Ylim',[ylims(1,:)]);
    ylabel(yl_topbot{1});
    
    axes(ha(2));plot(time_dn,u_surface_int,'LineWidth',1.2);
    set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
    ylabel(yl_topbot{2});
    
    axes(ha(3));plot(time_dn,v_surface_int,'LineWidth',1.2);
    set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
    ylabel(yl_topbot{3});
    
    axes(ha(4));plot(time_dn,u_bottom_int,'LineWidth',1.2);
    set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
    ylabel(yl_topbot{4});
    
    % reformat xlims,ylims, ylabels, and line color
    for h=1:length(ha)
        axes(ha(h));
        set(gca,'xlim',datenum(gregorian(xlims)))%,'xtick',datenum(gregorian(xt)))
        set(gca,'ylim',ylims(h,:))
    end
    
    % give it a title and export to .pdf
    switch filt
        case 'yes'
            tit1e=get(ha(1),'Title');
            set(tit1e,'String',['Low-Passed Water Flow Observations Fire Island 2012'])
            set(gcf,'PaperPositionMode','Auto','PaperOrientation','portrait')
            switch coord
                case 'en'
                    print(fullfile(plotDirectory,['FI12_en_sb_lp']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,['FI12_xl_sb_lp']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
        otherwise
            tit1e=get(ha(1),'Title');
            set(tit1e,'String',['Water Flow Observations Fire Island 2012'])
            set(gcf,'PaperPositionMode','Auto','PaperOrientation','portrait')
            
            switch coord
                case 'en'
                    print(fullfile(plotDirectory,['FI12adcp_en_sb']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,['FI12adcp_xl_sb']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
    end
    
    
    
