#!/usr/bin/python

from math import sqrt
from sys import argv

def segment(x):
    y = dy / dx * x + b
    return y

def dist(x0, y0, f, x):
    y = f(x)
    d = sqrt((x - x0) ** 2 + (y - y0) ** 2)
    return d

def quadratic_optimization(x1, x2, x3, f):
    for i in range(1, 5):
        y1 = f(x1) 
        y2 = f(x2) 
        y3 = f(x3) 
        denom = (x1 - x2) * (x1 - x3) * (x2 - x3)
        A = (x3 * (y2 - y1) + x2 * (y1 - y3) + x1 * (y3 - y2)) / denom
        B = (x3 ** 2 * (y1 - y2) + x2 ** 2
                * (y3 - y1) + x1 ** 2 * (y2 - y3)) / denom
        C = (x2 * x3 * (x2 - x3) * y1
                + x3 * x1 * (x3 - x1) * y2
                + x1 * x2 * (x1 - x2) * y3) / denom
        x0 = -B / (2 * A)
        y0 = (4 * A * C - B ** 2) / (4 * A)
        print 'p%d(x) = %f * x ** 2 + %f * x + %f' % (i, A, B, C)
        if y1 > y2 and y1 > y3:
            x1 = x0
        elif y2 > y1 and y2 > y3:
            x2 = x0
        else:
            x3 = x0

x0 = float(argv[1])
y0 = float(argv[2])
dy = float(argv[3])
dx = float(argv[4])
b = float(argv[5])

def f(x):
    return dist(x0, y0, segment, x)

print 'x0 = %f\ny0 = %f\n' % (x0, y0)
print 'set xrange [-5:15]'
print 'set yrange [-10:20]'
print 'set size square'
print 'set label at x0, y0 "" point pointtype 7 pointsize 1'
print 's(x) = %f / %f * x + %f' % (dy, dx, b)
print 'd(x) = sqrt((x - x0) ** 2 + (s(x) - y0) ** 2)'

quadratic_optimization(5, 10, 15, f)

print '''plot s(x) with lines, \\
    d(x) with lines, \\
    p1(x) with lines, \\
    p2(x) with lines, \\
    p3(x) with lines, \\
    p4(x) with lines
'''

