NB. Calculate the M_ij matrix, Equation (18) 
NB. Elitzur & Ramos (2006) MNRAS 365, 779.
NB. Usage: MMs tau (tau is array tau_i of boundary tau's).
MMs=: monad define
z=. - -:Alphas |y-/y  NB. |y-/y -> 2D array of delta tau's
(}.}."1 z)+(-}:}."1 z)+(-}.}:"1 z)+ }:}:"1 z
)

NB. read in data file:
'logt lalpha Dlal'=: readB '/your_path/log_alpha_table_161.dat'
load'/your_path/splines.ijs'
NB. calculate alpha's for array of tau's by spline interpolation.
Alphas=: monad : '10^ (logt;lalpha;Dlal) spltrp 10^. y'

NB. This is a good example of the compact power of J. The whole program to 
NB. compute the matrix of M_ij values can be written as a one line of J code!
NB. MMs=: monad : '(}.}."1 z)+(-}:}."1 z)+(-}.}:"1 z)+ }:}:"1 z=. - -:Alphas |y-/y'