*** Perturbations to the neutrino distribution function (massless) ***
May 13, 2011: E.Komatsu
Ref: Shoji & Komatsu, arXiv:1003.0942

The program, "massless.f90", solves the Boltzmann equations describing
the evolution of MASSLESS neutrino perturbations numerically. It also
evaluates the analytical solutions. 

Specifically, it solves Eq.[18-20] in Shoji&Komatsu (2010) using the 
Runge-Kutta method. The analytical solution is given by Eq.[B12].

- To compile and use the sample program (given below), edit Makefile
and simply "./make"
- It will generate an executable called "massless"

The output files contain:

 - solution_psi0.txt: l=0, needed for the density perturbation
 - solution_psi1.txt: l=1, needed for the velocity perturbation
 - solution_psi2.txt: l=2, needed for the anisotropic stress

with the following format:

- 1st column: x=k*tau (where tau is the conformal time)
- 2nd column: numerical solution with lmax=100
- 3rd column: exact solution 

We have made the following assumptions:
- Newtonian gauge.
- Spatially flat and matter-dominated universe
  1. Scale factor is proportional to tau^2, where tau is the 
     conformal time.
  2. Gravitational potential is constant and set to unity.
- Adiabatic initial conditions (Ma&Bertschinger 1995, Eq.[98])