1program answer
2 implicit none
3
4 integer :: i, nr, nbin
5 real(8), allocatable :: r(:), hist(:), binc(:)
6 real(8) :: xmin, xmax, lambda
7
8 lambda = 1.0_8
9 xmin = 0.0_8 / lambda
10 xmax = 12.0_8 / lambda
11
12 nr = 60000
13 nbin = 24
14 allocate(r(nr))
15 allocate(hist(nbin))
16 allocate(binc(nbin))
17
18 call random_seed_clock()
19 call rand_exp(r, lambda)
20 call histogram(r, xmin, xmax, nbin, binc, hist)
21
22 do i = 1, nbin
23 write(*, fmt='(e12.4, 1x, e12.4)') binc(i), hist(i)
24 enddo
25
26 deallocate(r)
27 deallocate(hist)
28 deallocate(binc)
29
30 stop
31contains
32
33 ! 指数乱数
34 subroutine rand_exp(r, lambda)
35 implicit none
36 real(8), intent(out) :: r(:)
37 real(8), intent(in) :: lambda
38
39 real(8) :: ur(size(r))
40
41 call random_number(ur)
42 r = -log(1 - ur) / lambda
43
44 endsubroutine rand_exp
45
46 ! 乱数のseedをシステムクロックに応じて変更
47 subroutine random_seed_clock()
48 implicit none
49 integer :: nseed, clock
50 integer, allocatable :: seed(:)
51
52 call system_clock(clock)
53
54 call random_seed(size=nseed)
55 allocate(seed(nseed))
56
57 seed = clock
58 call random_seed(put=seed)
59
60 deallocate(seed)
61 endsubroutine random_seed_clock
62
63 ! ヒストグラムを作成
64 subroutine histogram(x, xmin, xmax, nbin, binc, hist)
65 implicit none
66 real(8), intent(in) :: x(:)
67 real(8), intent(in) :: xmin, xmax
68 integer, intent(in) :: nbin
69 real(8), intent(out) :: binc(nbin)
70 real(8), intent(out) :: hist(nbin)
71
72 integer :: i, j
73 real(8) :: h, norm
74
75 h = (xmax - xmin) / nbin
76 norm = 1.0_8 / (size(x) * h)
77 hist = 0.0_8
78
79 do i = 1, size(x)
80 j = int(x(i) / h) + 1
81 if(j < 1 .or. j > nbin) cycle
82
83 hist(j) = hist(j) + 1.0_8 * norm
84 enddo
85
86 do j = 1, nbin
87 binc(j) = (j - 0.5_8) * h
88 enddo
89
90 endsubroutine histogram
91
92endprogram answer