% Routine to wave statistics computed for the entire cblast-2 period over
% successive increments.  Plots with three panels are generated for
% each increment.  The top panel has wind sea and dominant swell
% period superimposed on a 1-d spectral plot of wave energy.  The lower
% panels display wind speed and wind sea height, and wind direction and
% wind sea direction.  
% Required inputs are: a file containing a time series of the 1-d spectra
% and a file with the computed wave statistics.

clear
dday = 3;  % time increment (days) of each plot

load spec_1d
ydayall = year_day;
fr = frequency;
spec1dall = spec_1d;
clear year_day frequency spec_1d;
%load c:\cblast2\results\all\stats\stat_all_1g.mat   % load the statistics file; the last letter in the file name designates the parameter set.
load wave_statistics.mat   % load the statistics file; the last letter in the file name designates the parameter set.

ytickn = [0.03 0.05 0.1 0.2 0.4 0.6 0.8];
ytick = log10(ytickn); % frequency ticks on  the y (logarithmic) axis

%ipl = 1;
%for iday = 211:dday:365
for iday = 280:dday:300
    drange = [iday iday+dday];
    jj = find(ydayall>=iday&ydayall<=iday+dday);
    yday = ydayall(jj);
    spec1d = spec1dall(:,jj);

    % raw data is pretty noisy, try smoothing by 3 in freq and time
    if(exist('spec_lp1'));clear spec_lp1 spec_lp2;end
    yday_lp = rmfltr( yday, length(yday), 3, 1 );
    for n = 1:length(yday)
       spec_lp1(:,n) = rmfltr(spec1d(:,n), length(fr), 3, 1 )';
    end
    f_lp = rmfltr( fr, length(fr), 3, 1 );
    for n = 1:length(f_lp)
        spec_lp2(n,:) = rmfltr( spec_lp1(n,:), length(yday), 3, 1 );
    end
    clear n

    clf;
    subplot(311)
    [C,H] = contourf( yday_lp,log10(f_lp), log10(spec_lp2));
    set( H, 'EdgeColor', 'none' );
    title('ADCP Filtered 1-D Spectra') 
    hold on
    set(gca,'xlim',drange)
    set(gca,'ytick',ytick)
    set(gca,'yticklabel',num2str(ytickn'))
    ylabel('frequency - Hz','fontsize',15,'fontangle','italic')


    fwindsea = 1./windseaper;
    h0(1)=plot(ydstat,log10(fwindsea),'color',[0.8 0.8 0.8],'linewidth',2)

    fdomswell = 1./domswellper;
    h0(2)=plot(ydstat,log10(fdomswell),'color',[0.5 0.5 0.5],'linewidth',2)
    l1 = legend(h0,'wind sea','dom swell');
    lpos = get(l1,'position');
    lpos(1) = .72; lpos(2) = .91;
    set(l1,'position',lpos);
    set(gca,'xticklabels',' ')

    subplot(312)
    h(1)=plot(ydstat,windsp,'k');hold on
    h(2)=plot(ydstat,windseahgt*10,'r');
    title('Wind Speed and Wind Sea Hgt.')
    ylabel('m/s and m','fontsize',15,'fontangle','italic')
    l2 = legend(h,'wind speed','wind sea hgt*10');
    set(gca,'xlim',drange)
    lpos = get(l2,'position');
    lpos(1) = .72; lpos(2) = .616;
    set(l2,'position',lpos);
    set(gca,'xticklabels',' ')


    subplot(313)
    h2(1)=plot(ydstat,winddir,'k');hold on
    h2(2)=plot(ydstat,windseadir,'r');
    title('Wind and Wind Sea Direction')
    set(gca,'ylim',[0 360],'ytick',[0 90 180 270 360]);
    ylabel('degrees','fontsize',15,'fontangle','italic')
    xlabel('year-day-2003','fontsize',15,'fontangle','italic')
    grid on
    set(gca,'xlim',drange)
    l3 = legend(h2,'wind dir','wind sea dir');
    lpos = get(l3,'position');
    lpos(1) = .72; lpos(2) = .316;
    set(l3,'position',lpos);
    pause
 %   saveas(gcf,['cont' int2str(ipl)],'tiff')
 %   ipl = ipl+1;
end