SUBROUTINE sample_distribution(params,eta,etao)
TYPE(KORC_PARAMS), INTENT(IN) :: params
REAL(rp), DIMENSION(:), ALLOCATABLE, INTENT(INOUT) :: eta
REAL(rp), INTENT(OUT) :: etao
REAL(rp) :: go_root
REAL(rp) :: etao_root
REAL(rp) :: eta_buffer, eta_test
REAL(rp) :: ratio, rand_unif
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_samples
REAL(rp), DIMENSION(:), ALLOCATABLE :: eta_tmp
REAL(rp) :: minmax, min_pitch_angle, max_pitch_angle
REAL(rp) :: deta
LOGICAL :: leta
INTEGER :: num_accepted
INTEGER :: ii,jj,ppp,nsamples
INTEGER :: mpierr
ppp = SIZE(eta)
nsamples = ppp*params%mpi_params%nmpi
deta = (pdf_params%max_pitch_angle - pdf_params%min_pitch_angle)/100.0_rp
if (pdf_params%min_pitch_angle.GE.korc_zero) then
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = pdf_params%min_pitch_angle - REAL(jj,rp)*deta
if (minmax.GT.0.0_rp) then
min_pitch_angle = minmax
end if
end do
else
min_pitch_angle = pdf_params%min_pitch_angle
end if
do jj=1_idef,INT(minmax_buffer_size,idef)
minmax = pdf_params%max_pitch_angle + REAL(jj,rp)*deta
if (minmax.LE.90.0_rp) then
max_pitch_angle = minmax
end if
end do
! write(6,*) min_pitch_angle,max_pitch_angle
if (params%mpi_params%rank.EQ.0_idef) then
ALLOCATE(eta_samples(nsamples))! Number of samples to distribute among all MPI processes
ALLOCATE(eta_tmp(nsamples))! Number of samples to distribute among all MPI processes
!* * * Transient * * *!
call RANDOM_SEED()
call RANDOM_NUMBER(rand_unif)
eta_buffer = pdf_params%min_pitch_angle + (pdf_params%max_pitch_angle - pdf_params%min_pitch_angle)*rand_unif
ii=2_idef
do while (ii .LE. 1000_idef)
eta_test = eta_buffer + random_norm(0.0_rp,deta)
do while ((ABS(eta_test) .GT. pdf_params%max_pitch_angle).OR.(ABS(eta_test) .LT. pdf_params%min_pitch_angle))
eta_test = eta_buffer + random_norm(0.0_rp,deta)
end do
ratio = fRE(eta_test,pdf_params%po)/fRE(eta_buffer,pdf_params%po)
if (ratio .GE. 1.0_rp) then
eta_buffer = eta_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
eta_buffer = eta_test
ii = ii + 1_idef
end if
end if
end do
!* * * Transient * * *!
eta_tmp(1) = eta_buffer
call RANDOM_SEED()
call RANDOM_NUMBER(rand_unif)
num_accepted = 0_idef
do while(num_accepted.LT.nsamples)
ii=2_idef
do while (ii .LE. nsamples)
eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,deta)
do while ((ABS(eta_test) .GT. max_pitch_angle).OR.(ABS(eta_test) .LT. min_pitch_angle))
eta_test = eta_tmp(ii-1) + random_norm(0.0_rp,deta)
end do
ratio = fRE(eta_test,pdf_params%po)/fRE(eta_tmp(ii-1),pdf_params%po)
if (ratio .GE. 1.0_rp) then
eta_tmp(ii) = eta_test
ii = ii + 1_idef
else
call RANDOM_NUMBER(rand_unif)
if (rand_unif .LT. ratio) then
eta_tmp(ii) = eta_test
ii = ii + 1_idef
end if
end if
end do
eta_tmp = ABS(eta_tmp)
ii = 1_idef
do while ( (ii.LT.nsamples).AND.(num_accepted.LT.nsamples) )
leta = (eta_tmp(ii).LE.pdf_params%max_pitch_angle).AND.(eta_tmp(ii).GE.pdf_params%min_pitch_angle)
if (leta) then
num_accepted = num_accepted + 1_idef
eta_samples(num_accepted) = eta_tmp(ii)
end if
ii = ii + 1_idef
end do
end do
etao = SUM(eta_samples)/nsamples
end if
CALL MPI_SCATTER(eta_samples,ppp,MPI_REAL8,eta,ppp,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
CALL MPI_BCAST(etao,1,MPI_REAL8,0,MPI_COMM_WORLD,mpierr)
call MPI_BARRIER(MPI_COMM_WORLD,mpierr)
if (params%mpi_params%rank.EQ.0_idef) then
DEALLOCATE(eta_samples)
DEALLOCATE(eta_tmp)
end if
END SUBROUTINE sample_distribution