# Locate minimum in X,Y data set, fit to quadratic
# polynomial P(x) = A*x**2 + B*x + C, and find the minimum
# xmin of P(x).
def QuadraticMin(X,Y):
    n = len(X)
    h = X[1] - X[0]

    x0 = None

    # Locate possible minimum
    for i in range(1,n-1):
        if ( (Y[i-1] > Y[i]) and (Y[i] < Y[i+1]) ):
            x0 = X[i]
            (ym,y0,yp) = (Y[i-1], Y[i], Y[i+1])
            break

    # No minimum
    if (x0 == None):
        print '%s' % ('No minimum found')
        return

    # Fit polynomial to minimum interval
    A = (yp + ym - 2.0*y0)/(2.0*h*h)
    B = (yp - ym)/(2.0*h)
    xmin = x0 - B/(2.0*A)
    print 'Quadratic fit: xmin = %g' % (xmin)

    return


# Read two-column data from file
f = open('energy.dat', 'r')
data = f.readlines()
f.close()

# Partition data into separate lists 
X = [float(i.split()[0]) for i in data]
Y = [float(i.split()[1]) for i in data]

QuadraticMin(X,Y)