714.sagews

var('x,y,z,t')

(x, y, z, t)

def plot_vector_field3d(vec, xrange, yrange, zrange, plot_points=[5,5,5], center_arrows=False,**kwds):
    xvar, xmin, xmax = xrange
    yvar, ymin, ymax = yrange
    zvar, zmin, zmax = zrange
    if not isinstance(plot_points, (list, tuple)):
        plot_points = [plot_points]*3
    ff, gg, hh = fast_float(vec, xvar, yvar, zvar)
    xpoints = [xmin..xmax, step=float(xmax-xmin)/(plot_points[0]-1)][0:plot_points[0]]
    ypoints = [ymin..ymax, step=float(ymax-ymin)/(plot_points[1]-1)][0:plot_points[1]]
    zpoints = [zmin..zmax, step=float(zmax-zmin)/(plot_points[2]-1)][0:plot_points[2]]
    points = [vector((i,j,k)) for i in xpoints for j in ypoints for k in zpoints]
    vectors = [vector((ff(*point), gg(*point), hh(*point))) for point in points]
    # scale the vectors
    max_len = max(v.norm() for v in vectors)
    scaled_vectors = [v/max_len for v in vectors]
    if center_arrows:
        return sum([plot(v,**kwds).translate(p-v/2) for v,p in zip(scaled_vectors, points)])
    else:
        return sum([plot(v,**kwds).translate(p) for v,p in zip(scaled_vectors, points)])

r=sqrt(x^2+y^2);Ex=(x/(sqrt((z-1)^2+r^2))^3 -x/(sqrt((z+1)^2+r^2))^3 ) ;Ey=(y/(sqrt((z-1)^2+r^2))^3 -y/(sqrt((z+1)^2+r^2))^3 ) ; Ez=((z-1)/(sqrt((z-1)^2+r^2))^3 -(z+1)/(sqrt((z+1)^2+r^2))^3);

sage: plot_vector_field3d((Ex,Ey,Ez), (x,-3,3), (y,-3,3), (z,-3,3), plot_points=8)