Solving the equation for 'r' as a function of 'time' is tricky because one really does need an off-the-shelf computational 'widget' such as Brent's Method. Enjo's widget does the job well here, I'm sure.
It is easier, though, to drive the equation in reverse: instead of solving for 'r' as a function of 'time', we can easily solve for 'time' as function of 'r'. In fact, as written, that's exactly what the equation does. The left-hand-side is just 'time'; and the right-hand-side is just a function of 'r' - albeit a messy one.
For example, let's suppose that one wants to calculate the time it will take after the retrograde burn to hit the ground. And let's use the same scenario as before - i.e., starting a 20 x 20 km circular orbit around the moon, the craft executes a periapsis lowering retrograde burn that lowers the orbital periapsis altitude to -20 km. When does the spacecraft impact the surface?
To answer this question, we basically just need to plug the mean radius of the Moon is the right-hand-side as 'r', crank the handle, and calculate the answer. But let's break the process down into its component tasks. First, we note that the equation can be written as:
[MATH]\Delta \tau =A \times\left(\pi -2 \left(B(r)-C(r)\right)\right)[/MATH]
where
[MATH]A = \sqrt{\frac{\left(r_a+r_p\right){}^3}{8 \mu }}[/MATH]
[MATH]B(r) = \tan ^{-1}\left(\sqrt{\frac{r-r_p}{r_a-r}}\right)[/MATH]
[MATH]C(r) = \frac{\sqrt{\left(r_a-r\right) \left(r-r_p\right)}}{r_a+r_p}[/MATH]
Here, 'A' is just a constant. For the specific example that we are using with [MATH]r_a=1758100\,m[/MATH]; [MATH]r_p=1718100\,m[/MATH]; and [MATH]\mu=4.90487\times 10^{12}[/MATH] so that [MATH]A=1034.66[/MATH]. You might recognise that 'A' is closely related to the orbital period such that [MATH]T=2\,\pi\,A[/MATH] and all that the factor 'A' is doing here is converting the rest of the right-hand-side to units of time - i.e., its a conversion factor.
The rest of the right-hand-side - the [MATH]\pi -2 \left(B(r)-C(r)\right)[/MATH] bit - is just Kepler's Equation in disguise. Rather than deal with the eccentric anomaly and orbital eccentricity, it has just converted everything to expressions that just involve [MATH]r_a[/MATH], [MATH]r_p[/MATH] and 'r' - which are the things that we probably already know.
[MATH]B(r)[/MATH] and [MATH]C(r)[/MATH] are functions of the transfer radius, 'r'. If we want to know the time to impact, then, we just have to plug in values of [MATH]r_a[/MATH], [MATH]r_p[/MATH] and [MATH]r[/MATH]. In this case, the value of 'r' that we want to use is [MATH]1738100\,m[/MATH] - i.e., the Moon's mean radius. To calculate [MATH]B(r)[/MATH] we need to use the inverse tan function. This is found on most calculators and spreadsheets. One rarely calculates it 'by hand'. Sometimes it is called 'atan' or 'arctan'.
Anyway, if we plug in the values into the expression for [MATH]B(r)[/MATH], we get the value [MATH]0.785398[/MATH]. And if we plug in the values for [MATH]C(r)[/MATH] we get the value [MATH]0.00575341[/MATH]. So, then, the value of [MATH]\pi -2 \left(B(r)-C(r)\right)[/MATH] is [MATH]3.14159 - 2 \times (0.785398 -0.00575341) = 1.5823[/MATH]. And, so, to calculate [MATH]\Delta \tau[/MATH], the time to impact with the surface, we multiply this number by the value of 'A' (i.e., 1034.66) and get: [MATH]\Delta \tau = 1034.66 \times 1.5823 = 1637.14[/MATH]. So, the time to impact after the periapsis lowering burn is 1,637.14 seconds.
Now, there is nothing special about the value 'r' that we have chosen We could have chosen any other value of 'r' and thereby calculated the time taken to pass through the corresponding orbital altitude. The inversion of this process, calculating 'r' as a function pf 'time', is the tricky bit. Basically a root finding algorithm is needed and all these algorithms work in more or less the same way. They take an initial guess at 'r' and work out the corresponding value for [MATH]\Delta \tau[/MATH]. If that value doesn't equal the target value of [MATH]\Delta \tau[/MATH] - e.g. 1000 seconds, say - then the algorithm makes an intelligent update of the guess for 'r' and tries again. And the algorithm keeps on going through this process until it guesses right. Rather than invent their own guessing mechanism, most people most of the time just use other people's previously developed guessing mechanisms - hence Enjo's suggestion that you might want to try his which is based on a popular method, Brent's Method.
Feel free to ask questions about this. Everything is open for discussion.
---------- Post added at 10:11 AM ---------- Previous post was at 03:26 AM ----------
Actually, there is an easier way of calculating 'r' as a function of 'time'. The technique is to calculate a series expansion for 'r' in terms of [MATH]\Delta \tau[/MATH]. Such that:
[MATH]r = \alpha_0 + \alpha_2\,\Delta \tau^2 + \alpha_4\,\Delta \tau^4 + \alpha_6\,\Delta \tau^6 + \alpha_8\,\Delta \tau^8 + \alpha_{10}\,\Delta \tau^{10} + \dots[/MATH]
If one goes through the maths (and it helps to use a tool such as Mathematica or Matlab here), one can show that [MATH]\alpha_0[/MATH] is just your starting altitude when you execute the retrograde burn to lower periapsis and:
[MATH]\alpha_2 = \frac{\mu \left(-r_a+r_p\right)}{2 \,r_a^2 \left(r_a+r_p\right)} [/MATH]
[MATH]\alpha_4 = -\frac{\mu ^2 \left(r_a-2 \,r_p\right) \left(r_a-r_p\right)}{12 \,r_a^5
\left(r_a+r_p\right){}^2}[/MATH]
[MATH]\alpha_6 = -\frac{\mu ^3 \left(r_a-r_p\right) \left(11 \,r_a^2-44 \,r_a r_p+35 \,r_p^2\right)}{360 \,r_a^8
\left(r_a+r_p\right){}^3}[/MATH]
[MATH]\alpha_8 = -\frac{\mu ^4 \left(r_a-r_p\right) \left(73 \,r_a^3-438 \,r_a^2 r_p+714 \,r_a r_p^2-350
\,r_p^3\right)}{5040 \,r_a^{11} \left(r_a+r_p\right){}^4}[/MATH]
[MATH]\alpha_{10} = -\frac{\mu ^5 \left(r_a-r_p\right) \left(3548\, r_a^4-28384 \,r_a^3 r_p+70653 \,r_a^2
r_p^2-70840 \,r_a r_p^3+25025\, r_p^4\right)}{453600 \,r_a^{14} \left(r_a+r_p\right){}^5}[/MATH]
These 5 easy to calculate terms approximate 'r' to millimetre accuracy for the example used earlier. Additional terms can be calculated, but I doubt that they would be needed in most circumstances.
