Хабрахабр

Визуализация линий напряженности и движений электростатических зарядов, симулирование движения планет солнечной системы

Привет, сегодня я хочу вам предложить наглядное пособие по моделированию некоторых физических процессов и показать как получить красивые изображения и анимации. Осторожно много картинок.
Весь код можно найти в Google colab.

Теория

Для начала нам понадобится небольшой теоретический минимум по этой теме. Начнем с понимания что такое линии напряженности и как их считать. По сути данные линии являются слиянием множества векторов напряжённости, которую можно посчитать так:$E=\frac{k*|q|}{r^2}$. image

Метод вычисления E

Я рассчитывал вектор напряженности через подобие треугольников, получая тем самым проекции на оси x и y dx и dy соответственно.
Из подобия следует, что радиуса вектора от заряда до точки в пространстве r и длинны вектора напряженности E равно отношению проекций этих векторов (x1 и dx соответственно) $\frac{R}{E}=\frac{x}{dx}=\frac{y}{dy}$. Формула результирующего вектора $\vec{E}=\left(\begin{array}{c} \sum_{i=0}^{N}\frac{x_i*|\vec{E}_i|}{R_i} \\ \sum_{i=0}^{N}\frac{y_i*|\vec{E}_i|}{R_i} \end{array}\right) \qquad $
с этими знаниями получаем первый результат.

Функция расчета проекций

def E(q_prop, xs, ys, nq): #q_prop=[[xq1, yq1, q1, mq1, vxq1, vyq1], [xq2, yq2, q2, mq2, vxq2, vyq2] ... ] l=1 k=9*10**9 Ex=0 Ey=0 c=0 for c in range(len(q_prop)):#проходимся по всем зарядам в массиве вычисляем проекции напряженности в заданной точке и обновляем значение результирующей напряженности q=q_prop[c] r=((xs-q[0])**2+(ys-q[1])**2)**0.5 dEv=(k*q[2])/r**2 dEx=(xs-q[0])*(dEv/r)*l dEy=(ys-q[1])*(dEv/r)*l Ex+=dEx Ey+=dEy return Ex, Ey

Метод построения линий

Для начала нужно определиться с начальной и конечной точкой, откуда будет идти линия и докуда. Началом являются точки на окружности с радиусом r вокруг заряда, а концом точки отдаленные от зарядов не более чем на r.

код для определения начальных точек

theta = np.linspace(0, 2*np.pi, n)mask=q_prop[ q_prop[:,2]>0 ]#оставляем только положительные заряды for cq in range(len(mask)): qmask=mask[cq] xr = r_q*np.cos(theta)+qmask[0]#определение х-ов точек окружности вокруг заряда yr = r_q*np.sin(theta)+qmask[1]#аналогично

Так-же стоит сказать, что линии строятся только из положительных зарядов.
И наконец построение линий. Для этого мы из начальной точки строим линию вектора напряженности в ней, обновляем начальную точку на конец построенной линии и повторяем пока не будет достигнуто условия окончания, названные выше.

функция вычисления координат линий

def Draw(size, q_prop,r_q, n): linen=np.empty((np.count_nonzero(q_prop[:,2]>0),n, 2000000), dtype=np.float64) linen[:] = np.nan theta = np.linspace(0, 2*np.pi, n) mask=q_prop[ q_prop[:,2]>0 ][ q_prop[q_prop[:,2]>0][:,3]==1 ] for cq in range(len(mask)): qmask=mask[cq] x11 = r_q*np.cos(theta)+qmask[0] x22 = r_q*np.sin(theta)+qmask[1] for c in range(len(x11)): xs=x11[c] ys=x22[c] lines=np.empty((2,1000000), dtype=np.float64) lines[:]=np.nan stop=0 nnn=0 lines[0][nnn]=xs lines[1][nnn]=ys while abs(xs)<size+2 and abs(ys)<size+2: nnn+=1 for cq1 in range(len(q_prop)): q=q_prop[cq1] if ((ys-q[1])**2+(xs-q[0])**2)**0.5<r_q/2 : stop=1 break if stop==1: break dx, dy = E1(q_prop,xs,ys) xs+=dx ys+=dy lines[0][nnn]=xs lines[1][nnn]=ys linen[cq,c,:]=lines.reshape(-1) return linen 

Взаимодействие между зарядами

Чтобы отразить их взаимодействие, нужно изменять его координаты и скорость через каждое маленькое время dt.

$x+=\frac{E_x \cdot q \cdot dt^2}{2\cdot m} + v_x\cdot dt$

$y+=\frac{E_y \cdot q \cdot dt^2}{2\cdot m} + v_y \cdot dt$

$v_x+= \frac{E_x \cdot q \cdot dt}{ m} $

$v_y+= \frac{E_y \cdot q \cdot dt}{ m} $

Функция обновления координат и проекций скоростей зарядов

def Update_all(q_prop): vx=0 vy=0 x=0 y=0 q_prop_1=np.copy(q_prop) for c in range(len(q_prop)):#проход по зарядам и обновление их координат и проекций скоростей xs=q_prop[c][0] ys=q_prop[c][1] q =q_prop[c][2] m =q_prop[c][3] vx=q_prop[c][4] vy=q_prop[c][5] Ex, Ey= E(q_prop, xs, ys, c) x=(((Ex*q)/m)*dt**2)/2+vx*dt+xs y=(((Ey*q)/m)*dt**2)/2+vy*dt+ys vx+=((Ex*q)/m)*dt vy+=((Ey*q)/m)*dt #print(q_prop[c]-[x,y,q,m,vx,vy]) q_prop_1[c]=[x,y,q,m,vx,vy] return q_prop_1#возвращение обновлённого массива характеристик зарядов

Гравитация

На основе имеющегося кода я написал симулятор, отражающий движения тел под действием гравитации. Изменения в коде в основном для функции напряженности т.к. теперь будет считаться ускорение по похожей формуле.

$g=\frac{G*m}{r^2}$

$\vec{g}=\left(\begin{array}{c} \sum_{i=0}^{N}\frac{x_i*|\vec{g}_i|}{R_i} \\ \sum_{i=0}^{N}\frac{y_i*|\vec{g}_i|}{R_i} \end{array}\right) \qquad $

Планеты стартуют с оси х в перигелийном расстояние и с перигелийной скоростью. Все значения планет и солнца (массы, расстояния, экстренциситеты) из справочника.
Анимация для первых 4 планет + cолнца.

Жду критику и предложений. До свидания.

Показать больше

Похожие публикации

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»