[nuclear] ~/for > cat all/random_walk.for
subroutine random_walk
REAL VARIA(23)
real dummy(19)
COMMON /urandom_walk/ rit, dist, rNum, r_end, dummy
EQUIVALENCE (VARIA, rit)
character*8 NAME(23)
DATA NAME /'rit','dist','rNum',20*' '/
DATA VARIA/ 1, 1, 10, 20*0. /
cjr real*8 co,thetao,si,alo,beo,gao,alb,beb,gab,al,be,ga,thetap, xb,yb
real co,thetao,si,alo,beo,gao,alb,beb,gab,al,be,ga,thetap, xb,yb
character*80 inpLIN ! Buffer for lines input
data twopi / 6.2830/
C *******************
999 imess = MESSmc( NAME,VARIA,
1'iterate this number; -1:Description of Program',
2'distance between collisions (any unit)',
3'Number of displacements (in 3-dim)',
4'Radius at the end of the 3-dim walk ',
4 0)
if (imess.eq.-1) return !^Z (ctrl Z) returns
C
if ( rit .lt. 0) then
WRITE(*,3)
3 FORMAT(//' RANDOM_WALK: ',/
1 ' This program was written to check subroutine DIRECT_COS, ',/
2 ' used in several programs such as PROTON_IMAGING. ',/
3 ' RANDOM_WALK deals with the isotropic random walk of an ',/
4 ' object (molecule) in 3-dimensions. That is, a molecule in a ',/
5 ' gas moves equal distances D between collisions with equal ',/
6 ' probability in any direction. After a total of N of such ',/
7 ' displacements, the mean square displacement from the ORIGEN ',/
8 ' is SUM(r^2)/N, and the RMS displacement is ',/
9 ' rms(r) = sqrt( SUM(r^2)/N ) ',/
9 ' This problem can be solved analytically (see e.g. Rice, prob',/
1 ' 1.18), and the answwer is ',/
1 ' rms(r) = D * sqrt(N) ',/
1 ' So by running a Monte carlo in this program, one can test ',/
2 ' that the transport of the particle is done correctly, by ',/
3 ' generating a distribution and calculating rms(r) ',/
4 ' (from the ORIGEN, not from the average of the distribution).' )
rit = 1
cjr accept '(A)',I
read (*,'(a)') inpLIN ! read line
go to 999
end if
c pick initial phi, theta (i.e., phio,thetao):
Rn = random() !random number between 0 and 1
co = 1. - 2.*Rn !isotropic dist. in thetap
cjr thetao = Dacos (co)
thetao = acos (co)
cjr si = Dsqrt(1.d0-co**2)
si = sqrt(1.d0-co**2)
Rn=random() !random number between 0 and 1
phio = twopi*Rn
alb = cos(phio)*si !init direction cosine along X-axis for unit vector
beb = sin(phio)*si ! Y-axis
gab = co ! Z-axis
xb = 0
yb = 0
zb = 0
cjr thetaf = Dacos(co) * 1000
thetaf = acos(co) * 1000
Number_steps = rNum
rb = dist
do i=1, Number_steps
c propagate along direction of unit vector k':
xb = xb + alb*rb !New position in XYZ coord system
yb = yb + beb*rb
zb = zb + gab*rb
c Set for iteration algorithm
al = alb
be = beb
ga = gab
c pick phip, thetap (thetap is theta', phip is phi') (in system X'Y'Z', wich
c has its Z' axis along last direction k')
Rn=random() !random number between 0 and 1
co = 1. - 2.*Rn !isotropic dist in thetap for 3-d random walk
cjr thetap = Dacos(co)
thetap = acos(co)
Rn=random() !random number between 0 and 1
phip = twopi*Rn
call direct_cos (thetap, phip, al, be, ga, 0., 0,
1 alb, beb, gab)
end do
2 continue
cjr r_end = Dsqrt (xb**2 + yb**2 + zb**2 )
r_end = sqrt (xb**2 + yb**2 + zb**2 )
GO TO 999
END !of s.r. random_walk
c*******************************************************************
subroutine direct_cos (thetap,phip,al,be,ga,deg,iflg, alb,beb,gab)
c iflg = 0, angles in rad; = 1, angles in deg ( deg = 180./3.1415)
thetapr = thetap
phipr = phip
if (iflg.eq.1) then !thetap and phip in deg
thetapr = thetap/deg !in rad
phipr = phip/deg !in rad
end if
s = sin (thetapr)
co = cos (phipr)
si = sin (phipr)
alp = co*s
bep = si*s
gap = cos (thetapr)
c type *, alp,bep,gap
den = sqrt(al*al + be*be)
c What follows depends on axis convention
c Take z' axis along incident unit vector k'
c Take x' axis along unit vector j', defined as : j' = (k x k')//k x k'/
c Take y' axis as a right-hand system: i' = j'x k'
c Unit vector k is taken along Z axis in the fixed cartesian system X Y Z.
gab = gap*ga - alp*den
if (al.eq.0 .and. be.eq.0) then
alb = 0
beb = 0
return
end if
alb = gap*al + (-bep*be + alp*al*ga)/den !notice convention
beb = gap*be + ( bep*al + alp*be*ga)/den
return
end !of subroutine direct_cos
c*****************************************************************
c This is not used. It uses an incovenient convention.
subroutine direct_cost (thetap,phip,al,be,ga,deg,iflg, als,bes,gas)
cjr real*8 thetap,phip,al,be,ga,als,bes,gas, thetapr,phipr
real thetap,phip,al,be,ga,als,bes,gas, thetapr,phipr
cjr real*8 s, co, si, alp, bep, gap, den
real s, co, si, alp, bep, gap, den
thetapr = thetap
phipr = phip
if (iflg.eq.1) then !thetap and phip in deg
thetapr = thetap*deg !in rad
phipr = phip*deg !in rad
end if
cjr s = Dsin (thetapr)
s = sin (thetapr)
cjr co = Dcos (phipr)
co = cos (phipr)
cjr si = Dsin (phipr)
si = sin (phipr)
alp = co*s
bep = si*s
cjr gap = Dcos (thetapr)
gap = cos (thetapr)
c type *, alp,bep,gap
cjr den = Dsqrt(al*al + be*be)
den = sqrt(al*al + be*be)
c What follows depends on axis convention
c Take z' axis along incident unit vector k'
c Take x' axis along unit vector i', defined as : i' = (k'x k)//k'x k/
c Take y' axis as a right-hand system: j' = k'x i'
c Unit vector k is taken along Z axis in the fixed cartesian system X Y Z.
gas = gap*ga - bep*den
if (al.eq.0 .and. be.eq.0) then
als = 0
bes = 0
return
end if
als = gap*al + (alp*be + bep*al*ga)/den !notice convention
bes = gap*be + (-alp*al + bep*be*ga)/den
return
end !of subroutine direct_cost
C---------------------------------------------------------------------------