function [pdf, cdf] = pdf_convolution(time_steps,par1,par2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Convolution integral of two Weibull pdfs with parameters par1 = (a1,b1) and
%%% par2 = (a2,b2), evaluated within the given time domain (time_steps)
%%%
%%% Outputs: convolution integral (pdf) and its cumulative (cdf)
%%%
%%% Dmitry Yumashev, 20/11/2019
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

total_time_steps = length(time_steps);

a1 = par1(1);
b1 = par1(2);
a2 = par2(1);
b2 = par2(2);

%%%

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

% first 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
   
        if i_min==1
        
            delta = time_steps(j+1) - time_steps(j);
        
        elseif i_min==2 
            
            delta = time_steps(j) - time_steps(j-1);
            
        end
        
        sum = sum + delta * wblpdf(time_steps(j),a1,b1) * wblpdf(time_steps(i-j+i_min),a2,b2);
        
    end
    
    pdf(i-i_min+1) = sum;
    
end
    
pdf_left = pdf;

% second approximation: right-hand sides of integration intervals
i_min = 2; % NOTE: one step shift to the right compared with previous case
i_max = total_time_steps; % NOTE: one step shift to the right compared with previous case

for i=i_min:i_max

    sum = 0;
    
    for j = i_min:i
   
        if i_min==1
        
            delta = time_steps(j+1) - time_steps(j);
        
        elseif i_min==2 
            
            delta = time_steps(j) - time_steps(j-1);
            
        end
        
        sum = sum + delta * wblpdf(time_steps(j),a1,b1) * wblpdf(time_steps(i-j+i_min),a2,b2);
        
    end
    
    pdf(i-i_min+1) = sum;
    
end

pdf_right = pdf;

% take average of the two

pdf = 0.5 * (pdf_left + pdf_right);


%%% 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