Skip to content

Really tiny bug fix in interp_table()

John Helly requested to merge interp_table_fix into master

I think I found a tiny bug (which probably has no effect on anything!) in interp_table() in cosmology.c:

  if (ii <= 1)
    return table[0] * xx;
  else
    return table[ii - 1] + (table[ii] - table[ii - 1]) * (xx - ii);

My interpretation of this is that table[0] is the value at the function at the end of the first interval in dloga and we have an implicit table[-1] which is assumed to be always zero. So returning table[0]*xx for ii=0 is correct.

For ii=1 the interpolation expression in the else reduces to

    table[0] + (table[1] - table[0]) * (xx - 1);

which looks like the right answer to me, but we're returning table[0]*xx instead, where 1 <= xx < 2. The effect of this is just that for ii=1 we're extrapolating from the first interval when we could interpolate over the second.

(I found this while adding a comoving distance function which didn't work because I wasn't aware of the assumption that the table starts at zero).

Merge request reports