Drawing elliptical orbit in Python (using numpy, matplotlib)

I wonder how can I draw elliptical orbit by using the equation ay² + bxy + cx + dy + e = x² ?

I have first determined the a,b,c,d,e constants and now I assume by giving x values I will obtain y and this will give me the graph I want but I couldn't do it by using matplotlib.

I would really appreciate, if you could help me!

EDIT: I added the code here.

from numpy import linalg
from numpy import linspace
import numpy as np
from numpy import meshgrid
import random
import matplotlib.pyplot as plt
from scipy import optimize

x = [1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.30, 0.16, 0.01]
y = [0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.12, 0.15]

my_list = [] #It is the main list.
b = [0] * len(x) # That is the list that contains the results that are given as x^2 from the equation.

def fxn():  # That is the function that solves the given equation to find each parameter.
    global my_list
    global b
    for z in range(len(x)):
        w = [0] * 5
        w[0] = y[z] ** 2
        w[1] = x[z] * y[z]
        w[2] = x[z]
        w[3] = y[z]
        w[4] = 1
        my_list.append(w)
        b[z] = x[z] ** 2

    t = linalg.lstsq(my_list, b)[0]
    print 'List of list representation is', my_list
    print 'x^2, the result of the given equation is', b
    print '\nThe list that contains the parameters is', t

fxn()
t = linalg.lstsq(my_list, b)[0]

print '\nThe constant a is', t[0]
print 'The constant b is', t[1]
print 'The constant c is', t[2]
print 'The constant d is', t[3]
print 'The constant e is', t[4]

EDIT: Here are the constant values:

a = -4.10267300566
b = 1.10642410023
c = 0.39735696603
d = 3.05101004127
e = -0.370426134994

The problem can be solved for y as a function of x

The catch is that there are 2 values of y for every valid x, and no (or imaginary) y solutions outside the range of x the ellipse spans

below is 3.5 code, sympy 1.0 should be fine but print, list comps may not be backwards compatable to 2.x

from numpy import linalg
from numpy import linspace
import numpy as np
from numpy import meshgrid
import random
import matplotlib.pyplot as plt
from scipy import optimize
from sympy import *

xs = [1.02, 0.95, 0.87, 0.77, 0.67, 0.56, 0.44, 0.30, 0.16, 0.01]
ys = [0.39, 0.32, 0.27, 0.22, 0.18, 0.15, 0.13, 0.12, 0.12, 0.15]

b = [i ** 2 for i in xs] # That is the list that contains the results that are given as x^2 from the equation.

def fxn(x, y):  # That is the function that solves the given equation to find each parameter.
    my_list = [] #It is the main list.
    for z in range(len(x)):
        w = [0] * 5
        w[0] = y[z] ** 2
        w[1] = x[z] * y[z]
        w[2] = x[z]
        w[3] = y[z]
        w[4] = 1
        my_list.append(w)
    return my_list

t = linalg.lstsq(fxn(xs, ys), b)


def ysolv(coeffs):
    x,y,a,b,c,d,e = symbols('x y a b c d e')
    ellipse = a*y**2 + b*x*y + c*x + d*y + e - x**2
    y_sols = solve(ellipse, y)
    print(*y_sols, sep='\n')

    num_coefs = [(a, f) for a, f in (zip([a,b,c,d,e], coeffs))]
    y_solsf0 = y_sols[0].subs(num_coefs)
    y_solsf1 = y_sols[1].subs(num_coefs)

    f0 = lambdify([x], y_solsf0)
    f1 = lambdify([x], y_solsf1)
    return f0, f1

f0, f1 = ysolv(t[0])

y0 = [f0(x) for x in xs]
y1 = [f1(x) for x in xs]

plt.scatter(xs, ys)
plt.scatter(xs, y0, s=100, color = 'red', marker='+')
plt.scatter(xs, y1, s=100, color = 'green', marker='+')
plt.show()  

when the above is ran in Spyder:

runfile('C:/Users/john/mypy/mySE_answers/ellipse.py', wdir='C:/Users/john/mypy/mySE_answers')
(-b*x - d + sqrt(-4*a*c*x - 4*a*e + 4*a*x**2 + b**2*x**2 + 2*b*d*x + d**2))/(2*a)
-(b*x + d + sqrt(-4*a*c*x - 4*a*e + 4*a*x**2 + b**2*x**2 + 2*b*d*x + d**2))/(2*a)

 The generated functions for the y values aren't valid everywhere:

f0(0.1), f1(0.1)
Out[5]: (0.12952825130864626, 0.6411040771593166)

f0(2)
Traceback (most recent call last):

  File "<ipython-input-6-9ce260237dcd>", line 1, in <module>
    f0(2)

  File "<string>", line 1, in <lambda>

ValueError: math domain error


In [7]:

The domain error would require a try/execpt to "feel out" the valid x range or some more math

like the try/except below: (Edited to "close" drawing re comment )

def feeloutXrange(f, midx, endx):
    fxs = []
    x = midx
    while True:
        try: f(x)
        except:
            break
        fxs.append(x)
        x += (endx - midx)/100
    return fxs

midx = (min(xs) + max(xs))/2    

xpos = feeloutXrange(f0, midx, max(xs))
xnegs = feeloutXrange(f0, midx, min(xs))
xs_ellipse = xnegs[::-1] + xpos[1:]

y0s = [f0(x) for x in xs_ellipse]
y1s = [f1(x) for x in xs_ellipse]

ys_ellipse = y0s + y1s[::-1] + [y0s[0]] # add y start point to end to close drawing

xs_ellipse = xs_ellipse + xs_ellipse[::-1] + [xs_ellipse[0]] # added x start point


plt.scatter(xs, ys)
plt.scatter(xs, y0, s=100, color = 'red', marker='+')
plt.scatter(xs, y1, s=100, color = 'green', marker='+')
plt.plot(xs_ellipse, ys_ellipse)
plt.show()

Edit: added duplicate start point to the end of ellipse point lists to close the plot figure

ys_ellipse = y0s + y1s[::-1] + [y0s[0]] # add y start point to end to close drawing

xs_ellipse = xs_ellipse + xs_ellipse[::-1] + [xs_ellipse[0]] # added x start point