Homec4science
Diffusion JuNoTE_git_repository (master) Browse

R9442/notebooks/ElectrostaticsCode.py
630809e856b0master

Last Change
Raw File

ElectrostaticsCode.py
View Options

import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from ipywidgets import *
import matplotlib
from numpy import *
from pylab import *
from scipy.integrate import ode
import math
from itertools import combinations
import matplotlib.patches as patches
import copy
#%matplotlib notebook
def update_two_point_charges(wx, wy, q1, q2):
k = 9*10^9
f = lambda r : (k*q1*q2)/r**2
r = np.linspace(0.01, 0.25, 64)
force = [abs(f(ri)) for ri in list(r)]
fig1 = plt.figure(figsize=(15, 8))
ax = plt.subplot(121)
ax.set_title('Force vs Separation')
ax.set_xlabel('Separation')
ax.set_ylabel('Force')
ax.plot(r, force)
ax2 = plt.subplot(122)
ax2.add_artist(Circle((0.25,0.5), 0.005, color="red"))
#plt.figure(figsize=(10,8))
plt.figure(figsize=(10, 8))
ax2.add_artist(Circle((wx,wy), 0.005, color="blue"))
f_mod = f(math.sqrt((wx-0.25)**2 + (wy-0.5)**2))*1e-4
vec_norm = math.sqrt((-wx+0.25)**2 + ( -wy+0.5)**2)
vec1 = f_mod*(wx-0.25)/vec_norm, f_mod*( wy-0.5)/vec_norm
vec2 = f_mod*(-wx+0.25)/vec_norm, f_mod*( -wy+0.5)/vec_norm
ax2.arrow(0.25, 0.5, *vec2, head_width=0.005, head_length=0.01)
ax2.arrow(wx, wy, *vec1, head_width=0.005, head_length=0.01)
plt.show()
output_slider_variable = widgets.Text()
def set_n_charges(nq):
output_slider_variable.value = str(nq)
return output_slider_variable
def return_widget_list(nc):
w_list_x = []
w_list_y = []
w_dic = {}
val_pos = []
for i in range(nc):
w_dic[str(i)+"xxxx"] = widgets.FloatSlider(min=-2, max=2, step=.2, value=np.cos(2*np.pi*i/nc))
w_dic[str(i)+"y"] = widgets.FloatSlider(min=-2, max=2, step=.2, value=np.sin(2*np.pi*i/nc))
w_dic[str(i)+"c"] = widgets.IntSlider(min=-2, max=2, step=1, value=1)
val_pos.append((np.cos(2*np.pi*i/nc), np.sin(2*np.pi*i/nc)))
return w_dic
def multiple_point_charges(**w_dic):
k = 9*1e9
fig = plt.figure(figsize=(10, 8))
ax = plt.subplot(111)
charges =[]
for i in range(int(output_slider_variable.value)):
q = w_dic[str(i)+"c"]
charges.append((q, (w_dic[str(i)+"x"], w_dic[str(i)+"y"])))
#calculating and displaying force on the 0th charge
charge_colors = {True: '#aa0000', False: '#0000aa'}
force_on = charges[0]
vec = []
force_ = []
ax.set_xlim(-2, 2)
ax.set_ylim(-1.5, 1.5)
for q, pos in charges:
ax.add_artist(Circle(pos, 0.05, color=charge_colors[q>0]))
vec_= {}
force_ = {}
indx = [i for i in range(int(output_slider_variable.value))]
pairs = combinations(indx, 2)
for p1, p2 in pairs:
x_v = w_dic[str(p1)+"x"]-w_dic[str(p2)+"x"]
y_v = w_dic[str(p1)+"y"]-w_dic[str(p2)+"y"]
v_norm = math.sqrt(x_v**2 + y_v**2)
u_x, u_y = x_v/v_norm, y_v/v_norm
f_norm = k*w_dic[str(p1)+"c"]*w_dic[str(p2)+"c"]/(v_norm**2)
vec_.setdefault(p1, []).append(( f_norm*u_x, f_norm*u_y))
vec_.setdefault(p2, []).append(( -f_norm*u_x, -f_norm*u_y))
vec_each_charge = {}
for i in indx:
vec = vec_[i]
x_vec = sum(i for i, j in vec)
y_vec = sum(j for i, j in vec)
vec_each_charge[i] = (1e-10*x_vec, 1e-10*y_vec)
for i in indx:
ax.arrow(w_dic[str(i)+"x"], w_dic[str(i)+"y"], vec_each_charge[i][0], vec_each_charge[i][1], head_width=0.1, head_length=0.1)
plt.show()
class charge:
def __init__(self, q, pos):
self.q=q
self.pos=pos
def E_point_charge(q, a, x, y):
return q*(x-a[0])/((x-a[0])**2+(y-a[1])**2)**(1.5), \
q*(y-a[1])/((x-a[0])**2+(y-a[1])**2)**(1.5)
def E_total(x, y, charges):
Ex, Ey=0, 0
for C in charges:
E=E_point_charge(C.q, C.pos, x, y)
Ex=Ex+E[0]
Ey=Ey+E[1]
return [ Ex, Ey ]
def E_dir(t, y, charges):
Ex, Ey=E_total(y[0], y[1], charges)
n=sqrt(Ex**2+Ey*Ey)
return [Ex/n, Ey/n]
def update_dipole(wx):
plt.figure(figsize=(6, 4.5))
charges = []
# charges and positions
charges=[ charge(1, [-1, 0]), charge(-1, [wx, 0]) ]
# plot field lines
x0, x1=-2, 2
y0, y1=-1.5, 1.5
R=0.01
# loop over all charges
for C in charges:
# plot field lines starting in current charge
dt=0.8*R
if C.q<0:
dt=-dt
# loop over field lines starting in different directions
# around current charge
for alpha in np.linspace(0, 2*pi*15/16, 16):
r=ode(E_dir)
r.set_integrator('vode')
r.set_f_params(charges)
x=[ C.pos[0] + cos(alpha)*R ]
y=[ C.pos[1] + sin(alpha)*R ]
r.set_initial_value([x[0], y[0]], 0)
while r.successful():
r.integrate(r.t+dt)
x.append(r.y[0])
y.append(r.y[1])
hit_charge=False
# check if field line left drwaing area or ends in some charge
for C2 in charges:
if sqrt((r.y[0]-C2.pos[0])**2+(r.y[1]-C2.pos[1])**2)<R:
hit_charge=True
if hit_charge or (not (x0<r.y[0] and r.y[0]<x1)) or \
(not (y0<r.y[1] and r.y[1]<y1)):
break
ax = plt.gca()
as_xi = int(len(x)/2)
as_yi = int(len(y)/2)
a_s_x = x[as_xi]
a_s_y = y[as_yi]
a_e_x = x[as_xi+1]
a_e_y = y[as_yi+1]
if (C.q > 0):
ax.arrow(a_s_x, a_s_y, (a_e_x-a_s_x), (a_e_y-a_s_y) , head_width=0.1)
else:
ax.arrow(a_s_x, a_s_y, -(a_e_x-a_s_x), -(a_e_y-a_s_y) , head_width=0.1)
plot(x, y, '-k')
# plot point charges
for C in charges:
if C.q>0:
plot(C.pos[0], C.pos[1], 'bo', ms=8*sqrt(C.q))
if C.q<0:
plot(C.pos[0], C.pos[1], 'ro', ms=8*sqrt(-C.q))
output_slider_variable = widgets.Text()
def update(**w_dic):
plt.figure(figsize=(6, 4.5))
charges = []
# charges and positions
for i in range(int(output_slider_variable.value)):
q = w_dic[str(i)+"c"]
charges.append(charge(q, [w_dic[str(i)+"x"], w_dic[str(i)+"y"]]))
# plot field lines
x0, x1=-2, 2
y0, y1=-1.5, 1.5
R=0.01
# loop over all charges
xs = []
ys = []
for C in charges:
# plot field lines starting in current charge
dt=0.8*R
if C.q<0:
dt=-dt
# loop over field lines starting in different directions
# around current charge
for alpha in np.linspace(0, 2*pi*15/16, 16):
r=ode(E_dir)
r.set_integrator('vode')
r.set_f_params(charges)
x=[ C.pos[0] + cos(alpha)*R ]
y=[ C.pos[1] + sin(alpha)*R ]
r.set_initial_value([x[0], y[0]], 0)
while r.successful():
r.integrate(r.t+dt)
x.append(r.y[0])
y.append(r.y[1])
hit_charge=False
# check if field line left drwaing area or ends in some charge
for C2 in charges:
if sqrt((r.y[0]-C2.pos[0])**2+(r.y[1]-C2.pos[1])**2)<R:
hit_charge=True
if hit_charge or (not (x0<r.y[0] and r.y[0]<x1)) or \
(not (y0<r.y[1] and r.y[1]<y1)):
break
ax = plt.gca()
as_xi = int(len(x)/2)
as_yi = int(len(y)/2)
a_s_x = x[as_xi]
a_s_y = y[as_yi]
a_e_x = x[as_xi+1]
a_e_y = y[as_yi+1]
if (C.q > 0):
ax.arrow(a_s_x, a_s_y, (a_e_x-a_s_x), (a_e_y-a_s_y) , head_width=0.1)
else:
ax.arrow(a_s_x, a_s_y, -(a_e_x-a_s_x), -(a_e_y-a_s_y) , head_width=0.1)
plt.plot(x, y, color="k")
xs.append(x)
ys.append(y)
# plot point charges
for C in charges:
if C.q>0:
plot(C.pos[0], C.pos[1], 'bo', ms=8*sqrt(C.q))
if C.q<0:
plot(C.pos[0], C.pos[1], 'ro', ms=8*sqrt(-C.q))
c4science · Help