Визуализация линий напряженности и движений электростатических зарядов, симулирование движения планет солнечной системы
Привет, сегодня я хочу вам предложить наглядное пособие по моделированию некоторых физических процессов и показать как получить красивые изображения и анимации. Осторожно много картинок.
Весь код можно найти в Google colab.
Теория
Для начала нам понадобится небольшой теоретический минимум по этой теме. Начнем с понимания что такое линии напряженности и как их считать. По сути данные линии являются слиянием множества векторов напряжённости, которую можно посчитать так:.
Метод вычисления E
Я рассчитывал вектор напряженности через подобие треугольников, получая тем самым проекции на оси x и y dx и dy соответственно.
Из подобия следует, что радиуса вектора от заряда до точки в пространстве r и длинны вектора напряженности E равно отношению проекций этих векторов (x1 и dx соответственно) . Формула результирующего вектора
с этими знаниями получаем первый результат.
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.
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#возвращение обновлённого массива характеристик зарядов
Гравитация
На основе имеющегося кода я написал симулятор, отражающий движения тел под действием гравитации. Изменения в коде в основном для функции напряженности т.к. теперь будет считаться ускорение по похожей формуле.
Планеты стартуют с оси х в перигелийном расстояние и с перигелийной скоростью. Все значения планет и солнца (массы, расстояния, экстренциситеты) из справочника.
Анимация для первых 4 планет + cолнца.
Жду критику и предложений. До свидания.