NB. Isotropic scattering from a point source at x=y=0, z=z_s of slab.
NB. Usage: (T; z_s; hnmu; nxy; dxy) PaS k; nscat
NB. Here, T=slab thickness, z_s= location of source along z-axis
NB. hnmu= number of (positive) mu bins, nxy= number of positive x and y 
NB. bins, total array will be (2*nxy x 2*nxy), dxy= size of x & y bins
NB. k= number of Monte Carlo scattering photons, scatter "nscat" times
PaS=: 4 : 0
'T z_s 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 [-1,1] range of mu
fr2=: -:nx*dxy*(_1+2*int01 nx=. +:nxy) NB. divide up x-coord 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=0           
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    
aStep^:(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.
aStep=: 3 : 0      
'Int z C P'=. y
mu=. (_1+2*rand k)              NB. random set of mu's
r=. %:1- *:mu                   NB. sin(theta)
phi=. 2p1*rand k                NB. random set of phi's
C=. >|.(mu;(r*sin);(r*cos))phi  NB. direction cosines
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. Get array of direction cosines from array of (mu,phi) angles
NB. Angles in radians, 0=<phi<2*pi and -1<mu<1
DCs=: 3 : 0                     
dc=. _1+ 2*rand y  NB. random mu= cos theta = z direction cosine
st=. %: 1- *:dc    NB. sin theta(s) 
xc=. st*cos phi=. 2p1* rand y        
yc=. st*sin phi 
xc,yc,:dc
)

NB. Given the array A, which is a function of 3 dimensions, sort it into  
NB. bins where the bin edges are provided in the arrays fr1,fr2 & fr3, 
NB. and sum the contents of each bin. The arrays a1, a2, and a3 are the  
NB. coordinates of A that we are sorting on. The output Sp has dimensions 
NB. (1+#fr1,1+#fr2,1+#fr3). Usage: Sp=.  (fr1;fr2;fr3) Sort a1;a2;a3;A 
NB. (Thanks to Raul Miller and Roger Hui.)
Sort=: 4 : 0
b=. |:x I.&> }:y
idx=. <"1 ~. b 
(b +//. >{:y ) idx} 0$~>:#&>x        
)

trim=: 3 : 0  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 1st (y=0) bin to 2nd & trim off 1st
|: (|.a),a             NB. join to reflection in x-axis, undo transpose
)

L=: +/ &. (*:"_)        NB. Gives the lengths of vector arrays
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<tau<inf.
int01=: i.@>:%]         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. Run "PaS" ncycle times:
PaS_CYC=: 3 : 0
'T z_s hnmu nxy dxy nphoton nscat ncycle'=. y
z=. (T;z_s;hnmu;nxy;dxy) PaS nphoton;nscat 
'mubin xbins ybins SSortup SSortdn'=. z
count=. 1
while. count< ncycle do.
  z=. (T;z_s;hnmu;nxy;dxy) PaS nphoton;nscat 
  'mubin xbins ybins Sortup Sortdn'=. z
  SSortup =. SSortup  + Sortup   
  SSortdn =. SSortdn  + Sortdn   
  count=. >: count
  ((": count),LF) (1!:2)<Cycle_log
end.
SSortup=. count%~ SSortup NB. divide by number of cycles
SSortdn=. count%~ SSortdn 
mubin;xbins;ybins;SSortup;SSortdn
)
Cycle_log=: '/your_path/PaS_CYC.log'

NB.  An example consider a slab of total thickness T=3 with the source
NB. located at z=1 (0.5 below the upper surface). Lets scatter 40*4e4
NB. = 16e6 photons 35 times each, looking at a 80x80 grid of the surface
NB. with a width of (-5.6,5.6)=11.2 in tau. Let's consider 30 angles for
NB. the emerging radiation. Thus we will run:
NB.    z=. PaS_CYC 3;1;30;40;0.14;4e5;35;40
NB.    'mubin xbins ybins Su Sd'=. z
NB. To look at the results, consider mu index 21, which we see corresponds
NB. to an angle with respect to the perpendicular to the slab of
NB.    (180%1p1)*acos 21{mubin
NB. 44.22   
NB. degrees. We can then plot the log of the emergent radiation:
NB.    'surface' plot 10^. 1e_6+ 21{Su
NB.    'surface' plot 10^. 1e_6+ 21{Sd
NB. If we look at this as a contour plot:
NB.    'contour' plot xbins;ybins; 10^. 1e_6+ 21{Sd
NB. We see sharp peak at y=0, x=2.4. This due to the radiation that escapes
NB. directly from the source without scttering. The source is (T/2)+z_s =
NB. 1.5+1=2.5 units above the bottom surface, so from this angle, we see it
NB. displaced 
NB.    2.5*tan acos 21{mubin
NB. 2.43258 units. The scattered radiation, on the other hand, is shifted 
NB. toward the subsource point of x=y=0. If we average over the entire surface
NB.    Au=. +/"1 +/"2 Su
NB.    Ad=. +/"1 +'"2 Sd
NB. we can see the different I(mu) laws for the two surfaces.
NB.    'pensize 2' plot mubin;>Au;Ad  
NB. This looks different from the I(mu) plots for for a mid-plane source from
NB. the previous code, since the source is not an infinite plane, but a point.
NB. To recover the I(mu) law for a source at z=z_s but infinite in x and y, we
NB. need only divide the results by mu to account for projection:
NB.    'pensize 2' plot mubin;>(Au%mubin);(Ad%mubin)
NB. The lower curve is of course the radiation out the bottom. The total flux
NB. from top and bottom are
NB.    +/&.> Au; Ad; Au+Ad
NB. +--------+--------+--------+
NB. |0.730642|0.264995|0.995637|
NB. +--------+--------+--------+
NB. but about 0.4% have still not escaped after 35 scatterings:
NB.    remain, remain+ +/Au+Ad
NB. 0.00435893 0.999996
NB. where the last number should be unity and indicated the error.