function [pdf, cdf] = pdf_convolution(time_steps,pdf1,pdf2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Convolution integral of two pdfs (pdf1, pdf2) evaluated within the 
%%% given time domain (time_steps, one step more than in the pdf's domain)
%%%
%%% Outputs: convolution integral (pdf) and its cumulative (cdf)
%%%
%%% Dmitry Yumashev, 20/11/2019
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

total_time_steps = length(time_steps);

%%%

pdf = zeros(total_time_steps-1,1);

% NOTE: introduce the integration step 'delta' in the convolution sum to 
% ensure the correct numerical approximation of the convolution integral

% only possible approximation: left-hand sides of integration intervals
i_min = 1; % NOTE the value
i_max = total_time_steps-1; % NOTE the value 

for i=i_min:i_max

    sum = 0;
    
    for j = i_min:i
   
        delta = time_steps(j+1) - time_steps(j);
                
        sum = sum + delta * pdf1(j) * pdf2(i-j+i_min);
        
    end
    
    pdf(i-i_min+1) = sum;
    
end
    



%%% now work out the cdf

cdf = zeros(total_time_steps,1);

cdf(1) = 0; % initial value (just in case)

for i=2:total_time_steps % NOTE the starting value
    
    delta = time_steps(i) - time_steps(i-1);
      
    cdf(i) = cdf(i-1) + pdf(i-1) * delta;
   
end


%%% end of the function

end