Bielajew A.F. Fundamentals of the Monte Carlo method for neutral and charged particle transport
Подождите немного. Документ загружается.
18.1. UTILITY/GENERAL 301
stats.f
Statistical analysis routine.
Given:
N
X
i=1
x
i
,
N
X
i=1
x
2
i
and N,
calculates:
x =
1
N
N
X
i=1
x
i
and s
x
=
v
u
u
t
1
N
P
N
i=1
x
2
i
− x
2
N − 1
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine stats(sumx,sumx2,avgx,sigmax,N)
C Given a sum of values and the associated sum of squares, this routine
C computes the average and the estimated error of the mean.
implicit none
integer N ! Input: Number of samples
real
* sumx, ! Input: The sum of measurements
* sumx2 ! Input: The sum of the squares of measurements
real
* avgx, ! Ouput: The average
* sigmax ! Ouput: The estimated error of the mean
avgx = sumx/N
sigmax = sqrt(max(0e0,sumx2/N - avgx**2)/(N-1))
return
end
302 CHAPTER 18. CODE LIBRARY
18.2 Subroutines for random number generation
rngmc1.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 ≤ r ≤ 1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 69069, c =0,X
0
= 123456789.
The sequence length is 2
30
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple multiplicative congruential pseudo random number generator using
C one of Knuth’s magic multpliers, 69069
C Sequence length: 2**30
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 69069 ! The scrambling process
rng = 0.5 + ixx*0.23283064e-09 ! The conversion to real 0 <= rng <= 1
return
end
18.2. SUBROUTINES FOR RANDOM NUMBER GENERATION 303
rngmc2.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 <r<1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 69069, c =0,X
0
= 123456789.
The sequence length is 2
30
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple multiplicative congruential pseudo random number generator using
C one of Knuth’s magic multpliers, 69069
C Sequence length: 2**30
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 69069 ! The scrambling process
rng = 0.5 + ixx*0.2328306e-09 ! The conversion to real0<rng<1
return
end
304 CHAPTER 18. CODE LIBRARY
rngmc3.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 ≤ r ≤ 1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 663608941, c =0,X
0
= 123456789.
The sequence length is 2
30
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C13456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple multiplicative congruential pseudo random number generator using
C the default magic multiplier, 663608941
C Sequence length: 2**30
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 663608941 ! The scrambling process
rng = 0.5 + ixx*0.23283064e-09 ! The conversion to real 0 <= rng <= 1
return
end
18.2. SUBROUTINES FOR RANDOM NUMBER GENERATION 305
rngmc4.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 <r<1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 663608941, c =0,X
0
= 123456789.
The sequence length is 2
30
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C13456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple multiplicative congruential pseudo random number generator using
C the default magic multiplier, 663608941
C Sequence length: 2**30
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 663608941 ! The scrambling process
rng = 0.5 + ixx*0.2328306e-09 ! The conversion to real0<rng<1
return
end
306 CHAPTER 18. CODE LIBRARY
rnglc1.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 ≤ r ≤ 1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 69069, c = 987654321, X
0
= 123456789.
The sequence length is 2
32
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C13456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple congruential pseudo random number generator using
C one of Knuth’s magic multpliers, 69069
C Sequence length: 2**32
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 69069 + 987654321 ! The scrambling process
rng = 0.5 + ixx*0.23283064e-09 ! The conversion to real 0 <= rng <= 1
return
end
18.2. SUBROUTINES FOR RANDOM NUMBER GENERATION 307
rnglc2.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 <r<1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 69069, c = 987654321, X
0
= 123456789.
The sequence length is 2
32
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C13456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple congruential pseudo random number generator using
C one of Knuth’s magic multpliers, 69069
C Sequence length: 2**32
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 69069 + 987654321 ! The scrambling process
rng = 0.5 + ixx*0.2328306e-09 ! The conversion to real0<rng<1
return
end
308 CHAPTER 18. CODE LIBRARY
rnglc3.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 ≤ r ≤ 1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 663608941, c = 987654321, X
0
= 123456789.
The sequence length is 2
32
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C13456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple congruential pseudo random number generator using
C the default magic multiplier, 663608941
C Sequence length: 2**32
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 663608941 + 987654321 ! The scrambing process
rng = 0.5 + ixx*0.23283064e-09 ! The conversion to real 0 <= rng <= 1
return
end
18.2. SUBROUTINES FOR RANDOM NUMBER GENERATION 309
rnglc4.f
Calculates:
X
n+1
=mod(aX
n
+ c, 2
32
) ,
and returns random numbers uniform over the range 0 <r<1.:
r
n+1
=1/2+X
n+1
/2
32
,
where a = 663608941, c = 987654321, X
0
= 123456789.
The sequence length is 2
32
.
Restrictions: This routine only works on computer architectures that provide 2’s complement
32-bit integers.
C13456789|123456789|123456789|123456789|123456789|123456789|123456789|12
real function rng()
C A simple congruential pseudo random number generator using
C the default magic multiplier, 663608941
C Sequence length: 2**32
implicit none
integer ixx ! The pseudo-random integer
data ixx/123456789/ ! The intialization value
save ixx ! Must be saved for subsequent passes
ixx = ixx * 663608941 + 987654321 ! The scrambing process
rng = 0.5 + ixx*0.2328306e-09 ! The conversion to real 0 < rng < 1
return
end
310 CHAPTER 18. CODE LIBRARY
rngranlux.f
The original RANLUX subroutine coded by Fred James [Jam94] following the theoretical
development by L¨uscher [L
¨
94].
C The references are Lu:scher, Computer Physics Communications 79 (1994) 100
C and James, CPC 79 (1994) 111.
SUBROUTINE RANLUX(RVEC,LENV)
C Subtract-and-borrow random number generator proposed by
C Marsaglia and Zaman, implemented by F. James with the name
C RCARRY in 1991, and later improved by Martin Luescher
C in 1993 to produce "Luxury Pseudorandom Numbers".
C Fortran 77 coded by F. James, 1993
C
C references:
C M. Luscher, Computer Physics Communications 79 (1994) 100
C F. James, Computer Physics Communications 79 (1994) 111
C
C LUXURY LEVELS.
C ------ ------ The available luxury levels are:
C
C level 0 (p=24): equivalent to the original RCARRY of Marsaglia
C and Zaman, very long period, but fails many tests.
C level 1 (p=48): considerable improvement in quality over level 0,
C now passes the gap test, but still fails spectral test.
C level 2 (p=97): passes all known tests, but theoretically still
C defective.
C level 3 (p=223): DEFAULT VALUE. Any theoretically possible
C correlations have very small chance of being observed.
C level 4 (p=389): highest possible luxury, all 24 bits chaotic.
C
C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C!!! Calling sequences for RANLUX: ++
C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++
C!!! 32-bit random floating point numbers between ++
C!!! zero (not included) and one (also not incl.). ++
C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++
C!!! one 32-bit integer INT and sets Luxury Level LUX ++
C!!! which is integer between zero and MAXLEV, or if ++
C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++
C!!! should be set to zero unless restarting at a break++