Particle-in-cell simulacija trajektorije nabitega delca.
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

223 lines
6.1 KiB

#!/usr/bin/python3
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
#rc('font',**{'family':'serif','serif':['Computer Modern']})
#rc('text', usetex=True)
from scipy.special import ellipe,ellipk
e = 1.602176565e-19 #osnovni naboj [C]
m_pr = 1.672621777e-27 #masa protona [kg]
m_el = 9.10938291e-31 #masa elektrona [kg]
c = 299792458 #hitrost svetlobe [m/s]
def boris(x0, v0, E, B, dt, tdur, q=1, m=1):
''' Borisov algoritem zabjega skoka za
integracijo diferencialne enačbe
x0 = [x,y,z] : vektor zacetne pozicije
v0 = [vx,vy,vz] : vektor zacetne hitrosti
E : funkcija, ki daje jakost polja E(x)
B : funkcija, ki daje jakost polja B(x)
dt: casovni korak '''
duration = int(tdur/dt)
X = np.zeros((duration,3))
V = np.zeros((duration,3))
x = x0
v = v0
u = v0 / np.sqrt(1-(np.sum(v0**2)/c**2))
q_prime = dt * q * 0.5 / m
# https://en.wikipedia.org/wiki/Particle-in-cell
# https://warpx.readthedocs.io/en/latest/theory/picsar_theory.html
for time in range(duration):
u_min = u + q_prime*E(x)
t = q_prime*B(x)/np.sqrt(1+np.sum(u_min**2)/c**2)
u_prime = u_min + np.cross(u_min,t)
u_plu = u_min + np.cross(u_prime,2*t)/(1+np.sum(t**2))
u = u_plu + q_prime*E(x)
v = u/np.sqrt(1+np.sum(u_min**2)/c**2)
x = x + dt*v
V[time,:] = v
X[time,:] = x
#implementacija po wikipediji - nerelativisticno
#h = q_prime*B(x)
#s = 2*h / (1+np.dot(h,h))
#u = v + q_prime*E(x)
#u_ = u + np.cross(u+np.cross(u,h),s)
#
#v = u_ + q_prime * E(x)
#x = x + dt*v
#V[time,:] = v
#X[time,:] = x
return X, V
def vay(x0, v0, E, B, dt, tdur, q=1, m=1):
''' doi: 19,1963/1.2837054 '''
print("ni implementiran")
return 0
def B_zemlja(x):
''' model zemljinega magnentnega polja aproksimiran z
magnetnim dipolom '''
B0 = 3.07e-5 # [T]
Rz = 6378137 # radii zemlje [m]
r = np.sqrt(np.sum(x**2))
k = -B0*Rz**3/r**5
return k*np.array([3*x[0]*x[2],3*x[1]*x[2],2*x[2]**2-x[0]**2-x[1]**2])
def B_loop(p, a):
''' polje magnetne zanke polmera a v xy ravnini, normiran tok 1A
glej: Sinigoj, ELMG polje, stran 200'''
mu = 2e-7 # mu/2/pi
x,y,z = p
r = np.sqrt(x**2+y**2)
if r < 1e-15:
Bz = 2*np.pi*1e-7*a**2/(a**2+z**2)**1.5
return np.array([0.,0.,Bz])
m = np.sqrt(4*a*r/( (r+a)**2+z**2))
K = ellipk(m**2)
E = ellipe(m**2)
C = 1./np.sqrt((a+r)**2+z**2)
Br = -mu * z/r*C * (-K + (a**2+r**2+z**2)*E/((a-r)**2+z**2))
Bz = mu*C * (K - (a**2-r**2-z**2)*E/((a-r)**2+z**2))
return np.array([Br*x/r, Br*y/r, Bz])
def B_bottle(p, a, h, I):
''' polje magnetne steklenice visine 2h in polmera a
'''
B1 = B_loop(p-np.array([0.,0.,h]),a)*I
B2 = B_loop(p+np.array([0.,0.,h]),a)*I
return B1+B2
def tokamak(p, a, b, I, n):
''' polje v modelu tokamaka, zanke polmera a
nataknjene na toroid polmera b '''
Rx = np.array([[1.,0.,0.],[0.,0.,-1.],[0.,1.,0.]])
B = np.zeros(3)
premik = np.array([b,0.,0.])
for i in range(n):
theta = i*2.*np.pi/n
Rz = np.array([[np.cos(theta), np.sin(theta),0.],
[-np.sin(theta),np.cos(theta), 0.],
[0.,0.,1.]])
#print(Rz)
p_ = np.matmul(Rx.T, np.matmul(Rz.T,p)) + premik
b_ = I*B_loop(p_,a)
b = np.matmul(Rz, np.matmul(Rx,b_))
B = B+b
return B
def Wk(m,V):
return 0.5*m*np.sum(V**2,1)
def plot3(X,enote='m',zemlja=False,lim=[],bottle=[],tokamak=[]):
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
if zemlja:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x,y,z,zorder=-4)
#ax.set_box_aspect(np.ptp(X,axis=0))
ax.plot(X[:,0],X[:,1],X[:,2],zorder=4)
ax.set_xlabel(r'$x$ [{}]'.format(enote))
ax.set_ylabel(r'$y$ [{}]'.format(enote))
ax.set_zlabel(r'$z$ [{}]'.format(enote))
if lim:
ax.set_xlim(lim)
ax.set_ylim(lim)
ax.set_zlim(lim)
if bottle:
r = bottle[0]
h = bottle[1]
fi = np.linspace(0,2*np.pi,300)
a = r*np.cos(fi)
b = r*np.sin(fi)
ax.plot(a,b,+h,color='red',zorder=100)
ax.plot(a,b,-h,color='red',zorder=100)
if tokamak:
ax=plot_tokamak(tokamak[0],tokamak[1],tokamak[2],ax)
plt.show()
def plot_tokamak(a,b,n,ax=None):
fi = np.linspace(0,2.*np.pi,300)
x = a*np.cos(fi)
y = a*np.sin(fi)
z = np.zeros(300)
X=np.zeros([300,3,n])
x = x-b
Rx = np.array([[1.,0.,0.],[0.,0.,-1.],[0.,1.,0.]])
if ax == None:
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
for i in range(n):
theta = i*2.*np.pi/n
Rz = np.array([[np.cos(theta),-np.sin(theta),0.],
[np.sin(theta),np.cos(theta),0.],
[0.,0.,1]])
for j in range(300):
p = np.array([x[j],y[j],z[j]])
#X=np.zeros([300,3])
X[j,:,i] = np.matmul(Rz, np.matmul(Rx,p.T))
ax.plot(X[:,0,i],X[:,1,i],X[:,2,i],color='red',zorder=1)
plt.draw()
plt.show()
def rplot(X, enote='m', xy=[0,1],os='xy', fname=None):
fig = plt.figure()
ax = fig.subplots()
ax.plot(X[:,xy[0]],X[:,xy[1]], 'k')
ax.grid()
dx = X[-1,xy[0]]-X[-2,xy[0]]
dy = X[-1,xy[1]]-X[-2,xy[1]]
#ax.arrow(X[-1,xy[0]],X[-1,xy[1]],dx,dy,head_width=0.1,head_length=0.2,
# overhang=0.3, zorder=10,color='k')
ax.set_xlabel(r'${}$ [{}]'.format(os[0],enote))
ax.set_ylabel(r'${}$ [{}]'.format(os[1],enote))
if fname:
plt.savefig(fname,bbox_inches='tight')
plt.show()
if __name__ == "__main__": # testni del kode
E = lambda x: np.array([0.0,0.2,0.3])
B = lambda x: np.array([0.,0.0,1.*x[0]])
x0 = np.array([-1.,0.,0.])
v0 = np.array([0.,1.,0.])
dt = 1e-2
#X, V = boris(x0,v0,E,B,dt,20)
#plot3(X)
plot_tokamak(1,2,16)