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

% C. Sullivan 03/14/06, 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:
%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
% Dependencies:  
%       the netcdf toolbox
%       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\PCADP';

% define ADV data filenames 
dataFiles={'728pcvp-cal5.nc','748pcvp-cal6.nc'};
nFiles=length(dataFiles);

% define parameters for plots
filt='yes';                                         %filtering switch, 'yes' or 'no'
col=get(0,'defaultaxescolororder');                 %plot color order
ylab={ {'Site 2','Significant Wave','Height (m)'},...
       {'Sensor Range','to Boundary','(m)'},...
       {'Northward','Velocity','(cm/s)'},...
       {'Eastward','Velocity','(cm/s)'},...
       {'Northward','Velocity STD','(cm/s)'},...
       {'Eastward','Velocity STD','(cm/s)'} };
site='7A';
startTime=julian(2003,10,25,00);                            %plot start stime
stopTime=julian(2004,04,25,00);                             %plot stop time  
timediff=stopTime-startTime;                                %plot xlimits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %beyond start time and stop time
ylims=[0 2; 0.85 1.1; -12 12; -12 12; 0 35; 0 35;];      %u,v plot ylimits
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];   % xticks on the first of every month
        end
    end
end

% '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 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);

% combine wave data from deployments 1 and 2
jd_waves_all=[jd_waves{1}; nan; jd_waves{2}];
hsig_all=[hsig{1}; nan; hsig{2}];

% filter the waves (optional)
if strcmp(filt,'yes')
    hsig_all = PL66TN(hsig_all,waveDelta_t*24,33);
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=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'}(:);
        for n=1:nVars
            ncobj=nc{ncVars{n}};
            theFillValue=ncobj.FillValue_(:);
            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
        nc=close(nc);

end
        
% combine data from deployments 1 and 2
jd_all=[jd{f-1}; nan; jd{f}];
depth_all=[depth{f-1}; depth{f}];
u_all=[u{f-1}; nan(1,9); u{f}];
v_all=[v{f-1}; nan(1,9); v{f}];
ustd_all=[ustd{f-1}; nan(1,9); ustd{f}];
vstd_all=[vstd{f-1}; nan(1,9); vstd{f}];
brange_all=[brange{f-1}; nan(1,3); brange{f}];
press_all=[press{f-1}; nan; press{f}];
            
%%%%%%%%%%%%%%%%%%%%%%%%%% Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
dolandscape
fixpaper

Han(1)=subplot(6,1,1,'align');
plot(jd_waves_all, hsig_all)
hold on
Han(2)=subplot(6,1,2,'align');
plot(jd_all,brange_all(:,1))
hold on
plot(jd_all,brange_all(:,2),'r')
plot(jd_all,brange_all(:,3),'k')
LEG2=legend('Beam 1','Beam 2','Beam 3');
set(LEG2,'orientation','horizontal')
set(LEG2,'location','NorthEast')
set(LEG2,'box','off')

% to contour the velocity bins ...
x=repmat(jd_all,1,length(depth_all));
z=repmat(depth_all',length(jd_all),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},'{1};'])
    eval(['data2 = ',vars{n},'{2};'])
    y(1:length(jd{1}),1:length(depth{1})) = data1;
    y(length(jd{1})+2:length(jd_all),length(depth{1})+1:length(depth_all)) = data2;
    y2=y';
    y2_sort=y2(sortInd,:);
    y3=y2_sort';
    Han(n+2)=subplot(6,1,n+2,'align');
    for ii=1:length(depth_all)
        [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(Han)
    axes(Han(h))
    axpos=get(gca,'position');
    set(gca,'xlim',xlims,'xtick',xt)
    set(gca,'xticklabel',[])
    set(gca,'ylim',ylims(h,:))
    ylabhan=get(gca,'Ylabel');
    set(ylabhan,'String',ylab{h})
    pos=get(ylabhan,'position');
    set(get(gca,'Ylabel'),'position',[julian([2003,10,6,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(Han(1))
switch filt
    case 'yes'
        text(julian(2004,01,29),3,{'Low-Passed PCADP Bottom Water Flow Observations',...
            'Site 7A'},'horizontalalignment','center','fontsize',14);
        print(fullfile(plotDirectory,['pcadp_site',site,'_lp']),'-dpdf');
    otherwise
        text(julian(2004,01,29),3,{'PCADP Bottom Water Flow Observations',...
            'Site 7A'},'horizontalalignment','center','fontsize',14);
        print(fullfile(plotDirectory,['pcadp_site',site]),'-dpdf');
end