OPT = struct(...
    'id', 0,...
    'irun', [0 0 0],...
    'maindir', fullfile(fileparts(which(mfilename)), '..', 'simulations'),...
    'models', {{'durosplus' 'durosta' 'xbeach'}},...
    'markers', {{'o' 'd' 'x'}},...
    'variables', {{'WL' 'Hsig' 'Tp' 'D50'}},...
    'xlabels',  {{'Water level above NAP [m]' 'Significant wave height [m]' 'Peak wave period [s]' 'Grain size D_{50} [m]'}},...
    'erosion_measure', 'EV',...
    'ylabel', 'Erosion volume [m^3/m]',...
    'sqlfilename', 'results.sqlite',...
    'xrelative', false,...
    'yrelative', true);

if isscalar(OPT.irun)
    OPT.irun = repmat(OPT.irun, size(OPT.models));
end

% Open database
database = fullfile(OPT.maindir, '..', 'sql', OPT.sqlfilename);
table = 'data';
mksqlite('open', database);
for imodel = 1:length(OPT.models)
    model = OPT.models{imodel};
    irun = OPT.irun(imodel);
    temp = mksqlite(['select ' sprintf('%s, ', OPT.variables{:}) OPT.erosion_measure ' from ' table ' where model = "' model '" and run = ' num2str(irun)]);
    data.(model) = cell2mat(struct2cell(temp));
end
mksqlite('close', database);


for ivar = 1:length(OPT.variables)
    variable = OPT.variables{ivar};
    for imodel = 1:length(OPT.models)
        model = OPT.models{imodel};
        refcase = median(data.(model), 2);
        istart = find(data.(model)(ivar,:) ~= refcase(ivar), 1, 'first');
        iend   = find(data.(model)(ivar,:) ~= refcase(ivar), 1, 'last');
        xdata.(variable)(:,imodel) = data.(model)(ivar,istart:iend)';
        ydata.(variable)(:,imodel) = data.(model)(end,istart:iend)';
        if OPT.yrelative
            iref = ismember(xdata.(variable)(:,imodel), refcase(ivar));
            ydata.(variable)(:,imodel) = ydata.(variable)(:,imodel) / ydata.(variable)(iref,imodel);
        end
        if OPT.xrelative
            iref = ismember(xdata.(variable)(:,imodel), refcase(ivar));
            xdata.(variable)(:,imodel) = xdata.(variable)(:,imodel) / xdata.(variable)(iref,imodel);
        end
    end
    subplot(2,2,ivar)
    ph = plot(xdata.(variable), ydata.(variable), ':');
    xlabel(OPT.xlabels{ivar})
    ylabel(OPT.ylabel)
    set(ph,...
        {'displayname'}, OPT.models',...
        {'marker'}, OPT.markers')
    if isempty(findobj(gcf, 'tag', 'legend'))
        leg = legend('show');
        set(leg,...
            'location', 'northwest',...
            'box', 'off')
    end
    title(variable)
end

%% adjust axis limits
ah = findobj(gcf,...
    'type', 'axes',...
    '-and', '-not',...
    'tag', 'legend');

ph = findall(gcf,...
    'type', 'line');
ydata = get(ph, 'ydata');
ylen = cellfun(@length, ydata);
ydata(ylen==5) = cellfun(@(y) y/y(3), ydata(ylen==5),...
    'uniformoutput', false);
set(ph, {'ydata'}, ydata)

yall = cell2mat(struct2cell(ydata));
yyall = cell2mat(ydata(ylen==5));
ylim = minmax(yall(:)');
lim = ylim + diff(ylim) * .05 * [-1 1];
set(ah, 'ylim', ylim)
set(ah, {'xlim'}, get(ah, 'xlim'))

%%
figureStyle(@styleThesis, {'figure', 'Position', [1 1 15 12]})
fpath = fileparts(which(mfilename));
fname = 'des_sensitivity_overview';
figurepath = fullfile(fpath, '..', 'figures');
set(gcf, 'tag', '$HeadURL$ $Revision$')
print(gcf, '-depsc2', fullfile(figurepath, [fname '.eps']))
if ispc
    saveas(gcf, fullfile(figurepath, [fname '.emf']))
end
saveas(gcf, fullfile(figurepath, [fname '.fig']))