subroutine aptscan (kph, au, av, aw, ph, np, bu, bv, bw, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTSCAN c c call aptscan (kph, au, av, aw, ph, np, bu, bv, bw, nerr) c c Version: aptscan Updated 1990 May 2 11:00. c aptscan Originated 1990 May 1 14:00. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To randomly sample np unit direction vectors b = (bu, bv, bw), c each from a cone at an angle phi from a unit direction vector c a = (au, av, aw) in 3-D space. c The angle phi is specified by "ph", based on the option "kph". c Flag nerr indicates any input error (np not positive). c c Input: kph, au, av, aw, ph, np. 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 Size np. Coordinates u, v and w may be any three c orthogonal coordinates. c c bu Output The u component of a unit vector sampled randomly from c all directions in 3-D space that are a specified c angle (see kph, ph) from unit vector "a". Size np. c c bv Output The v component of a unit vector sampled randomly from c all directions in 3-D space that are a specified c angle (see kph, ph) from unit vector "a". Size np. c c bw Output The w component of a unit vector sampled randomly from c all directions in 3-D space that are a specified c angle (see kph, ph) from unit vector "a". Size np. c c kph Input Specifies the meaning of ph: c 0 if ph is the scattering angle in degrees. c 1 if ph is the scattering angle in radians. c 2 if ph is the cosine of the scattering angle. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c 2 if kph is not 0, 1 or 2. c c np Input Size of arrays au, av, aw, ph, bu, bv, bw. c c ph Input Specifies the scattering angle in degrees, from c 0.0 to 180.0 (kph = 0), in radians, from 0.0 to pi c (kph = 1), or its cosine, from -1.0 to 1.0 (kph = 2). c Size np. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Component u of unit vector "a". dimension au (1) c---- Component v of unit vector "a". dimension av (1) c---- Component w of unit vector "a". dimension aw (1) c---- Component u of random unit vector. dimension bu (1) c---- Component v of random unit vector. dimension bv (1) c---- Component w of random unit vector. dimension bw (1) c---- Scattering angle specifier. dimension ph (1) c.... Local variables. c---- Cosine of scattering angle. common /laptscan/ cosph (64) c---- Cosine of azimuth angle. common /laptscan/ costh c---- Component u of random unit vector. common /laptscan/ cu (64) c---- Component v of random unit vector. common /laptscan/ cv (64) c---- Component w of random unit vector. common /laptscan/ cw (64) c---- Function 1.0 / (1.0 + aw + fuz). common /laptscan/ faw c---- A very small number. common /laptscan/ fuz c---- Index in local array. common /laptscan/ n c---- First index of subset of data. common /laptscan/ n1 c---- Last index of subset of data. common /laptscan/ n2 c---- Index in external array. common /laptscan/ nn c---- Size of current subset of data. common /laptscan/ ns c---- Scattering angle in radians. common /laptscan/ phi c---- Mathematical constant, pi. common /laptscan/ pi c---- Random number (0.0 to 1.0). common /laptscan/ ranfp (64) c---- Rotation matrix. common /laptscan/ rotm (3,3) c---- Sine of scattering angle. common /laptscan/ sinph (64) c---- Sine of angle of rotation. common /laptscan/ sinr c---- Sine of angle of rotation + fuz. common /laptscan/ sinr1 c---- Sine of azimuth angle. common /laptscan/ sinth c---- Azimuth angle. common /laptscan/ theta cbugc***DEBUG begins. cbug common /laptscan/ avgu, avgv, avgw, devu, devv, devw cbug common /laptscan/ avgu2, avgv2, avgw2, sumsqs cbug 9901 format (/ 'aptscan scattering by PH around vector A:' / cbug & ' kph=',i2 / cbug & (i3,' ph= ',1pe22.14 / cbug & ' au,av,aw=',1p3e22.14)) cbug write (3, 9901) kph, (n, ph(n), au(n), av(n), aw(n), n = 1, np) cbugc***DEBUG ends. c.... initialize. c---- A very small number. fuz = 1.e-99 c---- Mathematical constant, pi. pi = 3.14159265358979323 nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 go to 210 endif if ((kph .lt. 0) .or. (kph .gt. 2)) then nerr = 2 go to 210 endif c.... Set up the indices of the first subset of data. n1 = 1 n2 = min (np, 64) c.... Loop over the data subsets. 110 ns = n2 - n1 + 1 c.... Find the scattering angle cosine and sine. c---- Unit of ph is degrees. if (kph .eq. 0) then c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 phi = ph(nn) * pi / 180.0 cosph(n) = cos (phi) sinph(n) = sin (phi) c---- End of loop over subset of data. 120 continue c---- Unit of ph is radians. elseif (kph .eq. 1) then c---- Loop over subset of data. do 130 n = 1, ns nn = n + n1 - 1 phi = ph(nn) cosph(n) = cos (phi) sinph(n) = sin (phi) c---- End of loop over subset of data. 130 continue c---- Unit of ph is cosine of angle. elseif (kph .eq. 2) then c---- Loop over subset of data. do 140 n = 1, ns nn = n + n1 - 1 cosph(n) = ph(nn) sinph(n) = sqrt (1.0 - cosph(n)**2) c---- End of loop over subset of data. 140 continue endif c.... Generate the needed random numbers. c---- Loop over subset of data. do 150 n = 1, ns ranfp(n) = ranf( ) c---- End of loop over subset of data. 150 continue c.... Sample randomly on a cone at angle phi from direction (0, 0, 1). c---- Loop over vectors. do 160 n = 1, ns theta = 2.0 * pi * ranfp(n) costh = cos (theta) sinth = sin (theta) cu(n) = costh * sinph(n) cv(n) = sinth * sinph(n) cw(n) = cosph(n) c---- End of loop over vectors. 160 continue c.... Rotate random vectors to a cone at angle phi from unit vector "a". c---- Loop over vectors. do 170 n = 1, ns nn = n + n1 - 1 sinr = sqrt (au(nn)**2 + av(nn)**2) sinr1 = sinr + fuz faw = (1.0 - aw(nn)) / sinr1**2 rotm(1,1) = av(nn)**2 * faw + aw(nn) rotm(1,2) = -au(nn) * av(nn) * faw rotm(1,3) = au(nn) * sinr / sinr1 rotm(2,1) = -au(nn) * av(nn) * faw rotm(2,2) = au(nn)**2 * faw + aw(nn) rotm(2,3) = av(nn) * sinr / sinr1 rotm(3,1) = -au(nn) * sinr / sinr1 rotm(3,2) = -av(nn) * sinr / sinr1 rotm(3,3) = aw(nn) bu(nn) = cu(n) * rotm(1,1) + cv(n) * rotm(1,2) + & cw(n) * rotm(1,3) bv(nn) = cu(n) * rotm(2,1) + cv(n) * rotm(2,2) + & cw(n) * rotm(2,3) bw(nn) = cu(n) * rotm(3,1) + cv(n) * rotm(3,2) + & cw(n) * rotm(3,3) c---- End of loop over vectors. 170 continue c.... See if all data subsets are done. c---- Do another subset of data. if (n2 .lt. np) then n1 = n2 + 1 n2 = min (np, n1 + 63) go to 110 endif cbugc***DEBUG begins. cbug cbug 9902 format (/ 'aptscan results:') cbug 9903 format (i3,' bu,bv,bw=',1p3e22.14 / cbug & ' sumsq= ',1pe22.14) cbug 9904 format (/ 'aptscan mean and deviation:' / cbug & ' avg(u,v,w)=',1p3e22.14 / cbug & ' dev(u,v,w)=',1p3e22.14) cbug cbug write ( 3, 9902) cbug cbug avgu = 0.0 cbug avgv = 0.0 cbug avgw = 0.0 cbug avgu2 = 0.0 cbug avgv2 = 0.0 cbug avgw2 = 0.0 cbug cbug do 180 n = 1, np cbug avgu = avgu + bu(n) cbug avgv = avgv + bv(n) cbug avgw = avgw + bw(n) cbug avgu2 = avgu2 + bu(n)**2 cbug avgv2 = avgv2 + bv(n)**2 cbug avgw2 = avgw2 + bw(n)**2 cbug sumsqs = bu(n)**2 + bv(n)**2 + bw(n)**2 cbug write ( 3, 9903) n, bu(n), bv(n), bw(n), sumsqs cbug 180 continue cbug cbug avgu = avgu / np cbug avgv = avgv / np cbug avgw = avgw / np cbug avgu2 = avgu2 / np cbug avgv2 = avgv2 / np cbug avgw2 = avgw2 / np cbug devu = sqrt (avgu2 - avgu**2) cbug devv = sqrt (avgv2 - avgv**2) cbug devw = sqrt (avgw2 - avgw**2) cbug cbug write ( 3, 9904) avgu, avgv, avgw, devu, devv, devw cbugc***DEBUG ends. 210 return c.... End of subroutine aptscan. (+1 line.) end UCRL-WEB-209832