function [fig1 fig2] = xb_compare_profiles(x_xb, z_xb,x_data, z_data, tsglobal,zs0, name, xlim, ylim,skip,volume,bermslope,retreat)

% --- store iriginal data for plotting
z_xb_orig       = z_xb;
z_data_orig     = z_data;
tsglobal_orig   = tsglobal;

% --- skip first hour except for first datapoint
%I           = tsglobal>=3600 | tsglobal==tsglobal(1);
% --- select all datapoints
I           = tsglobal>=0;
z_xb        = z_xb(:,I);
z_data      = z_data(:,I);
tsglobal    = tsglobal(I);

N_divide = 0;

%% 1 score volume
V_xb            = NaN(1,length(tsglobal)-1);
V_data          = NaN(1,length(tsglobal)-1);
relBiasVolume   = NaN(1,length(tsglobal)-1);
relBiasVolume2  = NaN(1,length(tsglobal)-1);
if volume
    N_divide = N_divide+1;
    V_data_last = compute_Vero(x_data, z_data(:,1),z_data(:,end), zs0  );
    for ii = 1:length(tsglobal)-1% skip first output time
        % --- xb
        V_xb(ii)        = compute_Vero(x_xb, z_xb(:,1),z_xb(:,ii+1), zs0  );
        % --- data
        V_data(ii)      = compute_Vero(x_data, z_data(:,1),z_data(:,ii+1), zs0  );
    end
    % -- rel bias
    relBiasVolume       = (V_xb-V_data)./V_data_last;
    relBiasVolume2      = (V_xb-V_data)./V_data;
end
%% 2 score berm slope
slope_xb            = NaN(1,length(tsglobal)-1);
slope_data          = NaN(1,length(tsglobal)-1);
relSlope            = NaN(1,length(tsglobal)-1);
if bermslope
    N_divide = N_divide+1;
    for ii = 1:length(tsglobal)-1% skip first output time
        % --- xb
        slope_xb(ii)    = compute_bedslope_berm_based_on_accretion_zone(x_xb,z_xb(:,1),z_xb(:,ii+1),zs0);
        % --- data
        zb_data_int0    = interp1( x_data, z_data(:,1),x_xb);
        zb_data_int     = interp1( x_data, z_data(:,ii+1), x_xb);
        slope_data(ii)  = compute_bedslope_berm_based_on_accretion_zone(x_xb,zb_data_int0,zb_data_int,zs0);

    end
    relSlope     = (slope_xb-slope_data)./(slope_data);
end
%% 3 migration
dx_xb                   = NaN(1,length(tsglobal)-1);
dx_data                 = NaN(1,length(tsglobal)-1);
relMigration            = NaN(1,length(tsglobal)-1);
relMigration2           = NaN(1,length(tsglobal)-1);
if retreat
    N_divide = N_divide+1;
    dxmin = min(diff(unique(x_xb)));

    zb_data_int0    = interp1( x_data, z_data(:,1),x_xb);
    zb_data_int     = interp1( x_data, z_data(:,end), x_xb);
    dx_data_last = migration(x_xb, zb_data_int,zb_data_int0,dxmin);
    for ii = 1:length(tsglobal)-1% skip first output time
        % --- xb
        dx_xb(ii) = migration(x_xb,z_xb(:,ii+1),z_xb(:,1),dxmin);
        % --- data
        zb_data_int0    = interp1( x_data, z_data(:,1),x_xb);
        zb_data_int     = interp1( x_data, z_data(:,ii+1), x_xb);
        dx_data(ii) = migration(x_xb, zb_data_int,zb_data_int0,dxmin);
    end
    relMigration     = (dx_xb-dx_data)/(dx_data_last);
    relMigration2    = (dx_xb-dx_data)./(dx_data);
end

% --- skip values (e.g. profile is not complete)
relBiasVolume(skip) = NaN;
relSlope(skip)      = NaN;
relMigration(skip)  = NaN;

% --- compute error only over datapoints outside first hour
I               = tsglobal(2:end)>=3600;
RMS_vol         = sqrt( mean(relBiasVolume(I).^2,'omitnan') );
RMS_slope       = sqrt( mean(relSlope(I).^2,'omitnan') );    
RMS_migratie    = sqrt( mean(relMigration(I).^2,'omitnan') ); 

skill = sum([RMS_vol RMS_slope RMS_migratie],'omitnan')/N_divide * sign(mean(relBiasVolume(I),'omitnan'));


%% plot

colors_experiments = jet(size(tsglobal_orig,1));

fig1 = figure; set(gcf, 'Position',  [50, 50, 750, 400])
subplot('position',[0.15 0.15 0.6 0.7]);
for it = 1:length(tsglobal_orig)
    %% plot
    h1(it) = plot(x_xb,z_xb_orig(:,it),'k-','LineWidth',2,'color',colors_experiments(it,:),'displayName',sprintf('xbeach t = %i min',tsglobal_orig(it)/60)); hold on
    hold on; h2(it) = plot(x_data, z_data_orig(:,it),'k:','LineWidth',2,'color',colors_experiments(it,:),'displayName',sprintf('meas t = %i min',tsglobal_orig(it)/60)); hold on
    legendInfo{it} = ['Time =' num2str(round(tsglobal_orig(it)/3600,2)) ' Hour'];
    

end
plot([x_data(1) x_data(end)],[zs0 zs0],'k--','LineWidth',2);

xlabel('x [m w.r.t. wavemaker]','FontSize',15); ylabel('z [m]','FontSize',15);
title(name)
legend(h1,legendInfo,'Location',[0.8 0.15 0.15 0.7]);
grid on; box on; 
set(gca,'FontSize',12,'xlim',xlim,'ylim',ylim);



%text(0.5, 0.1,['RMS(rel. error_{V}) = ' num2str(round(RMS_vol,2))] ,'Units','normalized','FontSize',15);
%text(0.5, 0.2,['RMS(rel. error_{Slope}) = ' num2str(round(RMS_slope,2))] ,'Units','normalized','FontSize',15);
%text(0.5, 0.3,['RMS(rel. error_{dx}) = ' num2str(round(RMS_migratie,2))] ,'Units','normalized','FontSize',15);

%% plot 2
fig2 = figure; set(gcf, 'Position',  [50, 50, 750, 400])
subplot('position',[0.15 0.15 0.6 0.7]);
if volume
plot(round(tsglobal(1:end-1)/3600,2), relBiasVolume*100,'b-o','LineWidth',2,'displayName','V_{rel}');hold on
end
if bermslope
plot(round(tsglobal(1:end-1)/3600,2), relSlope*100,'r-o','LineWidth',2,'displayName','S_{rel}');hold on
end
if retreat
plot(round(tsglobal(1:end-1)/3600,2), relMigration*100,'g-o','LineWidth',2,'displayName','dx_{rel}');hold on
end
xlabel('Time [hour]'); ylabel('Rel. error [%]'); grid on
legend('Location',[0.8 0.15 0.15 0.7]);



save('skill2.mat','skill','RMS_vol','RMS_slope','RMS_migratie','relBiasVolume','relSlope','relMigration','tsglobal','V_xb','V_data','slope_xb','slope_data','dx_data','dx_xb')

% --- get times
for mm=1:length(tsglobal)
    rowlabel{mm} = num2str(round(tsglobal(mm)/3600,2));
end
rowlabel{1} = 'Time \newline[hour]'

% --- get data
data    = [];
format = {'%d'};
collabel = {};
if volume
    data    = [data V_xb', V_data' relBiasVolume'*100 relBiasVolume2'*100];
    format  = [format '%4.2f', '%4.2f','%4.0f','%4.0f'];
    collabel= [collabel '$V_{xb}$\newline[$m^3/m$]','$V_{data}$\newline[$m^3/m$]','V_{rel}\newline [\%]','V_{rel,t}\newline [\%]'];
end
if bermslope
    data    = [data slope_xb' slope_data' relSlope'*100];
    format  = [format '%4.2f', '%4.2f','%4.0f'];
    collabel= [collabel '$S_{xb}\newline[-]$','$S_{data}$\newline[-]','S_{rel}\newline [\%]'];
end
if retreat
    data    = [data dx_xb' dx_data' relMigration'*100 relMigration2'*100];
    format  = [format '%4.2f', '%4.2f','%4.0f','%4.0f'];
    collabel= [collabel '$dx_{xb}\newline[m]$','$dx_{data}\newline[m]$','dx_{rel}\newline [\%]','dx_{rel,t}\newline [\%]'];
end  

if volume & bermslope
    caption = sprintf('The computed and observed volume ($V$) and berm slopes ($S$) for %s. Volumes (V) are given in $m^3/m$. The relative error is expressed in a percentage. Both the relative error in terms of the same time ($rel,t$) and the final time are shown ($rel$). ',name);
elseif volume & bermslope & retreat
    caption = sprintf('The computed and observed volume ($V$), berm slopes ($S$) and retreat ($dx$) for %s. Volumes (V) are given in $m^3/m$ and retreat in $m$. The relative error is expressed in a percentage. Both the relative error in terms of the same time ($rel,t$) and the final time are shown ($rel$). ',name);
else
    'Warning';
end
matrix2latex(data, 'filename', 'indicator.tex', ...
    'caption',caption  , ...
    'rowlabel', rowlabel, ...
	'rowjustification','p{0.75cm}',...
    'collabel', collabel, ...
    'format',	format);

end




function dVero = compute_Vero(x,z0,zend,zs0) %based on rijkswaterstaat conventions
    dz = max(zend, zs0)-z0;
    dA = [diff(x);0].*dz;
    dVero = sum(dA(dA<0)); %m3/m
    dVero = dVero*-1;
end

function dx = migration(x,z,z0,dxmin)
    elevation_to_check_propagation_front = max(z0)-dxmin*1.5;
    
    iF              = find(z0>elevation_to_check_propagation_front,1,'first');
    %x_inital_XBeach = interp1(z0(iF-1:iF,1,1),x(iF-1:iF,1,1),elevation_to_check_propagation_front);
    x_inital_XBeach = x(iF-1);
    % X at tend
    elevation_to_check_propagation_front = max(z)-dxmin*1.5;
    iF = find(z(:,1,end)>elevation_to_check_propagation_front,1,'first');
    %x_final_XBeach = interp1(z(iF-1:iF,1,end),x(iF-1:iF,1,end),elevation_to_check_propagation_front);
    
    x_final_XBeach =x(iF-1);
    
    dx = x_final_XBeach-x_inital_XBeach;
end


function [s] = compute_bedslope_berm_based_on_accretion_zone(x,zt0,zte,zs0)
    
    plot_result = 0;    % Only use this [1] while finetuning this function (as it will generate a figure and pause at the end of this function)
      
    % Step 1: determine x_coordinate SSL
    if zte(end)<zte(1)
        s = 0;
        return
        %error('Assuming the dunes to be at the right side of the domain');
    end
    last_point_below_zs0 = max(find(zte<zs0));
    x_SSL = interp1(zte(last_point_below_zs0:(last_point_below_zs0+1)),x(last_point_below_zs0:(last_point_below_zs0+1)),zs0);
    
    % Step 2: determine which cells to consider
    % Top point: SSL point
    z_to_consider = zte;
    z_to_consider(min(find(zte>zs0)):end) = NaN;    
    % Bottom point, end of the accretion zone
    % Here defined as the lowest bed level where still 50% of the maximum accretion occured (excluding first point)
    max_accretion = max(zte-zt0);
    I = find(zte-zt0>0.5*max_accretion,1,'first');
    z_to_consider(1:(I-1)) = NaN;
    
    % Step 3: determine the average slope of the remaining cells   
    % First interpolate to a uniform (tight) grid to ensure fair averaging
    dx = min(diff(x))/10;
    new_x = x(min(find(~isnan(z_to_consider))))+dx:dx:x(max(find(~isnan(z_to_consider))))-dx;
    z_new = interp1(x,z_to_consider,new_x);
    slopes = (z_new(2:end)-z_new(1:end-1))./(new_x(2:end)-new_x(1:end-1));
    s = mean(slopes);
    
    % ADDITIONAL: bed slope on initial profile from the bottom of the
    % accretion zone to the original SSL intersect --> s0
    z_to_consider_for_s0 = zt0;
    z_to_consider_for_s0(1:min(find(~isnan(z_to_consider)))-1) = NaN; % Remove cells left of deposition zone
    z_to_consider_for_s0(min(find(zt0>zs0)):end) = NaN;
    new_x_for_s0 = x(min(find(~isnan(z_to_consider_for_s0)))):min(diff(x))/10:x(max(find(~isnan(z_to_consider_for_s0))));
    z_new_for_s0 = interp1(x,z_to_consider_for_s0,new_x_for_s0);
    slopes_for_s0 = (z_new_for_s0(2:end)-z_new_for_s0(1:end-1))./(new_x_for_s0(2:end)-new_x_for_s0(1:end-1));
    s0 = mean(slopes_for_s0);

    if plot_result == 1 
        fig = figure;
        plot(x,zt0,'k.-');
        hold on
        plot(x,zte,'.-');
        hold on
        plot(x,z_to_consider,'.-r')
        plot(x,x.*0+zs0)
        plot(x_SSL,zs0,'ok')
        hold on
        plot([x_SSL x(max(find(~isnan(z_to_consider))))]  ,[zs0 z_to_consider(max(find(~isnan(z_to_consider))))],'k')
        title('Slope at t')
        
        
        fig = figure;
        plot(x,zt0,'k.-');
        hold on
        plot(x,zte,'.-');
        hold on
        plot(x,z_to_consider_for_s0,'.-r')
        plot(x,x.*0+zs0)
        plot(x_SSL,zs0,'ok')
        hold on
        plot([x_SSL x(max(find(~isnan(z_to_consider_for_s0))))]  ,[zs0 z_to_consider_for_s0(max(find(~isnan(z_to_consider_for_s0))))],'k')
        title('Slope at t0')
        
        pause
    end
end