import numpy as np
from mayavi import mlab
from sympy import *
from mayavi.mlab import *

q = 1 #電荷量
e0 = 4*3.14*8.85 #4×pi×ε0 放大10的12次方

def cacualte():
    x, y, z = np.mgrid[-1:1:20j, -1:1:20j,-1:1:20j]#x y z網格化，形成填充三維坐標數組

    r=np.sqrt(x**2+y**2+z**2)#與原點距離
   
    E_scale = q/(e0*(r**2))#電場大小
    Exoy = E_scale*(np.sqrt(x**2+y**2)/r)#電場在XOY投影大小
    Ex = Exoy*(x/np.sqrt(x**2+y**2))#電場在X分量大小
    Ey = Exoy*(y/np.sqrt(x**2+y**2))#電場在Y分量大小
    Ez = E_scale*(z/r)#電場在Z分量大小
    
    u = Ex #三坐標分量賦值
    v = Ey
    w = Ez
    
    scalars = x * x * 0.5 + y * y + z * z * 2.0#切片范圍
    ball = mlab.points3d(np.array(0), np.array(0.), np.array(0),color=(1,0,0),extent = [-0.01,0.01,-0.01,0.01,-0.01,0.01])#插入小球代表電荷
    obj1 = volume_slice(x,y,z,scalars, plane_orientation='x_axes',extent = [-1,1,-1,1,1,1])#插入切片以觀察矢量場
    obj = quiver3d(x, y, z, u, v, w, mask_points=1)#繪制3d矢量場
    return obj
    
cacualte()
            
          

?

?