#! /usr/bin/env python

import numpy as np

def f(x,y):
    return 3*x**3 + 4*y**3 - 10*x*y**2 - 2*x**2 + 0.5*y**2 + 2*x + 9*y - 2*x*y + 1

A = np.linspace(0, 20, 20)
M_X, M_Y = np.meshgrid(A, A)
M_Z = np.empty(M_X.shape)

P, Y = [], []
for i in xrange(M_X.shape[0]):
    for j in xrange(M_X.shape[1]):
        P.append([i,j])
        Y.append(f(i,j))

import professor2 as prof
import sys
# order = int(sys.argv[1]) if len(sys.argv) > 1 else 3
# ipol = prof.Ipol(P, Y, order)
ipol= prof.SimplexInterpolation( P, Y )

for i in xrange( M_X.shape[0] ):
    for j in xrange( M_X.shape[1] ):
        M_Z[i,j] = ipol.value( [ M_X[i,j], M_Y[i,j] ] )


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X1, Y1 = zip(*P)
Z1 = Y
ax.scatter(X1, Y1, Z1, c="r")
ax.plot_wireframe(M_X, M_Y, M_Z, rstride=1, cstride=1)
plt.show()
