function plot_ABS_stats(statCdfFile, vars, xax, pFlag)
% plot_ABS_stats.m  A function to plot ABS statistics data.
%
% usage:    plot_ABS_stats(statCdfFile, varFlag, pFlag);
%
%       where:  staCdfFile is the name of an ABS statistics NetCDF file, 
%                          surrounded by single quotes with the .cdf or .nc
%                          file extension
%               vars is a cell array of strings to specify which statistics
%                    to plot. Those statistics for all transducers will be
%                    shown. Choose from
%                           {}: empty cell array will plot all statistics
%                               in the NetCDF file.
%                           {'mean'}: mean amplitudes
%                           {'mean_sq'}: mean amplitudes squared                                   xducers
%                           {'std'}: standard deviation of amplitudes
%                           {'std_sq'}: standard deviation of amplitudes
%                                       squared
%                           {'pctl'}: amplitude percentiles
%                           {'pctl_sq'}: amplitude percentiles squared
%               xax is a string that specifies the xaxis units. 'burst'
%                   is for burst, 'time' is for time. 
%               pFlag is a flag to indicate whether plots are to be
%                       saved
%                           pFlag=0 --> no plots output
%                           pFlag=1 --> output .pdf's
%                           pFlag=2 --> output .png's
%                          
% Copyright 2005
% USGS Woods Hole Field Center
% Written by Charlene Sullivan
% csullivan@usgs.gov
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is self described.
% 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:
%   gregorian.m

% C. Sullivan   08/26/05,   version 1.2
% Changed name from plot_ABS_stat_cdf.m to plot_ABS_stats.m, because this
% function will work with any ABS statistics data file. Added ability
% to plot specific variables or all variables (default). Added ability to
% specify the x-axis as time or burst.
% C. Sullivan   07/15/05,   version 1.1
% Add option to output either .pdf's or .png's with the flag 'pFlag'
% C. Sullivan   07/05/05,   version 1.0
% This function pcolors all data in an ABS statistics netCDF file.  The
% x-axis is time and the y-axis is distance from the ABS. The user can
% optionally output these plots to pdf files by setting pdfFlag=1.


more off

version = '1.2';

tic

% Load the ABS statistics file
cdfs = netcdf(statCdfFile);

% Set up time and depth axis
jd = cdfs{'time'}(:) + cdfs{'time2'}(:)/(3600*1000*24);
dn = datenum(gregorian(jd));
nTime = size(dn,1);
burst = cdfs{'burst'}(:);
nBurst = size(burst,1);
pctiles = {'1st','16th','50th','84th','99th'}; 
depth = cdfs{'r'}(:);
nDepth = size(depth,1);
if size(depth,2)>1
    cols = size(depth,2);
    expr = [];
    for c = 1:cols
        eval(['depth_chk',num2str(c),' = depth(:,c);'])
        expr = [expr,' depth_chk',num2str(c),','];
        if c==cols; expr(end)=[]; end
    end
    eval(['tf = isequal(',expr,');'])    
    if tf == 1
        depth = depth(:,1);
    end
    clear S depth_chk* cols expr tf
end
nDepths = length(depth);
time = repmat(dn', nDepths, 1);
burst = repmat(burst', nDepths, 1);
pctiles = repmat(pctiles', nDepths, 1);
depth = repmat(depth, 1, nTime);

% Variables to plot
if isempty(vars)
    theVars = ncnames(var(cdfs));
else
    nVars = length(vars);
    theVars = cell(nVars*3,1);
    for v = 1:nVars
            theVars{v*3-2} = ['abs_trans1_',vars{v}];
            theVars{v*3-1} = ['abs_trans2_',vars{v}];
            theVars{v*3} = ['abs_trans3_',vars{v}];
    end
end

% X-axis units
switch xax
    case 'burst'
        xdata = burst;
        xlab = 'Burst number';
    case 'time'
        xdata = time;
        xlab = 'Time';
end

for i = 1:length(theVars)
    if ~strcmp(theVars{i},'time') && ...
       ~strcmp(theVars{i},'time2') && ...
       ~strcmp(theVars{i},'burst') && ...
       ~strcmp(theVars{i},'r') && ...
       ~strcmp(theVars{i},'pctile')
            
            theVar = theVars{i};
            eval(['data = cdfs{''',theVar,'''}(:);'])
            f = find(data >= 1e35);
            if ~isempty(f)
                data(f) = nan;
            end
            eval(['ttl = cdfs{''',theVar,'''}.long_name(:);'])
            eval(['clrttl = cdfs{''',theVar,'''}.units(:);'])
            [row, col, lvl] = size(data);
            if isequal(lvl,1)
                j = str2num(theVar(10)); %xducer #
                switch j
                    case 1
                        figure('pos',[100 100 1000 800])
                        orient landscape
                        subplot(311)
                    case 2
                        subplot(312)
                    case 3
                        subplot(313)
                end
                pcolor(xdata, flipud(depth), log(flipud(data')))
                shading flat
                if strcmp(xax,'time')
                    datetick('x',2,'keepticks','keeplimits')
                end
                ylabel('Distance from ABS (m)')
                set(gca,'YDir','reverse')
                set(gca,'Ytick',[0:.25:depth(end)])
                set(gca,'YtickLabel',num2str([0:.25:depth(end)]',2))
                set(gca,'Layer','top')
                axpos = get(gca,'position');
                if j ==1;
                    ca=caxis;
                else
                    caxis(ca);
                end
                L = axpos(1)-.07; B = axpos(2); W = axpos(3)+.09; H = .2;
                set(gca,'position',[L B W H])
                clrbr = colorbar;
                cpos = get(clrbr,'position');
                set(clrbr,'position',[L+W+.01 B 0.02 H])
                set(get(clrbr,'Ylabel'),'String',clrttl)
                title(ttl)
                if isequal(j,3)
                    xlabel(xlab)
                    if ~isequal(pFlag,0)
                        set(gcf,'renderer','painters')
                        switch pFlag
                            case 1
                                print('-dpdf',theVars{i})
                            case 2
                                print('-dpng',theVars{i})
                        end
                    end
                end
            else
                figure('pos',[100 100 1000 800])
                orient landscape
                for k = 1:lvl
                    subplot(lvl,1,k)
                    pcolor(xdata, flipud(depth), log(flipud(data(:,:,k)')))
                    shading flat
                    if strcmp(xax,'time')
                        datetick('x',2,'keepticks','keeplimits')
                    end
                    ylabel({'Distance from','ABS (m)'})
                    set(gca,'YDir','reverse')
                    set(gca,'Layer','top')
                    title([pctiles{k},' Percentile'])
                    axpos = get(gca,'position');
                    if k ==1;
                        ca=caxis;
                    else
                        caxis(ca);
                    end
                    L = axpos(1)-.05; B = axpos(2); W = axpos(3)+.07; H = .1;
                    set(gca,'position',[L B W H])
                    clrbr = colorbar;
                    cpos = get(clrbr,'position');
                    set(clrbr,'position',[L+W+.01 B 0.02 H])
                    set(get(clrbr,'Ylabel'),'String',clrttl)
                    if isequal(k,1)
                        if strcmp(xax,'time')
                            text(((time(1)+time(end))/2),-0.05,ttl,'horizontalalignment',...
                                'center','fontsize',14)
                        else
                            text(((burst(1)+burst(end))/2),-0.05,ttl,'horizontalalignment',...
                                'center','fontsize',14)
                        end
                    elseif isequal(k,5)
                        xlabel(xlab)
                        if ~isequal(pFlag,0)
                            set(gcf,'renderer','painters')
                            switch pFlag
                                case 1
                                    print('-dpdf',theVar)
                                case 2
                                    print('-dpng',theVar)
                            end
                        end
                    end
                end

            end
    end
end
    
close(cdfs)

toc