% adcp_magdir_ts2.m  An mfile to plot the (optionally filtered) timeseries
%                 of current velocities from North Carolina aquadopp and ADCP data.    

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use of this program is described in:
%
% Armstrong, Brandy N., Warner, John C., List, Jeff H., Martini, 
% Marinna A., Montgomery, Ellyn T., Traykovski, Peter, Voulgaris, George 
% 2014, Coastal Change Processes Project data report for oceanographic 
% observations near Fire Island, New York, February through May  
% 2014 : U.S. Geological Survey Open-File Report 2014-xxxx, available 
% online at http://pubs.usgs.gov/of/2014/xxxx/.
%
% Program written in Matlab v8.0.0.755
% Program ran on PC with Windows 7 64bit
%
%
% "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."
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% B. Armstrong 6/20/14, Raw ADCP 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.

%
% Dependencies:  
%       rot.m
%       smart_interp.m
%       PL66TN.m
%       gmean.m
%       draw.m
%       julian.m
%       gregorian.m
%       replace.m
%       tight_subplot

clear

%addpath('c:Users\emontgomery\mtl\m_cmg\branches\barmstrong\BNA\t_tide_v1.3beta','-end');
more off

% define data and output directories
dataDirectory = 'c:\home\data\validation\FI14\candidate_final_files';
plotDirectory = 'C:\Users\emontgomery\doc\OFR_FI14\FI_14_WEBOFR\plotsem';

% define ADCP data filenames - precursors of final published names used
%dataFiles={'9861wh_prx2_trm.nc','9891whall.nc','9892aw-cal.nc','9921wh3.nc','9922V.nc',...
 dataFiles={'9861wh_prx2_trm.nc','9891whall.nc','9921wh3.nc',...
           '9981wh.nc','10001wh_prx2_trm.nc','10021wh2.nc','10041wh_trm.nc'};
%waveFiles={'9861whp-calBS-trm.nc','9891p-cal.nc','9892awWvsall-trm.nc','9921p-cal404-trm.nc','9922Vp-cal.nc',...
waveFiles={'9861whp-calBS-trm.nc','9891p-cal.nc','9921p-cal404-trm.nc',...
            '9981p-cal.nc','10021whp-cal404-trm.nc'};
%           '9981whp-cal.nc','10001wh-p.nc','10021whp-cal.nc','10041wh-p.nc'};
names={'Site 1','Site 2','Site 3','Site 6','Site 7','Site 8','Site 9'};
wavenameFiles={'Site 1','Site 2','Site 3','Site 6','Site 8'};
%wavenameFiles={'Site 1','Site 2','Site 2awac','Site 3','Site 3V','Site 6','Site 7','Site 8','Site 9'};


nFiles=length(names);

% define switches
filt='no';   %filtering switch, 'yes' to filter data, or 'no'
coord='en';  %coordinates switch,'xl' for cross-shore, longshore; 'en' for
             %east,north
sticks='no';%sticks switch, 'yes' to create magnitude,direction plots and
             %u,v plots, 'no' to create just u,v plots

%% set x-limits, ticks; x-limits extend 5% beyond startTime and stopTime
startTime=julian(2014,2,8,00);                   
stopTime=julian(2014,05,5,00);         
timediff=stopTime-startTime;
xlims=[(startTime-timediff*0.05) (stopTime+timediff*0.05)]; 
xt=[];
for yy=2014:2014
    for mm=1:5
%         for dd=1:2:28
            tic=julian(yy,mm,1,0);
            if tic>=xlims(1) & tic<=xlims(2)
                xt=[xt; tic];
            end
%         end
    end
end

% set y-limits, labels
ylims=[0 4;0 360; 0 100; 0 360; 0 60];
yl_topbot={ {'Significant Wave','Height (m)'},{'Direction Surface','(degrees)'},...
    {'Magnitude Surface','(cm/s)'},...
    {'Direction Bottom','(degrees)'},...
    {'Magnitude Bottom','(cm/s)'},...
     };
    
switch sticks
    case 'yes'
        yl_sticks={ {'Velocity Bottom','(cm/s)'},...
             {'Velocity Surface','(cm/s)'},...
             {'Significant Wave','Height (m)'} };
end

%% Interpolate all data onto the array 'jd_int' for timeplt.m. Data gaps
% larger than 'maxGap' hours are not interpolated over
delta_t = 900/(3600*24); %for 15 min. data
jd_int = [startTime:delta_t:stopTime]';
maxGap = 1; %hr

% Get the color scheme to distinguish between data sites
col=get(0,'defaultaxescolororder');

% Initialize variables
u_bottom_int=zeros(length(jd_int),nFiles);
u_surface_int=zeros(length(jd_int),nFiles);
v_bottom_int=zeros(length(jd_int),nFiles);
v_surface_int=zeros(length(jd_int),nFiles);
hsig_int=zeros(length(jd_int),nFiles);
varList={'u','v'};
ncVarList={'u_1205','v_1206'};
nVars=length(ncVarList);
%%
for f=1:length(waveFiles) 
ncwaves=(fullfile(dataDirectory,waveFiles{f}));
    time=double(ncread(ncwaves,'time'));
    time2=double(ncread(ncwaves,'time2'));
    jd_waves = time + (time2/3600/1000/24);
    waveDelta_t=gmean(diff(jd_waves));
    hsig=squeeze(ncread(ncwaves,'wh_4061'));
    try
    theFillValue=ncreadatt(ncwaves,'wh_4061','_FillValue');
    catch
        theFillValue=1.0E35
    end
    hsig(hsig >= theFillValue)=nan;
    hsig(hsig <= 0)=nan;
    hsig_int(:,f) = smart_interp([jd_waves], [hsig],...
        jd_int, 3);     
end
for f=1:length(names)
    nc=(fullfile(dataDirectory,dataFiles{f}));
    time=double(ncread(nc,'time'));
    time2=double(ncread(nc,'time2'));
    jd= time + (time2/3600/1000/24);
    delta_t=gmean(diff(jd));
    depth=ncread(nc,'depth');
    site=names{f};
    for v=1:length(ncVarList)
        ncobj=squeeze(ncread(nc,ncVarList{v}));
        theFillValue=ncreadatt(nc,ncVarList{v},'_FillValue');
        data=ncobj;
        bads=find(data>=theFillValue);
        data(bads)=nan;
        %velocity in bottom bin
        eval([varList{v},'_bottom=data(1,:);'])
        %velocity in surface bin where surface is upward-most good bin
        eval([varList{v},'_surface=data(end-10,:);'])
        eval(['plot(',varList{v},'_surface);']);
        if v==length(ncVarList)
            % (optionally convert velocities from east, north coordinates
            % to longshore, cross-shore coordinates using rot.m. Assume
            % a shoreline orientation of 45 deg.
            switch coord
                case 'xl'
                    disp(['Converting from north,east coords '...
                        'to cross,long coords'])
                    [u_bottom,v_bottom]=rot(u_bottom,v_bottom,-45);
                    [u_surface,v_surface]=rot(u_surface,v_surface,-45);
                case 'en'
                    disp('Not converting coordinate systems')
                otherwise
                    disp('Not a valid coordinate system!')
            end
        end
    end


    % interpolate data onto 'jd_int'
    for var=1:nVars
        eval([varList{var},'_bottom_int(:,f)=smart_interp(jd, ',...
            varList{var},'_bottom, jd_int,3);'])
        eval([varList{var},'_surface_int(:,f)=smart_interp(jd, ',...
            varList{var},'_surface, jd_int,3);'])

    end
eval(['plot(',varList{v},'_surface_int(:,f));']);
    % (optionally filter the velocities)
    switch filt
        case 'yes'
            disp('Low-pass filtering u,v current velocities');
            for var=1:nVars
                eval([varList{var},'_bottom_int(:,f)=PL66TN(',...
                    varList{var},'_bottom_int(:,f), delta_t*24, 33);'])
                eval([varList{var},'_surface_int(:,f)=PL66TN(',...
                    varList{var},'_surface_int(:,f), delta_t*24, 33);'])
            end
        otherwise
            disp('Not filtering current velocities');
    end
end
    
    angle_r=zeros(size(u_bottom_int));
    velb=u_bottom_int  + v_bottom_int  * sqrt(-1);
    vels=u_surface_int  + v_surface_int  * sqrt(-1);
    velb=velb .* exp ( sqrt(-1) * angle_r);
    vels=vels .* exp ( sqrt(-1) * angle_r);
    qub=real(velb);
    qus=real(vels);
    qvb=imag(velb);
    qvs=imag(vels);
    dirb=(-atan2(qvb,qub)+(pi()/2))*180/pi();
    dirs=(-atan2(qvs,qus)+(pi()/2))*180/pi();
    dirb(dirb<0)=dirb(dirb<0)+360;
    dirs(dirs<0)=dirs(dirs<0)+360;
    magb=sqrt(qub.^2+qvb.^2);clear qub qvb velb;
    mags=sqrt(qus.^2+qvs.^2);clear qus qvs vels;


    %% %%%%%%%%%%%%%%%%%%%% U,V PLOTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % create U, V plots at the BOTTOM and TOP
    xdata=jd_int;
    % give more space to allow legend to be plotted         
    xlims(2)=xlims(2)+17;           
    clf
    set(gcf,'position',[100 25 850 1100])
    
    ha = tight_subplot(5,1,[.04 .3],[.1 .1],[.1 .1]) ;
    axes(ha(5));
    plot(jd_int,magb,'LineWidth',1.2);
    set(gca,'ylim',ylims(5,:))

      set(gca,'xlim',xlims,'xtick',xt);
  %   datetick('x',2,'keeplimits','keepticks');
       set(ha(5),'xticklabel',{[datestr(datenum(gregorian(xt)),'mmm dd yyyy')]})
   ylabel(yl_topbot{5});
    
    axes(ha(1));
    plot(jd_int,hsig_int,'LineWidth',1.2);
      [LEGH,OBJH,OUTH,OUTM]=legend(ha(1),wavenameFiles,'Location','NorthEast');%,'Orientation','Horizontal');
       set(LEGH,'Fontsize',8);
       legpos=get(LEGH,'position');
        set(LEGH,'position',[legpos(1) legpos(2)+.0056 legpos(3) legpos(4)])

    set(gca,'xlim',xlims,'xtick',xt,'XTickLabel','');
    set(gca,'Ylim',[ylims(1,:)]);
    ylabel(yl_topbot{1});
    
    axes(ha(2));plot(jd_int,dirs,'.');
   set(gca,'xlim',xlims,'xtick',xt,'XTickLabel','');
    ylabel(yl_topbot{2});
    
    axes(ha(3));plot(jd_int,mags,'LineWidth',1.2);
    ylabel(yl_topbot{3});
       set(gca,'xlim',xlims,'xtick',xt,'XTickLabel','');
 
    axes(ha(4));plot(jd_int,dirb,'.');
    set(gca,'xlim',xlims,'xtick',xt,'XTickLabel','');
    ylabel(yl_topbot{4});
   
% 
  [LEGH,OBJH,OUTH,OUTM]=legend(ha(3),names,'Location','SouthEast');%,'Orientation','Horizontal');
       set(LEGH,'Fontsize',8);
       legpos=get(LEGH,'position');
        set(LEGH,'position',[legpos(1) legpos(2)+.055 legpos(3) legpos(4)])

   
  [LEGH,OBJH,OUTH,OUTM]=legend(ha(5),names,'Location','SouthEast');%,'Orientation','Horizontal');
       set(LEGH,'Fontsize',8);
       legpos=get(LEGH,'position');
        set(LEGH,'position',[legpos(1) legpos(2)+.055 legpos(3) legpos(4)])


    % give it a title and export to .pdf
      switch filt
        case 'yes'
            tit1e=get(ha(5),'Title');
            set(tit1e,'String',['Low-Passed Water Flow Observations Fire Island 2014'])
            set(gcf,'PaperPositionMode','Auto')
            switch coord
                case 'en'
                    print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_adcp_en_tb_lp']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_adcp_xl_tb_lp']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
        otherwise
            tit1e=get(ha(1),'Title');
            set(tit1e,'String',['Water Flow Observations Fire Island 2014'])
            set(gcf,'PaperPositionMode','Auto')
            
            switch coord
                case 'en'
                    print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_adcp_magdir_sb']),'-dpdf');
                case 'xl'
                    print(fullfile(plotDirectory,['ofr2014-xxxx-FI14_adcp_magdirrot_sb']),'-dpdf');
                otherwise
                    disp('Not a valid coordinate system!')
            end
    end
