NB. Isotropic scattering of source distributed iniformly in slab.
NB. Usage: (T; 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. since this is a lambda-iteration, we must set nscat >> T^2
Uslab=: 4 : 0
'T 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# 1                      NB. Set Int=1 for the "k" initial photons   
z=. T2- T*rand k                NB. initially emitted photons have uniform 
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=: 3 : 0      
'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=: 4 : 0      
'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.

NB. Example: slab of optical depth 2, scattering 2e6 
NB. photons 30 times each: 
NB.    'mubin Su Sd'=. (2;10) Uslab 2e6; 30
NB.    remain
NB. 0.000628216     NB. 0.06% still trapped after 30 scatterings 
NB.    load'plot'
NB.    plot mubin;>Su;Sd   NB. plot emergent intensities 
NB. they should be equal -- difference indicates statistical error
NB. Let's compare to the "exact" solution at the mu's of bin centers:
NB. load '/your_path/Uniform.ijs'  NB. get matrix transform routine
NB.    tau=. tau_grid2 0.001 10 1  NB. define a 63-point tau grid over [0,2].
NB.    'Ftot Itop'=. Uniform tau;mubin     NB. run code
NB.    plot mubin;>(1p1*Itop);((1-remain)%~(-:Su+Sd))  NB. plot comparison
NB. (1-remain)%~ corrects for the photons still trapped in the slab
NB. Also, we multiply the matrix transform result by pi: physical vs. 
NB. "astrophysical flux" (i.e., per steradian) used in matrix equations.

NB. We sholudn't make nphoton so large that memory is filled and there
NB. are page faults. Then to get better statistics, do many runs:
NB. Run "Uslab" ncycle times:
UCYC=: 3 : 0
'T hnmu nphoton nscat ncycle'=. y
z=. (T;hnmu) Uslab nphoton;nscat 
'mubin SSu SSd'=. z
count=. 1
while. count< ncycle do.
  'mubin Su Sd'=. (T;hnmu)Uslab nphoton;nscat
  SSu=. SSu + Su
  SSd=. SSd + Sd
  count=. >: count
  NB. Look at 'UCYC_count.log' to monitor progress:
  ((": count),LF) (1!:2)<Cycle_log
end.
SSu=. count%~ SSu      NB. divide by number of cycles
SSd=. count%~ SSd
mubin;SSu;SSd
)
Cycle_log=: '/your_path/UCYC_count.log'