subroutine aptscaz (au, av, aw, cosph, bu, bv, bw, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTSCAZ c c call aptscaz (au, av, aw, cosph, bu, bv, bw, nerr) c c Version: aptscaz Updated 1991 August 13 17:30. c aptscaz Originated 1991 July 25 16:30. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample a unit direction vector b = (bu, bv, bw), c at an angle phi from a unit direction vector a = (au, av, aw) c in 3-D space, where cosph = cos (phi). c c Input: au, av, aw, cosph. c c Output: bu, bv, bw, nerr. c c Glossary: c c au,av,aw Input The u, v and w components of the unit vector "a". c Coordinates u, v and w may be any three c orthogonal coordinates, such as x, y and z. c c bu,bv,bw Output The u, v and w components of a unit vector randomly c sampled from a cone at angle phi from unit vector c "a", where cosph = cos (phi). c To get vector "b" isotropic in 3-D space, randomly c sample cosph between -1.0 and 1.0. c c cosph Input The cosine of the angle phi. c c nerr Output Indicates an input error, if not 0. c 1 if cosph is not between -1.0 and 1.0. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Local variables. c---- Absolute value of aw. common /laptscaz/ awabs c---- Component u of a random vector. common /laptscaz/ cu c---- Component v of a random vector. common /laptscaz/ cv c---- Has magnitude 1.0, same sign as aw. common /laptscaz/ faxis c---- Is sin (phi) / normalization factor. common /laptscaz/ fnorm c---- Rotation factor. common /laptscaz/ frota c---- Rotation factor. common /laptscaz/ frotb c---- Normalization factor for cu, cv. common /laptscaz/ sumsq cbugc***DEBUG begins. cbug common /laptscaz/ sumsqs cbug 9901 format (/ 'aptscaz scattering by cosph around vector A:' / cbug & ' cosph= ',1pe22.14 / cbug & ' au,av,aw=',1p3e22.14) cbug write (3, 9901) cosph, au, av, aw cbugc***DEBUG ends. c.... Test input for errors. nerr = 0 if ((cosph .lt. -1.0) .or. (cosph .gt. 1.0)) then nerr = 1 go to 210 endif c.... Randomly sample a unit vector at angle phi from the (+/-)w axis, c.... and do rotation that moves (+/-)w axis onto vector "a". 100 cu = 1.0 - 2.0 * ranf( ) cv = 1.0 - 2.0 * ranf( ) sumsq = cu**2 + cv**2 if (sumsq .gt. 1.0) go to 100 fnorm = sqrt ((1.0 - cosph**2) / sumsq) cu = cu * fnorm cv = cv * fnorm faxis = sign (1.0, aw) awabs = faxis * aw frota = cu * au + cv * av frotb = cosph - frota / (1.0 + awabs) bu = cu + au * frotb bv = cv + av * frotb bw = cosph * aw - faxis * frota cbugc***DEBUG begins. cbug 9902 format (/ 'aptscaz results:') cbug 9903 format (' bu,bv,bw=',1p3e22.14 / cbug & ' sumsq= ',1pe22.14) cbug write ( 3, 9902) cbug sumsqs = bu**2 + bv**2 + bw**2 cbug write ( 3, 9903) bu, bv, bw, sumsqs cbugc***DEBUG ends. 210 return c.... End of subroutine aptscaz. (+1 line.) end UCRL-WEB-209832