% abs_ts_test.m  An mfile to plot the timeseries of statistics calculated from
%           South Carolina ABS data.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Charlene Sullivan
% USGS Woods Hole Science Center
% csullivan@usgs.gov
%
% Dependencies:  
%       the netcdf toolbox
%       gregorian.m
%       suptitle.m

% C. Sullivan   06/13/06,   version 1.1
% Update paths. Show Hs from site 2 in top panel.
% C. Sullivan   09/05/05,   version 1.0


more off

version = '1.0';

tic

% Define data and output directories
dataDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES';
plotDirectory = 'E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\ABS';

% Define ABS data filenames 
dataFiles={'7264abss-cal.nc','7464abss-cal.nc',...
           '7322abss-cal.nc','7522abss-cal.nc',...
           '','7484abss-cal.nc'};
nFiles=length(dataFiles);

% Define x-axis (time) limits, ticks,
% and labels
startTime=julian(2003,10,25,00);                            %plot start time
stopTime=julian(2004,04,25,00);                             %plot stop time
timediff=stopTime-startTime;                                %plot x-limits extend 5 percent
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; %beyond start time and stop time
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

% Significant wave height from site 2 is
% shown in the top panel of all plots
wvs1 = netcdf('7221adw-a.nc');
wvs2 = netcdf('7421adw-a.nc');
wave_time = [wvs1{'time'}(:) + wvs1{'time2'}(:)/3600/1000/24;...
             wvs2{'time'}(:) + wvs2{'time2'}(:)/3600/1000/24];
hsig = [wvs1{'wh_4061'}(:); wvs2{'wh_4061'}(:)];
theFillValue = wvs1{'wh_4061'}.FillValue_(:);
badInd = find(hsig >= theFillValue);
hsig(badInd) = nan;

% Initialize ABS index variables
time = cell(nFiles,1);
time2 = cell(nFiles,1);
burst = cell(nFiles,1);
r = cell(nFiles,1);
pctile = cell(nFiles,1);
pctiles = {'1st','16th','50th','84th','99th'}; 

for f = 5:6%1:length(dataFiles)
    
    % Define the site (b/c its not in the .nc files!)
    if f==1 || f==2
        site = '6';
    elseif f==3 || f==4
        site = '7B';
    elseif f==5 || f==6
        site = '7A';
    end
    
    if ~isempty(dataFiles{f})
        abs = netcdf(fullfile(dataDirectory,dataFiles{f}));
        absVars = var(abs);
        for v = 1:length(absVars)
            absVar = absVars{v};
            absVarName = name(absVar);
            if f==1
                eval([absVarName,' = cell(nFiles,1);']) %Initialize ABS variables
            end
            eval([absVarName,'{f} = absVar(:);']) %Assign ABS data to variables
        end
    end
       
    if mod(f,2)==0
        % Units on 'r' for USC instrument are
        % incorrect
        if f==2
            r{f} = r{f}*10;
        end
        % If the depth array, 'r', has more than 1
        % column and the columns are all the same,
        % keep just one column
        sz1 = size(r{f-1});
        if sz1(2)>1
            if diff(r{f-1},1,2) == 0
                r{f-1} = r{f-1}(:,1);
            end
        end
        sz2 = size(r{f});
        if sz2(2)>1
            if diff(r{f},1,2) == 0
                r{f} = r{f}(:,1);
            end
        end
        
        % Get depth (y) axis
        fall_r = r{f-1};
        spring_r = r{f};
        max_r = max([fall_r;spring_r]);
        fall_nR = size(fall_r,1);
        spring_nR = size(spring_r,1);
        
        % Get time (x) axis
        fall_abs_time = time{f-1}(:) + time2{f-1}(:)/3600/1000/24;
        spring_abs_time = time{f}(:) + time2{f}(:)/3600/1000/24;
        fall_nTime = size(fall_abs_time,1);
        spring_nTime = size(spring_abs_time,1);
        
        % Make depth and time axis for pcolor
        fall_r = repmat(fall_r, 1, fall_nTime);
        spring_r = repmat(spring_r, 1, spring_nTime);
        fall_abs_time = repmat(fall_abs_time', fall_nR, 1);
        spring_abs_time = repmat(spring_abs_time', spring_nR, 1);
        
        
        for v = 1:length(absVars)
            absVar = absVars{v};
            absVarName = name(absVar);
            if  ~strcmp(absVarName,'time') && ...
                ~strcmp(absVarName,'time2') && ...
                ~strcmp(absVarName,'burst') && ...
                ~strcmp(absVarName,'r') && ...
                ~strcmp(absVarName,'pctile')
                    ttl = absVar.long_name(:);
                    clrttl = absVar.units(:);
                    theFillValue = absVar.FillValue_(:);
                    xducer = str2num(absVarName(10));
                    eval(['fall_abs_data = [',absVarName,'{f-1}];'])
                    eval(['spring_abs_data = [',absVarName,'{f}];'])
                    badInd = find(fall_abs_data >= theFillValue);
                    fall_abs_data(badInd) = nan;
                    badInd = find(spring_abs_data >= theFillValue);
                    spring_abs_data(badInd) = nan;
                    [row, col, lvl] = size(spring_abs_data);
                    if isequal(lvl,1) %2D data
                        switch xducer
                            case 1
                                figure('pos',[100 100 1000 800])
                                orient landscape
                                supttl = '';
                                % wave height
                                P(1) = subplot(411);
                                plot(wave_time,hsig)
                                H = .18;
                                ylabel({'Significant wave','height (m)'})
                                title('Significant Wave Height from Site 2')
                                % abs data
                                P(2) = subplot(412);
                            case 2
                                P(3) = subplot(413);
                            case 3
                                P(4) = subplot(414);
                        end
                        if ~isempty(fall_abs_data)
                            pcolor(fall_abs_time(:,1:2:end), flipud(fall_r(:,1:2:end)), log(flipud(fall_abs_data(1:2:end,:)')))
                        else
                            text(datenum(2003,11,15),0.75,'No Data')
                        end
                        hold on
                        pcolor(spring_abs_time(:,1:2:end), flipud(spring_r(:,1:2:end)), log(flipud(spring_abs_data(1:2:end,:)')))
                        shading flat
                        title(ttl)
                        fname = absVarName(12:end);
                    else %3D data
                       figure('pos',[100 100 1000 800])
                       orient landscape
                       supttl = ttl;
                       % wave height
                       P(1) = subplot(611);
                       plot(wave_time,hsig)
                       H = .1;
                       ylabel({'Significant wave','height (m)'})
                       title('Significant Wave Height from Site 2')
                       % abs data
                       for k = 1:lvl
                            P(k+1) = subplot(6,1,k+1);
                            if ~isempty(fall_abs_data)
                                pcolor(fall_abs_time(:,1:2:end), flipud(fall_r(:,1:2:end)), log(flipud(fall_abs_data(1:2:end,:,k)')))
                            else
                                text(datenum(2003,11,15),0.75,'No Data')
                            end
                            hold on
                            pcolor(spring_abs_time(:,1:2:end), flipud(spring_r(:,1:2:end)), log(flipud(spring_abs_data(1:2:end,:,k)')))
                            shading flat
                            title([pctiles{k},' Percentile'])
                       end
                       fname = ['trans',num2str(xducer),'_',absVarName(12:end)];
                    end
                    if xducer == 3 && ~exist('k') || exist('k') && k == lvl
                    % format axis
                    for p = 1:length(P);
                        ax = P(p);
                        % format axis position
                        axpos = get(ax,'position');
                        L = axpos(1)-.05; B = axpos(2); W = axpos(3)+.07;
                        set(ax,'position',[L B W H])
                        % format time (x) axis
                        set(ax,'xlim',xlims,'xtick',xt)
                        set(ax,'xticklabel',[])
                        if p>1
                            % format depth (y) axis
                            ylabel(ax,{'Distance from','ABS (m)'})
                            set(ax,'YDir','reverse')
                            set(ax,'Ytick',[0:.25:1.25])
                            set(ax,'YtickLabel',num2str([0:.25:1.25]'))
                            set(ax,'Ylim',[0 1.35])
                            set(ax,'box','on')
                            set(ax,'Layer','top')
                            if p >= 2;
                                if strfind(absVarName,'sq')
                                    caxis(ax,[8 22]);
                                else
                                    caxis(ax,[2 10])
                                end
                                % format colorbar
                                clrbr = colorbar('peer',ax);
                                cpos = get(clrbr,'position');
                                set(clrbr,'position',[L+W+.01 B 0.02 H])
                                set(get(clrbr,'Ylabel'),'String',clrttl)
                                if p == length(P)
                                    for t=1:length(xt)
                                        xtg=gregorian(xt(t));
                                        xtd=datenum(xtg);
                                        xtl={[datestr(xtd,3),' 1'];datestr(xtd,10)}; %2 line label
                                        xtlpos=max_r+0.15;
                                        text(xt(t),xtlpos,xtl,'horizontalalignment','center',...
                                            'verticalalignment','top')
                                    end
                                    suptitle({['Site ',site,' ABS Observations'],supttl})
                                    set(gcf,'renderer','painters')
                                    % print the .pdf --> you will need to
                                    % manually rotate these
                                    print('-dpdf',fullfile(plotDirectory,['site',site,'_abs_',fname]))
                                    % print the .png.  manual rotation of
                                    % these is not required
                                    orient portrait
                                    set(gcf,'papersize',[11 8.5])
                                    set(gcf,'paperposition',[.25 .25 10.5 8])
                                    print('-dpng',fullfile(plotDirectory,['site',site,'_abs_',fname]))
                                    close
                                    clear k
                                end
                            end
                        end
                    end
                    end
               end
            end
        end
    end

    
toc

ncclose