Computing the distance between two points in polar coordinates using floating-point

Let b = (a1-a2)/2 
then using
cos( a1-a2) = 1 - 2*sin(b)*sin(b)
D = sqrt( (r1-r2)*(r1-r2) + 4*r1*r2*sin(b)*sin(b))

This at least gets rid of square roots of negative numbers but would still problems with overflow.

Perhaps

x = 2*sqrt(r1)*sqrt(r2)*sin(b)
D = hypot( r1-r2, x)

would resolve that?