function [integral, flow_extrap] = probability_integral_advanced(flow,timespan_flow,...
    pdf,timespan_pdf,pdf_scale_factor,timesteps_pdf_scale_factor,timespan_out)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Integral of a flow back in time weighted with a given pdf (POM -> WG
%%% through residence time, etc), with the pdf itself varying with time:
%%% pdf = rho(t-t';t'); all the entries in the pdf matrix are assumed to
%%% have been foramted to the same timespan, while t' covers the entire 
%%% timespan_out
%%%
%%% Dmitry Yumashev, 13/12/2019 
%%% Dmitry Yumashev, 22/12/2019 
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% first, work out scaled pdfs in all years corresponding to timespan_out

[pdf_matrix] = time_distrib_historic_extrap(pdf,timespan_pdf,pdf_scale_factor,timesteps_pdf_scale_factor,timespan_out);

timespan_pdf_matrix = timespan_out;

%%% 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 sum(pdf_matrix(:,1)) < 1 && sum(pdf_matrix(1,:)) == 1 % if the distributions are placed along the 2nd dimension (sum = 1)
    pdf_matrix = transpose(pdf_matrix); % make sure the distributons are place along the 1st dimension
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_matrix);

%%% 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_matrix(i-j+1,j); % working out the integral; note the 2nd argument in the pdf matrix: each year j in the past could have a different pdf!!!
        
    end

    integral(i) = sss;
    
end