% all_wvs_ts_mtg29Mar05.m  An mfile to plot the timeseries of wave
%                          parameters from South Carolina ADCP, Argonaut,
%                          and Aquadoppdata.    

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 10/21/04, Raw ADCP data is first processed through WavesMon
%                       software.  The mcode adcpWvs2nc.m transforms the
%                       data output by WavesMon into a netCDF file.  It is
%                       assumed that the netCDF files are stored in a 
%                       directory named according to mooring number and 
%                       instrument type.  For ADCP data on my machine this 
%                       'E:\CSULLIVAN\SOUTH_CAROLINA\DATAFILES\###\###1wh\'
%                       where ### is mooring number.
%             02/07/05, Datafiles now stored in 'dataDirectory'.  
%                       Interpolate data onto a standard time axis.  Name
%                       pdf's waves_site#.pdf. 
%             03/03/05, Reformatting of labels and ticks
%             03/15/05, Create 'all observations' version for OFR.  This 
%                       version does not show Vspec.
%             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 
dataFiles={'7201adwp-cal.nc','7401adwp-cal.nc',...
           '7221adwp-cal.nc','7421adwp-cal.nc',...
           '','7441adwp-cal.nc',...
           '7331adwp-cal.nc','7531adwp-cal.nc',...
           '7263arwp-cal.nc','7463arwp-cal.nc',...
           '7321adwp-cal.nc','7521adwp-cal.nc',...
           '7311aqwp-cal.nc','7511aqwp-cal.nc'};
nFiles=length(dataFiles);

% define parameters for plots
col=get(0,'defaultaxescolororder');                 %plot color order
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=[-2 2; 180 360; 0 20; 0 2];      %stick 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'}(:);                    %peak wave direction, degrees from
            depth{f}=nc{'hght_18'}(:);                  %water depth, meters
            freq=nc{'freq'}(:);                         %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);
            
            % 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 site 1, or 7 don't show depth (pressure) b/c its
            % funky. this data should be filled in the .nc files
            if ismember(f,[1 11])
                depth{f} = nan(size(jd{f}));
            end
        end
        
        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}];
            vspec_all(vspec_all>=1e30)=nan;
            
            % interpolate data onto 'jd_int'
            hs_int(:,f/2) = smart_interp(jd_all, hs_all, jd_int, 3);
            tp_int(:,f/2) = smart_interp(jd_all, tp_all, jd_int, 3);
            wvdir_int(:,f/2) = smart_interp(jd_all, wvdir_all, jd_int, 3);
            depth_int(:,f/2) = smart_interp(jd_all, depth_all, jd_int, 3);

            % remove the mean from the water depth
            depth_int(:,f/2)=depth_int(:,f/2)-gmean(depth_int(:,f/2));
        end
        
end


% call timeplt.m
xdata=jd_int;
ydata=[depth_int wvdir_int tp_int hs_int];
panelOrd=[repmat(1,1,nFiles/2),...
          repmat(2,1,nFiles/2),...
          repmat(3,1,nFiles/2),...
          repmat(4,1,nFiles/2)];
figure
han=timeplt(xdata,ydata,panelOrd,ylims);
clear global timeplt*
set(gcf,'ResizeFcn',[])
set(gcf,'Tag','')
LEG=legend(han(4),'Site 1','Site 2','Site 3','Site 4',...
                  'Site 6','Site 7B','Site 8');
set(LEG,'Location','NorthOutside')   
set(LEG,'Orientation','Horizontal')
figCol=get(gcf,'color');
set(LEG,'color',figCol)
set(LEG,'edgecolor',figCol)
dolandscape
fixpaper
legpos=get(LEG,'position');
pos=get(han(4),'position');
set(LEG,'position',[pos(1) .95 pos(3) legpos(4)])

% 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',ylabels{h})
    pos=get(ylabhan,'position');
    set(get(gca,'Ylabel'),'position',[julian(2003,10,8) pos(2) pos(3)])
    if h==1
        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

% output to pdf
set(gcf,'PaperPositionMode','Auto')
suptitle([{'Wave Observations','All Sites'}])
print(fullfile(plotDirectory,['waves_all']),'-dpdf')

