NB. Isotropic scattering from a point source at x=y=z=0 of slab.
NB. Usage: (T; hnmu; nxy; dxy) PiS k; nscat
NB. Here, T=slab thickness, hnmu= number of (positive) mu bins
NB. nxy= number of positive x-axis and y-axis bins, dxy= size of xy bins
NB. k= number of Monte Carlo scattering photons, scatter "nscat" times
PiS=: 4 : 0
'T 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=: int01 hnmu                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: positive only by symmetry
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#0                         NB. initially emitted photons at z=0           
P=. >0;0;z                      NB. starting x,y,z coordinates of all photons 
NB. if run with small k, PP can store scattering points for each photon 
NB. PP=: >P;P                   NB. un-NB. to capture all scattering points    
C=. DCs k                       NB. get direction cosines of random directions    
Step^:(nscat) Int;z;C;P         NB. do "nscat" escape/scatter cycles      
mubin=. (%nmu)+ }: fr1          NB. the value of mu at the bin centers
ybins=. xbins=.(-:dxy)+ }:fr2   NB. x & y coords of bin centers (reflection makes y same)
Sorted=: Sorted% 2*k            NB. normalize so +/+/+/Sorted = nmu (less "remain")
Sorted=: trim Sorted            NB. trim empty borders & reflect in x-axis
NB.  PP=: }. PP
mubin;xbins;ybins;Sorted
)

NB. One cycle of escape and re-scatter.
Step=: 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'

NB. Run "PiS" ncycle times:
PiS_CYC=: 3 : 0
'T hnmu nxy dxy nphoton nscat ncycle'=. y
z=. (T;hnmu;nxy;dxy) PiS nphoton;nscat 
'mubin xbins ybins SSorted'=. z
count=. 1
while. count< ncycle do.
  z=. (T;hnmu;nxy;dxy) PiS nphoton;nscat 
  'mubin xbins ybins Sorted'=. z
  SSorted =. SSorted  + Sorted   
  count=. >: count
  ((": count),LF) (1!:2)<Cycle_log
end.
SSorted=. count%~ SSorted NB. divide by number of cycles
mubin;xbins;ybins;SSorted
)
Cycle_log=: '/your_path/PiS_CYC.log'

NB.  An example run
NB.  z=. PiS_CYC 2;20;20;0.25;2e5;30;20
NB.  'mubin xbins ybins SS'=. z
NB.  (180%1p1)*acos 8{mubin
NB. 64.8493
NB. index mubin=8 holds radiation emerging near 65 deg wrt normal to slab
NB. 'surface' plot 8{SS  NB. look radiation leaving surface of slab
NB. 'surface' plot 10^. 0.001+ 8{SS   NB. look at log to see fainter regions
NB. We may plot contours - see how peak is shifted along x-axis, which is
NB. the axis along which the observer is looking:
NB. 'contour; pensize 2' plot xbins;xbins;10^. 0.001+ 8{SS
NB. mubin=19 --> 13 deg, looking almost straight down:
NB. 'contour; pensize 2' plot xbins;xbins;10^. 0.001+ 19{SS
NB. mubin=1 --> 85.7 deg, line of sight just skims surgace:
NB. 'contour; pensize 2' plot xbins;xbins;10^. 0.001+ 1{SS
NB. 'surface; viewpoint _1.7 2.2 0.5' plot %: 1{SS
NB. (square root shows low levels better, but not as extreme as log)