d = [(-70,2000,330),(-70,3000,303),(-70,4000,280),(-70,4500,265),(-70,5000,245),(-70,5500,220),(-70,6000,195),(-70,6500,160),(-70,7000,120),(-70,7500,80),(-60,1500,345),(-60,2000,344),
(-60,2500,343),(-60,3000,341),(-60,3500,341),(-60,4000,340),(-60,4500,335),(-60,5000,330),(-60,5500,322),(-60,6000,312),(-60,6500,302),(-60,7000,290),(-60,7500,274),(-60,8000,250),
(-60,8500,225),(-60,9000,198),(-60,9500,158),(-60,10000,120),(-50,1500,375),(-50,2000,376),(-50,2500,377),(-50,3000,378),(-50,3500,378),(-50,4000,379),(-50,4500,378.5),(-50,5000,378),
(-50,5500,377),(-50,6000,375),(-50,6500,370),(-50,7000,365),(-50,7500,360),(-50,8000,352),(-50,8500,345),(-50,9000,340),(-50,9500,330),(-50,10000,320),(-50,11000,280),(-50,11500,255),
(-50,12000,220),(-50,13000,145),(-50,13500,100),(-40,1500,380),(-40,2000,380),(-40,2500,381),(-40,3000,383),(-40,3500,385),(-40,4000,390),(-40,4500,394),(-40,5000,398),(-40,5500,400),
(-40,6000,407),(-40,6500,410),(-40,7000,415),(-40,7500,417),(-40,8000,419),(-40,8500,418),(-40,9000,417),(-40,9500,415),(-40,10000,410),(-40,10500,405),(-40,11000,396),(-40,12000,379),
(-40,12500,370),(-40,13000,355),(-40,13500,338),(-40,14000,320),(-40,14500,300),(-40,15000,275),(-40,16000,215),(-30,1500,400),(-30,2000,403),(-30,2500,405),(-30,3000,410),(-30,3500,415),
(-30,4000,418),(-30,4500,420),(-30,5000,425),(-30,5500,430),(-30,6000,435),(-30,6500,438),(-30,7000,440),(-30,7500,442),(-30,8000,445),(-30,8500,447),(-30,9000,450),(-30,9500,452),
(-30,10000,453),(-30,10500,455),(-30,11000,457),(-30,11500,455),(-30,12000,450),(-30,12500,447),(-30,13000,445),(-30,13500,442),(-30,14000,440),(-30,14500,435),(-30,15000,427),
(-30,15500,419),(-30,16000,400),(-30,16500,380),(-30,17000,360),(-30,17500,330)]


import numpy
import itertools
import matplotlib.pyplot as plt
from matplotlib import cm
import mpl_toolkits.mplot3d.axes3d as axes3d

#least square fitting function taken from https://stackoverflow.com/questions/7997152/python-3d-polynomial-surface-fit-order-dependent
def polyfit2d(x, y, z, order=3):
    ncols = (order + 1)**2
    G = numpy.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), range(order+1))
    for k, (i,j) in enumerate(ij):
        G[:,k] = x**i * y**j
    m, _, _, _ = numpy.linalg.lstsq(G, z)
    return m

def polyval2d(x, y, m):
    order = int(numpy.sqrt(len(m))) - 1
    ij = itertools.product(range(order+1), range(order+1))
    z = numpy.zeros_like(x)
    for a, (i,j) in zip(m, ij):
        z += a * x**i * y**j
    return z

def equation(m):
    order = int(numpy.sqrt(len(m))) - 1
    ij = itertools.product(range(order+1), range(order+1))
    print " z = ",
    for a, (i,j) in zip(m, ij):
        print "%1.8f*(x^%i*y^%i) +"%(a,i,j),

#SCALING
d = numpy.array(d)
dOffset = d.min(axis=0)
d = d - dOffset #offset the d closer to zero
dScale = d.max(axis=0)
d = d/dScale #scale the d so it ranges from 0-1 on every axis

#FITTING
m = polyfit2d(d[:,0],d[:,1],d[:,2],order = 3) # fit d, polynms returnded

#CALCULATION OF FIT DATA
nx, ny = 50, 50 # gridsize
xx, yy = numpy.meshgrid(numpy.linspace(d[:,0].min(), d[:,0].max(), nx), numpy.linspace(d[:,1].min(), d[:,1].max(), ny))
zz = polyval2d(xx, yy, m) # calculate gridpoints
fx = polyval2d(d[:,0], d[:,1], m) # calculate values at the original locations
error = fx-d[:,2] #so we can get the error between actual and fit d

#SCALING
xx = xx*dScale[0]+dOffset[0]  #scale the results
yy = yy*dScale[1]+dOffset[1]
zz = zz*dScale[2]+dOffset[2]
d = d*dScale+dOffset  # scale the source d to original
error = error*dOffset[2] #scale the errors too

#PRINT FIT FUNC
equation(m)

#PLOTTING
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
p = ax.plot_wireframe(xx,yy,zz, rstride=4, cstride=4, alpha=0.4) # plot a mesh
p = ax.scatter3D(d[:,0],d[:,1],d[:,2], c=error, cmap='hot', s=100) # plot the original d with error coloring
fig.colorbar(p, ax=ax)
plt.show()