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