NB. Non-isotropic scattering of radiation incident on a slab. The radiation NB. illumintes an infinite plane, so we compute only the angular distribution NB. of reflected and transmitted radiation (I vs. mu(=cos theta) and phi). NB. Slab surface is at z=T, bottom at z=0. Beam enters at (x,y,z)=(0,0,T). NB. Scattering follows the Henyey-Greenstein phase function (parameter g NB. is gHG). Here, T= slab thickness, hnmu= number of (positive) mu bins NB. k= number of Monte Carlo scattering photons, each scatter "nscat" times. NB. mu_i= cos(theta_i), theta_i is the angle of incident beam to the normal. NB. Since this is in effect a lambda-iteration, we must set nscat >> T^2 NB. Usage: (T; mu_i; gHG; hnmu; nang) Beam_HG k; nscat NB. hnmu= number of (positive) mu bins, nang = number of phi angle bins. NB. The albedo=1. For any other albedo, we multiply the series of outputs NB. following each scattering (SNu & SNd) by "albvec" (see below). NB. 16 Jun 2020: fixed error in cross product - vectors are (3,n) not (n,3). Illum_HG=: 4 : 0 'T mu_i gHG hnmu nang'=: x NB. T is total optical depth through slab '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] bfr=: trim &. > fr1;fr2 NB. boxed trimmed frets 'fr1 fr2'=. 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 point of scattering 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 rmn=: +/Int % k NB. -------------------------------- Step^:(nscat) Int;C;P NB. do "nscat" escape/scatter cycles NB. -------------------------------- SS=: k%~ }.SS 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) Iesc=. Int* extn=. ^-DD NB. flux of escaping photons Sp=. bfr Sort mu;phi;Iesc SS=: SS, Sp Int=. Int*(1-extn) NB. intensities of scattered fraction rmn=: rmn, (k%~ +/Int) NB. integrated intensity still in slab (global) 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, which is a function of 2 dimensions, sort it into NB. bins where the bin edges are provided in the arrays fr1 & fr2, NB. and sum the contents of each bin. The arrays a1 and a2 are the NB. coordinates of A that we are sorting on. The output Sp has dimensions NB. (1+#fr1,1+#fr2). Usage: Sp=. (fr1;fr2) Sort a1;a2;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' Pr=: [: }. [ QM (0"_ , ]) QM [: Qb [ reflect=: 3 : '-: y + |."1 y' NB. symmetry about phi=0, average +,- phi's 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 "Illum_HG" ncycle times: ICY=: 3 : 0 'T mu_i gHG hnmu nang nphoton nsct ncycle'=. y R0=. (T;mu_i;gHG;hnmu;nang) Illum_HG nphoton;nsct 'mubins abins SSu SSd SSNu SSNd'=. R0 remain=: rmn count=. 1 while. count<: ncycle do. count=. >: count Rn=. (T;mu_i;gHG;hnmu;nang) Illum_HG nphoton;nsct 'mubins abins Su Sd SNu SNd'=. Rn SSu=. SSu+ Su SSd=. SSd+ Sd SSNu=. SSNu+ SNu SSNd=. SSNd+ SNd remain=: remain + rmn ((": count),LF) (1!:2)0) & (phi<0) results to SAd=. reflect SSd NB. exploit +/- symetery in phi scattering mubins;abins;SSu;SSd;SSNu;SSNd;SAu;SAd ) Cycle_log=: '/your_path/ICY.log' NB. ============================================================================= NB. load'Illum_HG.ijs' NB. OUT=. ICY 5;0.6;0.9;20;30;4e5;50;50 NB. 20 million photons, 1e9 scatters NB. 'mubins abins Su Sd SNu SNd SAu SAd'=. OUT NB. +/+/Su NB. 0.322784 NB. +/+/Sd NB. 0.674425 NB. (+/+/Su)+(+/+/Sd)+ {:remain NB. 0.99976 NB. p1=. 'contour;pensize 3' NB. p1 plot abins; mubins; |: Su NB. p1 plot abins; mubins; |: SAu NB. p1 plot abins; mubins; |: Sd NB. p1 plot abins; mubins; |: +/(i.5){SNd NB. later escapes are NB. p1 plot abins; mubins; |: +/(10+i.25){SNd NB. more diffuse NB. 'pensize 3' plot >(+/+/|:SNu);(+/+/|:SNd) NB. escaping intensity vs. scattering