Merge reservoir 2 into reservoir 1
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(NeutralParticleReservoir), | intent(inout) | :: | res1 | Reservoir 1 |
||
type(NeutralParticleReservoir), | intent(inout) | :: | res2 | Reservoir 2 |
subroutine merge_reservoirs(res1,res2)
!+ Merge reservoir 2 into reservoir 1
type(NeutralParticleReservoir), intent(inout) :: res1
!+ Reservoir 1
type(NeutralParticleReservoir), intent(inout) :: res2
!+ Reservoir 2
integer :: i
type(NeutralParticleReservoir) :: res3
real(Float64) :: p, randomu(1)
if(res1%n.le.0) then
call init_reservoir(res1) !res2%k)
endif
if(res2%n.le.0) return
!! If neither reservoir is full just fill res1 from res2 normally
if ((res1%n.le.res1%k).and.(res2%n.le.res2%k)) then
do i=1,res2%n
call update_reservoir(res1, res2%R(i)%v, res2%R(i)%w)
enddo
res2 = res1
return
endif
!! Switch res1 with res2 if res1 is not full but res2 is
if (res1%n.le.res1%k) then
res3 = res1
res1 = res2
res2 = res3
endif
!! if res2 is not full then fill res1 from res2 normally
if (res2%n.le.res2%k) then
do i=1,res2%n
call update_reservoir(res1, res2%R(i)%v, res2%R(i)%w)
enddo
else !! if both are full select particles from res2 according to # particles seen
p = float(res1%n)/(res1%n + res2%n)
do i=1,res1%k
call randu(randomu)
if(randomu(1).gt.p) then
res1%R(i) = res2%R(i)
endif
enddo
endif
!! set res1 to res2
res2 = res1
end subroutine merge_reservoirs