Как мне использовать колчан в Python для полярной?

Question

Во-первых, да, я читал предыдущие темы и документацию по этой проблеме, например, Как построить колчан в полярных координатах .Это не помогло мне полностью.Позвольте мне показать вам, с чем я работаю, а затем немного кода. Это сходящийся канал, он показывает поле скорости / вектора.Ясно, что у меня есть только радиальная составляющая, но она меняется в зависимости от угла тета.Этот образец стрелок повторяется, когда мы спускаемся (поток) к альфе.Так что сюжет должен быть простым, верно.Вот уравнение для радиальной составляющей скорости:

Теперь, прежде чем я покажу свой код, я сохранил значения f (тета) для ряда тэт.Эта функция, f, должна быть численно решена, и я сохранил ее как вектор u [0].Вот что я сейчас делаю в своем коде:

radii = np.linspace(0.1,1,11)
thetas = np.linspace(-alpha,alpha,20)
theta, r = np.meshgrid(thetas, radii)

q = 0.0001


dr = [-q/x for x in radii]*u_sol[0]
dt = 0

f = plt.figure()
ax = f.add_subplot(111, polar=True)

ax.quiver(theta, r, dr * cos(theta) - dt * sin (theta), dr * sin(theta) +     
dt* cos(theta))

Пятое выражение для переменной dr было отчаянной попыткой умножить все r фиксированной длины в сетке на u [0], но онине имеют одинаковые размеры, поэтому он не работает.Так что я застрял.

Мой вопрос: как получить векторное поле для сходящегося канала?Я не могу на самом деле поставить последние части вместе, могу ли я манипулировать сеткой?

Результаты на данный момент в MATLAB:

Редактировать Код выше был взят по ссылке в начале моего текста.Я сделал некоторые изменения в dr и dt, но в остальном ничего.

Answer 1

Единственная реальная проблема с вашим кодом - это проблема numpy, т. Е. В вашем dr неправильные размеры.С небольшими изменениями в вашем коде:

from matplotlib import pyplot as plt
import numpy as np

#to make the code runnable
u_sol = [1]
alpha0 = 5*np.pi/180
alpha = 10*np.pi/180

radii = np.linspace(0.2,1,10)
print(radii)
thetas = np.linspace(alpha0-alpha,alpha0+alpha,20)
print(thetas)
theta, r = np.meshgrid(thetas, radii)

q = 0.0001


#dr = [-q/x for x in radii]*u_sol[0]
dr = -q/r
dt = 0

f = plt.figure()
ax = f.add_subplot(111, polar=True)

ax.quiver(
    theta, r,
    dr * np.cos(theta) - dt * np.sin(theta),
    dr * np.sin(theta) + dt * np.cos(theta),
)

plt.show()

я получаю следующее изображение:

Обратите внимание, что в radiiопределение, я переместил нижний предел с 0,1 до 0,2, так как в противном случае стрелки становятся настолько длинными, что они указывают на другую сторону начала координат, что выглядит довольно странно.