Removed input array "Lambda_star" from argument list of "calc_Lambda_star".
Removed extra input array "source_fn" from "calc_source_fn."

I forgot it was only there for the space and was abusing it terribly.
ALI looks exactly the way it does in OAB86 and OK87, but still doesn't give the correct answer...
Changed matrix elements of Lstar from individual arrays to a full matrix.

This is much less efficient storage-wise, but it makes the construction
of the full Lambda operator (or ALOs of variable bandwidth) much easier.
Changed Lstar to be the full Lambda operator so we can see if the matrix elements are correct.

Result: they are not correct. Hmm...
Added parameter "n_ray_pts" to distinguish it from "n_depth_pts" for clarity.

In the plane-parallel case the two are equal, but in, e.g., spherical
symmetry they are not.
Fixed full Lambda operator (almost).

Now when I operate Lstar (which is actually the full Lambda operator) on
the source function, it looks like regular Lambda iteration, except that
the depth boundary condition is wonky. Need to fix that...
Fixed some missing "elif"s
Changed "value" of the Planck function from 1 to 5.

This is to distinguish it from things which really do need to be one
(unity).
FINALLY FIXED APPROXIMATE LAMBDA OPERATOR!

The boundary value at depth was screwing up. It turns out for
upward-pointing rays, the coefficient for the matrix element at depth
needs to be 1, not 0. (This is confusing to me, given that Olson &
Kunasz (1987) says to "set the incident boundary intensities to zero".)

I'll need to think a lot more about why this is the right answer. But at
least for now it works. So I should be able to strip out any bandwidth
of this matrix to get a correct Lambda_star.