Как построить поверхность с разными функциями на разных интервалах?

Question

У меня есть некоторые числовые данные для части поверхности потенциальной энергии, для которых я использовал бикубическую интерполяцию, как объяснено на этой странице Википедии . Теперь я хочу построить полученную поверхность, возможно, таким образом, чтобы цвета указывали значения функции (вид тепловой карты).

У меня есть регулярное смешение (m x m) точек и (m-1) x (m-1) функций. У меня есть значения x, y (при p (0,0) в нотации Википедии) и функциональные коэффициенты, хранящиеся в файле с именем coefficients_ci.txt .

Вот моя попытка сделать это. Я получаю некоторую поверхность , но это не только поверхность - пространство между плоскостью xy (при z = 0) и моей поверхностью заполнено цветами.
Вот мой код:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt


with open ('coefficients_ci.txt', 'r') as all_data:
    lines3 = np.array([line.split() for line in all_data])
all_input_data3 = lines3.astype(float)



def fun(u, v, method, i):


    a00 = all_input_data3[i,3]
    a10 = all_input_data3[i,4]
    a20 = all_input_data3[i,5]
    a30 = all_input_data3[i,6]
    a01 = all_input_data3[i,7]
    a11 = all_input_data3[i,8]
    a21 = all_input_data3[i,9]
    a31 = all_input_data3[i,10]
    a02 = all_input_data3[i,11]    
    a12 = all_input_data3[i,12]
    a22 = all_input_data3[i,13]
    a32 = all_input_data3[i,14]
    a03 = all_input_data3[i,15]    
    a13 = all_input_data3[i,16]    
    a23 = all_input_data3[i,17]    
    a33 = all_input_data3[i,18]   

    x0 = all_input_data3[i, 1]
    y0 = all_input_data3[i, 2]
#    transform coordinates so that starting point is 0,0    
    x = (u-x0)*100
    y = (v-y0)*100


    return a00 + a01*y + a02*y**2 + a03*y**3\
        + a10*x + a11*x*y + a12*x*y**2 + a13*x*y**3\
        + a20*x**2 + a21*x**2*y + a22*x**2*y**2 + a23*x**2*y**3\
        + a30*x**3 + a31*x**3*y + a32*x**3*y**2 + a33*x**3*y**3   


def plotData(method):


    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')


    for i in range(0, len(all_input_data3)):
        print(i)
        u = np.arange(all_input_data3[i, 1], all_input_data3[i, 1]+0.01, 0.0001)
        v = np.arange(all_input_data3[i, 2], all_input_data3[i, 2]+0.01, 0.0001)

        X, Y = np.meshgrid(u, v)
        zs = np.array([fun(u,v, method, i) for u,v in zip(np.ravel(X), np.ravel(Y))])
        print('Calculated function values')        

        Z = zs.reshape(X.shape)
        ax.plot_surface(X, Y, Z)

        continue
#    plt.show()
    plt.savefig('plot_%s.pdf' % method)
    return

Я надеюсь, что моя проблема описана достаточно ясно и что кто-то сможет мне помочь.