Matplotlib Streamplot error 'Столбцы' y 'должны быть равны' только для одной итерации пространственных переменных? - PullRequest
0 голосов
/ 04 мая 2020

Я попытался написать код для генерации визуализации магнитных c полей z-образных проводов. Все работает нормально, функция "magi c" генерирует трехмерный вектор поля magneti c для каждой точки P моей 3D-сетки.

Недавно я увидел простой способ (на первый взгляд) генерировать двухмерную стримплот из моих ранее трехмерных данных, используя срезы моих трехмерных данных, определяя плоскость 1D.

Это работает со следующим кодом: ax1.streamplot(X[:,:,zplane],Y[:,:,zplane],t_X[:,:,zplane],t_Y[:,:,zplane]) Здесь я фиксирую плоскость z, чтобы увидеть xy-streamplot, в то время как я обрезаю ось z в середине максимального значения z. Интересно, что если я делаю это для других измерений (например, обрезаю ось y или ось x), я либо получаю ошибку «Столбцы у должны быть равны», либо получаю пустой Streamplot. Работает только приведенный выше код, я также попытался изменить входные переменные x и y функции streamplot. Затем он переключается между пустым графиком и ошибкой «Столбцы у должны быть равны».

Ниже вы видите код, здесь я прокомментировал работающий пример «ax1», в то время как график ax2 показывает обсуждаемые ошибки .

Я надеюсь, что кто-нибудь может мне помочь с этим. Если это проблема стримплота, возможно, у кого-то есть подсказка использовать тип графика плотности для построения моего поля c magneti.

import numpy as np 
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


I = 10 #Amps
B0 = 0.05
mu0 = 4* np.pi*1e-7

#Einzelne Vektorpunkte des Z-Wires in 3D
A = np.array([10,10,10])
B = np.array([10,10,20])
C = np.array([10,20,20])
D = np.array([10,20,30])


nx, ny, nz= 20,20,20
XMAX, YMAX, ZMAX = 40,40,40

x = np.linspace(0, XMAX, nx)
y = np.linspace(0, YMAX, ny)
z = np.linspace(0, ZMAX, nz)
X,Y,Z = np.meshgrid(x,y,z) 
t_X,t_Y,t_Z = np.meshgrid(np.zeros(nx),np.zeros(ny),np.zeros(nz)) 


def magic(A,B):
    for i in range(len(x)):
        #print(i)    
        for j in range(len(y)):
            for k in range(len(z)):

                x_value=X[i][j][k]
                y_value=Y[i][j][k]
                z_value=Z[i][j][k]


                P = np.array([x_value,y_value,z_value])


                r1 = np.linalg.norm(P-A)
                r2 = np.linalg.norm(P-B)
                r3 = np.linalg.norm(B-A)


                absH = I/(4*np.pi) * (r1+r2)*(r3**2-r1**2-r2**2+2*r1*r2)/(r1*r2*np.sqrt((2*r1**2*r2**2 + 2*r1**2 * r3**2 +2 *r2**2 * r3**2 - r1**4 - r2**4- r3**4))) 

                vx = (B[1]-A[1])*(P[2]-A[2])-(P[1]-A[1])*(B[2]-A[2])
                vy = (P[0]-A[0])*(B[2]-A[2])-(B[0]-A[0])*(P[2]-A[2])
                vz = (B[0]-A[0])*(P[1]-A[1])-(P[0]-A[0])*(B[1]-A[1])

                Hx = absH * vx /(np.sqrt(vx**2+vy**2+vz**2))
                Hy = absH * vy /(np.sqrt(vx**2+vy**2+vz**2))
                Hz = absH * vz /(np.sqrt(vx**2+vy**2+vz**2))


                t_X[i][j][k]+=Hx
                t_Y[i][j][k]+=Hy
                t_Z[i][j][k]+=Hz+B0


magic(A,B)
magic(B,C)
magic(C,D)


#fig = plt.figure(figsize=(40,20))
#ax = fig.gca(projection='3d')
#Q = ax.quiver(X, Y, Z,t_X,t_Y, t_Z, length = 1)
#qk = ax.quiverkey(Q, 0.9, 0.9, 1,label='Quiver key, length = 10', labelpos='E',coordinates='figure')

#print(Z[yplane,:,:].shape)
#print(X[:,:,zplane].shape)
#print(X[:,yplane,:].shape)
#print(Z[:,yplane,:].shape)
#print(t_X[:,yplane,:].shape)
#print(t_Z[:,yplane,:].shape)
#print(t_Y[:,:,zplane].shape)

xplane = round(nx/2)
yplane = round(ny/2)
zplane = round(nz/2)

#print(yplane)

#fig, (ax1, ax2) = plt.subplots(1, 2)

fig = plt.figure()
#ax1 = fig.gca()

#THIS ONE WORKS, I DONT KNOW WHY JUST THIS ONE...

#ax1.streamplot(X[:,:,zplane],Y[:,:,zplane],t_X[:,:,zplane],t_Y[:,:,zplane])
#ax1.set_xlabel('$x$')
#ax1.set_ylabel('$y$')
#ax1.set_xlim(0, XMAX)
#ax1.set_ylim(0, YMAX)
#plt.show()


#fig2 = plt.figure()
ax2 = fig.gca()

ax2.streamplot(Z[:,yplane,:],X[:,yplane,:],t_Z[:,yplane,:],t_X[:,yplane,:])
ax2.set_xlabel('$x$')
ax2.set_ylabel('$y$')
ax2.set_xlim(0, XMAX)
ax2.set_ylim(0, YMAX)
plt.show()


Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...