NB. Isotropic scattering of a beam incident on a slab.
NB. Usage: (T; mu_i; hnmu) Uslab k; nscat
NB. 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
MC_beam=: dyad define
'T mu_i hnmu'=. x    
T2=. -: T                       NB. optical depth to midplane
nmu=. +: hnmu                   NB. total (positive and negative) mu bins
'k nscat'=. y                   NB. number of initial photons and scatterings
Int=. k# Fscat=: F T%mu_i       NB. Set Int to fraction of beam scattered
z=. T2- mu_i*Dt Fscat*rand k    NB. z of initially scattered photons
fret=: _1+ 2*int01 nmu          NB. frets to divide -1<mu<1 into nmu bins
Sorted=: nmu#0                  NB. placeholder for the escaping intensities
Step^:(nscat) T2;nmu;k;Int;z    NB. do "nscat" escape/scatter cycles      
mubin=. (%nmu)+ }: int01 hnmu   NB. the value of mu at the bin centers
Sorted=: (nmu%(4*k))*Sorted     NB. Scale so integral 2*Imu*mu*d mu = 1
Su=. (hnmu+ i. hnmu){Sorted     NB. intensities escaping upward (positive mu)
mubin;Su;Sd=.|.(i.hnmu){Sorted  NB. and Sd = escaping downward (negative mu's)
)

NB. One cycle of escape and re-scatter.
Step=: monad define
'T2 nmu k Int z'=. y
mu=. (_1+2*rand k)              NB. random set of mu's
dst=. (T2%|mu)-(z%mu)           NB. dst= optical depth to edge; extn= extinction
Iout=. (|mu)%~Int*extn=. ^-dst  NB. divide I by |mu| to account for projection
Std=. fret Sort mu; Iout        NB. sort escaping Imu into nmu angle bins
Sorted=: Sorted+ Std            NB. add to previously escaped (global)
Int=. Int*(1-extn)              NB. intensities of scattered fraction
remain=: k%~ +/Int              NB. integrated intensity still in slab (global)
z=. z+ mu*Dt (1-extn)*rand k    NB. new z of scattered (fractional) photons 
T2;nmu;k;Int;z
)

NB. Given the array A, where each element is tagged by an element of a,
NB. sort A into bins, where the bin edges are provided in the array
NB. "fret", and sum the contents of each bin. We assume no elements of
NB. a are outside the limits of fret. The output Sp has dimension <:#fret.
NB. Usage: Sp=.  fret Sort a;A
Sort=: dyad define
'a A'=. y
b=. x I. a
(b +//. A) (<: ~.b)} 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<tau<inf.
F=: 1- [: ^-               NB. 1-exp{-y} = fraction scattered by distance y

NB. Example: slab of optical depth 1.5, scatter 2e6 photons
NB. 30 times each, beam incident at mu_i=0.5 (60 deg):: 
NB.    'mubin Su Sd'=. (1.5;0.5;10) MC_beam 4e6; 30
NB.    remain
NB. 4.32197e_5      NB. 4e-5 still trapped after 30 scatterings
NB.    Fscat
NB. 0.950213        NB. 95% scattered by slab, 5% pass through
NB.    load'plot'
NB.  'pensize 2' plot mubin;>Su;Sd  NB. plot emergent intensities