% 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 = 'E:\Hatteras10\candidate_final_files'; 
dataDirectory ='\\IGSAGIEGGS-FILE\Stg-proc\data\HatterasNearshore10\BBV_reviewed';
% plotDirectory = 'E:\Hatteras10\Hatteras2010DVD\PLOTS';
plotDirectory = 'D:\BNA\data\Carolinas\webpage\dataDVD\Hatteras2010DVD\PLOTS';

% define ADCP data filenames 
% Define adcp data filenames.
plot_location=2;%1=north,2=south,3=offshore
if plot_location==1
    %current file names
    datafile={'8771aqd-a.nc','8781aqd-a.nc','8791aqd-a.nc',...
        '8801wh-a.nc','8641wh-a.nc','8651aqd-cal.nc'};
    
    %wave file names
    waveFiles={'8771aqwp-cal.cdf','8781aqdwp-cal.nc','8791aqwp-cal.cdf',...
        '8801adwp-cal.cdf','8641adwp-cal.cdf','8651aqwp-cal.cdf'};
    
    %names of stations
    names={'N1','N2','N3','N4','N5','N6'};
    wavenames=names;
    %location
    location='North Sites';
elseif plot_location==2
    %current file names
    datafile={'8751aqd-a.nc','8711wh-a.nc','8661aqd-cal.nc',...
        '8671Aaqd-cal.nc','8671Baqd-cal.nc'};
    
    %wave file names
    waveFiles={'8751aqwp-cal.cdf','8711whp-cal.cdf','8661aqwp-cal.cdf',...
        '8671Aaqwp-cal.cdf','8671Baqwp-cal.cdf'};
    
    %names of stations
    names={'N8','N9','N11','N12A','N12B'};
    wavenames=names;
    %location
    location='South Sites';
else %plot location ==3
    %current file names
    datafile={'8721wh-a.nc','8731wh-a.nc','8741wh-a.nc'};
    
    %wave file names
    waveFiles={'8721adwp-cal.nc','8731adwp-cal.cdf','8741Aadwp-cal.cdf','8741Badwp-cal.cdf'};
    
    %names of stations
    names={'O1','O2','O3'};
    wavenames={'O1','O2','O3'};
    %location
    location='Offshore Sites';
end

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
if plot_location==3
    startTime=julian(2010,2,5,00);
    stopTime=julian(2010,3,24,00);
else
    startTime=julian(2010,2,5,00);
    stopTime=julian(2010,02,23,00);
end
timediff=stopTime-startTime;                               
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; 
xt=[];
for yy=2010:2010
    for mm=2:3
        tic=julian(yy,mm,01,00); 
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic]; 
        end
    end
end

% set y-limits, labels
ylims=[-180 180; 0 150; -180 180; 0 150; 0 4];
yl_topbot={ {'Direction Surface','(degrees)'},...
    {'Magnitude Surface','(cm/s)'},...
    {'Direction Bottom','(degrees)'},...
    {'Magnitude Bottom','(cm/s)'},...
    {'Significant Wave','Height (m)'} };
    
switch sticks
    case 'yes'
        yl_sticks={ {'Velocity Bottom','(cm/s)'},...
             {'Velocity 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
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'));
    theFillValue=ncreadatt(ncwaves,'wh_4061','_FillValue');
    hsig(hsig >= theFillValue)=nan;
    hsig(hsig <= 0)=nan;
    hsig_int(:,f) = interp1([jd_waves], [hsig],...
        jd_int);
%         hsig_int(:,f) = smart_interp([jd_waves], [hsig],...
%         jd_int, 3);

end

for f=1:length(names)
    
    
    nc=(fullfile(dataDirectory,datafile{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-4,:);'])
        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, maxGap);'])
%         eval([varList{var},'_surface_int(:,f)=smart_interp(jd, ',...
%             varList{var},'_surface, jd_int, maxGap);'])
        eval([varList{var},'_bottom_int(:,f)=interp1(jd, ',...
            varList{var},'_bottom, jd_int);'])
        eval([varList{var},'_surface_int(:,f)=interp1(jd, ',...
            varList{var},'_surface, jd_int);'])

    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_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



    %%%%%%%%%%%%%%%%%%%%%%MAG and DIR Plots%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % interpolate data onto 'jd_int'
    
    angle_r=zeros(size(u_bottom_int));
    velb=u_bottom_int  + v_bottom_int  * sqrt(-1);
    vels=u_surface_int  + v_surface_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();
    dirs=(-atan2(qvs,qus)+(pi()/2))*180/pi();
    dirb(dirb<0)=dirb(dirb<0)+360;
    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;
end
time_dn=datenum(gregorian(jd_int));

    

 %%   %%%%%%%%%%%%%%%%%%%%%% Mag,Dir PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % call timeplt.m to create individual U, V plots at the BOTTOM
    % and TOP
    xdata=jd_int;
    
clf
set(gcf,'position',[100 25 850 1100])
if plot_location==3
corder=get(gca,'ColorOrder');
corder2=zeros(11,3);
corder2(1,:)=corder(1,:);
corder2(2,:)=corder(2,:);
corder2(3,:)=corder(3,:);
corder2(4,:)=corder(3,:);
end
clf

ha = tight_subplot(5,1,[.04 .3],[.1 .1],[.1 .1]) ;
if plot_location==3
    set(ha(:),'ColorOrder',corder2);
end

axes(ha(5));
if plot_location==3
set(gca,'NextPlot','replacechildren');
end
plot(time_dn,hsig_int,'LineWidth',1.2);
datetick('x',2,'keeplimits','keepticks');
if plot_location==3
set(gca,'YTickLabel',[0 1 2 3 4]);
end
% legend(wavenames,'FontSize',7)
ylabel(yl_topbot{5});

axes(ha(1));
plot(time_dn,dirs,'LineWidth',1.2);
legend(names,'FontSize',7);set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
ylabel(yl_topbot{1});

axes(ha(2));
plot(time_dn,mags,'LineWidth',1.2);
set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
ylabel(yl_topbot{2});

axes(ha(3));
plot(time_dn,dirb,'LineWidth',1.2);
set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
ylabel(yl_topbot{3});

axes(ha(4));
plot(time_dn,magb,'LineWidth',1.2);
set(gca,'xlim',xlim,'xtick',get(ha(5),'xtick'),'XTickLabel','');
ylabel(yl_topbot{4});

 
set(gcf,'PaperPositionMode','Auto','PaperOrientation','portrait')

    % give it a title and export to .pdf
    switch filt
        case 'yes'
            tit1e=get(ha(5),'Title');
            set(tit1e,'String',['Low-Passed Water Flow Observations ',location])
            set(gcf,'PaperPositionMode','Auto')
            switch coord
                case 'en'
                    print(fullfile(plotDirectory,[location,'_en_tb_lp']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,[location,'_xl_tb_lp']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
        otherwise
            tit1e=get(ha(1),'Title');
            set(tit1e,'String',['Water Flow Observations ',location])
            set(gcf,'PaperPositionMode','Auto')

            switch coord
                case 'en'
                    print(fullfile(plotDirectory,[location,'_magdir_sb']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,[location,'_magdirrot_sb']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
    end
