function [integral, flow_extrap] = probability_integral(flow,timespan_flow,pdf,timespan_pdf,timespan_out)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Integral of a flow back in time weighted with a given pdf (POM -> WG
%%% through residence time, etc)
%%%
%%% Dmitry Yumashev, 13/12/2019 
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% make sure the time dimensions are correct to allow for matrix
%%% operations
if size(timespan_out,1) > 1 && size(timespan_out,2) == 1 % if the vector variable is a column
    timespan_out = transpose(timespan_out); % make it a row
end
%
if size(pdf,1) == 1 && size(pdf,2) > 1 % if the vector variable is a row
    pdf = transpose(pdf); % make it a column
end


%%% general settings

total_timesteps_out = numel(timespan_out);
integral = zeros(1,total_timesteps_out); % make sure it's a row

total_timesteps_pdf = numel(timespan_pdf);

%%% extrapolate the input flow according to the specified timespan_out

flow_extrap = interp1(timespan_flow,flow,timespan_out,'linear','extrap');

flow_extrap(flow_extrap<0) = 0; % failsafe

%%% work out the integral

for i=1:total_timesteps_out % current time index
        
    % historic time loop
    j_start = max([(i - total_timesteps_pdf + 1) 1]); % make sure we don't go below 1
    j_end = i; 
    
    sss = 0;
    
    for j = j_end:-1:j_start % going back in time

        sss = sss + flow_extrap(j) * pdf(i-j+1); % working out the integral
        
    end

    integral(i) = sss;
    
end