import numpy as np
import matplotlib.pyplot as plt
import math

def normal(v):
    n = np.array([8 * v[0], 2 * v[1]])
    return n / np.linalg.norm(n)

def plot_normal(v):
    u = v - normal(v)
    plt.plot([u[0], v[0]], [u[1], v[1]])

def draw_ellipse():
    plt.figure(figsize=(4,5))

    X = np.linspace(-5, 5, 256)

    Y = [math.sqrt(100-4*x*x) for x in X]
    plt.plot(X, Y, color="red")

    Y = [-y for y in Y]
    plt.plot(X, Y, color="black")

def veclen(u):
    return (u[0] * u[0] + u[1] * u[1]) ** 0.5

def angle(u, v):
    dot = u[0] * v[1] - u[1] * v[0]
    sin = dot / veclen(u) / veclen(v)
    return math.asin(sin)

def rotate_vector(v, angle, origin):
    x, y = v
    x = x - origin[0]
    y = y - origin[1]
    cos_theta = math.cos(angle)
    sin_theta = math.sin(angle)
    nx = x*cos_theta - y*sin_theta
    ny = x*sin_theta + y*cos_theta
    nx = nx + origin[0]
    ny = ny + origin[1]
    return np.array([nx, ny])

def solve_quad(a, b, c):
    D = (b * b - 4 * a * c) ** 0.5
    x1 = (-b + D) / (2 * a)
    x2 = (-b - D) / (2 * a)
    return x1 + x2

def reflect(p, q):
    alpha = 2*angle(q-p,normal(q))
    p = rotate_vector(p, alpha, q)
    dx, dy = p - q
    x0, y0 = q

    a = 4 * dx * dx + dy * dy
    b = 8 * x0 * dx + 2 * y0 * dy
    c = -100 + 4 * x0 * x0 + y0 * y0;
    t = solve_quad(a, b, c)
    p = np.array([x0+dx*t, y0+dy*t])
    return q, p

draw_ellipse()
p = np.array([0, 10.1])
q = np.array([1.4, -9.6])
for i in xrange(6):
    plt.plot([p[0], q[0]], [p[1], q[1]])
    p, q = reflect(p, q)
plt.show()