function analyze_this_in_depth(info, dirs)

% read measurements
M = load(fullfile(dirs.data, [info.run 'data.mat']));
t = [0 .1 .3 1 2.04 6];

% read dimensions
d       = xb_read_dims;
[m i]   = closest(1000+t*3600, d.t);

% read input
xbi     = xb_read_input(pwd);
zs0     = max([0 xs_get(xbi,'zs0') max(xs_get(xbi,'zs0file.tide'))]);
fp      = min([1 1/xs_get(xbi,'Trep')]);

% read output
xbo1    = xb_read_output(pwd, 'vars', {'zb'}, 'index', i-1);
xbo3    = xb_read_output(pwd, 'vars', {'zs' 'ccg' 'DR' 'Tbore' 'kb' 'zb' 'dzav' 'Susg' 'Subg' 'u' 'ue' 'ua'});
xbo4    = xb_read_output(pwd, 'vars', {'zb_mean' 'zs_mean' 'u_mean' 'H_mean' 'urms_mean' 'ue_mean' 'ua_mean' 'Sk_mean' 'As_mean' 'hh_var'});

% analyze output
xbb     = xb_get_morpho(xbo3, 'level', zs0);
xbh     = xb_get_hydro(xbo4, 'Trep', 1/fp);
xbs     = xb_get_sedtrans(xbo3, 'Trep', 1/fp);

A = [];
for i = 1:size(M.z,1)
    A(i) = sum(max(0,M.z(1,:))-max(0,M.z(i,:)))*1;
end

beta = atan(M.Asu./M.Sku);
B    = sqrt(M.Asu.^2+M.Sku.^2);

% plot output
xb_plot_profile(xbo1, ...
    'measured', [M.x' M.z']);

xb_write_plot(gcf, dirs.output, 'profiles');

id = find(strcmp({xbb.data.name}, 'ero'),1); % find which row contains 'ero'
xbb.data(id).value = abs(xbb.data(id).value); %take abs of erosion - seems to solve problem with pictures, quick & dirty fix

xb_plot_morpho(xbb, ...
    'dz', [M.x' M.z(end,:)'-M.z(1,:)'], ...
    'ero', [t' A']);

xb_write_plot(gcf, dirs.output, 'morpho');

xb_plot_hydro(xbh, ...
    'Hrms_hf', [M.xhrms' M.Hrms_hf'], ...
    'Hrms_lf', [M.xhrms' M.Hrms_lf'], ...
    'Hrms_t', [M.xhrms' M.Hrms'], ...
    's', [M.xsetup' M.setup'], ...
    'umean', [M.xum' M.Um'], ...
    'urms_lf', [M.xurms' M.Urms_lf'], ...
    'urms_hf', [M.xurms' M.Urms_hf'], ...
    'urms_t', [M.xurms' M.Urms']);

subplot(2,1,1)
ylim([0 3])
h = legend('Location','Northwest')
h.NumColumns=2;
subplot(2,1,2)
h = legend('Location','Northwest')
h.NumColumns=2;
ylim([-0.5 2])

xb_write_plot(gcf, dirs.output, 'hydro');

f1 = figure; ax1 = axes;
f2 = figure; ax2 = subplot(2,1,1); ax3 = subplot(2,1,2);
f3 = figure; ax4 = axes;

xb_plot_hydro2(xbh, 'handles', [ax1 ax2 ax3 ax4], ...
    'urms_lf', [M.xurms' M.Urms_lf'], ...
    'urms_hf', [M.xurms' M.Urms_hf'], ...
    'urms_t', [M.xurms' M.Urms'], ...
    'uavg', [M.xhrms' .15*M.ua'], ...
    'SK', [M.xhrms' M.Sku'], ...
    'AS', [M.xhrms' M.Asu'], ...
    'rho', [M.xhrms' M.rho'], ...
    'B', [M.xhrms' B'], ...
    'beta', [M.xhrms' beta']);

set(ax2,'ylim',[-2,3])
set(ax3,'ylim',[-2,3])
set(ax4,'ylim',[0,1])
xb_write_plot(f1, dirs.output, 'corr');
xb_write_plot(f2, dirs.output, 'waveshape');
xb_write_plot(f3, dirs.output, 'undertow');

f4 = figure; ax5 = axes;
f5 = figure; ax6 = axes;
f6 = figure; ax7 = subplot(2,1,1); ax8 = subplot(2,1,2);

xb_plot_sedtrans(xbs, 'handles', [ax5 ax6 ax7 ax8], 'showall', true, ...
    'c', [M.xum' M.Cm'], ...
    'S_dz', [M.x' M.Sdzb]);

xb_write_plot(f4, dirs.output, 'concentration');
xb_write_plot(f5, dirs.output, 'turbulence');
xb_write_plot(f6, dirs.output, 'transports');

xb_write_skilltable( ...
    xs_join(xbo3, xbo4, xbb, xbh, xbs), ...
    {[M.x' M.z(end,:)'] [M.xhrms' M.Hrms_lf'] [M.xhrms' M.Hrms_hf'] [M.xhrms' M.rho'] ...
        [M.xhrms' M.Sku'] [M.xhrms' M.Asu'] [M.xhrms' B'] [M.xhrms' beta'] ...
        [M.xurms' M.Urms_lf'] [M.xurms' M.Urms_hf'] [M.xum' M.Um'] [M.xum' M.Cm'] ...
        [M.x' M.Sdzb]}, ...
    'vars', {'zb' 'Hrms_lf' 'Hrms_hf' 'rho' 'SK' 'AS' 'B' 'beta' 'urms_lf' 'urms_hf' ...
                'umean' 'c' 'dz'}, ...
    'file', fullfile(dirs.output, 'skills.tex'));


% remove all but the first and last outputs
for iv = 2:length(xbo1.data)
    if length(size(xbo1.data(iv).value)) == 3
        xbo1.data(iv).value = xbo1.data(iv).value([1 end],:,:);
    end
end

xb_plot_profile(xbo1,'measured',[M.x' M.z(end,:)'])
title(['profiles ' info.test ' ' info.run], 'interpreter', 'none')
xb_write_plot(gcf, dirs.output, 'final_profiles')