subroutine aptmaxw (tgas, np, erest, beta, gamma, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTMAXW
c
c call aptmaxw (tgas, np, erest, beta, gamma, nerr)
c
c Version: aptmaxw Updated 1990 May 7 13:10.
c aptmaxw Originated 1990 January 31 16:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np values of the relativistic velocity
c functions beta and gamma from np relativistic Maxwellian
c distributions at temperatures tgas, for particles with a rest
c mass energy of erest (same units as tgas).
c Flag nerr indicates any input error.
c
c The distribution function is:
c
c p(beta, tgas) = a * gamma**5 * exp (-f * gamma)
c
c where a = f * beta**2 / K2 (f), f = erest / tgas,
c and K2 is the modified Bessel function of the second kind.
c
c Input: tgas, np, erest.
c
c Output: beta, gamma, nerr.
c
c Glossary:
c
c beta Output Ratio of particle velocity to the speed of light:
c beta = v / c (c = 2.99792458e+08 m / s).
c Range is from 0.0 to 1. Size np.
c
c erest Input The rest mass of the particles (same units as tgas).
c
c gamma Output Ratio of relativistic mass to particle rest mass:
c gamma = 1.0 / sqrt (1.0 - beta**2).
c Range is from 1.0 to infinity. Size np.
c
c nerr Output Indicates an input error, if not zero.
c 1 if np is not positive.
c 2 if erest is not positive.
c
c np Input Size of arrays tgas, beta, gamma.
c
c tgas Input Average particle temperature (same units as erest).
c Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Ratio of velocity to light speed.
dimension beta (1)
c---- Ratio of mass to rest mass.
dimension gamma (1)
c---- Average particle temperature.
dimension tgas (1)
c.... Local variables.
c---- Sqrt (besq).
common /laptmaxw/ be (64)
c---- 0.5 * tgas / erest.
common /laptmaxw/ besq (64)
c---- Working variable.
common /laptmaxw/ csasq (64)
c---- Working variable.
common /laptmaxw/ csbsq (64)
c---- Ratio of kinetic energy to tgas.
common /laptmaxw/ etsq
c---- A random fraction of psumd.
common /laptmaxw/ fsum
c---- Index in local array.
common /laptmaxw/ n
c---- First index of subset of data.
common /laptmaxw/ n1
c---- Last index of subset of data.
common /laptmaxw/ n2
c---- Index in external array.
common /laptmaxw/ nn
c---- Size of current subset of data.
common /laptmaxw/ ns
c---- Numerical constant, pi.
common /laptmaxw/ pi
c---- 0.25 * sqrt (pi).
common /laptmaxw/ psuma
c---- Partial sum of a PDF factor.
common /laptmaxw/ psumb (64)
c---- Partial sum of a PDF factor.
common /laptmaxw/ psumc (64)
c---- Partial sum of a PDF factor.
common /laptmaxw/ psumd (64)
c---- Log (1.0 / ranf( )).
common /laptmaxw/ snorm
c---- Sum csasq**2 + csbsq**2.
common /laptmaxw/ sumsq (64)
c---- Normal distribution factor.
common /laptmaxw/ xnsq
c---- Next value of xnsq.
common /laptmaxw/ xnsqs
c---- Normal distribution factor.
common /laptmaxw/ ynsq
c---- Next value of ynsq.
common /laptmaxw/ ynsqs
cbugc***DEBUG begins.
cbugc---- Standard deviation from mean.
cbug common /laptmaxw/ argdev
cbugc---- Mean value.
cbug common /laptmaxw/ argmean
cbug 9901 format (/ 'aptmaxw finding Maxwellian energies. np=',i3 /
cbug & ' erest= ',1pe22.14 /
cbug & (i3,' tgas=',1pe22.14))
cbug write ( 3, 9901) np, erest, (n, tgas(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
c---- Numerical constant, pi.
pi = 3.14159265358979323
c---- Normal distribution seed.
xnsqs = -1.0
c---- Normal distribution seed.
ynsqs = -1.0
c.... Test for input errors.
nerr = 0
if (np .le. 0) then
nerr = 1
go to 210
endif
if (erest .le. 0.0) then
nerr = 2
go to 210
endif
c.... Initialize.
c---- 0.44311346272637897.
psuma = 0.25 * sqrt (pi)
c.... Do data sets in groups of 64.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find the partial sums of one factor of the relativistic Maxwellian
c.... probability distribution function from which the kinetic energy
c.... is sampled.
c---- Loop over subset of data.
do 120 n = 1, ns
besq(n) = 0.5 * tgas(n+n1-1) / erest
be(n) = sqrt (besq(n))
psumb(n) = psuma + 0.5 * be(n)
psumc(n) = psumb(n) + 3.0 * psuma * besq(n)
psumd(n) = psumc(n) + 2.0 * be(n) * besq(n)
c---- End of loop over subset of data.
120 continue
c.... Sample a kinetic energy from a relativistic Maxwellian distribution,
c.... and find beta and gamma.
c---- Loop over subset of data.
do 160 n = 1, ns
130 fsum = psumd(n) * ranf( )
c---- Sample from x**2 * exp (-x**2).
if (fsum .le. psuma) then
c---- Find new xnsq and xnsqs.
if (xnsqs .le. 0.0) then
140 csasq(n) = ranf( )**2
csbsq(n) = ranf( )**2
sumsq(n) = csasq(n) + csbsq(n)
if (sumsq(n) .ge. 1.0) go to 140
snorm = -alog (ranf( ))
xnsq = snorm * csasq(n) / sumsq(n)
xnsqs = snorm * csbsq(n) / sumsq(n)
c---- Use last xnsqs.
else
xnsq = xnsqs
xnsqs = 0.0
c---- Tested xnsqs.
endif
etsq = xnsq - alog (ranf( ))
c++++ Sample from x**3 * exp (-x**2).
elseif (fsum .le. psumb(n)) then
etsq = -alog (ranf( ) * ranf( ))
c++++ Sample from x**4 * exp (-x**2).
elseif (fsum .le. psumc(n)) then
c---- Find new ynsq and ynsqs.
if (ynsqs .le. 0.0) then
150 csasq(n) = ranf( )**2
csbsq(n) = ranf( )**2
sumsq(n) = csasq(n) + csbsq(n)
if (sumsq(n) .ge. 1.0) go to 150
snorm = -alog (ranf( ))
ynsq = snorm * csasq(n) / sumsq(n)
ynsqs = snorm * csbsq(n) / sumsq(n)
c---- Use last ynsqs.
else
ynsq = ynsqs
ynsqs = 0.0
c---- Tested ynsqs.
endif
etsq = ynsq - alog (ranf( ) * ranf( ))
c---- Sample from x**5 * exp (-x**2).
else
etsq = -alog (ranf( ) * ranf( ) * ranf( ))
c---- Tested fsum.
endif
c.... Reject the kinetic energy, on the factor
c.... sqrt (1.0 + besq * etsq) / (1.0 + sqrt (besq * etsq)).
csasq(n) = ranf( )**2
if ((((1.0 - csasq(n)) * (1.0 + besq(n) * etsq))**2) .le.
& (4.0 * besq(n) * etsq * csasq(n)**2)) go to 130
nn = n + n1 - 1
gamma(nn) = 1.0 + tgas(nn) * etsq / erest
beta(nn) = sqrt ((gamma(nn) + 1.0) * (gamma(nn) - 1.0)) /
& gamma(nn)
c---- End of loop over subset of data.
160 continue
c.... See if all sets are done.
c---- Do another data set.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptmaxw results: beta, gamma')
cbug write ( 3, 9902)
cbug call aptmean (beta, np, 1.e-11, argmean, argdev, nerr)
cbug call aptmean (gamma, np, 1.e-11, argmean, argdev, nerr)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptmaxw. (+1 line.)
end
UCRL-WEB-209832