% script to explot Papa data (PAL + WR)
%
% J. Thomson, Apr 2021
%

clear all, close all

runmatch = false; % binary flag to avoid matching the data again
fixedfrangeflag = true; % calc mss from dynamic or fixed tail

flow = 0.5; % freq limits for plots (kHz)
fhigh = 15; % freq limits for plots (kHz)

findex = 26; % freq for examples, 18 is 4.6 kHz, 24 is 10.4 kHZ

%% load data
load('PAPA_PALdata_2010_to_2012.mat')

%% quality control
QCf = 26;  % weird spike at 11.4 kHz (index 25)
ScN(:,QCf) = [];
f64(QCf) = [];
QCf = 25;  % weird spike at 11.4 kHz (index 25)
ScN(:,QCf) = [];
f64(QCf) = [];

maxf = 20; % kHz

ScN(:,f64>maxf)= NaN;


if runmatch

    tic

    load('PAPA_waveriderspectra_Jun2010-Jun2012_SWIFTformat.mat')

    load('PAPA_PMELwinds_Jun2010-Sep2012.mat')  % is this already U10?


    %% plot time series

    figure(1), clf

    subplot(3,1,1)
    pcolor(tn+datenum(2010,0,0),f64,ScN')
    shading flat
    ylabel('f [kHz]')
    set(gca,'Xlim',[datenum(2010,6,1) datenum(2012,6,1)])
    set(gca,'Ylim',[0 maxf])
    datetick('x','mmm','keeplimits')
    caxis([0 80])
    set(gca,'fontsize',14,'fontweight','demi')
    cb = colorbar; cb.Label.String = 'S(f) [dB]';

    subplot(3,1,2)
    plot(windtime,windspd,'k.')
    ylabel('U_{10} [m/s]')
    set(gca,'Xlim',[datenum(2010,6,1) datenum(2012,6,1)])
    datetick('x','mmm','keeplimits')
    set(gca,'fontsize',14,'fontweight','demi')
    cb = colorbar; set(cb,'visible','off')

    subplot(3,1,3)
    plot([WR.time],[WR.sigwaveheight],'k.')
    ylabel('H_s [m]')
    set(gca,'Xlim',[datenum(2010,6,1) datenum(2012,6,1)])
    datetick('x','mmm','keeplimits')
    set(gca,'fontsize',14,'fontweight','demi')
    cb = colorbar; set(cb,'visible','off')
    xlabel('2010-2012')
    
    print -dpng papatimeseries.png


    %% calculate mean square slope (mss) 
    % from whole spectra (first), then tail, then tail normalized by f range
    [msstotal, mss, mssnorm] = SWIFTmss( WR, fixedfrangeflag, false ); 

    %% match times

    PALtime = tn+datenum(2010,0,0);

    for pi=1:length(PALtime),
        
        [tdiff matchindex] = min( abs( PALtime(pi) - [WR.time] ) );
        if tdiff<1/24,
            matchWRindices(pi) = matchindex;
            matchedHs(pi)= WR(matchindex).sigwaveheight;
            matchedmsstotal(pi) = msstotal(matchindex);
            matchedmss(pi) = mss(matchindex);
            matchedmssnorm(pi) = mssnorm(matchindex);
        else
            matchWRindices(pi) = 0;
            matchedHs(pi)= NaN;
            matchedmsstotal(pi) = NaN;
            matchedmss(pi) = NaN;
            matchedmssnorm(pi) = NaN;
        end
    end

    for pi=1:length(PALtime),

        [tdiff matchindex] = min( abs( PALtime(pi) - windtime ) );

        if tdiff<1/24,
            matchwindindices(pi) = matchindex;
            matchedwind(pi) = windspd(matchindex);
        else
            matchwindindices(pi) = 0;
            matchedwind(pi) = NaN;
        end
    end

    save matchedwindandwaves.mat matched*

    toc

end

%% wind and wave relations

load matchedwindandwaves_mssdynamictail.mat matched*
%matchedwind = 1.2*matchedwind;  % adjust to U10?

U10 = [3:3:21];
fdHhat = 0.15; % fully developed non-dim wave heights
fdHs = fdHhat * U10.^2 * 9.8^-1;  % fully developed wave heights

figure(2), clf
%subplot(1,2,1)
binscatter(matchedwind,matchedHs)
set(gca,'fontsize',16,'fontweight','demi')
xlabel('wind speed, U_{10} [m/s]')
ylabel('wave height, H_s [m]')
hold on
plot(U10,1.2*fdHs+0.5,'k-')
%legend('data','P-M limit')

% subplot(1,2,2)
% binscatter(matchedmss,matchedHs)
% set(gca,'fontsize',16,'fontweight','demi')
% xlabel('wind speed, U [m/s]')
% ylabel('Wave steepness, mss []')

%print -dpng StationP_windandwaves.png
print -dpng papa_wavedev.png


%% binned PAL spectra

Hhat = (9.8*matchedHs) ./ matchedwind.^2; 
%Hhatmss = matchedHs ./ (9.8*matchedmssnorm.^2); 

figure(5), clf
colormap jet
cmap = colormap;
%set(gca,'color',[0.8 .8 .8])
for ui = 1:length(U10)
    thesewinds = abs(matchedwind-U10(ui))< median(diff(U10)/2);
    %smallwaves = abs(matchedwind-U10(ui))< median(diff(U10)) & matchedHs < nanmean(matchedHs(thesewinds));
    smallwaves = abs(matchedwind-U10(ui))< median(diff(U10)/2) & Hhat<fdHhat; sum(smallwaves);
    %largewaves = abs(matchedwind-U10(ui))< median(diff(U10)) & matchedHs > nanmean(matchedHs(thesewinds));
    largewaves = abs(matchedwind-U10(ui))< median(diff(U10)/2) & Hhat>fdHhat; sum(largewaves);

    ScN_inU10bins(ui,:) = nanmean( ScN(thesewinds,:) );
    ScN_inU10bins_smallHhat(ui,:) = nanmean( ScN(smallwaves,:) );
    ScN_inU10bins_smallHhat_std(ui,:) = nanstd( ScN(smallwaves,:) ) ./ sum(smallwaves)^.5;  % standard error, not std deviation
    ScN_inU10bins_largeHhat(ui,:) = nanmean( ScN(largewaves,:) );
    ScN_inU10bins_largeHhat_std(ui,:) = nanstd( ScN(largewaves,:) ) ./ sum(largewaves)^.5;  % standard error, not std deviation
    mss_U10bins(ui) = nanmean(matchedmssnorm(thesewinds));
    stdmss_U10bins(ui) = nanstd(matchedmssnorm(thesewinds));
    pickf_U10bins(ui) = mean(ScN(thesewinds,findex));
    
    % plotting
    cindex = ceil( ui./length(U10) * 256);
    %loglog(f64,ScN_inU10bins(ui,:),'color',cmap(cindex,:),'linewidth',2), hold on
    loglog(f64,ScN_inU10bins_smallHhat(ui,:),'color',cmap(cindex,:),'linewidth',2), hold on
    %loglog(f64,ScN_inU10bins_largeHhat(ui,:),'--','color',cmap(cindex,:),'linewidth',2), hold on

    % fitting
    fforfit = f64>0.5 & f64<15;
    fit = polyfit(log(f64(fforfit)), log(ScN_inU10bins(ui,fforfit)),1);
    dBerror_wind(ui,:) =  ScN_inU10bins(ui,fforfit) -  exp(fit(2))*f64(fforfit).^fit(1);
end
set(gca,'fontsize',16,'fontweight','demi')
xlabel('f [kHz]')
ylabel('Ambient sound, S(f) [dB re 1uPa^2/Hz]')
title('Ensembles by wind speed, U_{10} [m/s]')
axis([flow fhigh 30 85])
for ui=1:length(U10), 
        cindex = ceil( ui./length(U10) * 256);
        loglog(f64,ScN_inU10bins_largeHhat(ui,:),'--','color',cmap(cindex,:),'linewidth',2), hold on
end
legend(num2str(U10'),'Location','NortheastOutside')


print -dpng papaspectra_windspdbins.png

%% Sound level vs wind speed in specific bands

selectfs = [8 21 27];

figure(11), clf

for ii=1:length(selectfs)
    subplot(  1, length(selectfs), ii)
    plot(U10,ScN_inU10bins_smallHhat(:,selectfs(ii)),'k-','linewidth',3), hold on
    plot(U10,ScN_inU10bins_largeHhat(:,selectfs(ii)),'--','color',[0.5 0.5 0.5],'linewidth',3), hold on
    grid
    errorbar(U10,ScN_inU10bins_smallHhat(:,selectfs(ii)),ScN_inU10bins_largeHhat_std(:,selectfs(ii)),'k','linewidth',1.5)
    errorbar(U10,ScN_inU10bins_largeHhat(:,selectfs(ii)),ScN_inU10bins_largeHhat_std(:,selectfs(ii)),'color',[0.5 0.5 0.5],'linewidth',1.5)
    %plot(matchedwind,ScN(:,findex),'.'), hold on
    %plot(U10,ScN_inU10bins(:,findex),'ks','linewidth',2), grid
    %fd = (9.8*matchedHs) ./ matchedwind.^2 > 0.15;
    %fd = matchedHs ./ (9.8*matchedmssnorm.^2) > 30;
    %plot(matchedwind(fd), ScN( fd,findex),'g.'), hold on % more than fully developed
    %plot(matchedwind(~fd),ScN( ~fd ,findex),'r.'), hold on % less than fully developed
    set(gca,'fontsize',14,'fontweight','demi')
    xlabel('wind speed [m/s]')
    title([num2str(f64(selectfs(ii)), 2) ' kHz'])
    ylabel('Ambient Sound [dB re 1uPa^2/Hz]')
    legend('developing','fully developed','Location','Southeast')
end

print('-dpng',['papa_U10bins_selectfkHzbands.png'])