Reorganise the matrix scheme

This is a bit of reorganising of @Higginbottom 's rate matrix implementation, as discussed also with @kslong 

* rename matrix_solv to matrix_ion_populations
* move code which populates the rate matrix with rates into a subroutine populate_ion_rate_matrix
    * populate_ion_rate_matrix(xplasma, rate_matrix, pi_rates, rr_rates, b_temp, xne, xelem)
    * this populates the array rate_matrix with rates rr_rates and pi_rates in cell xplasma
    * it also uses the di and dr rates compited
* solve_matrix is a new, general routine, which solves the equation A x = b for x
    * it places the solution in an array populations
    * it takes a_data and b_data arrays as 1d versions of the matrix A and vector b
    * this routine is now also used by macro_pops

This is merging from a branch in my forked repository into the matrix_test branch in the main dev branch. I've verified that this produces identical results in the matrix ionization scheme, but need to test thoroughly for macro-atoms.