% determine TR12vars
cd('D:\Data\LanderLT\ADCP\Data\mlab\transport')

load('TR12')
load('trfl')
load('vel_field')
load('D:\Data\LanderLT\MEEALL\mlab\dataset\landerDAY')
load('D:\Data\LanderLT\MEEALL\mlab\dataset\T12data')
load('D:\Data\LanderLT\MEEALL\mlab\dataset\burst2')
load('D:\Data\LanderLT\ADCP\Data\mlab\processing\ADCPD')
load('D:\Data\LanderLT\MEEALL\mlab\dataset\vektoravg')

%% determine gross, nett, duration per tidal period
TR12.ebb = NaN(1,length(TR12.ebst));
TR12.flood = NaN(1,length(TR12.ebst));
TR12.ebbs = NaN(1,length(TR12.ebst));
TR12.floods = NaN(1,length(TR12.ebst));
TR12.ebbb = NaN(1,length(TR12.ebst));
TR12.floodb = NaN(1,length(TR12.ebst));
TR12.net = NaN(1,length(TR12.ebst));
TR12.nets = NaN(1,length(TR12.ebst));
TR12.netb = NaN(1,length(TR12.ebst));
TR12.gr = NaN(1,length(TR12.ebst));
TR12.grs = NaN(1,length(TR12.ebst));
TR12.grb = NaN(1,length(TR12.ebst));
TR12.ebbtime = NaN(1,length(TR12.ebst));
TR12.floodtime = NaN(1,length(TR12.ebst));
TR12.nume = NaN(1,length(TR12.ebst));
TR12.numf = NaN(1,length(TR12.ebst));
TR12.sh = NaN(1,length(TR12.ebst));
timestep = 600; %seconds
nl = 41;
sh = ADCPD.dep/nl;

                                                            % the units of trfli are 10^3 mg m-2 s-1
for k = 1:length(TR12.ebst)
    if k == 801
        k = 801
    elseif ~isnan(TR12.ebst(k))
    eb = TR12.ebst(k):TR12.eben(k);
    fl = TR12.flst(k):TR12.flen(k);
        if TR12.ebst(k) == TR12.eben(k);
            TR12.ebb(k)      = 0;
            TR12.ebbs(k)     = 0;
            TR12.ebbb(k)     = 0;
            TR12.ebbtime(k)  = 0;
            TR12.nume(k)     = NaN;
        else
            TR12.ebb(k)     = -sum(sum(trfli3(:,eb)).*sh(eb))*timestep;  % 10^3 mg / m coastwith / tide
            TR12.ebbs(k)    = -sum(mean(trfli3(36:41,eb)))*timestep;  % 10^3 mg / m coastwith / tide         %.*sh(eb))
            TR12.ebbb(k)    = -sum(mean(trfli3(2:5,eb)))*timestep;  % 10^3 mg / m coastwith / tide           %.*sh(eb))
            TR12.ebbtime(k) = length(eb)*10/60;           %in hours
            TR12.nume(k)     = DAY.num(TR12.ebst(k)); 
        end
    TR12.flood(k)    = sum(sum(trfli3(:,fl)).*sh(fl))*timestep;
    TR12.floods(k)   = sum(mean(trfli3(36:41,fl)))*timestep;  % 10^3 mg / m coastwith / tide                %.*sh(fl)
    TR12.floodb(k)   = sum(mean(trfli3(2:5,fl)))*timestep;   % 10^3 mg / m coastwith / tide                %.*sh(fl)

    TR12.net(k)      = -TR12.ebb(k)+TR12.flood(k);
    TR12.nets(k)     = sum(mean(trfli3(36:41,[eb fl])))*timestep;
    TR12.netb(k)     = sum(mean(trfli3(2:5,[eb fl])))*timestep;

    TR12.gr(k)       = min(TR12.ebb(k),TR12.flood(k));
    TR12.grs(k)      = min(TR12.ebbs(k),TR12.floods(k));
    TR12.grb(k)      = min(TR12.ebbb(k),TR12.floodb(k));
    TR12.floodtime(k)= length(fl)*10/60;         %in hours
    TR12.tidetime(k) = TR12.ebbtime(k)+TR12.floodtime(k);
    TR12.numf(k)     = DAY.num(TR12.flst(k)); 
    TR12.sh(k)       = nanmean(sh([eb fl]));
    end
end

figure;     
plot(TR12.numf,TR12.ebb)
hold on
plot(TR12.numf,TR12.flood)
plot(TR12.numf,TR12.net)
hline(0,'k:')
legend('ebb','flood','residual')
plot_deploylines(1,0)
datetick


figure
hist(TR12.net,100)
vline(0,'k:')
xx = sort(TR12.net); xy = find(xx>0,1,'first') %--> 1/3 van de tijd is de residueel in ebb richting

%% other parameters
A = isfinite(VEK.wdir);
VEK.wdir(A == 0) = NaN;

TR12.tmax   = NaN(1,length(TR12.ebst));
TR12.tmean  = NaN(1,length(TR12.ebst));
TR12.tmed   = NaN(1,length(TR12.ebst));
TR12.wdir   = NaN(1,length(TR12.ebst));
TR12.wspd   = NaN(1,length(TR12.ebst));
TR12.wmod   = NaN(1,length(TR12.ebst));
TR12.we     = NaN(1,length(TR12.ebst));
TR12.sei    = NaN(1,length(TR12.ebst));
TR12.wspda  = NaN(1,length(TR12.ebst));
TR12.wdira  = NaN(1,length(TR12.ebst));
TR12.wmoda  = NaN(1,length(TR12.ebst));
TR12.Hdir   = NaN(1,length(TR12.ebst));
TR12.Hmn    = NaN(1,length(TR12.ebst));
TR12.Hmax   = NaN(1,length(TR12.ebst));
TR12.Hmod   = NaN(1,length(TR12.ebst));

for k = 1:length(TR12.ebst)
    if ~isnan(TR12.ebst(k))
    eb = TR12.ebst(k):TR12.eben(k);
    fl = TR12.flst(k):TR12.flen(k);
    TR12.tmean(k)    = nanmean(burst.tmax([eb fl]));
    TR12.tmed(k)    = nanmedian(burst.tmax([eb fl]));
    TR12.tmax(k)    = nanmax(burst.tmax([eb fl]));
    [~,IA,~]        = intersect(VEK.IJind,[eb fl]);
    [TR12.wdir(k),TR12.wspd(k)]    = mean_angle_speed(VEK.wdir(IA),VEK.spd(IA));     %average wind direction and speed, following vector summation
    [TR12.wdira(k),TR12.wspda(k)]    = mean_angle_speed(VEK.wdira(IA),VEK.spda(IA));     %average wind direction and speed, following vector summation
%     TR12.wspd(k)    = nanmean(VEK.spd(IA));         %mean wind speed
    [TR12.wmod(k),TR12.we(k)]    = wind_mode(VEK.wdir(IA),VEK.spd(IA)); %wind mode and wind energy
    [TR12.wmoda(k),~]    = wind_mode(VEK.wdira(IA),VEK.spda(IA)); %wind mode and wind energy
    [TR12.Hdir(k),~]    = mean_angle_speed(VEK.dir(IA),VEK.hgt(IA));     %average wave direction, following vector summation
    [TR12.Hmod(k),~]    = wind_mode(VEK.dir(IA),VEK.hgt(IA)); %wind mode and wind energy
    TR12.Hmn(k)     = nanmean(VEK.hgt(IA));         % mean wave height
    TR12.Hmax(k)    = nanmax(VEK.hgt(IA));          % max wave height per tide
    if ismember(DAY.M(TR12.ebst(k)),4:9)            %season definition; summer =1 (Apr-Sept), winter = 0 (Oct-Mar)
        TR12.sei(k) = 1;
    else
        TR12.sei(k) = 0;
    end
    end
end
TR12.wdir2cst = mod(TR12.wdir-96,360); % = E,Soff,S, Son, W, Non, N, Noff, E 
%% cumulative average or something
TR12.cum   = NaN(1,length(TR12.ebst));
TR12.cummdn = NaN(1,length(TR12.ebst));
TR12.cummn = NaN(1,length(TR12.ebst));

% index for 2012 == 521
kk = find(TR12.nume>DAY.num(find(DAY.Y==2012,1,'first')),1,'first');

% cummdn  = nanmedian(TR12.net);
cummdn11 = nanmedian(TR12.net(1:kk-1));
cummdn12 = nanmedian(TR12.net(kk:end));
% precum = TR12.net - nanmedian(TR12.net);
precum11 = TR12.net - cummdn11;
precum12 = TR12.net - cummdn12;
cummn11 = nanmean(TR12.net(1:kk-1));
cummn12 = nanmean(TR12.net(kk:end));

%2011
TR12.cum(1)    = precum11(1);
TR12.cummdn(1) = cummdn11; 
TR12.cummn(1)  = cummn11;
cum = precum11(1);
cummdn = cummdn11;
cummn = cummn11;
for k = 2:kk-1  %length(TR12.ebst)
    if ~isnan(precum11(k))
        cum = cum + precum11(k);
        cummdn = cummdn+cummdn11; 
        cummn = cummn+cummn11;
        TR12.cum(k)    = cum;
        TR12.cummdn(k) = cummdn; 
        TR12.cummn(k)  = cummn; 
%         display(['not NaN k = ' num2str(k)])
%         precum(k)
%     else % == NOT necessary to define, as matrices are NaNs by definition
% %         display(['NaN k = ' num2str(k)])
%         TR12.cum(k)    = NaN;
%         TR12.cummdn(k)  = NaN; 
%         TR12.cummn(k) = NaN;
        
    end
end

%2012
TR12.cum(kk)    = precum12(1);
TR12.cummdn(kk) = cummdn12; 
TR12.cummn(kk)  = cummn12; 
cum    = precum12(1);
cummdn = cummdn12;
cummn  = cummn12;
for k = kk+1:numel(TR12.ebst)
    if ~isnan(precum12(k))
        cum    = cum + precum12(k);
        cummdn = cummdn+cummdn12;
        cummn  = cummn+cummn12;
        TR12.cum(k)    = cum;
        TR12.cummdn(k)  = cummdn; 
        TR12.cummn(k)  = cummn; 
%         display(['not NaN k = ' num2str(k)])
%         precum(k)
%     else % == NOT necessary to define, as matrices are NaNs by definition
% %         display(['NaN k = ' num2str(k)])
%         TR12.cum(k)    = NaN;
%         TR12.cummdn(k)  = NaN; 
%         TR12.cummn(k) = NaN;
        
    end
end

% TR12.cum = cumsum(precum);

%% Determine type of tide --> tau ranges
% H2 = find(VEK.IJMHs > 200 & VEK.IJMHs < 300);
% H3 = find(VEK.IJMHs > 300 & VEK.IJMHs < 400);
% H4 = find(VEK.IJMHs > 400);
TR12.tcse = NaN(1,length(TR12.ebst));
TR12.hcse = NaN(1,length(TR12.ebst));


% Tc = find(TR12.tmed <= 0.7);
% TR12.tcse(Tc) = 1;
% Ti = find(TR12.tmed >= 0.7 & TR12.tmed < 1.7);
% TR12.tcse(Ti) = 2;
% Ts = find(TR12.tmed >= 1.7);
% TR12.tcse(Ts) = 3;
% Te = (TR12.tmax > 30);
% TR12.tcse(Te) = 4;
% 
% %Hs (for period 16 and the rest)
% Hc = find(TR12.Hmn <= 0.8);
% TR12.hcse(Hc) = 1;
% Hi = find(TR12.Hmn >= 0.8 & TR12.Hmn < 1.5);
% TR12.hcse(Hi) = 2;
% Hs = find(TR12.Hmn >= 1.5);
% TR12.hcse(Hs) = 3;
% He = find(TR12.Hmn >= 3.5); %3.5
% TR12.hcse(He) = 4;

Tc = find(TR12.tmed <= 1.7);
TR12.tcse(Tc) = 1;
Ts = find(TR12.tmed >= 1.7);
TR12.tcse(Ts) = 2;

% %Hs (for period 16 and the rest)
% Hc = find(TR12.Hmn < 0.8);  % this was 1.5 in script version August 2019
% TR12.hcse(Hc) = 1;
% Hs = find(TR12.Hmn >= 1.5);  % this was 1.5 in script version August 2019
% TR12.hcse(Hs) = 3;
% Hi = find(TR12.Hmn >= 0.8 & TR12.Hmn < 1.5);
% TR12.hcse(Hi) = 2;
%Hs (for period 16 and the rest)
Hc = find(TR12.Hmn < 0.8);  % this was 1.5 in script version August 2019
TR12.hcse(Hc) = 1;
Hs = find(TR12.Hmn >= 0.8);  % this was 1.5 in script version August 2019
TR12.hcse(Hs) = 2;

%% save
save('TR12','TR12')


%% TR12 back to minute interval or to hour (wind) interval
trfnet.min = NaN(1,length(DAY.num));
trfnet.hr  = NaN(1,numel(VEK.IJind));

for k = 1:length(TR12.ebst)
    if ~isnan(TR12.ebst(k))
        rng12 = TR12.ebst(k):TR12.flen(k);
        trfnet.min(rng12) = TR12.net(k);
        vst = find(VEK.IJind >= rng12(1),1,'first');
        ven = find(VEK.IJind <= rng12(end),1,'last');
        trfnet.hr(vst:ven) = TR12.net(k);
        
    end
end

save('trfnet','trfnet')

%% histogram residual

ind1 = find(TR12.tmax < 2);
ind2 = find(TR12.tmax > 4);
ind3 = setxor(1:length(TR12.tmax),[ind1 ind2]);
v = -10:2:100;
[n1,p1] = hist(TR12.net(ind1)/1e6,v); n1 = n1/length(ind1)*100;
[n2,p2] = hist(TR12.net(ind2)/1e6,v); n2 = n2/length(ind2)*100;
[n3,p3] = hist(TR12.net(ind3)/1e6,v); n3 = n3/length(ind3)*100;

figure; 
subplot(1,3,3)
bar([p1;p2;p3]',[n1;n3;n2]','grouped')
ylabel('%')
legend('\tau < 2','2 < \tau < 4','\tau > 4')
xlim([-10 40])
vline(0,'k--')
xlabel('Mass transport (g m^-^1 coast per tide)')
title('Residual mass transport')

nansum(TR12.net)/1e6
%% histogram ebb

% ind1 = find(TR12.tmax < 2);
% ind2 = find(TR12.tmax > 4);
% ind3 = setxor(1:length(TR12.tmax),[ind1 ind2]);
% v = 0:2:100;
[n1,p1] = hist(TR12.ebb(ind1)/1e6,v); n1 = n1/length(ind1)*100;
[n2,p2] = hist(TR12.ebb(ind2)/1e6,v); n2 = n2/length(ind2)*100;
[n3,p3] = hist(TR12.ebb(ind3)/1e6,v); n3 = n3/length(ind3)*100;

subplot(1,3,1); bar([p1;p2;p3]',[n1;n3;n2]','grouped')
ylabel('%')
legend('\tau < 2','2 < \tau < 4','\tau > 4')
xlim([-1 20])
% vline(0,'k--')
xlabel('Mass transport (g m^-^1 coast per tide)')
title('Ebb mass transport')
%% histogram flood

% ind1 = find(TR12.tmax < 2);
% ind2 = find(TR12.tmax > 4);
% ind3 = setxor(1:length(TR12.tmax),[ind1 ind2]);
% v = 0:100;
[n1,p1] = hist(TR12.flood(ind1)/1e6,v); n1 = n1/length(ind1)*100;
[n2,p2] = hist(TR12.flood(ind2)/1e6,v); n2 = n2/length(ind2)*100;
[n3,p3] = hist(TR12.flood(ind3)/1e6,v); n3 = n3/length(ind3)*100;

subplot(1,3,2); bar([p1;p2;p3]',[n1;n3;n2]','grouped')
ylabel('%')
legend('\tau < 2','2 < \tau < 4','\tau > 4')
xlim([-1 40])
% vline(0,'k--')
xlabel('Mass transport (g m^-^1 coast per tide)')
title('Flood mass transport')

%% per TAU
figure; hist(TR12.net(ind1)/1e6,100)
vline(0,'r--')
xlabel('Residual mass transport (g m^-^1 coast per tide)')
ylabel('count')
title('\tau < 2 Nm^-^1')

figure; hist(TR12.net(ind2)/1e6,50)
vline(0,'r--')
xlabel('Residual mass transport (g m^-^1 coast per tide)')
ylabel('count')
title('\tau > 4 Nm^-^1')

figure; hist(TR12.net(ind3)/1e6,50)
vline(0,'r--')
xlabel('Residual mass transport (g m^-^1 coast per tide)')
ylabel('count')
title('intermediate')


%% 
hh = hist(TR12.net(ind1)/1e6,100)
hh(hh == 0) = NaN;
mode(hh)
[hh,jj] = hist(TR12.net(ind1)/1e6,100)
[x,y] = max(hh)
jj(33)
nanmedian(TR12.net(ind1)/1e6)
find(jj>0,1,'first')
sum(hh(1:30))
sum(hh(31:end))