TOVSolver.solveTOV contains unexplained magic numbers

[rhaas80]

the function contains

    r = 4.441e-16 * unit.cm
    dr = 10.0 * unit.cm

without explaining why 4.441e-16 cm and 10cm are good initial locations and step sizes.

[tuhinbits]

Thank you for your review, @rhaas80. You're absolutely right about the unexplained values in our code.

We agree these were "magic numbers" without proper documentation. The value 4.441e-16 * unit.cm was chosen as a small enough initial radius to avoid the singularity at r=0 when starting the integration, while the 10.0 * unit.cm was simply an initial step size.

We've made the following improvements:

1. We now use a clearer value of 1e-18 * unit.cm as the initial radius in both solveTOV and solveTOV_tidal functions. This is small enough to avoid the singularity but not so small as to cause numerical issues.

2. We've kept the initial step size as 10.0 * unit.cm since it's just a starting point - the step size adapts dynamically inside the integration loop based on pressure and mass gradients.

3. We've added proper documentation to explain these choices in both functions:

# Start with a tiny but non-zero radius to avoid the singularity at r=0
r = 1e-18 * unit.cm
# Initial step size - will be adaptively adjusted during integration
dr = 10.0 * unit.cm

Our testing confirms that the integration is stable with these values across a wide range of equations of state. The adaptive step size algorithm quickly adjusts to an appropriate size based on the local gradients, making the initial value less critical as long as it's reasonably small.

Let us know if you have any other concerns or suggestions for further improvements!