Implement recombination calculation

I've implemented a function to return x_e and T_m as a function of z (there are options for the z range because NDSolve isn't too happy outside the range ~100 < z < ~1800, this range is good enough for z_drag calculations).

Implemented tau_b calculation.

Added helium fraction to cosmological parameters.

Added photon density function.

You should be able to test this out by trying:

createCosmology[default];
myXe = recombinationXe[default][[1]];
FindRoot[\[Tau]b[default, myXe, z] - 1, {z, 1050}]

Should give 

{z -> 1060.25}