clear all
close all
% Paths to the Lab data
root='\PhySeaCS\Lab_data\';
RepRes='\PhySeaCS\0_Results\';
Study = [{'Atkinson'},{'Beuzen'},{'Bayle'}];
% Experiments IDs
Runs.AtkRuns=[{'A1'},{'A2'},{'A3'},{'E1'},{'E2'},{'E3'},{'E4'},{'E5'},{'C1'},{'C2'},{'C3'},{'C4'},{'C5'},{'C6'}];
Runs.BeuRuns=[{'E1'},{'E2'}];
Runs.BayRuns={'E1'};

% Laboratory Data
SLR=[]; hs=[]; tp=[]; d=[]; L=[];
dPF=[];  WR_end=[]; WR_t=[];
WR_endRaw=[]; PF_Raw=[]; WR_tRaw=[];

for k=1:length(Study)
    stu=cell2mat(Study(k));
    Exp=Runs.(sprintf('%sRuns', stu(1:3)));

    load([root sprintf('%s/Analysis/WR.mat',stu)],'WR','WR_raw');
    load([root sprintf('%s/Analysis/PF.mat',stu)],'PF','PF_raw');
    WR_end=[WR_end; WR];       WR_endRaw=[WR_endRaw; WR_raw];
    dPF=[dPF; PF];             PF_Raw=[PF_Raw; PF_raw];
    load([root sprintf('%s/Analysis/WR_Part.mat',stu)],'WR', 'WR_raw');
    WR_t=[WR_t; WR];     WR_tRaw=[WR_tRaw; WR_raw];
   for i=1:length(Exp)
        iExp=cell2mat(Exp(i));
        load ([root sprintf('%s/Analysis/dataExp_%s.mat',stu,iExp)])
        if i>8
            hs=[hs;mean(Hydro.H)];
        else
            hs=[hs;Hydro.H];
        end
        tp=[tp;Hydro.T];
        SLR=[SLR;Hydro.SLR];
        d=[d;Hydro.d];
        L=[L;Hydro.L];
    end
end

ws=[0.037*ones(14,1); 0.044*ones(6,1); 0.041*ones(4,1)];   
R=[0.11*ones(8,1); 0.08*ones(6,1); 0.135*ones(6,1); 0.84*ones(4,1)]; % BSM Scaling Run up (mean for each set of experiments)

clearvars -except 'root' 'RepRes' 'd' 'tp' 'hs' 'SLR' 'dPF' 'WR_end'  'WR_t' 'WR_Hot' 'ws' 'L' 'WR_endRaw' 'PF_Raw' 'WR_tRaw' 'Runs' 'R' 'Study'

%% Calculate Hydrodynamic quantities
g=9.81; % Gravity acceleration
Lo=(g.*tp.^2)./(2*pi()); % deepwater wavelength [m]
d=d+SLR;
Ho= []; 
% Invert Linear wave shoaling to derive deep-water limit wave height
for i = 1: length(hs)
    kh=(d(i))*2*pi()/L(i);
    ks=((2*cosh(kh)^2)/(sinh(2*kh)+2*kh))^0.5;
    iHo=hs(i)/ks;
    Ho= [Ho;iHo] ; %Offshore wave height [m]
end

do=Lo./2; % deepwater wave steepness [-]
Om=Ho./(ws.*tp); % Dimensionless Fall Velocity % %NB for A1 Omega calculated for H=0.05 is 0.8333 (1 for H=0.06)

%% Plot data
rSLR=SLR./Ho;         
BSM_var=[dPF,WR_t];               
rwv=2*Ho.*Om./Lo;        
%%  WR Final vs WR Partial
pWR=rSLR\WR_t;   pWRWR=WR_end\WR_t;
pPF=[rSLR(2:6);rSLR(8:end)]\[dPF(2:6);dPF(8:end)];
pBr=SLR\(PF_Raw+WR_endRaw);
%% Scatter Plots by dataset
for j=1:2
    rWR=BSM_var(:,j); % Select PF or WR(part)
    iPSize=50*ones(length(rwv),1); % Marker sizes
    fig=figure;
    subplot(1,2,1); grid on; hold on;
    scatter(rSLR(1:8),rWR(1:8),iPSize(1:8),rwv(1:8),'^','filled');
    scatter(rSLR(9:14),rWR(9:14),iPSize(9:14),rwv(9:14),'o','filled');
    scatter(rSLR(15:20),rWR(15:20),iPSize(15:20),rwv(15:20),'s','filled');
    scatter(rSLR(21:end),rWR(21:end),iPSize(21:end),rwv(21:end),'d','filled');
    colorbar; caxis([0, 0.5]); colormap Jet
    xlabel('SLR* [-]'); ylabel('WR* [-]') %'Wave Power (H^2 T) [m^2 s]'

    iiPSize=50*ones(length(rwv),1);
    subplot(1,2,2); hold on;
    scatter(rwv(1:8),rWR(1:8),iiPSize(1:8),rSLR(1:8),'^','filled');
    scatter(rwv(9:14),rWR(9:14),iiPSize(9:14),rSLR(9:14),'o','filled');
    scatter(rwv(15:20),rWR(15:20),iiPSize(15:20),rSLR(15:20),'s','filled');
    scatter(rwv(21:end),rWR(21:end),iiPSize(21:end),rSLR(21:end),'d','filled');
    colorbar; caxis([-0.6,0.6]); colormap Jet
    xlabel('P* [-]'); ylabel('WR*_p_a_r_t [-]')
    grid on;
end
%% Part WR vs Final WR
fig=figure; grid on; hold on;
plot (WR_end(1:8) ,WR_t(1:8),'^b');    plot (WR_end(9:14) ,WR_t(9:14),'oc')
plot (WR_end(15:20) ,WR_t(15:20),'xr')
plot (WR_end(21:end) ,WR_t(21:end),'^g')
plot([-0.2;WR_end;0.5],pWRWR*[-0.2;WR_end;0.5],'-k')
xlabel('WR* (Final) [m]'); ylabel('WR* (Partial) [m]');