% aqdp_wvs_ts.m  An mfile to plot the timeseries of wave parameters from
%                South Carolina Aquadopp 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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% C. Sullivan 11/17/04, Raw Aquadopp data is loaded into Matlab where a
%                       PUV analysis of pressure and velocity data occurs.
%                       From this analysis we extract wave information.
%             02/08/05, Datafiles now stored in 'dataDirectory'.  Rename
%                       pdf's waves_site#.pdf. 
%             03/03/05, Reformatting of labels and ticks
%             05/24/06, Update paths to files and file names.  Update to
%                       show Hs in the top panel of the plot. 
% Dependencies:
%       the netcdf toolbox
%       smart_interp.m
%       timeplt.m
%       dolandscape.m
%       suptitle.m

dataDirectory='E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\DATAFILES\';
plotDirectory='E:\CSULLIVAN\SOUTH_CAROLINA\DATA_DVD\PLOTS\WAVES';

% define data filenames and figure numbers to plot data in
dataFiles={'7311aqwp-cal.nc','7511aqwp-cal.nc'};

% define parameters for plots
col=get(0,'defaultaxescolororder');                 %plot color order
labels={ {'Frequency (Hz)'},...                    %plot ylabels
          {'Mean-Removed','Water Depth (m)'},...
          {'Wave Direction','(deg towards)'},...
          {'Peak Wave','Period (s)'},...
          {'Significant Wave','Height (m)'}};
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
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
ylims=[0 1; -2 2; 180 360; 0 20; 0 2];      %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]';

% initialize variables
jd=cell(length(dataFiles),1);
hs=cell(length(dataFiles),1);
tp=cell(length(dataFiles),1);
wvdir=cell(length(dataFiles),1);
depth=cell(length(dataFiles),1);
vspec=cell(length(dataFiles),1);

for f=1:length(dataFiles)
      
        % load data
        if ~isempty(dataFiles{f})
            nc=netcdf(dataFiles{f}, 'nowrite');
            jd{f}=nc{'time'}(:) + (nc{'time2'}(:)/3600/1000/24);
            hs{f}=nc{'wh_4061'}(:);                     %significant wave height, meters
            tp{f}=nc{'wp_peak'}(:);                     %peak period, seconds
            wvdir{f}=nc{'wvdir'}(:);                    %wave direction, degrees
            depth{f}=nc{'hght_18'}(:);                  %water depth, meters
            vspec{f}=nc{'vspec'}(:)./1000;              %non-directional velocity-derived spectra, M/sqrt(Hertz)
            freq=nc{'frequency'}(:);                         %frequency, Hertz
            site=nc.DESCRIPT(6);                        %site number
            theFillValue=nc{'u_1205'}.FillValue_(:);    %netcdf fill value
            delta_t = gmean(diff(jd{f}));               %time step, julian days
            nc=close(nc);
        end
        
        bads = find(vspec{f}>=1e30);
        vspec{f}(bads) = nan;
        
        % remove the mean from the depth
        depth{f}=depth{f}-gmean(depth{f});

        % convert direction from direction from to direction towards
        ind=find(wvdir{f} >= 180);
        wvdir{f}(ind) = wvdir{f}(ind) - 180;
        ind=find(wvdir{f} < 180);
        wvdir{f}(ind) = wvdir{f}(ind) + 180;
        
        if mod(f,2)==0
            
            % combine data from deployment 1 and
            % deployment 2, and replace netcdf
            % fill values with nan's.
            jd_all=[jd{f-1}; jd{f}];
            hs_all=[hs{f-1}; hs{f}];
            hs_all(hs_all>=1e30)=nan;
            tp_all=[tp{f-1}; tp{f}];
            tp_all(tp_all>=1e30)=nan;
            wvdir_all=[wvdir{f-1}; wvdir{f}];
            wvdir_all(wvdir_all>=1e30)=nan;
            depth_all=[depth{f-1}; depth{f}];
            depth_all(depth_all>=1e30)=nan;
            vspec_all=[vspec{f-1}; vspec{f}];
            
            % interpolate data onto 'jd_int'
            hs_int = smart_interp(jd_all, hs_all, jd_int, 3);
            tp_int = smart_interp(jd_all, tp_all, jd_int, 3);
            wvdir_int = smart_interp(jd_all, wvdir_all, jd_int, 3);
            depth_int = smart_interp(jd_all, depth_all, jd_int, 3);
            for j=1:length(freq)
                vspec_int(:,j) = smart_interp(jd_all, vspec_all(:,j),...
                                              jd_int, 3);
            end

            % call timeplt.m
            xdata=jd_int;
            ydata=[NaN(size(jd_int)) depth_int wvdir_int tp_int hs_int];
            panelOrd=[1:min(size(ydata))];
            figure             
            han=timeplt(xdata,ydata,panelOrd,ylims);
            clear global timeplt*
            set(gcf,'ResizeFcn',[])
            set(gcf,'Tag','')
            dolandscape
            fixpaper
            set(gcf,'PaperPositionMode','Auto')
            set(gcf,'Renderer','Painters')
            
            % pcolor velocity spectra in the empty
            % bottom panel
            axes(han(1))
            pcolor(jd_int(1:2:end),freq,vspec_int(1:2:end,:)')
            shading flat
            caxis([-0.5 1.5])
            axpos=get(han(1),'Position');
            cbar=colorbar;
            cbarpos=get(cbar,'Position');
            set(cbar,'Position',[0.91 axpos(2) cbarpos(4)/4 cbarpos(4)])
            set(get(cbar,'Ylabel'),'String','M/sqrt(Hz)')
            
            % reformat xlims, and ylabels
            for h=1:length(han)
                axes(han(h))
                set(gca,'xlim',xlims,'xtick',xt)
                ylabhan=get(gca,'Ylabel');
                set(ylabhan,'String',labels{h})
                pos=get(ylabhan,'position');
                set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)])
                obj=findobj('type','line','color','b');
                set(obj,'color',col(7,:))
                if h == 2
                    obj=findobj('type','line','color','k');
                    delete(obj)
                end
            end

            % re-do xlabels
            axes(han(1))
            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)-0.1;
                text(xt(tic),xtlpos,xtl,'horizontalalignment','center',...
                    'verticalalignment','top')
            end

            % the ouput pdf filename
           suptitle({'Wave Observations',['Site ',site]})
            print(fullfile(plotDirectory,['waves_site',site]),'-dpdf')
            
        end
            
end