I would like to discuss also my method for finding a DE, without derivatives (that's the method that Knighty followed in its script)
Too bad, it is a bit too complicated for some surfaces, because it's needed to know the signs of every "polynomial component" of the surface.
Let's examine a very very simple example.
x2 = x*x;
y2 = y*y + z*z;
And examine the "lathed" Scarabaeus curve expression;
(x2 + y2) (x2 + y2 + R1 x)^2 - R2^2 (x2 - y2)^2 = 0
The problem for DE is that the field needs to be "approximately linear", on x,y,z. For example a sphere equation is;
x*x + y*y + z*z - R*R = 0
To make it "linear", we can take square root. It does not exist for negative values? Oh it's simple, bring R*R to the right.
x*x + y*y + z*z = R*R
sqrt(x*x + y*y + z*z) - abs(R) = 0
This gives a very good sphere DE.
Okay, applying the same method to the scarab curve is not so plain; (x2 + y2) (x2 + y2 + R1 x)^2 is always positive (a sum of squares) but (x2 - y2)^2 has a variable sign. If we expand we have (x2 - y2)^2 = (x^4 + y^4) - 2*(x*x*y*y)
Again we can separate positive and negative, because x*x*y*y is always >0 ;
ssqrt((x2 + y2) (x2 + y2 + R1 x)^2 + 2*R2*R2* (x*x*y*y))) - sqrt(R2)*ssqrt(x^4+y^4) where ssqrt(a) is sqrt(sqrt(a)) to make the vector field of power 1 because the highest power is 4.
That seems to works like a charm!
The problem is that many surfaces have lots mixed terms like x*y, or even power terms like x, x^3 ... , so separate those correctly is a trouble.
In particular Enneper implicit equation...