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

The program, "massive.f90", solves the Boltzmann equations describing
the evolution of MASSIVE 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.[B18].

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

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])

Note added (Oct 7, 2015): The calculation of line 71 and the subroutine "source0" is the result of the recurrence relation of the spherical Bessel function [ l j_{l-1}(x) - (l+1) j_{l+1}(x) = (2l+1) j'_{l}(x) ] and integration by parts from eq. (B19) of arXiv:1003.0942. The calculation of line 73 and the subroutine "souce1" as well as line 75 and subroutine "source2" use the same computation. [We thank Chi-Ting Chiang for clarifying this.]