function [ probabilities ] = probabilities( dist_north,dist_south,dist_east,dist_west,radii,n )

%% This function calculates the probability for one clast to hit a certain area on the summit. 
%For this purpose, a simple radial grid with four quadrants (east, west, north and south) is defined and it is counted how many of
%a large number of simulated clasts hit certain fields in this radial grid.

%% Input parameters
%radii is a vector that gives the distance thresholds of radial areas (see Strehlow et al., 2017, JAV)
%n is the total number of ejected (simulated) clasts
%dist_north, dist_south, dist_east and dist_west are column vectors with stored flight distances of simulated clasts in the respective direction (each row represents one clast)

%note: obviously, only one of the dist-vectors will have an entry >0 in each row, because each clast has only one direction.

%% Storage vectors for successes (i.e. clasts that hit certain areas)
success_east=zeros(length(radii)-1,1); 
success_west=zeros(length(radii)-1,1);
success_north=zeros(length(radii)-1,1);
success_south=zeros(length(radii)-1,1);
%Note: the ith entry is the number of successes for the ith area, which in turn is defined by entries in the input vector radii. 
%Areas are defined by their lower bound (area i defined by radii(i))- therefore the success vectors have one less entry than radii, because
%the last entry in radii is just an upper bound (needed later in the probabilistic study) and not a definition of another area.

%% Count successes
for i=1:(length(radii)-1)
    for j=1:length(dist_east) %note: all dist vectors are of same length (n)
        if dist_east(j)>=radii(i) && dist_east(j)<radii(i+1) %condition for a clast to hit area i in the eastern quadrant
            success_east(i)=success_east(i)+1;
        elseif dist_west(j)>=radii(i) && dist_west(j)<radii(i+1) %condition for a clast to hit area i in the western quadrant
            success_west(i)=success_west(i)+1;
        elseif dist_north(j)>=radii(i) && dist_north(j)<radii(i+1) %condition for a clast to hit area i in the northern quadrant
            success_north(i)=success_north(i)+1;
        elseif dist_south(j)>=radii(i) && dist_south(j)<radii(i+1) %condition for a clast to hit area i in the southern quadrant
            success_south(i)=success_south(i)+1;
        end
    end  
end

%% Calculate Probabilities by dividing successes through total number of simulated clasts
probabilities_east=success_east./n;
probabilities_west=success_west./n;
probabilities_north=success_north./n;
probabilities_south=success_south./n;

%% Create output
probabilities=[probabilities_north probabilities_south probabilities_east probabilities_west];

end