Drawing Elliptical Orbit In Python (using Numpy, Matplotlib)

October 30, 2022 Post a Comment

I wonder how can I draw elliptical orbit by using the equation ay2 + bxy + cx + dy + e = x2 ? I have first determined the a,b,c,d,e constants and now I assume by giving x values I

Solution 1:

Intro

Easiest thing would be to parametrize this equation. As @Escualo suggests, you could introduce a variable t and parametrize x and y along that. Parametrizing means separating your equation into two separate equations for x and y individually in terms of t. So you would have x = f(t) and y = g(t) for some values of t. You could then plot the x, y pairs that result for each value of t.

The catch here is that your ellipse is rotated (the x*y coupling term is an indication of that). To separate the equations, you have to first transform the equation to get rid of the coupling term. This is the same as finding a set of axes that are rotated by the same angle as the ellipse, parametrzing along those axes, then rotating the result back. Check out this forum post for a general overview.

Math

You first need to find the angle of rotation of the ellipse's axes with respect to the x-y coordinate plane.

$\theta = \frac{1}{2} tan^{-1}(\frac{b}{-1 - a})$

Your equation then transforms to

$a'x'^2+b'y'^2+c'x+d'y'+e'=0$

where

$a'=-cos^2(\theta)+b*sin(\theta)cos(\theta)+a*sin^2(\theta)$

$b'=-sin^2(\theta)-b*sin(\theta)cos(\theta)+a*cos^2(\theta)$

$c'=c*cos(\theta)+d*sin(\theta)$

$d'=-c*sin(\theta)+d*cos(\theta)$

$e'=e$

To find the (nearly) standard form of the ellipse, you can complete the squares for the $x'$ and $y'$ portions and rearrange the equation slightly:

$\frac{(x'-h)^2}{a''^2}+\frac{(y'-k)^2}{b''^2}=c''$

where

$a''=\frac{1}{\sqrt{a'}}$

$b''=\frac{1}{\sqrt{b'}}$

$c''=\frac{c'^2}{4a'}+\frac{d'^2}{4b'}-e'$

$h=-\frac{c'}{2a'}$

$k=-\frac{d'}{2b'}$

Since you know $\theta$ , you can now parametrize the equations for $x'$ and $y'$ :

$x'=h+a''\sqrt{c''}sin(t)$

$y'=k+b''\sqrt{c''}cos(t)$

You would then rotate back into normal x-y space using the formulas

$x=x'cos(\theta)-y'sin(\theta)$

and

$y=x'sin(\theta)+y'cos(\theta)$

Code

The code to get the x- and y- arrays to pass to plt.plot is now relatively straightforward:

def computeEllipse(a, b, c, d, e):
    """
    Returns x-y arrays for ellipse coordinates.
    Equation is of the form a*y**2 + b*x*y + c*x + d*y + e = x**2
    """
    # Convert x**2 coordinate to +1
    a = -a
    b = -b
    c = -c
    d = -d
    e = -e
    # Rotation angle
    theta = 0.5 * math.atan(b / (1 - a))
    # Rotated equation
    sin = math.sin(theta)
    cos = math.cos(theta)
    aa = cos**2 + b * sin * cos + a * sin**2
    bb = sin**2 - b * cos * sin + a * cos**2
    cc = c * cos + d * sin
    dd = -c * sin + d * cos
    ee = e
    # Standard Form
    axMaj = 1 / math.sqrt(aa)
    axMin = 1 / math.sqrt(bb)
    scale = math.sqrt(cc**2 / (4 * aa) + dd**2 / (4 * bb) - ee)
    h = -cc / (2 * aa)
    k = -dd / (2 * bb)
    # Parametrized Equation
    t = np.linspace(0, 2 * math.pi, 1000)
    xx = h + axMaj * scale * np.sin(t)
    yy = k + axMin * scale * np.cos(t)
    # Un-rotated coordinates
    x = xx * cos - yy * sin
    y = xx * sin + yy * cos

    return x, y

To actually use the code:

from matplotlib import pyplot as plt

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

lines = plt.plot(*computeEllipse(a, b, c, d, e))

To overplot your original points on the ellipse:

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]
ax = lines[0].axes
ax.plot(x, y, 'r.')

The result is the following image:

plot

Note

Keep in mind that the forum post I linked to uses a different notation than the one you do. Their equation is Ax + Bxy + Cy + Dx + Ey + F = 0. This is a bit more standard than your form of ay + bxy - x + cx + dy + e = 0. All of the math is in terms of your notation.

Solution 2:

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

Solution 3:

The easiest thing to do is rewrite in parametric form so that you end up with the expressions x = A cos(t); y = B sin(t). You then create a full ellipse by assigning t = [0, 2 * pi] and calculating the corresponding x and y.

Read this article which explains how to move from a general quadratic form into a parametric form.

Python Playground