% SBEACH toolbox - fsbeach_processresults
% Written by Joshua Simmons 11/2017

%   Syntax:
%   sbout = fsbeach_processresults(sblocs)
%
%   Input:
%   sbout    =  structure with input, computed and measured

%   Output:
%   sbout     = structure with sbeach data and statistics;

function sbout = fsbeach_processresults(sbout)
    for ii = 1:length(sbout.input)
        xin = sbout.input(ii).x;
        zin = sbout.input(ii).z;
        xcomp = sbout.computed(ii).x;
        zcomp = sbout.computed(ii).z;
        xmeas = sbout.measured(ii).x;
        zmeas = sbout.measured(ii).z;
        
        try
            [tmpbss, tmprmse, dVc, dVm, ddV] = fsbeach_modelskill([xmeas zmeas], [xcomp zcomp], [xin zin]);

            %calculate BSS above MSL
            %first get points above MSL - for measured and then apply that
            %to all - cant take these independently as may only get say one
            %point for modelled and then xb_skill interpolates.
            %BSS deals well if based on measured as the -10 and x
            %determinant - it only finds the interpolated values with
            %info so will deal if initpf doesnt stretch to all of
            %measured pf
            keepme1 = find(zmeas>0);
            if xmeas(keepme1(1))>xmeas(keepme1(end))
                keepmo1 = find(xcomp>xmeas(keepme1(end)) & xcomp<xmeas(keepme1(1)));
                keepin1 = find(xin>xmeas(keepme1(end)) & xin<xmeas(keepme1(1)));
            elseif xmeas(keepme1(1))<xmeas(keepme1(end))
                keepmo1 = find(xcomp<xmeas(keepme1(end)) & xcomp>xmeas(keepme1(1)));
                keepin1 = find(xin<xmeas(keepme1(end)) & xin>xmeas(keepme1(1)));
            end
            [tmpabMSLbss, tmpabMSLrmse, ~, ~, ~] = fsbeach_modelskill([xmeas(keepme1) zmeas(keepme1)], [xcomp(keepmo1) zcomp(keepmo1)], [xin(keepin1) zin(keepin1)]);

            %calculate BSS below MSL
            keepme2 = find(zmeas<0);
            if xmeas(keepme2(1))>xmeas(keepme2(end))
                keepmo2 = find(xcomp>xmeas(keepme2(end)) & xcomp<xmeas(keepme2(1)));
                keepin2 = find(xin>xmeas(keepme2(end)) & xin<xmeas(keepme2(1)));
            elseif xmeas(keepme2(1))<xmeas(keepme2(end))
                keepmo2 = find(xcomp<xmeas(keepme2(end)) & xcomp>xmeas(keepme2(1)));
                keepin2 = find(xin<xmeas(keepme2(end)) & xin>xmeas(keepme2(1)));
            end
            [tmpbeMSLbss,tmpbeMSLrmse, ~, ~, ~] = fsbeach_modelskill([xmeas(keepme2) zmeas(keepme2)], [xcomp(keepmo2) zcomp(keepmo2)], [xin(keepin2) zin(keepin2)]);

            %store results
            sbout.stats(ii).bss = tmpbss;
            sbout.stats(ii).rmse = tmprmse;
            sbout.stats(ii).dVc = dVc;
            sbout.stats(ii).dVm = dVm;
            sbout.stats(ii).ddV = ddV;

            sbout.stats(ii).aboveMSL.bss = tmpabMSLbss;
            sbout.stats(ii).aboveMSL.rmse = tmpabMSLrmse;

            sbout.stats(ii).belowMSL.bss = tmpbeMSLbss;
            sbout.stats(ii).belowMSL.rmse = tmpbeMSLrmse;
        catch
            disp(['Run number ' num2str(ii) ' failed.'])           
            continue
        end
    end
end