Yes, it's messy arc length calculus.
We can start with the relation between radius and true anomaly, which is rendered here in a form applicable to all conic section orbits by using the periapsis $p$:
$r=\dfrac{p(1+e)}{1+e\cos\theta}$
This is to be combined with the arc length differential
$(\dfrac{ds}{d\theta})^2=r^2+(\dfrac{dr}{d\theta})^2$
Direct substitution leads to
$\dfrac{ds}{d\theta}=\dfrac{p(1+e)\sqrt{(1+e\cos\theta)^2+e^2\sin^2\theta}}{(1+e\cos\theta)^2}$
$=\dfrac{p(1+e)\sqrt{1+e^2+2e\cos\theta}}{(1+e\cos\theta)^2}$
Unless $e=0$ (circular orbit) or $e=1$ (parabolic orbit), integrating this requires elliptic functions, and that is a bridge too far for me. More practical for your problem seems to be leaving the relation in differential form, updating the derivative with each timestep.