Crap.

Made a couple of mistakes.

If
(D/2)^2 = a^2 - b^2

then
a = sqrt((D/2)^2 + b^2)

The polar equation for an ellipse with centre at the origin is:
r(theta) = ab/sqrt((b.cos(theta))^2 + (a.sin(theta))^2)

If we substitute b = sqrt(a^2 - c^2) into the polar equation
then it should read:
r(theta) = a.sqrt(a^2 - c^2)/sqrt((sqrt(a^2 - c^2).cos(theta))^2 +
(a.sin(theta)^2))

(b.cos(theta))^2 rather than b.cos(theta)^2 !

Another potential option is to use parametric equations.

x = a.cos(t)
y = b.sin(t)

for 0 <= t <= 2.pi

substitute sqrt((D/2)^2 + b^2) for a.

Rob O'Connor