NB. Non-isotropic scattering of a beam incident on a slab. NB. Surface is at z=T, bottom at z=0. Beam enters at (x,y,z)=(0,0,T). NB. Scattering following the Henyey-Greenstein phase function (parameter NB. gHG). Here, T=slab thickness, hnmu= number of (positive) mu bins NB. k= number of Monte Carlo scattering photons, scatter "nscat" times NB. mu_i= cos(theta), theta is the angle of incident beam to normal NB. since this is a lambda-iteration, we must set nscat >> T^2 NB. Usage: (T; mu_i; gHG; hnmu; nxy; dxy; xoffset) Beam_HG k; nscat NB. hnmu= number of (positive) mu bins, nxy= number of *positive* x and y NB. bins, output size [nmu,2*nxy,2*nxy], dxy= size of x & y bins NB. 16 Jun 2020: fixed error in cross product - vectors are (3,n) not (n,3). Beam_Spot_HG=: 4 : 0 'T mu_i gHG hnmu nang nxy dxy xoffset'=: x NB. T = slab's optical thickness 'k nscat'=: y NB. number of initial photons and scatterings nmu=: +: hnmu NB. total (positive and negative) mu bins mubin=. BIN fr1=. _1+ 2*int01 nmu NB. frets to divide up the [-1,1] range of mu fr2=. 1p1*int01 nang NB. frets in angle over [0,pi] abins=. BIN fr2=. (-|. }.fr2),fr2 NB. frets in angle over [-pi,p1] fr4=. nxy*dxy*(_1+2*int01 +:nxy) NB. divide x-coordinates of emergent point ybins=. BIN fr4 NB. midpoints of bins xbins=. BIN fr3=. xoffset+ fr4 bfr=: trim &.> fr1;fr2;fr3;fr4 NB. boxed trimmed frets 'fr1 fr2 fr3 fr4'=. bfr NB. trimmed frets Sorted=. 0$~(>(1+#)&.>bfr) NB. place holder for sorted results SS=: (1,$Sorted)$ ,Sorted Int=. k# Fscat=: F T%mu_i NB. Int = fraction of beam scattered, F=1-^-tau DD=: Dt Fscat*rand k NB. distance along ray to scattering point mu_i=. - mu_i NB. beam is directed downwards z=. T+ mu_i*DD NB. z's of initially scattered photons xx=: DD*sin_th=. %:(1-*:mu_i) NB. x-distance along ray P=. >xx;0;z NB. starting x,y,z coordinates of all photons C=. k #"0 >sin_th;0;mu_i NB. (x,y,z) direction cosines of incoming beam NB. -------------------------------- Step^:(nscat) Int;C;P NB. do "nscat" escape/scatter cycles NB. -------------------------------- SS=: k%~ }.SS NB. normalize to total (+/+/+/Sorted)= 1. Sorted=. +/SS upind=. hnmu+ i. hnmu NB. mu indicies of upward radiation mubins=. upind {mubin NB. the seperated u & d just use 00)*mu%~(T-z=.2{P) Pout=. P+ DD*"1 C NB. coordinates of exiting photon's escape point Rout=. L }: Pout NB. radial distance from z-axis of escape point phi_out=. atan3 |: }: Pout NB. angular coordinate of escape point phi_ray=. atan3 |: }: C NB. angular coordinate of escape trajectory rphi=. phi_ray- phi_out NB. angle between position & trajectory xout=. Rout* cos rphi yout=. Rout* sin rphi Iesc=. Int* extn=. ^-DD NB. flux of escaping photons Sp=. bfr Sort mu;phi;xout;yout;Iesc SS=: SS, Sp NB. and append to the earlier escapes Int=. Int*(1-extn) NB. intensities of scattered fraction remain=: k%~ +/Int NB. integrated intensity still in slab Dtau=. Dt (1-extn)*rand k NB. travel distance of scattered (fractional) photons P=. P+ Dtau *"1 C NB. new position of next scattering Int;C;P ) NB. The inverse of the accumulated Henyey-Greenstein NB. phase function. "g iAHG r" for random numbers in NB. the interval [0,1] gives a distribution of mu= NB. cos(theta)'s over [-1,1] following the H-G function. iAHG=: 4 : 0 s=. _1+ +:y a=. 1+ g2=. *: g=. x if. 1e_5<|g do. b=. *:(1-g2)% >:g*s -: g%~ a-b else. s2=. 1- *:s s+ (1.5*g*s2)- +:g2*s*s2 end. ) NB. Compute new direction cosines from scattering angles phi and THeta. NB. Usage: dc' Swerve Theta,:phi -->dc (where dc' is incoming and dc is out) NB. dc' is [3 n] array and Theta,:phi a [2 n] array of angles NB. assumes Theta=0 ==> no deflection, Theta=pi ==> backscatter NB. phi is angle counterclockwise from z-axis plane along direction of propagation. Swerve=: 4 : 0 NB. x = DC's of incoming photons v1=. norm (1 0 2){ (_1 1 0)*x NB. [vectors perp to dc',z planes: dc' cross (0 0 1)] q1=. v1 Qrot 0{y NB. quaternion to rotate by Theta in (dc',z_axis) plane q2=. (-x) Qrot 1{y NB. and to rotate counter-clock by phi about dc' q2=. q2 QM q1 NB. quaternion for sequence of both rotations norm q2 Pr x NB. do rotation: norm fixes case of dc'=(0 0 1) ) NB. Given the array A, sort it into bins and sum bin contents. The NB. bin edges are provided in the arrays fr1,fr2,fr3 & fr4, and we` NB. and sum the contents of each bin. The arrays a1, a2, a3, and a4 NB. are the parameters of A that we are sorting on. The output Sp has NB. dimensions (1+#fr1,1+#fr2,1+#fr3,1+#fr4). NB. Usage: Sp=. (fr1;fr2;fr3;fr4) Sort a1;a2;a3;a4;A NB. (Thanks to Raul Miller and Roger Hui.) Sort=: 4 : 0 b=. |:x I.&> }:y idx=. <"1 ~. b (b +//. >{:y ) idx} 0$~>:#&>x ) Dt=: [: -[:^. 1-] NB. "Dt rand n" maps [0,1] into travel on 0:%] NB. int01 n --> divides [0 ... 1] into n intervals. atan3=: 13 : '(_3&o. %/"1 y)+1p1*(0>{:"1 y)*(*{."1 y)' NB. _1p1< out <1p1 BIN =: -:@(}. + }:) NB. midpoints within series F=: 1- [: ^- NB. 1-exp{-y} = fraction scattered by distance y trim=: [: }: }. NB. trim off both ends NB. Quaternion maths. NB. quaternions stored as scalars followed by 3-vectors NB. array of n quaternions stored as (4 n) array. L=: +/ &. (*:"_) NB. Gives the lengths of vector arrays dot=: [: +/ * NB. dot products of arrays of vectors norm=: ] %"1 L NB. normalize to unit length NB. vector cross product (right handed coordinate system): pm1=: 1&|. NB. 1st cyclic permutation pm2=: 2&|. NB. 2nd cyclic permutation pm12=: 4 : '(pm1 x)*pm2 y' cross=: pm12 - pm12~ NB. For vectors stored in (3 n) form NB. Multiply quaternions: q1 QM q2 --> q3 QM=: 4 : 0 s1=. 0{ x v1=. }. x s2=. 0{ y v2=. }. y s=. (s1*s2) - v1 dot v2 s,(s1*"1 v2)+(s2*"1 v1)+(v1 cross v2) ) NB. Quaternion for rotation of points by theta about unit vector u NB. (Rotates *points*; use -theta to rotate coordinate axes) NB. Usage: u Qrot theta --> q then q Pr P --> P' NB. u may be [3 n] array and theta [n] array. Qrot=: 4 : 0 s=. cos -: y f=. sin -: y s, (f*"1 x) ) NB. The conjuate of the quaternions: Qb q --> q_bar Qb=: 3 : '(1 _1 _1 _1)* y' NB. Rotate point using quaternion q NB. Usage: q Pr P --> P' where q is [4 n] array of quaternions, NB. P is [3 n] array of point coordinates, and P' also [3 n] NB. Pr=: 4 : '}. x QM (0,y) QM Qb x' NB. below is tacit form of this: Pr=: [: }. [ QM (0"_ , ]) QM [: Qb [ NB. Due to the symmetry of the problem, if we view the plane from -phi we will NB. see the same pattern as seen from +phi, except that the pattern is reflected NB. about the x-axis, so that y becomes (-y). Thus we can average the pattern NB. and its reflection to improve the statistics. E.g., Sd_av=. reflect Sd. reflect=: 3 : '-: y + |."1 |."3 y' NB. For an albedo other than 1, we form the albedo vector: NB. albvec=. albedo^ >: i. nscat NB. Then modify RNu=. albvec* SSNu and RNd=. albvec* SSNd NB. We then recompute SSu and SSd: Ru=. +/SSNu and Rd=. +/SSNd NB. Run "Beam_Spot_HG" ncycle times: BSCYC=: 3 : 0 'T mu_i gHG hnmu nang nxy dxy xoffset nphoton nsct ncycle'=. y R0=. (T;mu_i;gHG;hnmu;nang;nxy;dxy;xoffset) Beam_Spot_HG nphoton;nsct 'mubins abins xbins ybins SSu SSd SSNu SSNd'=. R0 count=. 1 remainS=. remain while. count< ncycle do. Rn=. (T;mu_i;gHG;hnmu;nang;nxy;dxy;xoffset) Beam_Spot_HG nphoton;nsct 'mubins abins xbins ybins Su Sd SNu SNd'=. Rn SSu=. SSu + Su SSd=. SSd + Sd SSNu=. SSNu+ SNu SSNd=. SSNd+ SNd remainS=. remainS + remain count=. >: count ((": count),LF) (1!:2)