%% fetch case study behind Jan Mayen
clear all, close all
%SWIFTpath = './L1_telemetry/';
%SWIFTpath = './L2_postprocessed/';
SWIFTpath = './L3_merged/';
SWIFTfiles = dir([ SWIFTpath '*.mat']);

JanMayen.lat = 71.0;
JanMayen.lon = -8.5;
JanMayen.rotation = 145;

U1toU10 = 1.25; % multiple SWIFT winds by this to get U10;

flow = 0.5; % lower freq limit (kHz) for spectral plots
fhigh = 12; % higher freq limit (kHz) for spectral plots

fetch_mintime = datenum(2021,9,12,9,0,0);
fetch_maxtime = datenum(2021,9,12,10,0,0);

figure(11), clf
subplot(3,2,5)
load('../jmBathy.mat')
pcolor(lon,lat,z), shading flat, hold on
%plot(JanMayen.lon,JanMayen.lat,'ko','linewidth',3,'markersize',16)
axis([ -9 -5 70 71.5])
colormap bone
cb = colorbar; cb.Label.String = 'depth [m]';
caxis([-3000 500])


figure(12), clf
for fi = 1:length(SWIFTfiles)
    
    load([ SWIFTpath SWIFTfiles(fi).name ])
    
    if SWIFT(1).ID == 'SWIFT_11' | SWIFT(1).ID == 'SWIFT_12'
        symbol = 's'; linesymbol = ':s';
    else
        symbol = 'x'; linesymbol = ':x';
    end
    
    SWIFT( [SWIFT.time] < fetch_mintime | [SWIFT.time] > fetch_maxtime ) = [];
    
    [ x, y ] = GenericCoordinateTransform([SWIFT.lat], [SWIFT.lon], JanMayen.lat, JanMayen.lon, JanMayen.rotation ) ;
    %U10 = 14;
    U10 = U1toU10 * smooth([SWIFT.windspd])';
    Xfetch = 9.8 * x./ U10.^2;
    Hfetch = 9.8 * [SWIFT.sigwaveheight] ./ U10.^2;
    
    figure(11),
    subplot(3,1,1)
    plot([SWIFT.time],smooth(U1toU10*[SWIFT.windspd]),linesymbol,'linewidth',2,'markersize',12), hold on
    datetick, ylabel('U_{10} [m/s]')
    set(gca,'YLim',[0 20])
    set(gca,'FontSize',14,'FontWeight','demi')
    xlabel('12 Sept 2021')
    
    subplot(3,1,2)
    plot([SWIFT.time],smooth([SWIFT.sigwaveheight]),linesymbol,'linewidth',2,'markersize',12), hold on
    datetick, ylabel('H_s [m]')
    set(gca,'YLim',[0 5])
    set(gca,'FontSize',14,'FontWeight','demi')
    xlabel('12 Sept 2021')
    
    subplot(3,2,5)
    plot([SWIFT.lon],[SWIFT.lat],symbol,'linewidth',2,'markersize',12), hold on, grid on
    xlabel('lon'), ylabel('lat')
    set(gca,'FontSize',14,'FontWeight','demi')
    %plot(JanMayen.lon,JanMayen.lat,'ko','linewidth',8,'markersize',16)
    
    
    subplot(3,2,6)
    %plot(x,y,'x','linewidth',3,'markersize',12), hold on, axis equal, grid on
    %xlabel('x [m]'), ylabel('y [m]')
    plot(x./1000,smooth([SWIFT.sigwaveheight]),symbol,'linewidth',2,'markersize',12), hold on, grid on
    xlabel('x [km]'), ylabel('H_s [m]')
    set(gca,'FontSize',14,'FontWeight','demi')
    %plot(0,0,'ko','linewidth',3,'markersize',16)
    
    
    figure(12),
    plot(Xfetch, Hfetch, symbol,'linewidth',2,'markersize',12), hold on
    %xlabel('gX/U_{10}^2'), ylabel('gH_s/U_{10}^2')
    ylabel('$\hat{H}$','interp','latex')
    xlabel('$\hat{X}$','interp','latex')
    set(gca,'FontSize',14,'FontWeight','demi')
    
end

figure(12),
Hoffset = 0.05;
Hpredict = 0.002 * [1e1 1e4].^0.5;
plot([1e1 1e4],Hpredict+Hoffset,'k--','linewidth',3)
plot([1e1 1e4],Hpredict+2*Hoffset,'-','linewidth',3,'color',[.5 .5 .5])
plot([1e1 1e4],Hpredict+.25*Hoffset,'-','linewidth',3,'color',[.5 .5 .5])
legend('SWIFT 09','SWIFT 11','SWIFT 12','SWIFT 13','fetch law','Location','NortheastOutside')


figure(11)
subplot(3,2,5)
[landlat landlon] = find(z>0); plot(lon(landlon),lat(landlat),'.','color',[0.5 0.5 0])
axis([ -9 -5 70 71.5])

figure(11)
subplot(3,2,6)
axis([0 140 0 5])
plot(0,0,'o','color',[0.5 0.5 0],'linewidth',3,'markersize',16)
legend('SWIFT 09','SWIFT 11','SWIFT 12','SWIFT 13','Jan Mayen','Location','NortheastOutside')



% fetch case sound spectra
figure(13), clf
loglog(0,0, '-','color',[0.85,0.33,0.10],'linewidth',3 ), hold on % just for legend
loglog(0,0, '-','color',[0.93,0.69,0.13] ,'linewidth',3 ), hold on % just for legend
load('SWIFT11_11-16Sep2021_hydrophone.mat')
badf = find(f_kHz > 0.9 & f_kHz < 2);
PSD(badf,:) = NaN;
PSD(:, time < fetch_mintime | time > fetch_maxtime ) = [];
time( time < fetch_mintime | time > fetch_maxtime ) = [];
mPSD = mean(PSD,2);
%loglog(f_kHz, PSD, '-','color',[0.85,0.33,0.10],'linewidth',0.5), hold on
loglog(f_kHz, mPSD, '-','color',[0.85,0.33,0.10],'linewidth',3 ), hold on
loglog([f_kHz(badf(1)-1) f_kHz(badf(end)+1)], [mPSD(badf(1)-1) mPSD(badf(end)+1)], ':','color',[0.85,0.33,0.10],'linewidth',3 ), hold on


load('SWIFT12_11-25Sep2021_hydrophone.mat')
badf = find(f_kHz > 0.9 & f_kHz < 2);
PSD(badf,:) = NaN;
PSD(:, time < fetch_mintime | time > fetch_maxtime ) = [];
time( time < fetch_mintime | time > fetch_maxtime ) = [];
mPSD = mean(PSD,2);
%loglog(f_kHz, PSD, '-','color',[0.93,0.69,0.13] ,'linewidth', 0.5 ), hold on
loglog(f_kHz, mPSD, '--','color',[0.93,0.69,0.13] ,'linewidth',3 ), hold on
loglog([f_kHz(badf(1)-1) f_kHz(badf(end)+1)], [mPSD(badf(1)-1) mPSD(badf(end)+1)], ':','color',[0.93,0.69,0.13],'linewidth',3 ), hold on


set(gca,'FontSize',14,'FontWeight','demi')
xlabel('f [kHz]')
ylabel('Ambient sound, S(f) [dB re 1uPa^2/Hz]')
legend('SWIFT 11 (developing)','SWIFT 12 (fully developed)','Location','Northeast')
axis([flow fhigh 40 85 ])

print -dpng fetchcase_soundspectra.png

%% wave spectra
figure(14), clf
load([ SWIFTpath 'SWIFT11_NORSEpilot2021.mat'])
SWIFT( [SWIFT.time] < fetch_mintime | [SWIFT.time] > fetch_maxtime ) = [];
for si=1:length(SWIFT)
    f = SWIFT(si).wavespectra.freq;
    E(:,si) = SWIFT(si).wavespectra.energy;
end
loglog(f, nanmean(E'),'color',[0.85,0.33,0.10],'linewidth',3 ), hold on

load([ SWIFTpath 'SWIFT12_NORSEpilot2021.mat'])
SWIFT( [SWIFT.time] < fetch_mintime | [SWIFT.time] > fetch_maxtime ) = [];
for si=1:length(SWIFT)
    f = SWIFT(si).wavespectra.freq;
    E(:,si) = SWIFT(si).wavespectra.energy;
end
loglog(f, nanmean(E'),'--','color',[0.93,0.69,0.13],'linewidth',3 ), hold on

set(gca,'FontSize',14,'FontWeight','demi')
xlabel('f [Hz]')
ylabel('Wave energy, E(f) [m^2/Hz]')
axis([1e-2 1e0 5e-3 1e2])

legend('SWIFT 11 (developing)','SWIFT 12 (fuly developed)','Location','NorthEast')
%text(0.25,1e-2,'mss')
text(0.2,7e-3,'<-- equil. range -->')
%text(0.1,1e-2,'peak')
text(0.08,7e-3,'<-- peak range -->')

print -dpng fetchcase_wavespectra.png