NB. Scattering of photons from a point source at x=y=0, z=z_s of slab. NB. Scattering following the Henyey-Greenstein phase function. NB. Usage: (T; z_s; gHG; hnmu; nxy; dxy) HG_slab k; nscat NB. Here, T=slab thickness, z_s= location of source along z-axis NB. gHG is the parameter of the H-G scattering, -1< gHG < 1 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. k= number of Monte Carlo scattering photons, scatter "nscat" times HG_slab=: dyad define 'T z_s gHG hnmu nxy dxy'=: x NB. T is total optical depth through slab T2=: -: T NB. optical depth to midplane nmu=: +: hnmu NB. total (positive and negative) mu bins fr1=: _1+ 2*int01 nmu NB. frets to divide up the [0,1] range of |mu| fr2=: nxy*dxy*(_1+2*int01 +:nxy) NB. divide up x-coordinates of emergent point fr3=: nxy*dxy*int01 nxy NB. y-coordinate bins: by symmetry, (y>=0) only Sorted=: 0$~(>:#&>fr1;fr2;fr3) NB. place holder for sorted results 'k nscat'=: y NB. number of initial photons and scatterings Int=. k# 1 NB. Set Int=1 for the "k" initial photons z=. k#z_s NB. initially emitted photons at z=z_s P=. >0;0;z NB. starting x,y,z coordinates of all photons NB. PP=: >P;P NB. detailed scattering history C=. DCs k NB. get direction cosines of random directions Step^:(nscat) Int;z;C;P NB. do "nscat" escape/scatter cycles Sorted=: Sorted%(2*k) NB. normalize to total (+/+/+/Sorted)= 1. Sortup=: trim (hnmu+ i. hnmu+2){Sorted NB. seperate upward and downward streams Sortdn=: trim |. (i. hnmu+2){Sorted NB. trim empty borders & reflect in x-axis mubin=. (hnmu+ i. hnmu){ (%nmu)+ }: fr1 NB. the value of mu at the bin centers xbins=. ({.fr2),(BIN fr2),{:fr2 NB. x-coordinates of bin centers + ends ybins=. xbins NB. after reflection, y-coordinates same as x NB. PP=: }. PP mubin;xbins;ybins;Sortup;Sortdn ) NB. One cycle of escape and re-scatter. Step=: monad define 'Int z C P'=. y Theta=. acos(gHG iAHG rand k) NB. get random sample of H-G scattering angles phi=. 2p1*rand k NB. random set of phi's mu=. 2{C=. C Swerve Theta,:phi NB. find new direction cosines after scattering dst=. (T2%|mu)-(z%mu) NB. dst= optical depth to edge; extn= extinction Pout=. P+ dst *"1 C NB. coordinates of exiting photon's escape point Rout=. L }: Pout NB. radial distance from z-axis of escape point phi_out=. atan2 |: }: Pout NB. angular coordinate of escape point phi_ray=. atan2 |: }: C NB. angular coordinate of escape trajectory rphi=. 1p1-(|1p1-(|phi_ray- phi_out)) NB. angle between position & trajectory xout=. Rout* cos rphi yout=. Rout* sin rphi Iesc=. Int* extn=. ^-dst NB. flux of escaping photons Sp=. (fr1;fr2;fr3) Sort mu;xout;yout;Iesc Sorted=: Sorted+ Sp NB. and add to the earlier escapes Int=. Int*(1-extn) NB. intensities of scattered fraction remain=: k%~ +/Int NB. integrated intensity still in slab (global) Dtau=. Dt (1-extn)*rand k NB. travel distance of scattered (fractional) photons z=. 2{ P=. P+ Dtau *"1 C NB. new position of next scattering NB. PP=: PP,"3 P NB. add new scattering points to history Int;z;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=: dyad define 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=: dyad define 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. Get array of direction cosines from array of (mu,phi) angles NB. Angles in radians, 0= }:y idx=. <"1 ~. b (b +//. >{:y ) idx} 0$~>:#&>x ) rand=: ?@# 0: NB. rand n -- provides "n" random numbers on [0,1]. Dt=: [: -[:^. 1-] NB. "Dt rand n" maps [0,1] into travel on 0:%] NB. int01 n --> divides [0 ... 1] into n intervals. atan2=: 13 : '2p1 | ((0>1{"1 y)*1p1)+2p1+_3 o. %/"1 y' BIN =: -:@(}. + }:) NB. midpoints within series 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): cross=: ((1: |.[)*(_1: |. ]))-((_1: |.[)*(1:|.])) NB. Multiply quaternions: q1 QM q2 --> q3 QM=: dyad define 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=: dyad define s=. cos -: y f=. sin -: y s, (f*"1 x) ) NB. The conjuate of the quaternions: Qb q --> q_bar Qb=: monad : '(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=: dyad : '}. x QM (0,y) QM Qb x' Pr=: [: }. [ QM (0"_ , ]) QM [: Qb [ trim=: monad define NB. fix edges & reflect to get whole image a=. |: }: }. y NB. trim off mu<0 & mu>1 and transpose a=. }. (+/><0 1{a) 1}a NB. add 0th (y=0) bin to 1st & trim off 0th |: (|.a),a NB. join to reflection in x-axis, undo transpose ) NB. Run "HG_slab" ncycle times: HG_CYC=: monad define 'T z_s gHG hnmu nxy dxy nphoton nscat ncycle'=. y z=. (T;z_s;gHG;hnmu;nxy;dxy) HG_slab nphoton;nscat 'mubin xbins ybins SSortup SSortdn'=. z count=. 1 while. count< ncycle do. z=. (T;z_s;gHG;hnmu;nxy;dxy) HG_slab nphoton;nscat 'mubin xbins ybins Sortup Sortdn'=. z SSortup =. SSortup + Sortup SSortdn =. SSortdn + Sortdn count=. >: count ((": count),LF) (1!:2)wu;wd NB. intensity as function of mu NB. remain NB. 0.000194846 NB. remain+ +/wu+wd NB. 1 NB. check on total flux conservation