> restart:
a:=1:

> theta:=(b,r)->Re(arccos((b^4-4*a*a*r*r-r**4)/(4*a*r**3))):
B:=NULL:
for b in [seq(0.1*k,k=1..10)] do
B:=B,polarplot([[r,theta(b,r),r=a-sqrt(a**2-b**2)..sqrt(a**2+b**2)-a],[r,-theta(b,r),r=a-sqrt(a**2-b**2)..sqrt(a**2+b**2)-a],[r,theta(b,r),r=a+sqrt(a**2-b**2)..sqrt(a**2+b**2)+a],[r,-theta(b,r),r=a+sqrt(a**2-b**2)..sqrt(a**2+b**2)+a]],color=red,numpoints=200)
od:
for b in [seq(1+0.05*k,k=1..10)] do
B:=B,polarplot([[r,theta(b,r),r=sqrt(a**2+b**2)-a..sqrt(a**2+b**2)+a],[r,-theta(b,r),r=-a+sqrt(a**2+b**2)..sqrt(a**2+b**2)+a]],color=red,numpoints=200)
od:

> display([plot([[-2*a,0],[0,0]],style=point,symbol=circle),B],scaling=constrained,axes=none);