梯度场表现实现-matplotlib.tri.CubicTriInterpolator

以前想画两个电场之间的势能场的时候，只会用Axes.quiver函数，最后只能表现势能场的方向，但是势能场因为距离呈现的疏密却没有表现出来，在参考了文档之后，偷了文档的一些代码总结一下：

1、首先画出不规则三角坐标：

import matplotlib.tri as tri
import numpy
import matplotlib.pyplot as plt
theta=numpy.linspace(0,2*numpy.pi,30)
r=numpy.linspace(0.25,0.95,10)
theta=theta[...,numpy.newaxis]
theta=numpy.repeat(theta,10,axis=1)
x=r*numpy.cos(theta)
y=r*numpy.sin(theta)
x=x.flatten()
y=y.flatten()
triang=tri.Triangulation(x,y)
fig,ax=plt.subplots(1)
ax.axis('equal')
ax.triplot(triang)

这样的结果是：

中间会有一些奇异点需要滤除，加上代码：

triang.set_mask(numpy.hypot(x[triang.triangles].mean(axis=1),
                         y[triang.triangles].mean(axis=1))
                < 0.25)

mask掉一些奇异点就可以：

2、先求电场强度，为后面求电势做铺垫。

r_sq = x**2 + y**2
theta = numpy.arctan2(y, x)
z = numpy.cos(theta)/r_sq
V=(numpy.max(z) - z) / (numpy.max(z) - numpy.min(z))

因为电场强度需要有一个参考值，这里就用相对值来表示电场的强度，也就是用所求的电场强度的最大值作为参考，用(max(z)-z)/(max(z)-min(z))来表示相对值。

3、在画好三角场和求得电场强度之后，就要用三角场来求电场的梯度，也就是电势。

tci = tri.CubicTriInterpolator(triang, -V)
(Ex, Ey) = tci.gradient(triang.x, triang.y)
E_norm = numpy.sqrt(Ex**2 + Ey**2)

4、最后用quiver画出梯度场：

ax.quiver(triang.x, triang.y, Ex/E_norm, Ey/E_norm,
          units='xy', scale=10., zorder=3, color='blue',
          width=0.007, headwidth=3., headlength=4.)

然后我惊奇地发现跟我当年画的差不多。。。。。。

回去翻看代码，发现了一句很有意思的代码：

theta[:,1::2]+=numpy.pi/30

这是加上去的效果，莫非这个才是梯度图好看的诀窍吗。