NB. This code evaluates the Chandrasekhar H-funcion.
load '/your_path/splines.ijs'
NB. Usage: H =. omega H_funct mu 
NB. "mu" is the array of values for which we want H(mu).
NB. Here, omega is the single-scattering albedo
H_funct=: dyad define
om=: x
C88=: %:1- om
NB. Define array of 48 mu values, emphasis near mu=0:
mu=. (0.0005*10^(10%~ i. 26)),0.2+0.05*i. 17
mu=. 0,0.00007,0.00017,0.00028,0.0004, mu
H=. mu& H_step ^:_ (#mu)#2  NB. iterate to convergence
NB. Interpolate H(mu) grid for value(s) of mu requested:
(mu;H) splint y
)
NB. Do the integrals over mu' :
H_step=: dyad define
HH=. |: y* x%x +/x         NB. HH is a 48 x 48 array
NB. we use "sintegA" to get all integrals in one step:
%C88+ (om%2)* x sintegA HH
)

NB. Here we are using the equation (TeX notation):
NB. H(\mu)~=~\left[\sqrt{1-\omega}~+~\frac{\omega}{2}~
NB. \int_0^1 H(\mu')\frac{\mu' d\mu'}{\mu'+\mu}\right]^{-1}
NB. to iteratively find H(mu) by integrating over H(mu').
NB. Thus for 48 values of mu, we need to evaluate 48
NB. integrals, each over 48 values of mu', for each
NB. iteration. We do the integrals based on cubic spline
NB. fits to the integrands.
NB. This routine does not seem to converge for values of
NB. of omega *very* close to 1. It will converge (rather  
NB. slowly) for omega=0.999999 (& all smaller values).
 
NB. For a thick, isotropically scattering atmosphere, with 
NB. a single-scattering albedo "omega" and radiation incident 
NB. at an angle whose cosine is mu0, this bit of code gives the
NB. intensity scattered into the direction(s) defined by direction
NB. cosine(s) mu. The equation is:  I(\mu)~=~\frac{F}{4}~\omega~
NB. \frac{\mu_0}{\mu+\mu_0}~H(\mu) H(\mu_0) taking F = 1.
NB. Usage: (omega;mu0) Iso_beam mu --> I(mu0,mu)
Iso_beam=: dyad define
'omega mu0'=. x
mu=. y
F=. 4%~ omega*mu0%mu0+ mu
F* (omega H_funct mu0)* omega H_funct mu
)