NB. Number of photons escaping from slab of optical depth     
NB. "tau", with point of emission ranging uniformly over
NB. tau. Usage: Alpha tau --> escapes 
Alpha=: monad define
if. y <: 0 do. 0 else.
 p=. Phi xx=. 5.5*int01 1e5  NB. integrate to 5.5 Doppler widths
 bb=. 1- +: 3&EIn y*p        NB. 1 - 2* E_3(tau*Phi(x))
 +/(BIN bb)*(DEL xx)         NB. 10^5 point trapezoidal rule integration 
end.
)
Phi=: monad : '(^-*:y)% %:1p1'  NB. phi(x)= [(1/sqrt(pi)] exp(-x^2)
NB. load the exponential integral function (see bottom of this J page):
load'/your_path/EIn.ijs'