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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math.h>
# GAMMLN -- Return natural log of gamma function.
# POIDEV -- Returns Poisson deviates for a given mean.
# GASDEV -- Return a normally distributed deviate of zero mean and unit var.
# GAMMLN -- Return natural log of gamma function.
# Argument must greater than 0. Full accuracy is obtained for values
# greater than 1. For 0<xx<1, the reflection formula can be used first.
#
# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling.
# Used by permission of the authors.
# Copyright(c) 1986 Numerical Recipes Software.
real procedure gammln (xx)
real xx # Value to be evaluated
int j
double cof[6], stp, x, tmp, ser
data cof, stp / 76.18009173D0, -86.50532033D0, 24.01409822D0,
-1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
begin
x = xx - 1.0D0
tmp = x + 5.5D0
tmp = (x + 0.5D0) * log (tmp) - tmp
ser = 1.0D0
do j = 1, 6 {
x = x + 1.0D0
ser = ser + cof[j] / x
}
return (tmp + log (stp * ser))
end
# POIDEV -- Returns Poisson deviates for a given mean.
# The real value returned is an integer.
#
# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling.
# Used by permission of the authors.
# Copyright(c) 1986 Numerical Recipes Software.
#
# Modified to return zero for input values less than or equal to zero.
real procedure poidev (xm, seed)
real xm # Poisson mean
long seed # Random number seed
real oldm, g, em, t, y, ymin, ymax, sq, alxm, gammln(), urand(), gasdev()
data oldm /-1./
begin
if (xm <= 0)
em = 0
else if (xm < 12) {
if (xm != oldm) {
oldm = xm
g = exp (-xm)
}
em = 0
for (t = urand (seed); t > g; t = t * urand (seed))
em = em + 1
} else if (xm < 100) {
if (xm != oldm) {
oldm = xm
sq = sqrt (2. * xm)
ymin = -xm / sq
ymax = (1000 - xm) / sq
alxm = log (xm)
g = xm * alxm - gammln (xm+1.)
}
repeat {
repeat {
y = tan (PI * urand(seed))
} until (y >= ymin)
em = int (sq * min (y, ymax) + xm)
t = 0.9 * (1 + y**2) * exp (em * alxm - gammln (em+1) - g)
} until (urand(seed) <= t)
} else
em = xm + sqrt (xm) * gasdev (seed)
return (em)
end
# GASDEV -- Return a normally distributed deviate with zero mean and unit
# variance. The method computes two deviates simultaneously.
#
# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling.
# Used by permission of the authors.
# Copyright(c) 1986 Numerical Recipes Software.
real procedure gasdev (seed)
long seed # Seed for random numbers
real urand()
double v1, v2, r, fac
int iset
data iset/0/
begin
if (iset == 0) {
repeat {
v1 = 2 * urand (seed) - 1.
v2 = 2 * urand (seed) - 1.
r = v1 ** 2 + v2 ** 2
} until ((r > 0) && (r < 1))
fac = sqrt (-2. * log (r) / r)
iset = 1
return (v1 * fac)
} else {
iset = 0
return (v2 * fac)
}
end
|
