傾きの計算

import numpy as np
import pylab as plt
from mpl_toolkits.mplot3d import Axes3D

def yama(x1, y1):
    h = 15 - np.sqrt((x1+10)**4 +(y1+6)**2)*0.015-np.sqrt((x1-20)**2+(y1-15)**2)*0.12
    return h
MAX_R=25
MIN_R=-25
MESH=2.5
fig = plt.figure()
ax = fig.add_subplot(projection='3d')

x = np.arange(MIN_R, MAX_R, MESH)    # 1次元のリスト
y = np.arange(MIN_R, MAX_R, MESH)

print('◆　データ数    x : ',len(x),'   y : ',len(y))

X, Y = np.meshgrid(x, y)             # 2次元のリスト
Z = yama(X,Y)                        # 2次元のリストからの変数なので Z は２次元
Z[np.less(Z, 0)] = 0
dZ_dX, dZ_dY = np.gradient(Z)  # Z が２次元あので、微分すると２次元変数が２つとなる

ax.contour(X, Y, Z, zdir='z', offset=0, cmap=plt.cm.Reds)
ax.contour(X, Y, Z, zdir='x', offset=MIN_R, cmap=plt.cm.OrRd)
ax.contour(X, Y, Z, zdir='y', offset=MAX_R, cmap=plt.cm.Greens)
ax.set_zlim(0, 16)
ax.plot_wireframe(X, Y, Z ,linewidth=0.3)
plt.show()

print('\n 傾斜のベクトルは下から上ではなく、上から下になるように向きを逆にしている\n')

fig, ax = plt.subplots(figsize=(5,5))
ax.contour(X, Y, Z,cmap=plt.cm.Reds)
ax.quiver(X, Y, -dZ_dY, -dZ_dX, angles="xy", scale_units='xy', scale=0.7,color='b')
plt.show()

print('\n 地形をグラフとして表わす\n')

fig, ax = plt.subplots(figsize=(5,5))
plt.plot(X[0],Z[0],linewidth=0.3,color='r',label='x=-20')
plt.plot(X[2],Z[2],linewidth=0.3,color='g',label='x=-16')
plt.plot(X[4],Z[4],linewidth=0.3,color='b',label='x=-12')
plt.plot(X[6],Z[6],linewidth=0.3,color='y',label='x=-8')
plt.plot(X[8],Z[8],linewidth=0.3,color='m',label='x=-4')
plt.plot(X[10],Z[10],linewidth=0.3,color='c',label='x=0')
plt.plot(X[15],Z[15],linewidth=0.3,color='k',label='x=10')
plt.legend()
plt.show()

print('\n 傾きをグラフとして表わす\n')

fig, ax = plt.subplots(figsize=(5,5))
plt.ylim(-0.6,0.6)
plt.plot(X[0],dZ_dX[0],linewidth=0.4,color='r',label='x=-20')
plt.plot(X[8],dZ_dX[8],linewidth=0.4,color='g',label='x=-4')
plt.plot(X[12],dZ_dX[12],linewidth=0.4,color='b',label='x=4')
plt.plot(X[15],dZ_dX[15],linewidth=0.4,color='c',label='x=10')
plt.plot(X[19],dZ_dX[19],linewidth=0.4,color='m',label='x=18')
plt.legend()
plt.show()