function [fig1 fig2] = xb_morpho_analyse(globalx, zb,x_data, z_data, tsglobal,zs0)

% --- erosion/sedimentation
fig1 = figure
subplot('position',[0.1 0.6 0.6 0.3]);
plot(x_data,z_data(:,end)-z_data(:,1),'ko' )
hold on
plot(globalx, zb(:,end)-zb(:,1),'r-','linewidth',2)
grid on
xlabel('Distance [m]');ylabel('Height [m]'); title('bed level change')
legend({'Observed','Computed'},'Location',[.75 .6 .2 .3])

subplot('position',[0.1 0.1 0.6 0.3]);
plot(x_data,z_data(:,1),'k--' ,'linewidth',2)
hold on
plot(x_data,z_data(:,end),'k' ,'linewidth',2)
plot(globalx, zb(:,end),'r-','linewidth',2)
legend('Initial profile','Observed final profile','Computed final profile','Location',[.75 .1 .2 .3])
grid on
xlabel('Distance [m]'); ylabel('Height [m]'); title('Profiles')

% --- erosion in time
for ii=1:length(tsglobal)
    V_xb(ii)    = compute_Vero(globalx,zb(:,1),zb(:,ii),zs0);
    V_obs(ii)   = compute_Vero(x_data,z_data(:,1),z_data(:,ii),zs0);
end

fig2 = figure;
subplot('position',[0.1 0.1 0.6 0.8]);
plot(tsglobal/3600,V_obs,'ko')
hold on
plot(tsglobal/3600,V_xb,'r-','linewidth',2)
title('Erosion volumes')
grid on
xlabel('Time [hr]'); ylabel('Volume [m^3/m]')
legend('Observed','Computed','Location',[.75 .1 .2 .3])
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