-
Notifications
You must be signed in to change notification settings - Fork 0
/
rands.f90
61 lines (51 loc) · 1.43 KB
/
rands.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
module rands
implicit none
private
public :: set_seed,random_integer,random_gauss
contains
!==========================
! Random Integer
!==========================
subroutine random_integer(min,max,randint)
integer,intent(in) :: min, max
integer,intent(out) :: randint
real(kind=8) :: rand
call random_number(rand)
randint = min+floor((max-min+1)*rand)
end subroutine random_integer
!==========================
! Random Integer
!==========================
subroutine random_gauss(rand,mu,sigma)
real(kind=8),intent(in) :: mu,sigma
real(kind=8),intent(out) :: rand
real(kind=8),parameter :: twopi = 8.0d0*atan(1.0d0)
real(kind=8) :: t,phi,x1,x2,x3
call random_number(t)
call random_number(phi)
x1 = sqrt(-2.0d0*log(t))*cos(twopi*phi)
x1 = sqrt(-2.0d0*log(t))*sin(twopi*phi)
call random_number(x3)
if (x3.lt.0.5d0) then
rand = x1*sigma+mu
else
rand = x2*sigma+mu
endif
end subroutine random_gauss
!==========================
! Set Seed
!==========================
subroutine set_seed
integer,allocatable :: seed(:)
integer :: n_seed,clock,i
real(kind=8) :: x
call random_seed(size=n_seed)
allocate(seed(n_seed))
call system_clock(count=clock)
do i=1,n_seed
seed(i)=clock/i
enddo
call random_seed(put=seed)
deallocate(seed)
end subroutine set_seed
end module rands