root / synthbench / euroben-shm / mod2as / ranmod.f
History | View | Annotate | Download (2.7 kB)
1 |
Module ran_module |
---|---|
2 |
! ---------------------------------------------------------------------- |
3 |
Use numerics |
4 |
|
5 |
Integer, Parameter :: m1 = 259200, a1 = 7141, c1 = 54773, |
6 |
& m2 = 134456, a2 = 8121, c2 = 28411 |
7 |
Real(l_), Parameter :: one = 1.0_l_, rm1 = one/m1, |
8 |
& rm2 = one/m2 |
9 |
Integer :: x1, x2 |
10 |
! ---------------------------------------------------------------------- |
11 |
Contains |
12 |
|
13 |
Function dran0() Result( ran ) |
14 |
! ----------------------------------------------------------------- |
15 |
Implicit None |
16 |
|
17 |
! ----------------------------------------------------------------- |
18 |
! --- dran0 returns a uniform deviate in [0,1). |
19 |
! |
20 |
! --- The algorithm is loosely based on an algorithm from |
21 |
! Press & Teukolsky et.al. and based on the linear congruential |
22 |
! method with choices for M, A, and C that are given by |
23 |
! D. Knuth in "Semi-numerical algorithms". |
24 |
! --- Input parameters: |
25 |
! Integer - a1, c1, m1, a2, c2, m2. The parameters of the two |
26 |
! linear congruent relations used. They are passed |
27 |
! via module 'ran_module'. |
28 |
! Integer - x1, x2. Seeds for the two linear congruences. |
29 |
! Passed via module 'ran_module'. |
30 |
! |
31 |
! --- Output-parameters: |
32 |
! Real(l_) - ran. Uniform deviate in [0,1) |
33 |
! ------------------------------------------------------------------ |
34 |
! |
35 |
Real(l_) :: ran |
36 |
! ------------------------------------------------------------------ |
37 |
x1 = Mod( a1*x1 + c1, m1 ) |
38 |
x2 = Mod( a2*x2 + c2, m2 ) |
39 |
ran = ( Real( x1, l_ ) + Real( x2, l_ )*rm2 )*rm1 |
40 |
! ----------------------------------------------------------------- |
41 |
End Function dran0 |
42 |
|
43 |
Subroutine ranfil( a, n ) |
44 |
! ---------------------------------------------------------------------- |
45 |
Implicit None |
46 |
|
47 |
Integer :: n |
48 |
Real(l_) :: a(n) |
49 |
|
50 |
Integer :: i |
51 |
! ---------------------------------------------------------------------- |
52 |
Do i = 1, n |
53 |
a(i) = dran0() |
54 |
End Do |
55 |
! ---------------------------------------------------------------------- |
56 |
End Subroutine ranfil |
57 |
|
58 |
Subroutine rinit |
59 |
! ---------------------------------------------------------------------- |
60 |
Implicit None |
61 |
|
62 |
Integer :: i |
63 |
! ---------------------------------------------------------------------- |
64 |
Do i = 1, 97 ! --- Warming up phase. |
65 |
x1 = Mod( a1*x1 + c1, m1 ) |
66 |
x2 = Mod( a2*x2 + c2, m2 ) |
67 |
End Do |
68 |
! ---------------------------------------------------------------------- |
69 |
End Subroutine rinit |
70 |
|
71 |
! ----------------------------------------------------------------- |
72 |
End Module ran_module |