subroutine geometry_definition(lattice,prob_distr,M,N,L) 
    implicit none
    integer(8) N
    integer s, p, M
    real, parameter :: Pi = 3.141593
    real L
    real(8) phi, Q
    real, dimension(4,N) :: lattice 
    real, dimension(M,2) :: prob_distr 
    call random_seed()
    call random_number(lattice(1:3,:))          
    lattice(1:3,:)=L*lattice(1:3,:)                                                         ! "Разбрасывание" координат пор                                              
    do s = 1,N                                                                              ! "Разбрасывание" радиусов пор
        call random_number(Q)         
        if (Q < prob_distr(1,2)) then
            lattice(4,s) = prob_distr(1,1)
        else 
                do p = 2, M
                if ((Q < sum(prob_distr(1:p,2))) .and. (Q > sum(prob_distr(1:(p-1),2)))) then
                lattice(4,s) = prob_distr(p,1)
                end if
            end do
        end if
    end do
    
    end subroutine geometry_definition
    
    subroutine findclusters(N,R_max,lattice)
    implicit none
    integer(8) N
    integer i, p, j, l, q
    integer :: k = 1
    real R_max, R
    real, dimension(5,N)::lattice(5,N), locations(N,3), distances(N), indx(N), radius_vectors(N,3)
    integer, dimension(:),allocatable::indx1
    
    indx = (/(q,q=1,N)/) 
    lattice(5,:) = 0
    do l = 1, N  
        do p = 1, 3
            locations(l,p) = lattice(p,l)                                                                  ! Задание матрицы координат пор              
        end do                                                                            
    end do
    do i=1, N-1                                                                                                                            
            do p = 1, 3                                                                                    
                radius_vectors(:,p) = locations(:,p) - locations(i,p)                                      ! Определение матрицы радиус-векторов для текущей поры
            end do
        distances = sqrt(sum((radius_vectors)**2, dim = 1))                               ! Определение матрицы расстояний от текущей поры до всех остальных (включая "текущую")                                                              ! Массив индексов всех пор
        allocate(indx1(count(distances < 2*R_max) - 1))                                   ! Выделение памяти для массива индексов потенциальных соседей текущей поры indx1
        indx1 = pack(indx, distances < 2*R_max .and. distances /= 0)
        do p = 1,size(indx1,1)
                j = indx1(p)
                R = sqrt(sum((lattice(1:3,i) - lattice(1:3,j))**2))
            if (R < (lattice(4,i) + lattice(4,j))) then
                if ((lattice(5,i) == 0) .and. (lattice(5,j) == 0)) then                   ! Если поры соединены и не принадлежат ни одному из кластеров, они образуют новый кластер 
                    lattice(5,i) = k                                                  
                    lattice(5,j) = k
                    k = k + 1
                else if ((lattice(5,i) /= 0) .and. (lattice(5,j) == 0)) then              ! Если одна из пор уже принадлежит какому-либо кластеру, а другая - нет, то
                    lattice(5,j) = lattice(5,i)                                           ! пора без кластера присоединяется к кластеру другой
                else if ((lattice(5,i) == 0) .and. (lattice(5,j) /= 0)) then
                    lattice(5,i) = lattice(5,j)
                else if ((lattice(5,j) /= 0) .and. (lattice(5,i) /= 0) .and. (lattice(5,j) /= lattice(5,i))) then                       ! Если же обе поры являются частью кластеров (разных),
                     where(lattice(5,:) == max(lattice(5,i),lattice(5,j))) lattice(5,:) = min(lattice(5,i),lattice(5,j))                ! то происходит объединение кластеров в тот, что имел меньший индекс                                                                                                            ! то происходит объединение кластеров в тот, что имел меньший индекс
                end if  
             end if
        end do
        deallocate(indx1)
    end do
    
    end subroutine findclusters
    
    
    program experiment 
    implicit none
    integer(8) N                                                                            ! Число пор
    integer k                                                                                                                                                 
    real, parameter :: Pi = 3.141593
    real L, phi, V_avg, R_max
    real, dimension(:,:), allocatable :: lattice 
    real, dimension(:,:), allocatable :: prob_distr                                         ! Распределение пор по радиусам
    print*, 'Enter the parameters of your system'
    print*, 'L ='                                                                           ! Размер системы    
    read*, L    
    print*, 'phi ='                                                                         ! Пористость системы
    read*, phi
    print*, 'How many points will be in your distribution?'
    read*, k
    allocate(prob_distr(k,2))
    print*, 'Enter the points of your distribution. First row is for radiuses, second is for probabilities'
    1 read*, prob_distr
        if (sum(prob_distr(:,2)) /= 1) then
            print*, 'The probabilities of all possible events must total to 1. Try again'
            goto 1
        end if
    R_max = maxval(prob_distr(:,1))                                                         ! Максимально возможный радиус поры
    V_avg = (sum(((prob_distr(:,1))**3)*prob_distr(:,2)))*4*Pi/3                            ! Средне значение объёма поры
    N = nint((L**3*log(1/(1-phi)))/V_avg)
    allocate(lattice(4,N))
    call geometry_definition(lattice,prob_distr,size(prob_distr(:,1)),N,L)
    call findclusters(N,R_max,lattice)
    pause
    end