The Equation
I have this equation:
$$\frac{1}{2} \left( x^2+y^2 \right) + \frac{2(1-m)}{r_1} + \frac{2m}{r_2} = C,$$
where $m$ and $C$ are constants and $x,y$ are two variables and where
$$r_1 = \sqrt {(x-x_1)^2+y^2},\\ r_2 = \sqrt {(x-x_2)^2+y^2}.$$
with $x_1$ and $x_2$ being constants, i.e. two fixed values of the variable $x$.
My Question
Is it possible to solve this equation analytically; if not, how to solve it numerically?
I know that, if there are two variables, one needs two different equations to solve them. Thus, I can only aim to find a set of solutions.
Basically my intention is to find all those points in the $x$–$y$ plane such that the LHS of this equation becomes the RHS, i.e. constant $C$.
Background
This equation is from Forest Ray Moulton classic and great book on celestial mechanics, which was published in 1914. That time there were no computers. Still author has drawn the contours (i.e. solution of this equation) in his book! Its really amazing. My second intention behind this post is that how Moulton did this?
