用gyrm3D/python绘制三维矢量场



我在用mayavi的函数jump3d绘制三维矢量场时遇到了一些问题。当我在终端中执行程序时,会出现以下内容:

Traceback (most recent call last):
  File "inter3d.py", line 80, in <module>
    mlab.quiver3d(x,y,z,Bx,By,Bz,scale_factor=1 ,mask_points =5)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 34, in the_function
    return pipeline(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 79, in __call__
    output = self.__call_internal__(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 175, in __call_internal__
    g = Pipeline.__call_internal__(self, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 90, in __call_internal__
    self.source = self._source_function(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/sources.py", line 991, in vector_scatter
    x, y, z, u, v, w = process_regular_vectors(*args)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/sources.py", line 904, in process_regular_vectors
    w.shape == v.shape ), "argument shape are not equal"
AssertionError: argument shape are not equal

这是我的程序:

# -*- coding: utf-8 -*-
from numpy import *
from math import *
import numpy as numpy
import mayavi.mlab as mlab 

#valeur initiale
a=6.371*10**6
b0=31200*10**-9
Bp=0.
nb_x=30
nb_y=30
nb_z=30
Lx=10*10**7
Ly=10*10**7
Lz=10*10**7
dx=Lx/nb_x
dy=Ly/nb_y
dz=Lz/nb_z
N = int(nb_x*nb_y*nb_z)
x=zeros(nb_x)
y=zeros(nb_y)
z=zeros(nb_z)
Bx=zeros((nb_x,nb_y,nb_z))
By=zeros((nb_x,nb_y,nb_z))
Bz=zeros((nb_x,nb_y,nb_z))
#xp=12*10**6
#yp=12*10**6
l=0
i=0
j=0
u=0
mu=4*pi*10**-7
k=mu/4*pi
M=8*10**22
#print N, dx, dy
########################
fichier = open('data.dat','w')
while i < nb_x :
    j=0
    while j < nb_y:
        l=0
        while l < nb_z:
            x[i]=-Lx/2+i*dx
            y[j]=-Ly/2+j*dy
            z[l]=-Lz/2+l*dz
            #print B[u,0],B[u,1], B[u,2]
            r=sqrt(x[i]**2+y[j]**2+z[l]**2)
            if r>=6.37*10**6:
                if y[j] >= 0:
                    teta=acos(x[i]/r)
                else:
                    teta=2*pi-acos(x[i]/r)
                phi=z[l]/r
                Br=2*k*M*cos(teta)/(r**3)
                Bteta=k*M*sin(teta)/(r**3)
                By[i,j,l]=sin(teta)*(Bteta*cos(phi)+Br*sin(phi))
                Bx[i,j,l]=-cos(teta)*(Br*sin(phi)+Bteta*cos(phi))
                Bz[i,j,l]=cos(phi)*Br-sin(teta)*Bteta
                #B[i,j]=sqrt(Bx[i,j]**2+By[i,j]**2)#calcul de la norme du champs
            else:
                Bx[i,j,l]=0.
                By[i,j,l]=0.
                Bz[i,j,l]=0.
            #fichier.write(str(B[u,0])+" "+str(B[u,1])+" "+str(B[u,2])+"n")    
            fichier.write(str(i)+" "+str(j)+" "+str(l)+" "+str(Bx[i,j,l])+" "+str(By[i,j,l])+" "+str(Bz[i,j,l])+"n")
            l=l+1
            u=u+1
            print u
        j=j+1
    i=i+1

fichier.close()
print dx, dy, dz
#####################
mlab.quiver3d(x,y,z,Bx,By,Bz)

我在谷歌上发现有人使用这个函数来绘制电场,但它不是数组,也许这就是我问题的原因,但当我使用matplot中的2D箭袋时,它的工作原理是这样写的:quiver(x,y,Bx,By)

编辑另一种尝试:

# -*- coding: utf-8 -*-
from numpy import *
from math import *
import numpy as np
import mayavi.mlab as mlab 
from pylab import *
import pylab as pl
#valeur initiale
a=6.371*10**6
b0=31200*10**-9
Bp=0.
nb_x=30
nb_y=30
nb_z=30
Lx=10*10**7
Ly=10*10**7
Lz=10*10**7
dx=Lx/nb_x
dy=Ly/nb_y
dz=Lz/nb_z
N = int(nb_x*nb_y*nb_z)
x=0.
y=0.
z=0.
Bx=0.
By=0.
Bz=0.
#xp=12*10**6
#yp=12*10**6
l=0
i=0
j=0
u=0
teta=0.
phi=0.
mu=4*pi*10**-7
k=mu/4*pi
M=8*10**22
#print N, dx, dy
########################
def u(x,y,z,k,M):
    r=sqrt(x**2+y**2+z**2)
    print r
    if r>=6.37*10**2:
        if y >= 0:
            teta=acos(x/r)
        else:
            teta=2*pi-acos(x/r)
        phi=z/r
        Br=2*k*M*cos(teta)/(r**3)
        Bteta=k*M*sin(teta)/(r**3)
        Bx=-cos(teta)*(Br*sin(phi)+Bteta*cos(phi))
        #B[i,j]=sqrt(Bx[i,j]**2+By[i,j]**2)#calcul de la norme du champs
    else:
        Bx=0.
    return Bx
def v(x,y,z,k,M):
    r=sqrt(x**2+y**2+z**2)
    if r>=6.37*10**2:
        if y >= 0:
            teta=acos(x/r)
        else:
            teta=2*pi-acos(x/r)
        phi=z/r
        Br=2*k*M*cos(teta)/(r**3)
        Bteta=k*M*sin(teta)/(r**3)
        By=sin(teta)*(Bteta*cos(phi)+Br*sin(phi))
        #B[i,j]=sqrt(Bx[i,j]**2+By[i,j]**2)#calcul de la norme du champs
    else:
        By=0.
    return By
def w(x,y,z,k,M):
    r=sqrt(x**2+y**2+z**2)
    if r>=6.37*10**2:
        if y >= 0:
            teta=acos(x/r)
        else:
            teta=2*pi-acos(x/r)
        phi=z/r
        Br=2*k*M*cos(teta)/(r**3)
        Bteta=k*M*sin(teta)/(r**3)
        Bz=cos(phi)*Br-sin(teta)*Bteta
        #B[i,j]=sqrt(Bx[i,j]**2+By[i,j]**2)#calcul de la norme du champs
    else:
        Bz=0.
    return Bz
#####################
x, y, z = np.mgrid[-5*10**7:5*10**7:400000, -5*10**7:5*10**7:400000, -5*10**7:5*10**7:400000]
mayavi.mlab.quiver3d(x, y, z, u, v, w, line_width=3, scale_factor=1)

终端内:

    Traceback (most recent call last):
  File "test2.py", line 97, in <module>
    mayavi.mlab.quiver3d(x, y, z, u, v, w, line_width=3, scale_factor=1)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 34, in the_function
    return pipeline(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 79, in __call__
    output = self.__call_internal__(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 175, in __call_internal__
    g = Pipeline.__call_internal__(self, *args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/helper_functions.py", line 90, in __call_internal__
    self.source = self._source_function(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/sources.py", line 991, in vector_scatter
    x, y, z, u, v, w = process_regular_vectors(*args)
  File "/usr/lib/python2.7/dist-packages/mayavi/tools/sources.py", line 902, in process_regular_vectors
    u.shape == z.shape and
AttributeError: 'function' object has no attribute 'shape'

我认为你的问题是u、v和w的形状,你试过读u、v、w的形状吗?它们的大小必须与x、y、z完全相同。

相关内容

  • 没有找到相关文章

最新更新