% pcadp_ts.m  An mfile to plot the (optionally filtered)
%             timeseries of current velocity data from South
%             Carolina PCADP data.

%                       Raw PCADP 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.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Dependencies:
%       smart_interp.m
%       PL66TN.m
%       suptitle.m
%       replace.m
%       gmean.m
%       gregorian.m
%       julian.m
%       draw.m
%       tight_subplot.m

clear
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 ADV data filenames
dataFiles={'9203pcs-cal.nc','9233pcs-cal.nc'};
datanames={'Site 1','Site 2'};
waveFiles={'9241whp-cal-trm.nc','9261awWvs-p-ll-trm.nc',...
    '9281whp-cal-trm.nc','9303sgdataWvs-p-trm.nc','9321awWvs-p-ll-trm.nc'};
wavenameFiles={'Site 2','Site 3','Site 4','Site 5','site 6'};
nFiles=length(dataFiles);

% define parameters for plots
filt='no';                                         %filtering switch, 'yes' or 'no'
col=get(0,'defaultaxescolororder');                 %plot color order
ylab={ {'Significant Wave Height All Stations(m)'},...
    {'Sensor Range to Boundary (m)'},...
    {'Northward Velocity (cm/s)'},...
    {'Eastward Velocity (cm/s)'},...
    {'Northward Velocity STD (cm/s)'},...
    {'Eastward Velocity STD (cm/s)'} };
startTime=julian(2012,1,27,00);
stopTime=julian(2012,04,19,00);
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=2012:2012
    for mm=1:12
        tic=julian(yy,mm,01,00);
        if tic>=xlims(1) & tic<=xlims(2)
            xt=[xt; tic];
        end
    end
end

ylims=[0 5; 0.9 1.1; -25 25; -25 25; 0 30; 0 30;];      %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)';
nTime = length(jd_int);
maxGap = 3;

%% Load site wave data
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

% initialize variables for PCADP data
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);
press=cell(nFiles,1);

% define PCADP variables
ncVars={'u_1205','v_1206','USTD_4097','VSTD_4098','brange','P_4023'};
vars={'u','v','ustd','vstd','brange','press'};
nVars=length(vars);

% load the PCADP data
for f=1:nFiles
    nc=(fullfile(dataDirectory,dataFiles{f}));
    jd{f}=double(ncread(nc,'time')) + (double(ncread(nc,'time2'))/3600/1000/24);
    delta_t = gmean(diff(jd{f}));
    depth{f}=ncread(nc,'depth');
    for n=1:nVars
        ncobj=squeeze(ncread(nc,ncVars{n}));
        try
            theFillValue=ncreadatt(nc,'u_1205','_FillValue');
        catch
            theFillValue=1.0E35
        end
        data=ncobj;
        badInd=find(data>=theFillValue);
        data(badInd)=nan;
        if strcmp(vars{n},'press')
            % de-mean the pressure signal
            data=(data-1000)/100;
            data=data-gmean(data);
        end
        % filter the data (optional)
        switch filt
            case 'yes'
                [nRows,nCols]=size(data);
                for col=1:nCols
                    data(:,col)=PL66TN(data(:,col), delta_t*24, 33);
                end
            otherwise
                disp('Not filtering current velocities');
        end
        eval([vars{n},'{f}=data;'])
    end
end


%%    figure;set(gcf,'position',[100 25 850 1100]);
%%%%%%%%%%%%%%%%%%%%%%%%%% Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for st=1:2
        figure;set(gcf,'position',[100 25 850 1100]);
        ha = tight_subplot(6,1,[.04 .4],[.05 .05],[.1 .15]) ;
        
        axes(ha(1));
        plot(jd_int,hsig_int,'LineWidth',1.2);ylim(ylims(1,:));
        
        axes(ha(2));
        plot(jd{st},brange{st},'LineWidth',1.2);ylim(ylims(2,:));
        legend('beam 1','beam 2','beam 3','FontSize',7,'Location','SouthWest','Orientation','Horizontal');
        
        axes(ha(3));
        
        % to contour the velocity bins ...
        x=repmat(jd{st},1,length(depth{st}));
        z=repmat(depth{st}',length(jd{st}),1);
        x2=x';
        z2=z';
        [z2_sort,sortInd]=sortrows(z2,1);
        x2_sort=x2(sortInd,:);
        x3=x2_sort';
        z3=z2_sort';
        y=nan(size(x));
        for n=1:4
            eval(['data1 = ',vars{n},'{st};'])
            y = data1;
            y2=y';
            y2_sort=y2(:,sortInd);
            y3=y2_sort;y3(y3>40)=nan;
            axes(ha(n+2));
            for ii=1:length(depth{st})
                [Ct,Ht]=contour(x3, y3, z3, [z3(1,ii) z3(1,ii)]);
                hold on
                C(ii)=Ct(1);
                H(ii)=Ht(1);
            end
        end
        
        % reformat xlims, ylabels, lines
        for h=1:length(ha)
            axes(ha(h))
            axpos=get(gca,'position');
            set(gca,'xlim',xlims,'xtick',xt)
            set(gca,'xticklabel',[])
            set(gca,'ylim',ylims(h,:))
            ylabhan=get(gca,'Ylabel');
            title(ha(h),ylab{h})
            pos=get(ylabhan,'position');
            set(get(gca,'Ylabel'),'position',[julian([2008,12,30,0,0,0]) pos(2) pos(3)])
        end
        
        % make a legend and reformat xlabels
        LEG=legend(H,num2str(round(C'*100)/100,5));
        figCol=get(gcf,'color');
        set(LEG,'color',figCol)
        set(LEG,'edgecolor',figCol)
        legpos=get(LEG,'position');
        axpos=get(gca,'position');
        set(LEG,'position',[axpos(1)+axpos(3)+.003 axpos(2) legpos(3) legpos(4)])
        set(LEG,'XAxisLocation','top')
        LEGLBL=get(LEG,'Xlabel');
        set(LEGLBL,'String','Depth (mbs)')
        set(LEGLBL,'Color',[0 0 0],'fontweight','bold')
        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
        set(gcf,'paperpositionmode','auto')
        axes(ha(1));
        switch filt
            case 'yes'
                suptitle(['PCADP Low-PAss Filtered Bottom Water Flow Observations Site ',num2str(st)])%,'horizontalalignment','center','fontsize',14);
                legend(wavenameFiles,'FontSize',10,'Location','NorthWest','Orientation','Horizontal');
                print(fullfile(plotDirectory,['pcadp_site',num2str(st),'_lp']),'-dpdf');
            otherwise
                suptitle(['PCADP Bottom Water Flow Observations Site ',num2str(st)])%,'horizontalalignment','center','fontsize',14);
                legend(wavenameFiles,'FontSize',10,'Location','NorthWest','Orientation','Horizontal');
                print(fullfile(plotDirectory,['pcadp_site',num2str(st)]),'-dpdf');
        end
    end
